{ ================================================================================================ }
{ unité contenant les procedures de calcul de reflectivité                                         }
{ ================================================================================================ }
{ 23.02.02 : FO : * the normalisation factor is included in the reflectivity procedure calculation }
{                   that does not make sense: or is at least dangerous                             }
{ 16.04.02 : FO : * sortie du facteur de normalisation                                             }
{ 02.07.02 : FO : * nouvelle procedure CalculateReflectivity(var MyExpPar:TExperimentalParameters;var MyML:TMultiLayer;var DataCalc:TDataReflectivity); }
{ ================================================================================================ }
{ A FAIRE :                                                                                }
{ * corriger le calcul de la reflectivité en champ fort                                    }
{ * simplifier encore les paramètres d'entree sortie des procédures de calcul              }
{   ne plus sortir R,T, Ru etc individuellement                                            }
{ * ACHTUNG : pour le moment la convolution ne fonctionne que pour des data lues avec des abscisses en degrés   }
{ si on a des abscisses en Qz, comment reconstruit-on le spectre en Qz a partir des vecteurs "lambda" et "theta"?}
{ * la rugosité n'est pas incluse dans CalculateReflectivitySF ; A FAIRE !!                }

unit CalculReflec;

interface

uses constants, dataReflec, ndspl, Math387, MtxVec, Dialogs;

{ procédures de calcul de base }
procedure CalculReflecNP(var Qz:Tvec;var z,V,Vi,sigma:Tvec;var R,T:Tvec);
procedure CalculAmplitudeNP(var Qz:Tvec;var z,V,Vi,sigma:Tvec;var AR,AT:Tvec);
procedure CalculReflecSF(var Qz:Tvec;var z,Vu,Vd,Vi,phi:Tvec;var Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu:Tvec);
procedure CalculReflecSF_HF_V10(var lambda:double;var angles:Tvec;var z,Vu,Vd,sigma,BeffMag,BeffTheta,BeffPhi:Tvec;var Ruu,Rdu,Rud,Rdd:Tvec);
procedure CalculReflecSF_HF_V11(var lambda:double;var angles:Tvec;var z,Vu,Vd,sigma,BeffMag,BeffTheta,BeffPhi:Tvec;var Ruu,Rdu,Rud,Rdd:Tvec);
procedure CalculReflecSF_HF_V12(var angles:Tvec;Qz:Tvec;var z,Vu,Vd,sigma,BeffMag,BeffTheta,BeffPhi:Tvec;var Ruu,Rdu,Rud,Rdd:Tvec);
procedure CalculReflecSF_HF_V13(var angles:Tvec;Qz:Tvec;var z,Vu,Vd,sigma,BeffMag,BeffTheta,BeffPhi:Tvec;var Ruu,Rdu,Rud,Rdd:Tvec);
procedure CalculReflecSF_HF(var angles:Tvec;Qz:Tvec;var z,Vu,Vd,Vi,sigma,BeffMag,BeffTheta,BeffPhi:Tvec;var Ruu,Rdu,Rud,Rdd:Tvec);
{ procédures de calcul tenant compte de la résolution }
procedure CalculateReflectivityNP(var DataCalc:TDataReflectivity;var z,V,Vi,sigma:Tvec;var R,T:Tvec);
procedure CalculateReflectivityP(var DataCalc:TDataReflectivity;var z,Vup,Vdo,Vi,sigma:Tvec;var Rup,Rdo,Tup,Tdo:Tvec);
procedure CalculateReflectivitySF(var DataCalc:TDataReflectivity;var z,Vup,Vdo,Vi,phi:Tvec;var Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu:Tvec);
procedure CalculateReflectivitySF_HF(var DataCalc:TDataReflectivity;var z,Vup,Vdo,Vi,BeffMag,BeffTheta,BeffPhi,sigma:Tvec;var Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu:Tvec);
procedure CalculateReflectivity(var MyExpPar:TExperimentalParameters;var MyML:TMultiLayer;var MyDC:TDataReflectivity);

procedure CalculReflecNPamplitudes(Qz:double;var z,V,Vi,sigma:Tvec;var kz,AR,AT:Tvec);

implementation



{ version 2.7                                                                   }
{ dernière version fonctionnelle et testée                                      }
{ procédure de calcul de la réflectivité (neut. non polarisés)}
{ NEW : au lieu d'utilser l'angle d'incidence,, on va rentrer un vecteur de vecteurs de diffusion }
{       c'est en vue des modifications pour le calcul des convolutions expérimentales             }
{ Cette procédure est fonctionnelle quelle que soit la situation NML ou ML }
{ INPUT: qz, vecteur des vecteurs d'onde incidents (en nm-1)  }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vetceur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ *** en principe cette procédure n'utilise aucun paramètre implicite (ou ne devrait pas au moins) *** }
procedure CalculReflecNP_V28(var Qz:Tvec;var z,V,Vi,sigma:Tvec;var R,T:Tvec);
var i,n,nn:integer;
    nbrLayers,nbrML:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    zz:double;          { norme du vecteur d'onde incident  }
    q0:double;         { composante perpendiculaire du vecteur d'onde pour les différents angles d'incidence }
    qq0,qq1,zz0,aux,auxC1,auxC2:TCplx;
    pp,mm:Tcplx;
    aux1,aux2,aux1DW,aux2DW:TCplx;
    sig:double;
    aux1r,aux2r:double;
    Eq,Eim,q:Tvec;
    DD,M,Maux:TMtx;
begin
   createIt(DD,M,Maux);
   createIt(Eq,Eim,q);
   DD.size(2,2,true);M.size(2,2,true);Maux.size(2,2,true);
   nbrLayers := z.length;
   nbrPoints := Qz.length;
   R.Size(Qz);T.Size(Qz); { est-ce bien utile? }
   Eq.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0); { a quoi ca sert ???? }
   q.Size(nbrLayers,true); { les composantes de q peuvent etre complexes }
   for i:=0 to nbrPoints-1 do
   begin
      M.eye(2,2,true); { initialisation de la récursion }
      q0 := Qz.values[i]/2; { Qz est le vecteur de diffusion; q0 est la projection de k0 }
      Eq.setVal(hb2over2mn);Eq.scale(q0*q0);Eq.sub(V);
      Eim.setval(0);Eim.sub(Vi); { *** NEW *** prise en compte de la partie complexe du potentiel *** }
      q.RealToCplx(Eq,Eim);q.scale(1/hb2over2mn);q.sqrt;q.scale(-1); { q peut-etre complexe, racine }
      { on applique le facteur -1 pour des raisons de convention de  l'axe Oz }
      { initialisation de variables pour la récursion }
      qq1:=q.Cvalues[nbrLayers-1];
      zz:=0;
      n:=nbrLayers-2;
      for n:=nbrLayers-2 downto 0 do   (* boucle du vide vers le substrat *)
      begin
         sig:=sigma.values[n]; { rugosité de l'interface j+1/j }
         qq0:=q.Cvalues[n];
         zz:=z.values[n]; { epaisseur de la couche considérée }
         zz0:=Cplx(zz,0);      { passage en complexe pour faire les multiplications }
         { détermination de la matrice de transfert }
         { coefficients p et m }
         pp:=Cdiv(Cadd(qq1,qq0),qq1);mm:=Cdiv(Csub(qq1,qq0),qq1);
         { TEMPORAIRE TANT QUE SCALE N'EST PAS REGLE }
         pp:=CMul(Cplx(0.5,0),pp);mm:=CMul(Cplx(0.5,0),mm);
         { calcul des coefficients de rugosité }
         aux1DW :=exp(Cmul(Csqr(Csub(qq0,qq1)),Cplx(-sig*sig/2,0)));
         aux2DW :=exp(Cmul(Csqr(Cadd(qq0,qq1)),Cplx(-sig*sig/2,0)));
         { initialisation of the D matrix}
         aux1 := Exp(Cmul(C_I,Cmul(qq0,zz0)));        { coefficients de la matrice de propagation }
         aux2 := CDiv(C_ONE,aux1);                    { le signe est ok parce que qq0 est negatif ! }
         DD.Cvalues[0,0]:=Cmul(Cmul(pp,aux1),aux1DW);
         DD.Cvalues[0,1]:=Cmul(Cmul(mm,aux2),aux2DW);
         DD.Cvalues[1,0]:=Cmul(Cmul(mm,aux1),aux2DW);
         DD.Cvalues[1,1]:=Cmul(Cmul(pp,aux2),aux1DW);
{         DD.scale(Cplx(0.5,0)); } { CA MARCHE PAS !!!! }
         Maux.Mul(M,DD);M.copy(Maux);
         qq1:=qq0;
      end;
      aux1r:=Sqr(Cabs(M.Cvalues[0,0]));
      aux2r:=Sqr(Cabs(M.Cvalues[1,0]));
      R.values[i]:=aux2r/aux1r;
      T.values[i]:=(1/aux1r);
   end;
{   R.scale(NormalisationFactor);
   T.scale(NormalisationFactor); }
   freeIt(DD,M,Maux);
   freeIt(Eq,Eim,q);
end;

{ version 2.8                                                 }
{ procédure de calcul de la réflectivité (neut. non polarisés)}
{ NEW : au lieu d'utilser l'angle d'incidence,, on va rentrer un vecteur de vecteurs de diffusion }
{       c'est en vue des modifications pour le calcul des convolutions expérimentales             }
{ Cette procédure est fonctionnelle quelle que soit la situation NML ou ML }
{ INPUT: qz, vecteur des vecteurs d'onde incidents (en nm-1)  }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vetceur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ *** en principe cette procédure n'utilise aucun paramètre implicite (ou ne devrait pas au moins) *** }
{ on va essayer d'accelerer cette procédure                                     }
{ * l'utilisation de C_half et scale ne permet pas de gagner de temps mesurable }
{ * aucun gain avec Caux:=Cplx(-sig*sig/2,0);                                   }
{ * creation de DD, M, Maux : n'utilise que 5% du temps                         }
{ * l'essentiel du temps est passé dans les calculs...                          }
{ * 20% du temps passe sur la rugosité                                          }
{ La seule solution pour accélérer le calcul c'est d'utiliser les formules de récurrences explicites }
{ A FAIRE au moins dans le cas non polarisé                                     }
{ utiliser le bouquin orange                                                    }
procedure CalculReflecNP(var Qz:Tvec;var z,V,Vi,sigma:Tvec;var R,T:Tvec);
var i,n,nn:integer;
    nbrLayers,nbrML:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    zz:double;          { norme du vecteur d'onde incident  }
    q0:double;         { composante perpendiculaire du vecteur d'onde pour les différents angles d'incidence }
    qq0,qq1,zz0,aux,auxC1,auxC2:TCplx;
    pp,mm:Tcplx;
    aux1,aux2,aux1DW,aux2DW:TCplx;
    sig:double;
    aux1r,aux2r:double;
    Eq,Eim,q:Tvec;
    DD,M,Maux:TMtx;
    C_half,Caux : TCplx;
begin
   C_half:=Cplx(0.5,0);

{ on va créer ces matrices par défaut au lancement de la procédure              }
{ plutôt que de les créer à chaque appel de CalculReflec                        }
   createIt(DD,M,Maux);
   DD.size(2,2,true);M.size(2,2,true);Maux.size(2,2,true);

   createIt(Eq,Eim,q);
   nbrLayers := z.length;
   nbrPoints := Qz.length;
   R.Size(Qz);T.Size(Qz); { est-ce bien utile? }
   Eq.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0); { a quoi ca sert ???? }
   q.Size(nbrLayers,true); { les composantes de q peuvent etre complexes }
   for i:=0 to nbrPoints-1 do
   begin
      M.eye(2,2,true); { initialisation de la récursion }
      q0 := Qz.values[i]/2; { Qz est le vecteur de diffusion; q0 est la projection de k0 }
      Eq.setVal(hb2over2mn);Eq.scale(q0*q0);Eq.sub(V);
      Eim.setval(0);Eim.sub(Vi); { *** NEW *** prise en compte de la partie complexe du potentiel *** }
      q.RealToCplx(Eq,Eim);q.scale(1/hb2over2mn);q.sqrt;q.scale(-1); { q peut-etre complexe, racine }
      { on applique le facteur -1 pour des raisons de convention de  l'axe Oz }
      { initialisation de variables pour la récursion }
      qq1:=q.Cvalues[nbrLayers-1];
      zz:=0;
      n:=nbrLayers-2;
      for n:=nbrLayers-2 downto 0 do   (* boucle du vide vers le substrat *)
      begin   { cette boucle représente 70% du temps de calcul }
         { le temps sur ces 4 lignes est non mesurable }
         sig:=sigma.values[n]; { rugosité de l'interface j+1/j }
         qq0:=q.Cvalues[n];
         zz:=z.values[n]; { epaisseur de la couche considérée }
         zz0:=Cplx(zz,0);      { passage en complexe pour faire les multiplications }
         { détermination de la matrice de transfert }
         { coefficients p et m }
         pp:=Cdiv(Cadd(qq1,qq0),qq1);mm:=Cdiv(Csub(qq1,qq0),qq1);
         { calcul des coefficients de rugosité }
         Caux:=Cplx(-sig*sig/2,0); { 0.4ms }
         { les 2 lignes ci-dessous consomment 2 ms sur 11.4 }
         { est-ce que ce calcul est juste                       }
         { cela a-t-il un sens d'utiliser des valeurs complexes pour calculer le facteur de DW ? }
         aux1DW :=exp(Cmul(Csqr(Csub(qq0,qq1)),Caux));
         aux2DW :=exp(Cmul(Csqr(Cadd(qq0,qq1)),Caux));
         { initialisation of the D matrix}
         { les 2 lignes ci-dessous consomment 1 ms sur 11.4 }
         aux1 := Exp(Cmul(C_I,Cmul(qq0,zz0)));        { coefficients de la matrice de propagation }
         aux2 := CDiv(C_ONE,aux1);                    { le signe est ok parce que qq0 est negatif ! }
         { les 5 lignes ci-dessous consomment 0.8 ms sur 11.4 }
         DD.Cvalues[0,0]:=Cmul(Cmul(pp,aux1),aux1DW);
         DD.Cvalues[0,1]:=Cmul(Cmul(mm,aux2),aux2DW);
         DD.Cvalues[1,0]:=Cmul(Cmul(mm,aux1),aux2DW);
         DD.Cvalues[1,1]:=Cmul(Cmul(pp,aux2),aux1DW);
         DD.scale(C_half);
         { les 2 lignes ci-dessous consomment 4.4 ms sur 11.4 }
         Maux.Mul(M,DD);M.copy(Maux);
         qq1:=qq0;
      end;
      aux1r:=Sqr(Cabs(M.Cvalues[0,0]));
      aux2r:=Sqr(Cabs(M.Cvalues[1,0]));
      R.values[i]:=aux2r/aux1r;
      T.values[i]:=(1/aux1r);
   end;
{   R.scale(NormalisationFactor);
   T.scale(NormalisationFactor); }
   freeIt(DD,M,Maux);
   freeIt(Eq,Eim,q);
end;

{ procédure de calcul des AMPLITUDES COMPLEXES de réflectivité}
{ neutrons non polarisés                                      }
{ INPUT: longueur d'onde en nm                                }
{        angles, vecteur des angles d'incidence (en  °)       }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vetceur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ REM: aucun paramètre implicite!!                            }
procedure CalculAmplitudeNP(var Qz:Tvec;var z,V,Vi,sigma:Tvec;var AR,AT:Tvec);
var i,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    k:double;          { norme du vecteur d'onde incident  }
    q0:double;         { composante perpendiculaire du vecteur d'onde pour les différents angles d'incidence }
    D,P,M,Maux:TMtx;
    Eq,Eim:Tvec;  { Energie cinétique perp. dans chaque couche }
    q:Tvec;            { vecteur d'onde perpendiculaire dans chaque couche }
    qq,aux:TCplx;
    Caux:TCplx;
begin
   createIt(P,D,M,Maux);
   createIt(Eq,Eim,q);
   D.size(2,2,true);P.size(2,2,true);M.size(2,2,true);Maux.size(2,2,true);
   nbrLayers := z.length;
   nbrPoints := Qz.length;
   AR.Size(Qz.length,true);AT.Size(Qz.length,true);
   Eq.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0);
   q.Size(nbrLayers,true); { les composantes de q peuvent etre complexes }
   for i:=0 to nbrPoints-1 do
   begin
      q0 := Qz.values[i]/2;
      Eq.setVal(hb2over2mn);Eq.scale(q0*q0);Eq.sub(V);
      q.RealToCplx(Eq,Eim);q.scale(1/hb2over2mn);q.sqrt;q.scale(-1); { q peut-etre complexe, racine }
      for n:=nbrLayers-1 downto 0 do   (* boucle du vide vers le substrat *)
      begin
      qq:=q.Cvalues[n];
      { initialisation of the D matrix}
      D.SetIt([1,0,1,0,qq.Re,qq.Im,-qq.Re,-qq.Im]);
      { initialisation of the P matrix }
      aux := Exp(Cmul(Cmul(Neg(C_I),qq),Cplx(z.values[n],0)));
      P.setZero;
      P.Cvalues[0,0]:=aux;
      P.Cvalues[1,1]:=CDiv(C_ONE,aux);
      If (n = nbrLayers-1)
      then begin
           M.copy(D);   { cas du vide, initialisation de M }
           M.inv;
           end
      else If (n = 0)
           then begin
                Maux.Mul(M,D); { on est arrivé au substrat }
                M.copy(Maux);
                end
           else begin
                Maux.Mul(M,D);
                M.Mul(Maux,P);
                Maux.Mul(M,D.Inv);
                M.copy(Maux);
                end;
      end;
{      viewValues(D);       ca ne semble pas marcher avec les matrices
      viewValues(P);      viewValues(M);      viewValues(Maux); }
      Caux:=Cdiv(C_ONE,M.Cvalues[0,0]);
      AR.Cvalues[i]:=Cmul(M.Cvalues[1,0],Caux);
      AT.Cvalues[i]:=Caux;
   end;
   freeIt(P,D,M,Maux);
   freeIt(Eq,Eim,q);
end;

{ version 0.1                                                 }
{ procédure de calcul de la réflectivité (neut. non polarisés)}
{ en tenant compte du spin-flip                               }
{ on s'inspire de Blundell PRB 46 (1992) 3391.                }
{ rugosité (DW) non prise en compte !!!                       }
{ on derive de la version calculReflecNP V1.0                 }
{ cad matrices locales et recursion 100% explicite            }
{     et pas de rugosité                                      }
{ INPUT: longueur d'onde en nm                                }
{        angles, vecteur des angles d'incidence (en  °)       }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vetceur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ si on ne met pas de résolution expérimentale, on a parfois des creux très marqués }
{ dans les courbes de réflectivité                                                  }
{ verifier que les resultats sont OK avec d'autres programmes de simulation         }
{ de manière générale, les suffices "u" et "d" designent up et down                 }
procedure CalculReflecSF_V01(var Qz:Tvec;var z,Vu,Vd,Vi,phi:Tvec;var Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu:Tvec);
var i,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    k:double;          { norme du vecteur d'onde incident  }
    q0:double;         { composante perpendiculaire du vecteur d'onde pour les différents angles d'incidence }
    D,P,Rot,M,Maux:TMtx;
    EqU,EqD,Eim:Tvec;       { Energie cinétique perp. dans chaque couche }
    qU,qD:Tvec;            { vecteur d'onde perpendiculaire dans chaque couche }
    qqU,qqD:TCplx;
    aux1,aux2:double;
    auxU,auxD,aux11,aux13,aux21,aux23,aux31,aux33,aux41,aux43:Tcplx;
begin
   createIt(P,D,Rot);createIt(M,Maux);
   createIt(EqU,EqD,Eim);
   createIt(qU,qD);
   { les matrices de la recursion deviennent des matrices 4x4 }
   D.size(4,4,true);P.size(4,4,true);Rot.size(4,4,true);M.size(4,4,true);Maux.size(4,4,true);
   nbrLayers := z.length;
   nbrPoints := Qz.length;
   Ruu.Size(Qz);Tuu.Size(Qz);Rdd.Size(Qz);Tdd.Size(Qz);
   Rud.Size(Qz);Tud.Size(Qz);Rdu.Size(Qz);Tdu.Size(Qz);
   { calcul des composantes du vecteur d'onde dans les différentes couches }
   EqU.Size(nbrLayers,false);EqD.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0);
   qu.Size(nbrLayers,true); { les composantes de q peuvent etre complexes }
   qd.Size(nbrLayers,true);
   for i:=0 to nbrPoints-1 do
   begin
      q0 := Qz.values[i]/2; { Qz est le vecteur de diffusion; q0 est la projection de k0 }
      EqU.setVal(hb2over2mn);EqU.scale(q0*q0);EqU.sub(Vu);
      EqD.setVal(hb2over2mn);EqD.scale(q0*q0);EqD.sub(Vd);
      Eim.setval(0);Eim.sub(Vi);
      qu.RealToCplx(EqU,Eim);qu.scale(1/hb2over2mn);qu.sqrt; { SIGNE A VERIFIER } { q peut-etre complexe, racine }
      qd.RealToCplx(EqD,Eim);qd.scale(1/hb2over2mn);qd.sqrt; { q peut-etre complexe, racine }
      for n:=nbrLayers-1 downto 0 do   (* boucle du vide vers le substrat *)
      { N.B. chez Blundell "1" is vacuum and "N" is substrate }
      { opposite to my convention                             }
      begin
         qqU:=qU.Cvalues[n];qqD:=qD.Cvalues[n];
         { initialisation of the D matrix}
         D.SetIt([1,0,1,0,0,0,0,0,qqU.Re,qqU.Im,-qqU.Re,-qqU.Im,0,0,0,0,
                  0,0,0,0,1,0,1,0,0,0,0,0,qqD.Re,qqD.Im,-qqD.Re,-qqD.Im]);
         { initialisation of the P matrix }
         { attention au signe etant donné les convention z+ }
         auxU := Exp(Cmul(Cmul(Neg(C_I),qqU),Cplx(z.values[n],0)));
         auxD := Exp(Cmul(Cmul(Neg(C_I),qqD),Cplx(z.values[n],0)));
         P.setZero;
         P.Cvalues[0,0]:=auxU;P.Cvalues[1,1]:=CDiv(C_ONE,auxU);
         P.Cvalues[2,2]:=auxD;P.Cvalues[3,3]:=CDiv(C_ONE,auxD);
         { initialisation de la matrice de rotation }
         { Attention au signe     }
         if n<>0 then { quand on est dans le substrat, il n'y a plus a definir de matrice de rotation }
         begin
         aux1:=cos((phi.values[n-1]-phi.values[n])/2);aux2:=sin((phi.values[n-1]-phi.values[n])/2);
         Rot.SetIt([aux1,0,0,0,aux2,0,0,0,
                    0,0,aux1,0,0,0,aux2,0,
                    -aux2,0,0,0,aux1,0,0,0,
                    0,0,-aux2,0,0,0,aux1,0]);
         end;
      If (n = nbrLayers-1)
      then begin
           M.copy(D);   { cas du vide, initialisation de M }
           M.inv;
           Maux.Mul(M,Rot);M.copy(Maux);
           end
      else If (n = 0)
           then begin
                Maux.Mul(M,D); { on est arrivé au substrat }
                M.copy(Maux);
                end
           else begin
                Maux.Mul(M,D);
                M.Mul(Maux,P);
                Maux.Mul(M,D.Inv);
                M.mul(Maux,Rot);
                end;
      end;
{      viewValues(D);       ca ne semble pas marcher avec les matrices }
      aux11:=M.Cvalues[0,0];aux13:=M.Cvalues[0,2];
      aux21:=M.Cvalues[1,0];aux23:=M.Cvalues[1,2];
      aux31:=M.Cvalues[2,0];aux33:=M.Cvalues[2,2];
      aux41:=M.Cvalues[3,0];aux43:=M.Cvalues[3,2];
      auxD:=Csub(Cmul(aux11,aux33),Cmul(aux13,aux31));
      Ruu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux21,aux33),Cmul(aux23,aux31)),auxD)));
      Rdd.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux43,aux11),Cmul(aux41,aux13)),auxD)));
      Rud.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux41,aux33),Cmul(aux43,aux31)),auxD)));
      Rdu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux23,aux11),Cmul(aux21,aux13)),auxD)));
      Tuu.values[i]:=Sqr(Cabs(Cdiv(aux33,auxD)));
      Tdd.values[i]:=Sqr(Cabs(Cdiv(aux11,auxD)));
      Tud.values[i]:=Sqr(Cabs(Cdiv(Csub(C_ZERO,aux31),auxD)));
      Tdu.values[i]:=Sqr(Cabs(Cdiv(Csub(C_ZERO,aux13),auxD)));
   end;
   freeIt(P,D,Rot);
   freeit(M,Maux);freeit(qU,qD);
   freeIt(EqU,EqD,Eim);

end;

{ version 0.2                                                 }
{ procédure de calcul de la réflectivité (neut. non polarisés)}
{ en tenant compte du spin-flip                               }
{ on s'inspire de Blundell PRB 46 (1992) 3391.                }
{ rugosité (DW) non prise en compte !!!                       }
{ on derive de la version calculReflecNP V1.0                 }
{ cad matrices locales et recursion 100% explicite            }
{     et pas de rugosité                                      }
{ V0.2: Suppression des inversions de matrices = gain d'un facteur 4 en vitesse % V01}
{ il reste une fuite de mémoire quelque part...               }
{ INPUT: longueur d'onde en nm                                }
{        angles, vecteur des angles d'incidence (en  °)       }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vetceur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ si on ne met pas de résolution expérimentale, on a parfois des creux très marqués }
{ dans les courbes de réflectivité                                                  }
{ verifier que les resultats sont OK avec d'autres programmes de simulation         }
{ de manière générale, les suffices "u" et "d" designent up et down                 }
procedure CalculReflecSF_V02(var Qz:Tvec;var z,Vu,Vd,Vi,phi:Tvec;var Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu:Tvec);
var i,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    k:double;          { norme du vecteur d'onde incident  }
    q0:double;         { composante perpendiculaire du vecteur d'onde pour les différents angles d'incidence }
    D,invD,P,Rot,M,Maux:TMtx;
    EqU,EqD,Eim:Tvec;       { Energie cinétique perp. dans chaque couche }
    qU,qD:Tvec;            { vecteur d'onde perpendiculaire dans chaque couche }
    qqU,qqD:TCplx;
    aux1,aux2:double;
    auxU,auxD,aux11,aux13,aux21,aux23,aux31,aux33,aux41,aux43:Tcplx;
    auxUr,auxUi,auxUm,auxDr,auxDi,auxDm:double;
begin
   createIt(P,D,invD,Rot);createIt(M,Maux);
   createIt(EqU,EqD,Eim);
   createIt(qU,qD);
   { les matrices de la recursion deviennent des matrices 4x4 }
   D.size(4,4,true);invD.size(4,4,true);P.size(4,4,true);Rot.size(4,4,true);M.size(4,4,true);Maux.size(4,4,true);
   nbrLayers := z.length;
   nbrPoints := Qz.length;
   Ruu.Size(Qz);Tuu.Size(Qz);Rdd.Size(Qz);Tdd.Size(Qz);
   Rud.Size(Qz);Tud.Size(Qz);Rdu.Size(Qz);Tdu.Size(Qz);
   { calcul des composantes du vecteur d'onde dans les différentes couches }
   EqU.Size(nbrLayers,false);EqD.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0);
   qu.Size(nbrLayers,true); { les composantes de q peuvent etre complexes }
   qd.Size(nbrLayers,true);
   for i:=0 to nbrPoints-1 do
   begin
      q0 := Qz.values[i]/2; { Qz est le vecteur de diffusion; q0 est la projection de k0 }
      EqU.setVal(hb2over2mn);EqU.scale(q0*q0);EqU.sub(Vu);
      EqD.setVal(hb2over2mn);EqD.scale(q0*q0);EqD.sub(Vd);
      Eim.setval(0);Eim.sub(Vi);
      qu.RealToCplx(EqU,Eim);qu.scale(1/hb2over2mn);qu.sqrt; { SIGNE A VERIFIER } { q peut-etre complexe, racine }
      qd.RealToCplx(EqD,Eim);qd.scale(1/hb2over2mn);qd.sqrt; { q peut-etre complexe, racine }
      for n:=nbrLayers-1 downto 0 do   (* boucle du vide vers le substrat *)
      { N.B. chez Blundell "1" is vacuum and "N" is substrate }
      { opposite to my convention                             }
      begin
         qqU:=qU.Cvalues[n];qqD:=qD.Cvalues[n];
         auxUr:=qqU.Re;auxUi:=qqU.Im;auxUm:=auxUr*auxUr+auxUi*auxUi;
         auxDr:=qqD.Re;auxDi:=qqD.Im;auxDm:=auxDr*auxDr+auxDi*auxDi;
         { initialisation of the D matrix}
         D.SetIt([1,0,1,0,0,0,0,0,auxUr,auxUi,-auxUr,-auxUi,0,0,0,0,
                  0,0,0,0,1,0,1,0,0,0,0,0,auxDr,auxDi,-auxDr,-auxDi]);
         invD.SetIt([0.5,0,0.5*auxUr/auxUm,-0.5*auxUi/auxUm,0,0,0,0,
                     0.5,0,-0.5*auxUr/auxUm,0.5*auxUi/auxUm,0,0,0,0,
                     0,0,0,0,0.5,0,0.5*auxDr/auxDm,-0.5*auxDi/auxDm,
                     0,0,0,0,0.5,0,-0.5*auxDr/auxDm,0.5*auxDi/auxDm]);
         { initialisation of the P matrix }
         { attention au signe etant donné les convention z+ }
         auxU := Exp(Cmul(Cmul(Neg(C_I),qqU),Cplx(z.values[n],0)));
         auxD := Exp(Cmul(Cmul(Neg(C_I),qqD),Cplx(z.values[n],0)));
         P.setZero;
         P.Cvalues[0,0]:=auxU;P.Cvalues[1,1]:=CDiv(C_ONE,auxU);
         P.Cvalues[2,2]:=auxD;P.Cvalues[3,3]:=CDiv(C_ONE,auxD);
         { initialisation de la matrice de rotation }
         { Attention au signe     }
         if n<>0 then { quand on est dans le substrat, il n'y a plus a definir de matrice de rotation }
         begin
         aux1:=cos((phi.values[n-1]-phi.values[n])/2);aux2:=sin((phi.values[n-1]-phi.values[n])/2);
         Rot.SetIt([aux1,0,0,0,aux2,0,0,0,
                    0,0,aux1,0,0,0,aux2,0,
                    -aux2,0,0,0,aux1,0,0,0,
                    0,0,-aux2,0,0,0,aux1,0]);
         end;
      If (n = nbrLayers-1)
      then begin
           M.copy(invD);   { cas du vide, initialisation de M }
{           M.inv; }
           Maux.Mul(M,Rot);M.copy(Maux);
           end
      else If (n = 0)
           then begin
                Maux.Mul(M,D); { on est arrivé au substrat }
                M.copy(Maux);
                end
           else begin
                Maux.Mul(M,D);
                M.Mul(Maux,P);
                Maux.Mul(M,invD);
                M.mul(Maux,Rot);
                end;
      end;
{      viewValues(D);       ca ne semble pas marcher avec les matrices }
      aux11:=M.Cvalues[0,0];aux13:=M.Cvalues[0,2];
      aux21:=M.Cvalues[1,0];aux23:=M.Cvalues[1,2];
      aux31:=M.Cvalues[2,0];aux33:=M.Cvalues[2,2];
      aux41:=M.Cvalues[3,0];aux43:=M.Cvalues[3,2];
      auxD:=Csub(Cmul(aux11,aux33),Cmul(aux13,aux31));
      Ruu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux21,aux33),Cmul(aux23,aux31)),auxD)));
      Rdd.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux43,aux11),Cmul(aux41,aux13)),auxD)));
      Rud.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux41,aux33),Cmul(aux43,aux31)),auxD)));
      Rdu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux23,aux11),Cmul(aux21,aux13)),auxD)));
      Tuu.values[i]:=Sqr(Cabs(Cdiv(aux33,auxD)));
      Tdd.values[i]:=Sqr(Cabs(Cdiv(aux11,auxD)));
      Tud.values[i]:=Sqr(Cabs(Cdiv(Csub(C_ZERO,aux31),auxD)));
      Tdu.values[i]:=Sqr(Cabs(Cdiv(Csub(C_ZERO,aux13),auxD)));
   end;
   freeIt(P,D,invD,Rot);
   freeit(M,Maux);freeit(qU,qD);
   freeIt(EqU,EqD,Eim);
end;

{ version 0.3                                                 }
{ procédure de calcul de la réflectivité (neut. non polarisés)}
{ en tenant compte du spin-flip                               }
{ on s'inspire de Blundell PRB 46 (1992) 3391.                }
{ rugosité (DW) non prise en compte !!!                       }
{ on derive de la version calculReflecNP V1.0                 }
{ cad matrices locales et recursion 100% explicite            }
{     et pas de rugosité                                      }
{ V0.2: Suppression des inversions de matrices = gain d'un facteur 4 en vitesse % V01}
{ V0.3: Suppression du max de multiplications de matrices = gain d'un facteur 2 supplémentaire }
{ il est probablement possible de gagner en vitesse en évitant les calculs sur les complexes autant que possible }
{ il y a quand meme encore au moins un facteur 10 % au cas polarisé simple }
{ on pourrait essayer de tendre vers le facteur 4 du aux dimensions }
{ INPUT: longueur d'onde en nm                                }
{        angles, vecteur des angles d'incidence (en  °)       }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vetceur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ si on ne met pas de résolution expérimentale, on a parfois des creux très marqués }
{ dans les courbes de réflectivité                                                  }
{ verifier que les resultats sont OK avec d'autres programmes de simulation         }
{ de manière générale, les suffices "u" et "d" designent up et down                 }
procedure CalculReflecSF(var Qz:Tvec;var z,Vu,Vd,Vi,phi:Tvec;var Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu:Tvec);
var i,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    k:double;          { norme du vecteur d'onde incident  }
    q0:double;         { composante perpendiculaire du vecteur d'onde pour les différents angles d'incidence }
    DD,M,Maux:TMtx;
    EqU,EqD,Eim:Tvec;  { Energie cinétique perp. dans chaque couche }
    qU,qD:Tvec;        { vecteur d'onde perpendiculaire dans chaque couche }
    qqU,qqD:TCplx;
    aux1,aux2:double;
    auxU,auxD,auxUp,auxUm,auxDp,auxDm,aux11,aux13,aux21,aux23,aux31,aux33,aux41,aux43:Tcplx;
    auxUr,auxUi,auxUmod,auxDr,auxDi,auxDmod:double;
    aux1CplxP,aux1CplxM,aux2CplxP,aux2CplxM:TCplx;
begin
   createIt(DD,M,Maux);createIt(EqU,EqD,Eim);createIt(qU,qD);
   { les matrices de la recursion deviennent des matrices 4x4 }
   DD.size(4,4,true);M.size(4,4,true);Maux.size(4,4,true);
   nbrLayers := z.length;
   nbrPoints := Qz.length;
   Ruu.Size(Qz);Tuu.Size(Qz);Rdd.Size(Qz);Tdd.Size(Qz);
   Rud.Size(Qz);Tud.Size(Qz);Rdu.Size(Qz);Tdu.Size(Qz);
   { calcul des composantes du vecteur d'onde dans les différentes couches }
   EqU.Size(nbrLayers,false);EqD.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0);
   qu.Size(nbrLayers,true); { les composantes de q peuvent etre complexes }
   qd.Size(nbrLayers,true);
   for i:=0 to nbrPoints-1 do
   begin
      q0 := Qz.values[i]/2; { Qz est le vecteur de diffusion; q0 est la projection de k0 }
      EqU.setVal(hb2over2mn);EqU.scale(q0*q0);EqU.sub(Vu);
      EqD.setVal(hb2over2mn);EqD.scale(q0*q0);EqD.sub(Vd);
      Eim.setval(0);Eim.sub(Vi);
      qu.RealToCplx(EqU,Eim);qu.scale(1/hb2over2mn);qu.sqrt; { SIGNE A VERIFIER } { q peut-etre complexe, racine }
      qd.RealToCplx(EqD,Eim);qd.scale(1/hb2over2mn);qd.sqrt; { q peut-etre complexe, racine }
      { cas du vide, initialisation de la récursion }
      n:=nbrLayers-1;
      qqU:=qU.Cvalues[n];qqD:=qD.Cvalues[n];
      { coefficients de la matrice de rotation } { nota : on inclut facteur 1/2 dans ces coefficients }
      aux1:=0.5*cos((phi.values[n-1]-phi.values[n])/2);aux2:=0.5*sin((phi.values[n-1]-phi.values[n])/2);
      auxUr:=qqU.Re;auxUi:=qqU.Im;auxUmod:=auxUr*auxUr+auxUi*auxUi;
      auxDr:=qqD.Re;auxDi:=qqD.Im;auxDmod:=auxDr*auxDr+auxDi*auxDi;
      DD.SetIt([aux1,0,  aux1*auxUr/auxUmod,-aux1*auxUi/auxUmod, aux2,0, aux2*auxUr/auxUmod,-aux2*auxUi/auxUmod,
                aux1,0,  -aux1*auxUr/auxUmod,aux1*auxUi/auxUmod, aux2,0, -aux2*auxUr/auxUmod,aux2*auxUi/auxUmod,
                -aux2,0, -aux2*auxDr/auxDmod,aux2*auxDi/auxDmod, aux1,0, aux1*auxDr/auxDmod,-aux1*auxDi/auxDmod,
                -aux2,0, aux2*auxDr/auxDmod,-aux2*auxDi/auxDmod, aux1,0, -aux1*auxDr/auxDmod,aux1*auxDi/auxDmod]);
      M.Copy(DD);
      for n:=nbrLayers-2 downto 1 do   (* boucle du vide vers le substrat *)
      { N.B. chez Blundell "1" is vacuum and "N" is substrate }{ opposite to my convention }
      begin
         qqU:=qU.Cvalues[n];qqD:=qD.Cvalues[n];
         { coefficients de la matrice de propagation P }
         { attention au signe etant donné les convention z+ }
         auxU := Exp(Cmul(Cmul(Neg(C_I),qqU),Cplx(z.values[n],0)));
         auxD := Exp(Cmul(Cmul(Neg(C_I),qqD),Cplx(z.values[n],0)));
         auxUp:=Cadd(auxU,CDiv(C_ONE,auxU));auxUm:=Csub(auxU,CDiv(C_ONE,auxU)); { termes (auxU+1/auxU) etc... }
         auxDp:=Cadd(auxD,CDiv(C_ONE,auxD));auxDm:=Csub(auxD,CDiv(C_ONE,auxD));
         { coefficients de la matrice de rotation }
         aux1:=0.5*cos((phi.values[n-1]-phi.values[n])/2);aux2:=0.5*sin((phi.values[n-1]-phi.values[n])/2);
         { ACHTUNG, ACHTUNG, le facteur 0.5 est la parce que .scale ne marche pas correctement sur les matrices complexes }
         { voir la version officielle des coefficients de Rot dans la version V02 }
         aux1CplxP:=Cplx(aux1,0);aux2CplxP:=Cplx(aux2,0);aux1CplxM:=Cplx(-aux1,0);aux2CplxM:=Cplx(-aux2,0);
         { coefficients diagonaux }
         DD.Cvalues[0,0]:=Cmul(aux1CplxP,auxUp);DD.Cvalues[1,1]:=DD.Cvalues[0,0];
         DD.Cvalues[2,2]:=Cmul(aux1CplxP,auxDp);DD.Cvalues[3,3]:=DD.Cvalues[2,2];
         DD.Cvalues[1,0]:=Cmul(qqU,Cmul(aux1CplxP,auxUm));DD.Cvalues[2,0]:=Cmul(aux2CplxM,auxDp);DD.Cvalues[3,0]:=Cmul(qqD,Cmul(aux2CplxM,auxDm));
         DD.Cvalues[0,1]:=Cdiv(Cmul(aux1CplxP,auxUm),qqU);DD.Cvalues[2,1]:=Cdiv(Cmul(aux2CplxM,auxDm),qqD);DD.Cvalues[3,1]:=Cmul(aux2CplxM,auxDp);
         DD.Cvalues[0,2]:=Cmul(aux2CplxP,auxUp);DD.Cvalues[1,2]:=Cmul(qqU,Cmul(aux2CplxP,auxUm));DD.Cvalues[3,2]:=Cmul(qqD,Cmul(aux1CplxP,auxDm));
         DD.Cvalues[0,3]:=Cdiv(Cmul(aux2CplxP,auxUm),qqU);DD.Cvalues[1,3]:=Cmul(aux2CplxP,auxUp);DD.Cvalues[2,3]:=Cdiv(Cmul(aux1CplxP,auxDm),qqD);
         Maux.Mul(M,DD);M.copy(Maux);
      end;
      { cas du substrat }
      qqU:=qU.Cvalues[0];qqD:=qD.Cvalues[0];
      auxUr:=qqU.Re;auxUi:=qqU.Im;auxDr:=qqD.Re;auxDi:=qqD.Im;
      DD.SetIt([1,0,1,0,0,0,0,0,auxUr,auxUi,-auxUr,-auxUi,0,0,0,0,
               0,0,0,0,1,0,1,0,0,0,0,0,auxDr,auxDi,-auxDr,-auxDi]);
      Maux.Mul(M,DD);M.copy(Maux);
{      viewValues(D);       ca ne semble pas marcher avec les matrices }
      aux11:=M.Cvalues[0,0];aux13:=M.Cvalues[0,2];
      aux21:=M.Cvalues[1,0];aux23:=M.Cvalues[1,2];
      aux31:=M.Cvalues[2,0];aux33:=M.Cvalues[2,2];
      aux41:=M.Cvalues[3,0];aux43:=M.Cvalues[3,2];
      auxD:=Csub(Cmul(aux11,aux33),Cmul(aux13,aux31));
      Ruu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux21,aux33),Cmul(aux23,aux31)),auxD)));
      Rdd.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux43,aux11),Cmul(aux41,aux13)),auxD)));
      Rud.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux41,aux33),Cmul(aux43,aux31)),auxD)));
      Rdu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux23,aux11),Cmul(aux21,aux13)),auxD)));
      Tuu.values[i]:=Sqr(Cabs(Cdiv(aux33,auxD)));
      Tdd.values[i]:=Sqr(Cabs(Cdiv(aux11,auxD)));
      Tud.values[i]:=Sqr(Cabs(Cdiv(Csub(C_ZERO,aux31),auxD)));
      Tdu.values[i]:=Sqr(Cabs(Cdiv(Csub(C_ZERO,aux13),auxD)));
   end;
   freeit(DD,M,Maux);freeit(qU,qD);
   freeIt(EqU,EqD,Eim);
end;

{version 0.3    P.Humbert                                    }
{ procédure de calcul de la réflectivité (neut. polarisés)    }
{ on ne tient pas compte des k parallleles                    }
{prise en compte de la rugosite                               }
{ INPUT: longueur d'onde en nm                                }
{        angles, vecteur des angles d'incidence (en  °)       }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vetceur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ REM: aucun paramètre implicite!!                            }
{      dans cette version, les matrices auxilaires D,P,M et Maux       }
{      sont créées à l'initialisation (gain de 10% seulement)          }
procedure CalculReflecSF_HF_V10(var lambda:double;var angles:Tvec;var z,Vu,Vd,sigma,BeffMag,BeffTheta,BeffPhi:Tvec;var Ruu,Rdu,Rud,Rdd:Tvec);
var i,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    k0,zz:double;          { norme du vecteur d'onde incident  }
    E0,q0:double;
    qquu0,qqdu0,qqud0,qqdd0,qquu1,qqdu1,qqud1,qqdd1,zz0:TCplx;
    quu,qdu,qud,qdd:Tvec;            { vecteur d'onde perpendiculaire dans chaque couche }
    Equu,Eqdd,Equd,Eqdu,Eim:Tvec;    { Energie cinétique perp. dans chaque couche }
    Vaux:Tvec;
    ppuu,ppud,ppdu,ppdd,mmuu,mmud,mmdu,mmdd:Tcplx; {coefficient p et m (Gibaud p.101)}
    rruu,rrud,rrdu,rrdd,nnuu,nnud,nndu,nndd:Tcplx;
    auxupuu1,auxupuu2,auxupud1,auxupud2,auxupdu1,auxupdu2,auxupdd1,auxupdd2:Tcplx; {coefficient de propagation}
    auxdouu1,auxdouu2,auxdoud1,auxdoud2,auxdodu1,auxdodu2,auxdodd1,auxdodd2:Tcplx;
    sig:double;
    auxupuu1DW,auxupuu2DW,auxupud1DW,auxupud2DW:Tcplx; {coefficient de rugosité}
    auxupdu1DW,auxupdu2DW,auxupdd1DW,auxupdd2DW:Tcplx;
    auxdouu1DW,auxdouu2DW,auxdoud1Dw,auxdoud2Dw:Tcplx;
    auxdodu1DW,auxdodu2DW,auxdodd1DW,auxdodd2DW:Tcplx;
    DD,M,Maux:TMtx;
    auxupD,aux11,aux13,aux21,aux23,aux31,aux33,aux41,aux43:TCplx;
    auxdoD,aux55,aux57,aux65,aux67,aux75,aux77,aux85,aux87:TCplx;
    tt1,tt2,phi1,phi2:TCplx;
    trigup1,trigup2,trigup3,trigup4:TCplx;
    Equu1,Equd1:Tvec;
    quu1,qud1:Tvec;

begin
   createIt(Equu,Eqdd,Eim);
   createIt(Equd,Eqdu);
   createIt(Vaux);
   createIt(DD);DD.size(8,8,true);
   createIt(M);M.size(8,8,true);
   createIt(Maux);Maux.size(8,8,true);
   nbrLayers := z.length;
   nbrPoints := angles.length;
   k0:=2*Pi/lambda;
   Vaux.Size(nbrLayers,false);
   Equu.Size(nbrLayers,false);Equd.Size(nbrLayers,false);
   Eqdu.Size(nbrLayers,false);Eqdd.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0);
   createIt(quu,qdu,qud,qdd);
   quu.Size(nbrLayers,true);
   qud.Size(nbrLayers,true);
   qdu.Size(nbrLayers,true);
   qdd.Size(nbrLayers,true);

   createIt(quu1,qud1);
   createIt(Equu1,Equd1);
   Equu1.Size(nbrLayers,false);
   Equd1.Size(nbrLayers,false);
   quu1.Size(nbrLayers,true); { les composantes de q peuvent etre complexes }
   qud1.Size(nbrLayers,true);

   Ruu.Size(angles);Rdu.Size(angles);
   Rud.Size(angles);Rdd.Size(angles);
   for i:=0 to nbrPoints-1 do
   begin
      DD.eye(8,8,true);
      M.eye(8,8,true); { initialisation de la récursion }
      q0 := k0*Sin(angles.values[i]);   { on a supprimé *DEG }

      Vaux.setVal(Cplx(mun*BeffMag[nbrlayers-1]*Cos(angles.values[i])*Cos(angles.values[i]),0));

      {calcul des q++,q+-,....}
      Equu.setVal(hb2over2mn);Eqdd.setVal(hb2over2mn);
      Equd.setVal(hb2over2mn);Eqdu.setVal(hb2over2mn);
      Equu.scale(q0*q0);Eqdd.scale(q0*q0);
      Equd.scale(q0*q0);Eqdu.scale(q0*q0);

      Equu.sub(Vu);Eqdd.sub(Vd);
      Equd.sub(Vu);Eqdu.sub(Vd);
      Equu.sub(Vaux);Eqdd.add(Vaux);
      Equd.add(Vaux);Eqdu.sub(Vaux);
      quu.RealToCplx(Equu,Eim);
      quu.scale(1/hb2over2mn);
      quu.sqrt; { q peut-etre complexe, racine }
      qdu.RealToCplx(Eqdu,Eim);qdu.scale(1/hb2over2mn);qdu.sqrt; { q peut-etre complexe, racine }
      qdd.RealToCplx(Eqdd,Eim);qdd.scale(1/hb2over2mn);qdd.sqrt;
      qud.RealToCplx(Equd,Eim);qud.scale(1/hb2over2mn);qud.sqrt;
      { initialisation de variables pour la récursion } { beaucoup de lignes inutiles, voir V11 }
      { + soupcon d'une erreur : le sin est appliqué à la fin du calcul des qquu1 et qqud1 ?!?  }
      { MAIS avec la version upgradée V11, on obtient les mêmes résultats !                     }
       {essai avec initialisation des q0 + et -}
      Equu1.setVal(hb2over2mn*k0*k0-mun*BeffMag[nbrlayers-1]);
      Equd1.setVal(hb2over2mn*k0*k0+mun*BeffMag[nbrlayers-1]);
      {Equu1.scale(q0*q0);
      Equd1.scale(q0*q0);       }
     { Equu1.add(-mun*BeffMag[nbrlayers-1]);
      Equd1.sub(Cplx(-mun*BeffMag[nbrlayers-1],0));
      }quu1.RealToCplx(Equu1,Eim);
      quu1.scale(1/hb2over2mn);
      quu1.sqrt;              { q peut-etre complexe, racine }
      quu1.Scale(Sin(angles.values[i])); { on a supprimé *DEG }
      qud1.RealToCplx(Equd1,Eim);qud1.scale(1/hb2over2mn);qud1.sqrt;
      qud1.Scale(Sin(angles.values[i])); { on a supprimé *DEG }


      qquu1:=quu1.Cvalues[nbrLayers-1]; qqdu1:=qud1.Cvalues[nbrLayers-1];
      qqud1:=quu1.Cvalues[nbrLayers-1]; qqdd1:=qud1.Cvalues[nbrLayers-1];
     {qquu1:=quu.Cvalues[nbrLayers-1]; qqdu1:=qdd.Cvalues[nbrLayers-1];
      qqud1:=quu.Cvalues[nbrLayers-1]; qqdd1:=qdd.Cvalues[nbrLayers-1];  }
      zz:=0;

      for n:=nbrLayers-2 downto 0 do   (* boucle du vide vers le substrat *)
      begin
         qquu0:=quu.Cvalues[n];
         qqdu0:=qdu.Cvalues[n];
         qqud0:=qud.Cvalues[n];
         qqdd0:=qdd.Cvalues[n];

         sig:=sigma.values[n]; { rugosité de l'interface j+1/j }
         zz:=zz+z.values[n+1]; { on fixe le séro à la surface de la couche }
         zz0:=Cplx(zz,0);      { passage en complexe pour faire les multiplications }

         {détermination des coefficients de rugosité}
         auxupuu1DW :=exp(Cmul(Csqr(Csub(qquu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupuu2DW :=exp(Cmul(Csqr(Cadd(qquu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupud1DW :=exp(Cmul(Csqr(Csub(qqdu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupud2DW :=exp(Cmul(Csqr(Cadd(qqdu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupdu1DW :=exp(Cmul(Csqr(Csub(qquu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdu2DW :=exp(Cmul(Csqr(Cadd(qquu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdd1DW :=exp(Cmul(Csqr(Csub(qqdu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdd2DW :=exp(Cmul(Csqr(Cadd(qqdu0,qqdu1)),Cplx(-sig*sig/2,0)));
         { détermination de la matrice de transfert }
         { coefficient p ++}
         ppuu:=Cdiv(Cadd(qquu1,qquu0),qquu1);
         { coefficient p +-}
         ppud:=Cdiv(Cadd(qquu1,qqdu0),qquu1);
         { coefficient p -+}
         ppdu:=Cdiv(Cadd(qqdu1,qquu0),qqdu1);
         { coefficient p --}
         ppdd:=Cdiv(Cadd(qqdu1,qqdu0),qqdu1);
         { coefficient m ++}
         mmuu:=Cdiv(Csub(qquu1,qquu0),qquu1);
         { coefficient m +-}
         mmud:=Cdiv(Csub(qquu1,qqdu0),qquu1);
         { coefficient m -+}
         mmdu:=Cdiv(Csub(qqdu1,qquu0),qqdu1);
         { coefficient m --}
         mmdd:=Cdiv(Csub(qqdu1,qqdu0),qqdu1);
         { paramètres de "propagation"}
         auxupuu1 := Exp(Cmul(C_I,Cmul(Csub(qquu0,qquu1),zz0)));
         auxupuu2 := Exp(Cmul(C_I,Cmul(Cadd(qquu0,qquu1),zz0)));
         auxupud1 := Exp(Cmul(C_I,Cmul(Csub(qqdu0,qquu1),zz0)));
         auxupud2 := Exp(Cmul(C_I,Cmul(Cadd(qqdu0,qquu1),zz0)));
         auxupdu1 := Exp(Cmul(C_I,Cmul(Csub(qquu0,qqdu1),zz0)));
         auxupdu2 := Exp(Cmul(C_I,Cmul(Cadd(qquu0,qqdu1),zz0)));
         auxupdd1 := Exp(Cmul(C_I,Cmul(Csub(qqdu0,qqdu1),zz0)));
         auxupdd2 := Exp(Cmul(C_I,Cmul(Cadd(qqdu0,qqdu1),zz0)));
         {termes "trigonomètriques"}
         tt1:= Cplx(BeffTheta[n+1]/2,0);
         tt2:= Cplx(BeffTheta[n]/2,0);
         phi1:= Cplx(BeffPhi[n+1],0);
         phi2:= Cplx(BeffPhi[n],0);
         trigup1:= Cadd(Cmul(cos(tt2),cos(tt1)),Cmul(Exp(Cmul(C_I,Csub(phi2,phi1))),Cmul(sin(tt1),sin(tt2))));
         trigup4:= Cadd(Cmul(cos(tt2),cos(tt1)),Cmul(Exp(Cmul(C_I,Csub(phi1,phi2))),Cmul(sin(tt1),sin(tt2))));
         trigup2:= Csub(Cmul(Exp(Cmul(Neg(C_I),phi1)),Cmul(cos(tt2),sin(tt1))),Cmul(Exp(Cmul(Neg(C_I),phi2)),Cmul(cos(tt1),sin(tt2))));
         trigup3:=Csub(Cmul(Exp(Cmul(C_I,phi2)),Cmul(cos(tt1),sin(tt2))),Cmul(Exp(Cmul(C_I,phi1)),Cmul(cos(tt2),sin(tt1))));
         {bloc matrice ++}
         DD.Cvalues[0,0]:=Cmul(Cmul(Cmul(ppuu,auxupuu1),auxupuu1DW),trigup1);
         DD.Cvalues[1,1]:=Cmul(Cmul(Cmul(ppuu,CDiv(C_ONE,auxupuu1)),auxupuu1DW),trigup1);
         DD.Cvalues[1,0]:=Cmul(Cmul(Cmul(mmuu,auxupuu2),auxupuu2DW),trigup1);
         DD.Cvalues[0,1]:=Cmul(Cmul(Cmul(mmuu,CDiv(C_ONE,auxupuu2)),auxupuu2DW),trigup1);
         {bloc matrice +-}
         DD.Cvalues[0,2]:=Cmul(Cmul(Cmul(ppud,auxupud1),auxupud1DW),trigup2);
         DD.Cvalues[1,3]:=Cmul(Cmul(Cmul(ppud,CDiv(C_ONE,auxupud1)),auxupud1DW),trigup2);
         DD.Cvalues[1,2]:=Cmul(Cmul(Cmul(mmud,auxupud2),auxupud2DW),trigup2);
         DD.Cvalues[0,3]:=Cmul(Cmul(Cmul(mmud,CDiv(C_ONE,auxupud2)),auxupud2DW),trigup2);
         {bloc matrice -+}
         DD.Cvalues[2,0]:=Cmul(Cmul(Cmul(ppdu,auxupdu1),auxupdu1DW),trigup3);
         DD.Cvalues[3,1]:=Cmul(Cmul(Cmul(ppdu,CDiv(C_ONE,auxupdu1)),auxupdu1DW),trigup3);
         DD.Cvalues[3,0]:=Cmul(Cmul(Cmul(mmdu,auxupdu2),auxupdu2DW),trigup3);
         DD.Cvalues[2,1]:=Cmul(Cmul(Cmul(mmdu,CDiv(C_ONE,auxupdu2)),auxupdu2DW),trigup3);
         {bloc matrice --}
         DD.Cvalues[2,2]:=Cmul(Cmul(Cmul(ppdd,auxupdd1),auxupdd1DW),trigup4);
         DD.Cvalues[3,3]:=Cmul(Cmul(Cmul(ppdd,CDiv(C_ONE,auxupdd1)),auxupdd1DW),trigup4);
         DD.Cvalues[3,2]:=Cmul(Cmul(Cmul(mmdd,auxupdd2),auxupdd2DW),trigup4);
         DD.Cvalues[2,3]:=Cmul(Cmul(Cmul(mmdd,CDiv(C_ONE,auxupdd2)),auxupdd2DW),trigup4);
    {            on remplace pp par rr et mm par nn.....  }
         auxdouu1DW :=Exp(Cmul(Csqr(Csub(qqud0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdouu2DW :=Exp(Cmul(Csqr(Cadd(qqud0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdoud1DW :=Exp(Cmul(Csqr(Csub(qqdd0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdoud2DW :=Exp(Cmul(Csqr(Cadd(qqdd0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdodu1DW :=Exp(Cmul(Csqr(Csub(qqud0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodu2DW :=Exp(Cmul(Csqr(Cadd(qqud0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodd1DW :=Exp(Cmul(Csqr(Csub(qqdd0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodd2DW :=Exp(Cmul(Csqr(Cadd(qqdd0,qqdd1)),Cplx(-sig*sig/2,0)));

         { coefficient p ++}
         rruu:=Cdiv(Cadd(qqud1,qqud0),qqud1);rruu:=CMul(Cplx(0.5,0),rruu);
         { coefficient p +-}
         rrud:=Cdiv(Cadd(qqud1,qqdd0),qqud1);rrud:=CMul(Cplx(0.5,0),rrud);
         { coefficient p -+}
         rrdu:=Cdiv(Cadd(qqdd1,qqud0),qqdd1);rrdu:=CMul(Cplx(0.5,0),rrdu);
         { coefficient p --}
         rrdd:=Cdiv(Cadd(qqdd1,qqdd0),qqdd1);rrdd:=CMul(Cplx(0.5,0),rrdd);
         { coefficient m ++}
         nnuu:=Cdiv(Csub(qqud1,qqud0),qqud1);nnuu:=CMul(Cplx(0.5,0),nnuu);
         { coefficient m +-}
         nnud:=Cdiv(Csub(qqud1,qqdd0),qqud1);nnud:=CMul(Cplx(0.5,0),nnud);
         { coefficient m -+}
         nndu:=Cdiv(Csub(qqdd1,qqud0),qqdd1);nndu:=CMul(Cplx(0.5,0),nndu);
         { coefficient m --}
         nndd:=Cdiv(Csub(qqdd1,qqdd0),qqdd1);nndd:=CMul(Cplx(0.5,0),nndd);
         { initialisation of the D matrix}
         auxdouu1 := Exp(Cmul(C_I,Cmul(Csub(qqud0,qqud1),zz0)));
         auxdouu2 := Exp(Cmul(C_I,Cmul(Cadd(qqud0,qqud1),zz0)));
         auxdoud1 := Exp(Cmul(C_I,Cmul(Csub(qqdd0,qqud1),zz0)));
         auxdoud2 := Exp(Cmul(C_I,Cmul(Cadd(qqdd0,qqud1),zz0)));
         auxdodu1 := Exp(Cmul(C_I,Cmul(Csub(qqud0,qqdd1),zz0)));
         auxdodu2 := Exp(Cmul(C_I,Cmul(Cadd(qqud0,qqdd1),zz0)));
         auxdodd1 := Exp(Cmul(C_I,Cmul(Csub(qqdd0,qqdd1),zz0)));
         auxdodd2 := Exp(Cmul(C_I,Cmul(Cadd(qqdd0,qqdd1),zz0)));

         {bloc matrice ++}
         DD.Cvalues[4,4]:=Cmul(Cmul(Cmul(rrdd,auxdodd1),auxdodd1DW),Trigup4);
         DD.Cvalues[5,5]:=Cmul(Cmul(Cmul(rrdd,CDiv(C_ONE,auxdodd1)),auxdodd1DW),Trigup4);
         DD.Cvalues[5,4]:=Cmul(Cmul(Cmul(nndd,auxdodd2),auxdodd2DW),Trigup4);
         DD.Cvalues[4,5]:=Cmul(Cmul(Cmul(nndd,CDiv(C_ONE,auxdodd2)),auxdodd2DW),Trigup4);
         {bloc matrice +-}
         DD.Cvalues[4,6]:=Cmul(Cmul(Cmul(rrdu,auxdodu1),auxdodu1DW),Trigup3);
         DD.Cvalues[5,7]:=Cmul(Cmul(Cmul(rrdu,CDiv(C_ONE,auxdodu1)),auxdodu1DW),Trigup3);
         DD.Cvalues[5,6]:=Cmul(Cmul(Cmul(nndu,auxdodu2),auxdodu2DW),Trigup3);
         DD.Cvalues[4,7]:=Cmul(Cmul(Cmul(nndu,CDiv(C_ONE,auxdodu2)),auxdodu2DW),Trigup3);
         {bloc matrice -+}
         DD.Cvalues[6,4]:=Cmul(Cmul(Cmul(rrud,auxdoud1),auxdoud1DW),Trigup2);
         DD.Cvalues[7,5]:=Cmul(Cmul(Cmul(rrud,CDiv(C_ONE,auxdoud1)),auxdoud1DW),Trigup2);
         DD.Cvalues[7,4]:=Cmul(Cmul(Cmul(nnud,auxdoud2),auxdoud2DW),Trigup2);
         DD.Cvalues[6,5]:=Cmul(Cmul(Cmul(nnud,CDiv(C_ONE,auxdoud2)),auxdoud2DW),Trigup2);
         {bloc matrice --}
         DD.Cvalues[6,6]:=Cmul(Cmul(Cmul(rruu,auxdouu1),auxdouu1DW),Trigup1);
         DD.Cvalues[7,7]:=Cmul(Cmul(Cmul(rruu,CDiv(C_ONE,auxdouu1)),auxdouu1DW),Trigup1);
         DD.Cvalues[7,6]:=Cmul(Cmul(Cmul(nnuu,auxdouu2),auxdouu2DW),Trigup1);
         DD.Cvalues[6,7]:=Cmul(Cmul(Cmul(nnuu,CDiv(C_ONE,auxdouu2)),auxdouu2DW),Trigup1);
        { DD.Cvalues[0,4]:=0;
         DD.Cvalues[0,5]:=0;  }
{         DD.scale(0.5); } { ne marche pas correctement }
         Maux.Mul(M,DD);
         M.copy(Maux);
         qquu1:=qquu0;
         qqdu1:=qqdu0;
         qqud1:=qqud0;
         qqdd1:=qqdd0;

      end;
      aux11:=M.Cvalues[0,0];aux13:=M.Cvalues[0,2];
      aux21:=M.Cvalues[1,0];aux23:=M.Cvalues[1,2];
      aux31:=M.Cvalues[2,0];aux33:=M.Cvalues[2,2];
      aux41:=M.Cvalues[3,0];aux43:=M.Cvalues[3,2];
      aux55:=M.Cvalues[4,4];aux57:=M.Cvalues[4,6];
      aux65:=M.Cvalues[5,4];aux67:=M.Cvalues[5,6];
      aux75:=M.Cvalues[6,4];aux77:=M.Cvalues[6,6];
      aux85:=M.Cvalues[7,4];aux87:=M.Cvalues[7,6];
      auxupD:=Csub(Cmul(aux11,aux33),Cmul(aux13,aux31));
      auxdoD:=Csub(Cmul(aux55,aux77),Cmul(aux57,aux75));
      Ruu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux21,aux33),Cmul(aux23,aux31)),auxupD)));
      Rdu.values[i]:=Sqr(Cabs(CDiv(Csub(Cmul(aux41,aux33),Cmul(aux43,aux31)),auxupD)));
      Rdd.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux65,aux77),Cmul(aux67,aux75)),auxdoD)));
      Rud.values[i]:=Sqr(Cabs(CDiv(Csub(Cmul(aux85,aux77),Cmul(aux87,aux75)),auxdoD)));

   end;
   freeit(DD,M,Maux);
   freeit(quu,qdu,qud,qdd);
   freeIt(Equu,Eqdd,Eim);
   freeIt(Equd,Eqdu);
end;

{version 0.3    P.Humbert                                    }
{ procédure de calcul de la réflectivité (neut. polarisés)    }
{ on ne tient pas compte des k parallleles                    }
{prise en compte de la rugosite                               }
{ INPUT: longueur d'onde en nm                                }
{        angles, vecteur des angles d'incidence (en  °)       }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vetceur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ REM: aucun paramètre implicite!!                            }
{      dans cette version, les matrices auxilaires D,P,M et Maux       }
{      sont créées à l'initialisation (gain de 10% seulement)          }
{ * V10 to V11 : simplification du calcul des variables qquu1 et qqdd1                         }
{                les noms de ces variables sont vaseux qquu1 et qqdd1 seraient mieux!          }
procedure CalculReflecSF_HF_V11(var lambda:double;var angles:Tvec;var z,Vu,Vd,sigma,BeffMag,BeffTheta,BeffPhi:Tvec;var Ruu,Rdu,Rud,Rdd:Tvec);
var i,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    k0,zz:double;          { norme du vecteur d'onde incident  }
    E0,q0:double;
    qquu0,qqdu0,qqud0,qqdd0,qquu1,qqdu1,qqud1,qqdd1,zz0:TCplx;
    quu,qdu,qud,qdd:Tvec;            { vecteur d'onde perpendiculaire dans chaque couche }
    Equu,Eqdd,Equd,Eqdu,Eim:Tvec;    { Energie cinétique perp. dans chaque couche }
    Vaux:Tvec;
    ppuu,ppud,ppdu,ppdd,mmuu,mmud,mmdu,mmdd:Tcplx; {coefficient p et m (Gibaud p.101)}
    rruu,rrud,rrdu,rrdd,nnuu,nnud,nndu,nndd:Tcplx;
    auxupuu1,auxupuu2,auxupud1,auxupud2,auxupdu1,auxupdu2,auxupdd1,auxupdd2:Tcplx; {coefficient de propagation}
    auxdouu1,auxdouu2,auxdoud1,auxdoud2,auxdodu1,auxdodu2,auxdodd1,auxdodd2:Tcplx;
    sig:double;
    auxupuu1DW,auxupuu2DW,auxupud1DW,auxupud2DW:Tcplx; {coefficient de rugosité}
    auxupdu1DW,auxupdu2DW,auxupdd1DW,auxupdd2DW:Tcplx;
    auxdouu1DW,auxdouu2DW,auxdoud1Dw,auxdoud2Dw:Tcplx;
    auxdodu1DW,auxdodu2DW,auxdodd1DW,auxdodd2DW:Tcplx;
    DD,M,Maux:TMtx;
    auxupD,aux11,aux13,aux21,aux23,aux31,aux33,aux41,aux43:TCplx;
    auxdoD,aux55,aux57,aux65,aux67,aux75,aux77,aux85,aux87:TCplx;
    tt1,tt2,phi1,phi2:TCplx;
    trigup1,trigup2,trigup3,trigup4:TCplx;
    Equu1,Equd1:Tvec;
    quu1,qud1:Tvec;

begin
   createIt(Equu,Eqdd,Eim);
   createIt(Equd,Eqdu);
   createIt(Vaux);
   createIt(DD);DD.size(8,8,true);
   createIt(M);M.size(8,8,true);
   createIt(Maux);Maux.size(8,8,true);
   nbrLayers := z.length;
   nbrPoints := angles.length;
   k0:=2*Pi/lambda;
   Vaux.Size(nbrLayers,false);
   Equu.Size(nbrLayers,false);Equd.Size(nbrLayers,false);
   Eqdu.Size(nbrLayers,false);Eqdd.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0);
   createIt(quu,qdu,qud,qdd);
   quu.Size(nbrLayers,true);
   qud.Size(nbrLayers,true);
   qdu.Size(nbrLayers,true);
   qdd.Size(nbrLayers,true);

   createIt(quu1,qud1);
   createIt(Equu1,Equd1);
   Equu1.Size(nbrLayers,false);
   Equd1.Size(nbrLayers,false);
   quu1.Size(nbrLayers,true); { les composantes de q peuvent etre complexes }
   qud1.Size(nbrLayers,true);

   Ruu.Size(angles);Rdu.Size(angles);
   Rud.Size(angles);Rdd.Size(angles);
   for i:=0 to nbrPoints-1 do
   begin
      DD.eye(8,8,true);
      M.eye(8,8,true); { initialisation de la récursion }
      q0 := k0*Sin(angles.values[i]);   { on a supprimé *DEG }

      Vaux.setVal(Cplx(mun*BeffMag[nbrlayers-1]*Cos(angles.values[i])*Cos(angles.values[i]),0));
      { a quoi correspond ce potentiel ?? }

      {calcul des q++,q+-,....}
      Equu.setVal(hb2over2mn);Eqdd.setVal(hb2over2mn);
      Equd.setVal(hb2over2mn);Eqdu.setVal(hb2over2mn);
      Equu.scale(q0*q0);Eqdd.scale(q0*q0);
      Equd.scale(q0*q0);Eqdu.scale(q0*q0);

      Equu.sub(Vu);Eqdd.sub(Vd);
      Equd.sub(Vu);Eqdu.sub(Vd);
      Equu.sub(Vaux);Eqdd.add(Vaux);
      Equd.add(Vaux);Eqdu.sub(Vaux);
      quu.RealToCplx(Equu,Eim);
      quu.scale(1/hb2over2mn);
      quu.sqrt; { q peut-etre complexe, racine }
      qdu.RealToCplx(Eqdu,Eim);qdu.scale(1/hb2over2mn);qdu.sqrt; { q peut-etre complexe, racine }
      qdd.RealToCplx(Eqdd,Eim);qdd.scale(1/hb2over2mn);qdd.sqrt;
      qud.RealToCplx(Equd,Eim);qud.scale(1/hb2over2mn);qud.sqrt;

      { initialisation de variables pour la récursion }
      qquu1:=quu.Cvalues[nbrLayers-1]; qqdu1:=qdd.Cvalues[nbrLayers-1]; { ces lignes avaient été mises en commentaires }
      qqud1:=quu.Cvalues[nbrLayers-1]; qqdd1:=qdd.Cvalues[nbrLayers-1];
      zz:=0;

      for n:=nbrLayers-2 downto 0 do   (* boucle du vide vers le substrat *)
      begin
         qquu0:=quu.Cvalues[n];
         qqdu0:=qdu.Cvalues[n];
         qqud0:=qud.Cvalues[n];
         qqdd0:=qdd.Cvalues[n];

         sig:=sigma.values[n]; { rugosité de l'interface j+1/j }
         zz:=zz+z.values[n+1]; { on fixe le séro à la surface de la couche }
         zz0:=Cplx(zz,0);      { passage en complexe pour faire les multiplications }

         {détermination des coefficients de rugosité}
         auxupuu1DW :=exp(Cmul(Csqr(Csub(qquu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupuu2DW :=exp(Cmul(Csqr(Cadd(qquu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupud1DW :=exp(Cmul(Csqr(Csub(qqdu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupud2DW :=exp(Cmul(Csqr(Cadd(qqdu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupdu1DW :=exp(Cmul(Csqr(Csub(qquu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdu2DW :=exp(Cmul(Csqr(Cadd(qquu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdd1DW :=exp(Cmul(Csqr(Csub(qqdu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdd2DW :=exp(Cmul(Csqr(Cadd(qqdu0,qqdu1)),Cplx(-sig*sig/2,0)));
         { détermination de la matrice de transfert }
         { coefficient p ++}
         ppuu:=Cdiv(Cadd(qquu1,qquu0),qquu1);
         { coefficient p +-}
         ppud:=Cdiv(Cadd(qquu1,qqdu0),qquu1);
         { coefficient p -+}
         ppdu:=Cdiv(Cadd(qqdu1,qquu0),qqdu1);
         { coefficient p --}
         ppdd:=Cdiv(Cadd(qqdu1,qqdu0),qqdu1);
         { coefficient m ++}
         mmuu:=Cdiv(Csub(qquu1,qquu0),qquu1);
         { coefficient m +-}
         mmud:=Cdiv(Csub(qquu1,qqdu0),qquu1);
         { coefficient m -+}
         mmdu:=Cdiv(Csub(qqdu1,qquu0),qqdu1);
         { coefficient m --}
         mmdd:=Cdiv(Csub(qqdu1,qqdu0),qqdu1);
         { paramètres de "propagation"}
         auxupuu1 := Exp(Cmul(C_I,Cmul(Csub(qquu0,qquu1),zz0)));
         auxupuu2 := Exp(Cmul(C_I,Cmul(Cadd(qquu0,qquu1),zz0)));
         auxupud1 := Exp(Cmul(C_I,Cmul(Csub(qqdu0,qquu1),zz0)));
         auxupud2 := Exp(Cmul(C_I,Cmul(Cadd(qqdu0,qquu1),zz0)));
         auxupdu1 := Exp(Cmul(C_I,Cmul(Csub(qquu0,qqdu1),zz0)));
         auxupdu2 := Exp(Cmul(C_I,Cmul(Cadd(qquu0,qqdu1),zz0)));
         auxupdd1 := Exp(Cmul(C_I,Cmul(Csub(qqdu0,qqdu1),zz0)));
         auxupdd2 := Exp(Cmul(C_I,Cmul(Cadd(qqdu0,qqdu1),zz0)));
         {termes "trigonomètriques"}
         tt1:= Cplx(BeffTheta[n+1]/2,0);
         tt2:= Cplx(BeffTheta[n]/2,0);
         phi1:= Cplx(BeffPhi[n+1],0);
         phi2:= Cplx(BeffPhi[n],0);
         trigup1:= Cadd(Cmul(cos(tt2),cos(tt1)),Cmul(Exp(Cmul(C_I,Csub(phi2,phi1))),Cmul(sin(tt1),sin(tt2))));
         trigup4:= Cadd(Cmul(cos(tt2),cos(tt1)),Cmul(Exp(Cmul(C_I,Csub(phi1,phi2))),Cmul(sin(tt1),sin(tt2))));
         trigup2:= Csub(Cmul(Exp(Cmul(Neg(C_I),phi1)),Cmul(cos(tt2),sin(tt1))),Cmul(Exp(Cmul(Neg(C_I),phi2)),Cmul(cos(tt1),sin(tt2))));
         trigup3:=Csub(Cmul(Exp(Cmul(C_I,phi2)),Cmul(cos(tt1),sin(tt2))),Cmul(Exp(Cmul(C_I,phi1)),Cmul(cos(tt2),sin(tt1))));
         {bloc matrice ++}
         DD.Cvalues[0,0]:=Cmul(Cmul(Cmul(ppuu,auxupuu1),auxupuu1DW),trigup1);
         DD.Cvalues[1,1]:=Cmul(Cmul(Cmul(ppuu,CDiv(C_ONE,auxupuu1)),auxupuu1DW),trigup1);
         DD.Cvalues[1,0]:=Cmul(Cmul(Cmul(mmuu,auxupuu2),auxupuu2DW),trigup1);
         DD.Cvalues[0,1]:=Cmul(Cmul(Cmul(mmuu,CDiv(C_ONE,auxupuu2)),auxupuu2DW),trigup1);
         {bloc matrice +-}
         DD.Cvalues[0,2]:=Cmul(Cmul(Cmul(ppud,auxupud1),auxupud1DW),trigup2);
         DD.Cvalues[1,3]:=Cmul(Cmul(Cmul(ppud,CDiv(C_ONE,auxupud1)),auxupud1DW),trigup2);
         DD.Cvalues[1,2]:=Cmul(Cmul(Cmul(mmud,auxupud2),auxupud2DW),trigup2);
         DD.Cvalues[0,3]:=Cmul(Cmul(Cmul(mmud,CDiv(C_ONE,auxupud2)),auxupud2DW),trigup2);
         {bloc matrice -+}
         DD.Cvalues[2,0]:=Cmul(Cmul(Cmul(ppdu,auxupdu1),auxupdu1DW),trigup3);
         DD.Cvalues[3,1]:=Cmul(Cmul(Cmul(ppdu,CDiv(C_ONE,auxupdu1)),auxupdu1DW),trigup3);
         DD.Cvalues[3,0]:=Cmul(Cmul(Cmul(mmdu,auxupdu2),auxupdu2DW),trigup3);
         DD.Cvalues[2,1]:=Cmul(Cmul(Cmul(mmdu,CDiv(C_ONE,auxupdu2)),auxupdu2DW),trigup3);
         {bloc matrice --}
         DD.Cvalues[2,2]:=Cmul(Cmul(Cmul(ppdd,auxupdd1),auxupdd1DW),trigup4);
         DD.Cvalues[3,3]:=Cmul(Cmul(Cmul(ppdd,CDiv(C_ONE,auxupdd1)),auxupdd1DW),trigup4);
         DD.Cvalues[3,2]:=Cmul(Cmul(Cmul(mmdd,auxupdd2),auxupdd2DW),trigup4);
         DD.Cvalues[2,3]:=Cmul(Cmul(Cmul(mmdd,CDiv(C_ONE,auxupdd2)),auxupdd2DW),trigup4);
    {            on remplace pp par rr et mm par nn.....  }
         auxdouu1DW :=Exp(Cmul(Csqr(Csub(qqud0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdouu2DW :=Exp(Cmul(Csqr(Cadd(qqud0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdoud1DW :=Exp(Cmul(Csqr(Csub(qqdd0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdoud2DW :=Exp(Cmul(Csqr(Cadd(qqdd0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdodu1DW :=Exp(Cmul(Csqr(Csub(qqud0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodu2DW :=Exp(Cmul(Csqr(Cadd(qqud0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodd1DW :=Exp(Cmul(Csqr(Csub(qqdd0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodd2DW :=Exp(Cmul(Csqr(Cadd(qqdd0,qqdd1)),Cplx(-sig*sig/2,0)));

         { coefficient p ++}
         rruu:=Cdiv(Cadd(qqud1,qqud0),qqud1);rruu:=CMul(Cplx(0.5,0),rruu);
         { coefficient p +-}
         rrud:=Cdiv(Cadd(qqud1,qqdd0),qqud1);rrud:=CMul(Cplx(0.5,0),rrud);
         { coefficient p -+}
         rrdu:=Cdiv(Cadd(qqdd1,qqud0),qqdd1);rrdu:=CMul(Cplx(0.5,0),rrdu);
         { coefficient p --}
         rrdd:=Cdiv(Cadd(qqdd1,qqdd0),qqdd1);rrdd:=CMul(Cplx(0.5,0),rrdd);
         { coefficient m ++}
         nnuu:=Cdiv(Csub(qqud1,qqud0),qqud1);nnuu:=CMul(Cplx(0.5,0),nnuu);
         { coefficient m +-}
         nnud:=Cdiv(Csub(qqud1,qqdd0),qqud1);nnud:=CMul(Cplx(0.5,0),nnud);
         { coefficient m -+}
         nndu:=Cdiv(Csub(qqdd1,qqud0),qqdd1);nndu:=CMul(Cplx(0.5,0),nndu);
         { coefficient m --}
         nndd:=Cdiv(Csub(qqdd1,qqdd0),qqdd1);nndd:=CMul(Cplx(0.5,0),nndd);
         { initialisation of the D matrix}
         auxdouu1 := Exp(Cmul(C_I,Cmul(Csub(qqud0,qqud1),zz0)));
         auxdouu2 := Exp(Cmul(C_I,Cmul(Cadd(qqud0,qqud1),zz0)));
         auxdoud1 := Exp(Cmul(C_I,Cmul(Csub(qqdd0,qqud1),zz0)));
         auxdoud2 := Exp(Cmul(C_I,Cmul(Cadd(qqdd0,qqud1),zz0)));
         auxdodu1 := Exp(Cmul(C_I,Cmul(Csub(qqud0,qqdd1),zz0)));
         auxdodu2 := Exp(Cmul(C_I,Cmul(Cadd(qqud0,qqdd1),zz0)));
         auxdodd1 := Exp(Cmul(C_I,Cmul(Csub(qqdd0,qqdd1),zz0)));
         auxdodd2 := Exp(Cmul(C_I,Cmul(Cadd(qqdd0,qqdd1),zz0)));

         {bloc matrice ++}
         DD.Cvalues[4,4]:=Cmul(Cmul(Cmul(rrdd,auxdodd1),auxdodd1DW),Trigup4);
         DD.Cvalues[5,5]:=Cmul(Cmul(Cmul(rrdd,CDiv(C_ONE,auxdodd1)),auxdodd1DW),Trigup4);
         DD.Cvalues[5,4]:=Cmul(Cmul(Cmul(nndd,auxdodd2),auxdodd2DW),Trigup4);
         DD.Cvalues[4,5]:=Cmul(Cmul(Cmul(nndd,CDiv(C_ONE,auxdodd2)),auxdodd2DW),Trigup4);
         {bloc matrice +-}
         DD.Cvalues[4,6]:=Cmul(Cmul(Cmul(rrdu,auxdodu1),auxdodu1DW),Trigup3);
         DD.Cvalues[5,7]:=Cmul(Cmul(Cmul(rrdu,CDiv(C_ONE,auxdodu1)),auxdodu1DW),Trigup3);
         DD.Cvalues[5,6]:=Cmul(Cmul(Cmul(nndu,auxdodu2),auxdodu2DW),Trigup3);
         DD.Cvalues[4,7]:=Cmul(Cmul(Cmul(nndu,CDiv(C_ONE,auxdodu2)),auxdodu2DW),Trigup3);
         {bloc matrice -+}
         DD.Cvalues[6,4]:=Cmul(Cmul(Cmul(rrud,auxdoud1),auxdoud1DW),Trigup2);
         DD.Cvalues[7,5]:=Cmul(Cmul(Cmul(rrud,CDiv(C_ONE,auxdoud1)),auxdoud1DW),Trigup2);
         DD.Cvalues[7,4]:=Cmul(Cmul(Cmul(nnud,auxdoud2),auxdoud2DW),Trigup2);
         DD.Cvalues[6,5]:=Cmul(Cmul(Cmul(nnud,CDiv(C_ONE,auxdoud2)),auxdoud2DW),Trigup2);
         {bloc matrice --}
         DD.Cvalues[6,6]:=Cmul(Cmul(Cmul(rruu,auxdouu1),auxdouu1DW),Trigup1);
         DD.Cvalues[7,7]:=Cmul(Cmul(Cmul(rruu,CDiv(C_ONE,auxdouu1)),auxdouu1DW),Trigup1);
         DD.Cvalues[7,6]:=Cmul(Cmul(Cmul(nnuu,auxdouu2),auxdouu2DW),Trigup1);
         DD.Cvalues[6,7]:=Cmul(Cmul(Cmul(nnuu,CDiv(C_ONE,auxdouu2)),auxdouu2DW),Trigup1);
        { DD.Cvalues[0,4]:=0;
         DD.Cvalues[0,5]:=0;  }
{         DD.scale(0.5); } { ne marche pas correctement }
         Maux.Mul(M,DD);
         M.copy(Maux);
         qquu1:=qquu0;
         qqdu1:=qqdu0;
         qqud1:=qqud0;
         qqdd1:=qqdd0;

      end;
      aux11:=M.Cvalues[0,0];aux13:=M.Cvalues[0,2];
      aux21:=M.Cvalues[1,0];aux23:=M.Cvalues[1,2];
      aux31:=M.Cvalues[2,0];aux33:=M.Cvalues[2,2];
      aux41:=M.Cvalues[3,0];aux43:=M.Cvalues[3,2];
      aux55:=M.Cvalues[4,4];aux57:=M.Cvalues[4,6];
      aux65:=M.Cvalues[5,4];aux67:=M.Cvalues[5,6];
      aux75:=M.Cvalues[6,4];aux77:=M.Cvalues[6,6];
      aux85:=M.Cvalues[7,4];aux87:=M.Cvalues[7,6];
      auxupD:=Csub(Cmul(aux11,aux33),Cmul(aux13,aux31));
      auxdoD:=Csub(Cmul(aux55,aux77),Cmul(aux57,aux75));
      Ruu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux21,aux33),Cmul(aux23,aux31)),auxupD)));
      Rdu.values[i]:=Sqr(Cabs(CDiv(Csub(Cmul(aux41,aux33),Cmul(aux43,aux31)),auxupD)));
      Rdd.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux65,aux77),Cmul(aux67,aux75)),auxdoD)));
      Rud.values[i]:=Sqr(Cabs(CDiv(Csub(Cmul(aux85,aux77),Cmul(aux87,aux75)),auxdoD)));

   end;
   freeit(DD,M,Maux);
   freeit(quu,qdu,qud,qdd);
   freeIt(Equu,Eqdd,Eim);
   freeIt(Equd,Eqdu);
end;

{ version 1.2    P.Humbert                                    }
{ procédure de calcul de la réflectivité (neut. polarisés)    }
{ on ne tient pas compte des k parallleles                    }
{prise en compte de la rugosite                               }
{ INPUT: longueur d'onde en nm                                }
{        angles, vecteur des angles d'incidence (en  °)       }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vetceur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ REM: aucun paramètre implicite!!                            }
{      dans cette version, les matrices auxilaires D,P,M et Maux       }
{      sont créées à l'initialisation (gain de 10% seulement)          }
{ * V10 to V11 : simplification du calcul des variables qquu1 et qqdd1                         }
{                les noms de ces variables sont vaseux qquu1 et qqdd1 seraient mieux!          }
{ * V11 to V12 : lambda n'est plus passé comme paramètre        }
{                = prmemière version utilisable avec la V1.4    }
procedure CalculReflecSF_HF_V12(var angles:Tvec;Qz:Tvec;var z,Vu,Vd,sigma,BeffMag,BeffTheta,BeffPhi:Tvec;var Ruu,Rdu,Rud,Rdd:Tvec);
var i,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    zz:double;
    E0,q0:double;
    qquu0,qqdu0,qqud0,qqdd0,qquu1,qqdu1,qqud1,qqdd1,zz0:TCplx;
    quu,qdu,qud,qdd:Tvec;            { vecteur d'onde perpendiculaire dans chaque couche }
    Equu,Eqdd,Equd,Eqdu,Eim:Tvec;    { Energie cinétique perp. dans chaque couche }
    Vaux:Tvec;
    ppuu,ppud,ppdu,ppdd,mmuu,mmud,mmdu,mmdd:Tcplx; {coefficient p et m (Gibaud p.101)}
    rruu,rrud,rrdu,rrdd,nnuu,nnud,nndu,nndd:Tcplx;
    auxupuu1,auxupuu2,auxupud1,auxupud2,auxupdu1,auxupdu2,auxupdd1,auxupdd2:Tcplx; {coefficient de propagation}
    auxdouu1,auxdouu2,auxdoud1,auxdoud2,auxdodu1,auxdodu2,auxdodd1,auxdodd2:Tcplx;
    sig:double;
    auxupuu1DW,auxupuu2DW,auxupud1DW,auxupud2DW:Tcplx; {coefficient de rugosité}
    auxupdu1DW,auxupdu2DW,auxupdd1DW,auxupdd2DW:Tcplx;
    auxdouu1DW,auxdouu2DW,auxdoud1Dw,auxdoud2Dw:Tcplx;
    auxdodu1DW,auxdodu2DW,auxdodd1DW,auxdodd2DW:Tcplx;
    DD,M,Maux:TMtx;
    auxupD,aux11,aux13,aux21,aux23,aux31,aux33,aux41,aux43:TCplx;
    auxdoD,aux55,aux57,aux65,aux67,aux75,aux77,aux85,aux87:TCplx;
    tt1,tt2,phi1,phi2:TCplx;
    trigup1,trigup2,trigup3,trigup4:TCplx;

begin
   createIt(Equu,Eqdd,Eim);
   createIt(Equd,Eqdu);
   createIt(Vaux);
   createIt(DD);DD.size(8,8,true);
   createIt(M);M.size(8,8,true);
   createIt(Maux);Maux.size(8,8,true);
   nbrLayers := z.length;
   nbrPoints := Qz.length;
   Vaux.Size(nbrLayers,false);
   Equu.Size(nbrLayers,false);Equd.Size(nbrLayers,false);
   Eqdu.Size(nbrLayers,false);Eqdd.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0);
   createIt(quu,qdu,qud,qdd);
   quu.Size(nbrLayers,true);
   qud.Size(nbrLayers,true);
   qdu.Size(nbrLayers,true);
   qdd.Size(nbrLayers,true);

   Ruu.Size(Qz);Rdu.Size(Qz);Rud.Size(Qz);Rdd.Size(Qz);
   for i:=0 to nbrPoints-1 do
   begin
      DD.eye(8,8,true);
      M.eye(8,8,true); { initialisation de la récursion }
      q0 := Qz.values[i]/2;   { on a supprimé *DEG }

      Vaux.setVal(Cplx(mun*BeffMag[nbrlayers-1]*Cos(angles.values[i])*Cos(angles.values[i]),0));
      { a quoi correspond ce potentiel ?? }

      {calcul des q++,q+-,....}
      Equu.setVal(hb2over2mn);Eqdd.setVal(hb2over2mn);
      Equd.setVal(hb2over2mn);Eqdu.setVal(hb2over2mn);
      Equu.scale(q0*q0);Eqdd.scale(q0*q0);
      Equd.scale(q0*q0);Eqdu.scale(q0*q0);

      Equu.sub(Vu);Eqdd.sub(Vd);
      Equd.sub(Vu);Eqdu.sub(Vd);
      Equu.sub(Vaux);Eqdd.add(Vaux);
      Equd.add(Vaux);Eqdu.sub(Vaux);
      quu.RealToCplx(Equu,Eim);
      quu.scale(1/hb2over2mn);
      quu.sqrt; { q peut-etre complexe, racine }
      qdu.RealToCplx(Eqdu,Eim);qdu.scale(1/hb2over2mn);qdu.sqrt; { q peut-etre complexe, racine }
      qdd.RealToCplx(Eqdd,Eim);qdd.scale(1/hb2over2mn);qdd.sqrt;
      qud.RealToCplx(Equd,Eim);qud.scale(1/hb2over2mn);qud.sqrt;

      { initialisation de variables pour la récursion }
      qquu1:=quu.Cvalues[nbrLayers-1]; qqdu1:=qdd.Cvalues[nbrLayers-1]; { ces lignes avaient été mises en commentaires }
      qqud1:=quu.Cvalues[nbrLayers-1]; qqdd1:=qdd.Cvalues[nbrLayers-1];
      zz:=0;

      for n:=nbrLayers-2 downto 0 do   (* boucle du vide vers le substrat *)
      begin
         qquu0:=quu.Cvalues[n];
         qqdu0:=qdu.Cvalues[n];
         qqud0:=qud.Cvalues[n];
         qqdd0:=qdd.Cvalues[n];

         sig:=sigma.values[n]; { rugosité de l'interface j+1/j }
         zz:=zz+z.values[n+1]; { on fixe le séro à la surface de la couche }
         zz0:=Cplx(zz,0);      { passage en complexe pour faire les multiplications }

         {détermination des coefficients de rugosité}
         auxupuu1DW :=exp(Cmul(Csqr(Csub(qquu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupuu2DW :=exp(Cmul(Csqr(Cadd(qquu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupud1DW :=exp(Cmul(Csqr(Csub(qqdu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupud2DW :=exp(Cmul(Csqr(Cadd(qqdu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupdu1DW :=exp(Cmul(Csqr(Csub(qquu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdu2DW :=exp(Cmul(Csqr(Cadd(qquu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdd1DW :=exp(Cmul(Csqr(Csub(qqdu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdd2DW :=exp(Cmul(Csqr(Cadd(qqdu0,qqdu1)),Cplx(-sig*sig/2,0)));
         { détermination de la matrice de transfert }
         { coefficient p ++}
         ppuu:=Cdiv(Cadd(qquu1,qquu0),qquu1);
         { coefficient p +-}
         ppud:=Cdiv(Cadd(qquu1,qqdu0),qquu1);
         { coefficient p -+}
         ppdu:=Cdiv(Cadd(qqdu1,qquu0),qqdu1);
         { coefficient p --}
         ppdd:=Cdiv(Cadd(qqdu1,qqdu0),qqdu1);
         { coefficient m ++}
         mmuu:=Cdiv(Csub(qquu1,qquu0),qquu1);
         { coefficient m +-}
         mmud:=Cdiv(Csub(qquu1,qqdu0),qquu1);
         { coefficient m -+}
         mmdu:=Cdiv(Csub(qqdu1,qquu0),qqdu1);
         { coefficient m --}
         mmdd:=Cdiv(Csub(qqdu1,qqdu0),qqdu1);
         { paramètres de "propagation"}
         auxupuu1 := Exp(Cmul(C_I,Cmul(Csub(qquu0,qquu1),zz0)));
         auxupuu2 := Exp(Cmul(C_I,Cmul(Cadd(qquu0,qquu1),zz0)));
         auxupud1 := Exp(Cmul(C_I,Cmul(Csub(qqdu0,qquu1),zz0)));
         auxupud2 := Exp(Cmul(C_I,Cmul(Cadd(qqdu0,qquu1),zz0)));
         auxupdu1 := Exp(Cmul(C_I,Cmul(Csub(qquu0,qqdu1),zz0)));
         auxupdu2 := Exp(Cmul(C_I,Cmul(Cadd(qquu0,qqdu1),zz0)));
         auxupdd1 := Exp(Cmul(C_I,Cmul(Csub(qqdu0,qqdu1),zz0)));
         auxupdd2 := Exp(Cmul(C_I,Cmul(Cadd(qqdu0,qqdu1),zz0)));
         {termes "trigonomètriques"}
         tt1:= Cplx(BeffTheta[n+1]/2,0);
         tt2:= Cplx(BeffTheta[n]/2,0);
         phi1:= Cplx(BeffPhi[n+1],0);
         phi2:= Cplx(BeffPhi[n],0);
         trigup1:= Cadd(Cmul(cos(tt2),cos(tt1)),Cmul(Exp(Cmul(C_I,Csub(phi2,phi1))),Cmul(sin(tt1),sin(tt2))));
         trigup4:= Cadd(Cmul(cos(tt2),cos(tt1)),Cmul(Exp(Cmul(C_I,Csub(phi1,phi2))),Cmul(sin(tt1),sin(tt2))));
         trigup2:= Csub(Cmul(Exp(Cmul(Neg(C_I),phi1)),Cmul(cos(tt2),sin(tt1))),Cmul(Exp(Cmul(Neg(C_I),phi2)),Cmul(cos(tt1),sin(tt2))));
         trigup3:=Csub(Cmul(Exp(Cmul(C_I,phi2)),Cmul(cos(tt1),sin(tt2))),Cmul(Exp(Cmul(C_I,phi1)),Cmul(cos(tt2),sin(tt1))));
         {bloc matrice ++}
         DD.Cvalues[0,0]:=Cmul(Cmul(Cmul(ppuu,auxupuu1),auxupuu1DW),trigup1);
         DD.Cvalues[1,1]:=Cmul(Cmul(Cmul(ppuu,CDiv(C_ONE,auxupuu1)),auxupuu1DW),trigup1);
         DD.Cvalues[1,0]:=Cmul(Cmul(Cmul(mmuu,auxupuu2),auxupuu2DW),trigup1);
         DD.Cvalues[0,1]:=Cmul(Cmul(Cmul(mmuu,CDiv(C_ONE,auxupuu2)),auxupuu2DW),trigup1);
         {bloc matrice +-}
         DD.Cvalues[0,2]:=Cmul(Cmul(Cmul(ppud,auxupud1),auxupud1DW),trigup2);
         DD.Cvalues[1,3]:=Cmul(Cmul(Cmul(ppud,CDiv(C_ONE,auxupud1)),auxupud1DW),trigup2);
         DD.Cvalues[1,2]:=Cmul(Cmul(Cmul(mmud,auxupud2),auxupud2DW),trigup2);
         DD.Cvalues[0,3]:=Cmul(Cmul(Cmul(mmud,CDiv(C_ONE,auxupud2)),auxupud2DW),trigup2);
         {bloc matrice -+}
         DD.Cvalues[2,0]:=Cmul(Cmul(Cmul(ppdu,auxupdu1),auxupdu1DW),trigup3);
         DD.Cvalues[3,1]:=Cmul(Cmul(Cmul(ppdu,CDiv(C_ONE,auxupdu1)),auxupdu1DW),trigup3);
         DD.Cvalues[3,0]:=Cmul(Cmul(Cmul(mmdu,auxupdu2),auxupdu2DW),trigup3);
         DD.Cvalues[2,1]:=Cmul(Cmul(Cmul(mmdu,CDiv(C_ONE,auxupdu2)),auxupdu2DW),trigup3);
         {bloc matrice --}
         DD.Cvalues[2,2]:=Cmul(Cmul(Cmul(ppdd,auxupdd1),auxupdd1DW),trigup4);
         DD.Cvalues[3,3]:=Cmul(Cmul(Cmul(ppdd,CDiv(C_ONE,auxupdd1)),auxupdd1DW),trigup4);
         DD.Cvalues[3,2]:=Cmul(Cmul(Cmul(mmdd,auxupdd2),auxupdd2DW),trigup4);
         DD.Cvalues[2,3]:=Cmul(Cmul(Cmul(mmdd,CDiv(C_ONE,auxupdd2)),auxupdd2DW),trigup4);
    {            on remplace pp par rr et mm par nn.....  }
         auxdouu1DW :=Exp(Cmul(Csqr(Csub(qqud0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdouu2DW :=Exp(Cmul(Csqr(Cadd(qqud0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdoud1DW :=Exp(Cmul(Csqr(Csub(qqdd0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdoud2DW :=Exp(Cmul(Csqr(Cadd(qqdd0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdodu1DW :=Exp(Cmul(Csqr(Csub(qqud0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodu2DW :=Exp(Cmul(Csqr(Cadd(qqud0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodd1DW :=Exp(Cmul(Csqr(Csub(qqdd0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodd2DW :=Exp(Cmul(Csqr(Cadd(qqdd0,qqdd1)),Cplx(-sig*sig/2,0)));

         { coefficient p ++}
         rruu:=Cdiv(Cadd(qqud1,qqud0),qqud1);rruu:=CMul(Cplx(0.5,0),rruu);
         { coefficient p +-}
         rrud:=Cdiv(Cadd(qqud1,qqdd0),qqud1);rrud:=CMul(Cplx(0.5,0),rrud);
         { coefficient p -+}
         rrdu:=Cdiv(Cadd(qqdd1,qqud0),qqdd1);rrdu:=CMul(Cplx(0.5,0),rrdu);
         { coefficient p --}
         rrdd:=Cdiv(Cadd(qqdd1,qqdd0),qqdd1);rrdd:=CMul(Cplx(0.5,0),rrdd);
         { coefficient m ++}
         nnuu:=Cdiv(Csub(qqud1,qqud0),qqud1);nnuu:=CMul(Cplx(0.5,0),nnuu);
         { coefficient m +-}
         nnud:=Cdiv(Csub(qqud1,qqdd0),qqud1);nnud:=CMul(Cplx(0.5,0),nnud);
         { coefficient m -+}
         nndu:=Cdiv(Csub(qqdd1,qqud0),qqdd1);nndu:=CMul(Cplx(0.5,0),nndu);
         { coefficient m --}
         nndd:=Cdiv(Csub(qqdd1,qqdd0),qqdd1);nndd:=CMul(Cplx(0.5,0),nndd);
         { initialisation of the D matrix}
         auxdouu1 := Exp(Cmul(C_I,Cmul(Csub(qqud0,qqud1),zz0)));
         auxdouu2 := Exp(Cmul(C_I,Cmul(Cadd(qqud0,qqud1),zz0)));
         auxdoud1 := Exp(Cmul(C_I,Cmul(Csub(qqdd0,qqud1),zz0)));
         auxdoud2 := Exp(Cmul(C_I,Cmul(Cadd(qqdd0,qqud1),zz0)));
         auxdodu1 := Exp(Cmul(C_I,Cmul(Csub(qqud0,qqdd1),zz0)));
         auxdodu2 := Exp(Cmul(C_I,Cmul(Cadd(qqud0,qqdd1),zz0)));
         auxdodd1 := Exp(Cmul(C_I,Cmul(Csub(qqdd0,qqdd1),zz0)));
         auxdodd2 := Exp(Cmul(C_I,Cmul(Cadd(qqdd0,qqdd1),zz0)));

         {bloc matrice ++}
         DD.Cvalues[4,4]:=Cmul(Cmul(Cmul(rrdd,auxdodd1),auxdodd1DW),Trigup4);
         DD.Cvalues[5,5]:=Cmul(Cmul(Cmul(rrdd,CDiv(C_ONE,auxdodd1)),auxdodd1DW),Trigup4);
         DD.Cvalues[5,4]:=Cmul(Cmul(Cmul(nndd,auxdodd2),auxdodd2DW),Trigup4);
         DD.Cvalues[4,5]:=Cmul(Cmul(Cmul(nndd,CDiv(C_ONE,auxdodd2)),auxdodd2DW),Trigup4);
         {bloc matrice +-}
         DD.Cvalues[4,6]:=Cmul(Cmul(Cmul(rrdu,auxdodu1),auxdodu1DW),Trigup3);
         DD.Cvalues[5,7]:=Cmul(Cmul(Cmul(rrdu,CDiv(C_ONE,auxdodu1)),auxdodu1DW),Trigup3);
         DD.Cvalues[5,6]:=Cmul(Cmul(Cmul(nndu,auxdodu2),auxdodu2DW),Trigup3);
         DD.Cvalues[4,7]:=Cmul(Cmul(Cmul(nndu,CDiv(C_ONE,auxdodu2)),auxdodu2DW),Trigup3);
         {bloc matrice -+}
         DD.Cvalues[6,4]:=Cmul(Cmul(Cmul(rrud,auxdoud1),auxdoud1DW),Trigup2);
         DD.Cvalues[7,5]:=Cmul(Cmul(Cmul(rrud,CDiv(C_ONE,auxdoud1)),auxdoud1DW),Trigup2);
         DD.Cvalues[7,4]:=Cmul(Cmul(Cmul(nnud,auxdoud2),auxdoud2DW),Trigup2);
         DD.Cvalues[6,5]:=Cmul(Cmul(Cmul(nnud,CDiv(C_ONE,auxdoud2)),auxdoud2DW),Trigup2);
         {bloc matrice --}
         DD.Cvalues[6,6]:=Cmul(Cmul(Cmul(rruu,auxdouu1),auxdouu1DW),Trigup1);
         DD.Cvalues[7,7]:=Cmul(Cmul(Cmul(rruu,CDiv(C_ONE,auxdouu1)),auxdouu1DW),Trigup1);
         DD.Cvalues[7,6]:=Cmul(Cmul(Cmul(nnuu,auxdouu2),auxdouu2DW),Trigup1);
         DD.Cvalues[6,7]:=Cmul(Cmul(Cmul(nnuu,CDiv(C_ONE,auxdouu2)),auxdouu2DW),Trigup1);
        { DD.Cvalues[0,4]:=0;
         DD.Cvalues[0,5]:=0;  }
{         DD.scale(0.5); } { ne marche pas correctement }
         Maux.Mul(M,DD);
         M.copy(Maux);
         qquu1:=qquu0;
         qqdu1:=qqdu0;
         qqud1:=qqud0;
         qqdd1:=qqdd0;

      end;
      aux11:=M.Cvalues[0,0];aux13:=M.Cvalues[0,2];
      aux21:=M.Cvalues[1,0];aux23:=M.Cvalues[1,2];
      aux31:=M.Cvalues[2,0];aux33:=M.Cvalues[2,2];
      aux41:=M.Cvalues[3,0];aux43:=M.Cvalues[3,2];
      aux55:=M.Cvalues[4,4];aux57:=M.Cvalues[4,6];
      aux65:=M.Cvalues[5,4];aux67:=M.Cvalues[5,6];
      aux75:=M.Cvalues[6,4];aux77:=M.Cvalues[6,6];
      aux85:=M.Cvalues[7,4];aux87:=M.Cvalues[7,6];
      auxupD:=Csub(Cmul(aux11,aux33),Cmul(aux13,aux31));
      auxdoD:=Csub(Cmul(aux55,aux77),Cmul(aux57,aux75));
      Ruu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux21,aux33),Cmul(aux23,aux31)),auxupD)));
      Rdu.values[i]:=Sqr(Cabs(CDiv(Csub(Cmul(aux41,aux33),Cmul(aux43,aux31)),auxupD)));
      Rdd.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux65,aux77),Cmul(aux67,aux75)),auxdoD)));
      Rud.values[i]:=Sqr(Cabs(CDiv(Csub(Cmul(aux85,aux77),Cmul(aux87,aux75)),auxdoD)));

   end;
   freeit(DD,M,Maux);
   freeit(quu,qdu,qud,qdd);
   freeIt(Equu,Eqdd,Eim);
   freeIt(Equd,Eqdu);
end;

{ version 1.3    F. Ott                                       }
{ procédure de calcul de la réflectivité (neut. polarisés)    }
{ on ne tient pas compte des k parallleles                    }
{prise en compte de la rugosite                               }
{ INPUT: longueur d'onde en nm                                }
{        angles, vecteur des angles d'incidence (en  °)       }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vecteur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ REM: aucun paramètre implicite!!                            }
{      dans cette version, les matrices auxilaires D,P,M et Maux       }
{      sont créées à l'initialisation (gain de 10% seulement)          }
{ * V10 to V11 : simplification du calcul des variables qquu1 et qqdd1                         }
{                les noms de ces variables sont vaseux qquu1 et qqdd1 seraient mieux!          }
{ * V11 to V12 : lambda n'est plus passé comme paramètre        }
{                = prmemière version utilisable avec la V1.4    }
{ * V12 to V13 : optimisation de la vitesse de calcul           }
{                actuellement, l'utilisation de variables complexes auxilaires ralentit le calcul : Cplx(-sig*sig/2,0), Cplx(0.5,0) }
{                ca doit être du à la librairie  Cplx   }
{                * D.scale(0.5) gain de 0.5s seulement! }
{                * sig2=-sig*sig/2 -> perte de 10% en vitesse !! }
{               POUR LE MOMENT PAS D'EXPLICATION A CE COMPORTEMENT BIZARRE      }
procedure CalculReflecSF_HF_V13(var angles:Tvec;Qz:Tvec;var z,Vu,Vd,sigma,BeffMag,BeffTheta,BeffPhi:Tvec;var Ruu,Rdu,Rud,Rdd:Tvec);
var i,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    zz:double;
    E0,q0:double;
    qquu0,qqdu0,qqud0,qqdd0,qquu1,qqdu1,qqud1,qqdd1,zz0:TCplx;
    quu,qdu,qud,qdd:Tvec;            { vecteur d'onde perpendiculaire dans chaque couche }
    Equu,Eqdd,Equd,Eqdu,Eim:Tvec;    { Energie cinétique perp. dans chaque couche }
    Vaux:Tvec;
    ppuu,ppud,ppdu,ppdd,mmuu,mmud,mmdu,mmdd:Tcplx; {coefficient p et m (Gibaud p.101)}
    rruu,rrud,rrdu,rrdd,nnuu,nnud,nndu,nndd:Tcplx;
    auxupuu1,auxupuu2,auxupud1,auxupud2,auxupdu1,auxupdu2,auxupdd1,auxupdd2:Tcplx; {coefficient de propagation}
    auxdouu1,auxdouu2,auxdoud1,auxdoud2,auxdodu1,auxdodu2,auxdodd1,auxdodd2:Tcplx;
    sig:double;
    auxupuu1DW,auxupuu2DW,auxupud1DW,auxupud2DW:Tcplx; {coefficient de rugosité}
    auxupdu1DW,auxupdu2DW,auxupdd1DW,auxupdd2DW:Tcplx;
    auxdouu1DW,auxdouu2DW,auxdoud1Dw,auxdoud2Dw:Tcplx;
    auxdodu1DW,auxdodu2DW,auxdodd1DW,auxdodd2DW:Tcplx;
    DD,M,Maux:TMtx;
    auxupD,aux11,aux13,aux21,aux23,aux31,aux33,aux41,aux43:TCplx;
    auxdoD,aux55,aux57,aux65,aux67,aux75,aux77,aux85,aux87:TCplx;
    tt1,tt2,phi1,phi2:TCplx;
    trigup1,trigup2,trigup3,trigup4:TCplx;
{    sig2:double; }
{    half:TCplx; }
begin
{   half:=Cplx(0.5,0); }
   createIt(Equu,Eqdd,Eim);
   createIt(Equd,Eqdu);
   createIt(Vaux);
   createIt(DD);DD.size(8,8,true);
   createIt(M);M.size(8,8,true);
   createIt(Maux);Maux.size(8,8,true);
   nbrLayers := z.length;
   nbrPoints := Qz.length;
   Vaux.Size(nbrLayers,false);
   Equu.Size(nbrLayers,false);Equd.Size(nbrLayers,false);
   Eqdu.Size(nbrLayers,false);Eqdd.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0);
   createIt(quu,qdu,qud,qdd);
   quu.Size(nbrLayers,true);
   qud.Size(nbrLayers,true);
   qdu.Size(nbrLayers,true);
   qdd.Size(nbrLayers,true);

   Ruu.Size(Qz);Rdu.Size(Qz);Rud.Size(Qz);Rdd.Size(Qz);
   for i:=0 to nbrPoints-1 do
   begin
      DD.eye(8,8,true);
      M.eye(8,8,true); { initialisation de la récursion }
      q0 := Qz.values[i]/2;   { on a supprimé *DEG }

      Vaux.setVal(Cplx(mun*BeffMag[nbrlayers-1]*Cos(angles.values[i])*Cos(angles.values[i]),0));
      { a quoi correspond ce potentiel ?? }

      {calcul des q++,q+-,....}
      Equu.setVal(hb2over2mn);Eqdd.setVal(hb2over2mn);
      Equd.setVal(hb2over2mn);Eqdu.setVal(hb2over2mn);
      Equu.scale(q0*q0);Eqdd.scale(q0*q0);
      Equd.scale(q0*q0);Eqdu.scale(q0*q0);

      Equu.sub(Vu);Eqdd.sub(Vd);
      Equd.sub(Vu);Eqdu.sub(Vd);
      Equu.sub(Vaux);Eqdd.add(Vaux);
      Equd.add(Vaux);Eqdu.sub(Vaux);
      quu.RealToCplx(Equu,Eim);
      quu.scale(1/hb2over2mn);
      quu.sqrt; { q peut-etre complexe, racine }
      qdu.RealToCplx(Eqdu,Eim);qdu.scale(1/hb2over2mn);qdu.sqrt; { q peut-etre complexe, racine }
      qdd.RealToCplx(Eqdd,Eim);qdd.scale(1/hb2over2mn);qdd.sqrt;
      qud.RealToCplx(Equd,Eim);qud.scale(1/hb2over2mn);qud.sqrt;

      { initialisation de variables pour la récursion }
      qquu1:=quu.Cvalues[nbrLayers-1]; qqdu1:=qdd.Cvalues[nbrLayers-1]; { ces lignes avaient été mises en commentaires }
      qqud1:=quu.Cvalues[nbrLayers-1]; qqdd1:=qdd.Cvalues[nbrLayers-1];
      zz:=0;

      for n:=nbrLayers-2 downto 0 do   (* boucle du vide vers le substrat *)
      begin
         qquu0:=quu.Cvalues[n];
         qqdu0:=qdu.Cvalues[n];
         qqud0:=qud.Cvalues[n];
         qqdd0:=qdd.Cvalues[n];

         sig:=sigma.values[n]; { rugosité de l'interface j+1/j }
{         sig2:=-sig*sig/2; }
         zz:=zz+z.values[n+1]; { on fixe le séro à la surface de la couche }
         zz0:=Cplx(zz,0);      { passage en complexe pour faire les multiplications }

         {détermination des coefficients de rugosité}
         auxupuu1DW :=exp(Cmul(Csqr(Csub(qquu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupuu2DW :=exp(Cmul(Csqr(Cadd(qquu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupud1DW :=exp(Cmul(Csqr(Csub(qqdu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupud2DW :=exp(Cmul(Csqr(Cadd(qqdu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupdu1DW :=exp(Cmul(Csqr(Csub(qquu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdu2DW :=exp(Cmul(Csqr(Cadd(qquu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdd1DW :=exp(Cmul(Csqr(Csub(qqdu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdd2DW :=exp(Cmul(Csqr(Cadd(qqdu0,qqdu1)),Cplx(-sig*sig/2,0)));
         { détermination de la matrice de transfert }
         { coefficient p ++} ppuu:=Cdiv(Cadd(qquu1,qquu0),qquu1);{ ppuu:=CMul(Cplx(0.5,0),ppuu); }
         { coefficient p +-} ppud:=Cdiv(Cadd(qquu1,qqdu0),qquu1);{ ppud:=CMul(Cplx(0.5,0),ppud); }
         { coefficient p -+} ppdu:=Cdiv(Cadd(qqdu1,qquu0),qqdu1);{ ppdu:=CMul(Cplx(0.5,0),ppdu); }
         { coefficient p --} ppdd:=Cdiv(Cadd(qqdu1,qqdu0),qqdu1);{ ppdd:=CMul(Cplx(0.5,0),ppdd); }
         { coefficient m ++} mmuu:=Cdiv(Csub(qquu1,qquu0),qquu1);{ mmuu:=CMul(Cplx(0.5,0),mmuu); }
         { coefficient m +-} mmud:=Cdiv(Csub(qquu1,qqdu0),qquu1);{ mmud:=CMul(Cplx(0.5,0),mmud); }
         { coefficient m -+} mmdu:=Cdiv(Csub(qqdu1,qquu0),qqdu1);{ mmdu:=CMul(Cplx(0.5,0),mmdu); }
         { coefficient m --} mmdd:=Cdiv(Csub(qqdu1,qqdu0),qqdu1);{ mmdd:=CMul(Cplx(0.5,0),mmdd); }
         { paramètres de "propagation"}
         auxupuu1 := Exp(Cmul(C_I,Cmul(Csub(qquu0,qquu1),zz0)));
         auxupuu2 := Exp(Cmul(C_I,Cmul(Cadd(qquu0,qquu1),zz0)));
         auxupud1 := Exp(Cmul(C_I,Cmul(Csub(qqdu0,qquu1),zz0)));
         auxupud2 := Exp(Cmul(C_I,Cmul(Cadd(qqdu0,qquu1),zz0)));
         auxupdu1 := Exp(Cmul(C_I,Cmul(Csub(qquu0,qqdu1),zz0)));
         auxupdu2 := Exp(Cmul(C_I,Cmul(Cadd(qquu0,qqdu1),zz0)));
         auxupdd1 := Exp(Cmul(C_I,Cmul(Csub(qqdu0,qqdu1),zz0)));
         auxupdd2 := Exp(Cmul(C_I,Cmul(Cadd(qqdu0,qqdu1),zz0)));
         {termes "trigonomètriques"}
         tt1:= Cplx(BeffTheta[n+1]/2,0);
         tt2:= Cplx(BeffTheta[n]/2,0);
         phi1:= Cplx(BeffPhi[n+1],0);
         phi2:= Cplx(BeffPhi[n],0);
         trigup1:= Cadd(Cmul(cos(tt2),cos(tt1)),Cmul(Exp(Cmul(C_I,Csub(phi2,phi1))),Cmul(sin(tt1),sin(tt2))));
         trigup4:= Cadd(Cmul(cos(tt2),cos(tt1)),Cmul(Exp(Cmul(C_I,Csub(phi1,phi2))),Cmul(sin(tt1),sin(tt2))));
         trigup2:= Csub(Cmul(Exp(Cmul(Neg(C_I),phi1)),Cmul(cos(tt2),sin(tt1))),Cmul(Exp(Cmul(Neg(C_I),phi2)),Cmul(cos(tt1),sin(tt2))));
         trigup3:=Csub(Cmul(Exp(Cmul(C_I,phi2)),Cmul(cos(tt1),sin(tt2))),Cmul(Exp(Cmul(C_I,phi1)),Cmul(cos(tt2),sin(tt1))));
         {bloc matrice ++}
         DD.Cvalues[0,0]:=Cmul(Cmul(Cmul(ppuu,auxupuu1),auxupuu1DW),trigup1);
         DD.Cvalues[1,1]:=Cmul(Cmul(Cmul(ppuu,CDiv(C_ONE,auxupuu1)),auxupuu1DW),trigup1);
         DD.Cvalues[1,0]:=Cmul(Cmul(Cmul(mmuu,auxupuu2),auxupuu2DW),trigup1);
         DD.Cvalues[0,1]:=Cmul(Cmul(Cmul(mmuu,CDiv(C_ONE,auxupuu2)),auxupuu2DW),trigup1);
         {bloc matrice +-}
         DD.Cvalues[0,2]:=Cmul(Cmul(Cmul(ppud,auxupud1),auxupud1DW),trigup2);
         DD.Cvalues[1,3]:=Cmul(Cmul(Cmul(ppud,CDiv(C_ONE,auxupud1)),auxupud1DW),trigup2);
         DD.Cvalues[1,2]:=Cmul(Cmul(Cmul(mmud,auxupud2),auxupud2DW),trigup2);
         DD.Cvalues[0,3]:=Cmul(Cmul(Cmul(mmud,CDiv(C_ONE,auxupud2)),auxupud2DW),trigup2);
         {bloc matrice -+}
         DD.Cvalues[2,0]:=Cmul(Cmul(Cmul(ppdu,auxupdu1),auxupdu1DW),trigup3);
         DD.Cvalues[3,1]:=Cmul(Cmul(Cmul(ppdu,CDiv(C_ONE,auxupdu1)),auxupdu1DW),trigup3);
         DD.Cvalues[3,0]:=Cmul(Cmul(Cmul(mmdu,auxupdu2),auxupdu2DW),trigup3);
         DD.Cvalues[2,1]:=Cmul(Cmul(Cmul(mmdu,CDiv(C_ONE,auxupdu2)),auxupdu2DW),trigup3);
         {bloc matrice --}
         DD.Cvalues[2,2]:=Cmul(Cmul(Cmul(ppdd,auxupdd1),auxupdd1DW),trigup4);
         DD.Cvalues[3,3]:=Cmul(Cmul(Cmul(ppdd,CDiv(C_ONE,auxupdd1)),auxupdd1DW),trigup4);
         DD.Cvalues[3,2]:=Cmul(Cmul(Cmul(mmdd,auxupdd2),auxupdd2DW),trigup4);
         DD.Cvalues[2,3]:=Cmul(Cmul(Cmul(mmdd,CDiv(C_ONE,auxupdd2)),auxupdd2DW),trigup4);
    {            on remplace pp par rr et mm par nn.....  }
         auxdouu1DW :=Exp(Cmul(Csqr(Csub(qqud0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdouu2DW :=Exp(Cmul(Csqr(Cadd(qqud0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdoud1DW :=Exp(Cmul(Csqr(Csub(qqdd0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdoud2DW :=Exp(Cmul(Csqr(Cadd(qqdd0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdodu1DW :=Exp(Cmul(Csqr(Csub(qqud0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodu2DW :=Exp(Cmul(Csqr(Cadd(qqud0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodd1DW :=Exp(Cmul(Csqr(Csub(qqdd0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodd2DW :=Exp(Cmul(Csqr(Cadd(qqdd0,qqdd1)),Cplx(-sig*sig/2,0)));

         { coefficient p ++} rruu:=Cdiv(Cadd(qqud1,qqud0),qqud1);{rruu:=CMul(Cplx(0.5,0),rruu); }
         { coefficient p +-} rrud:=Cdiv(Cadd(qqud1,qqdd0),qqud1);{rrud:=CMul(Cplx(0.5,0),rrud); }
         { coefficient p -+} rrdu:=Cdiv(Cadd(qqdd1,qqud0),qqdd1);{rrdu:=CMul(Cplx(0.5,0),rrdu); }
         { coefficient p --} rrdd:=Cdiv(Cadd(qqdd1,qqdd0),qqdd1);{rrdd:=CMul(Cplx(0.5,0),rrdd); }
         { coefficient m ++} nnuu:=Cdiv(Csub(qqud1,qqud0),qqud1);{nnuu:=CMul(Cplx(0.5,0),nnuu); }
         { coefficient m +-} nnud:=Cdiv(Csub(qqud1,qqdd0),qqud1);{nnud:=CMul(Cplx(0.5,0),nnud); }
         { coefficient m -+} nndu:=Cdiv(Csub(qqdd1,qqud0),qqdd1);{nndu:=CMul(Cplx(0.5,0),nndu); }
         { coefficient m --} nndd:=Cdiv(Csub(qqdd1,qqdd0),qqdd1);{nndd:=CMul(Cplx(0.5,0),nndd); }
         { initialisation of the D matrix}
         auxdouu1 := Exp(Cmul(C_I,Cmul(Csub(qqud0,qqud1),zz0)));
         auxdouu2 := Exp(Cmul(C_I,Cmul(Cadd(qqud0,qqud1),zz0)));
         auxdoud1 := Exp(Cmul(C_I,Cmul(Csub(qqdd0,qqud1),zz0)));
         auxdoud2 := Exp(Cmul(C_I,Cmul(Cadd(qqdd0,qqud1),zz0)));
         auxdodu1 := Exp(Cmul(C_I,Cmul(Csub(qqud0,qqdd1),zz0)));
         auxdodu2 := Exp(Cmul(C_I,Cmul(Cadd(qqud0,qqdd1),zz0)));
         auxdodd1 := Exp(Cmul(C_I,Cmul(Csub(qqdd0,qqdd1),zz0)));
         auxdodd2 := Exp(Cmul(C_I,Cmul(Cadd(qqdd0,qqdd1),zz0)));

         {bloc matrice ++}
         DD.Cvalues[4,4]:=Cmul(Cmul(Cmul(rrdd,auxdodd1),auxdodd1DW),Trigup4);
         DD.Cvalues[5,5]:=Cmul(Cmul(Cmul(rrdd,CDiv(C_ONE,auxdodd1)),auxdodd1DW),Trigup4);
         DD.Cvalues[5,4]:=Cmul(Cmul(Cmul(nndd,auxdodd2),auxdodd2DW),Trigup4);
         DD.Cvalues[4,5]:=Cmul(Cmul(Cmul(nndd,CDiv(C_ONE,auxdodd2)),auxdodd2DW),Trigup4);
         {bloc matrice +-}
         DD.Cvalues[4,6]:=Cmul(Cmul(Cmul(rrdu,auxdodu1),auxdodu1DW),Trigup3);
         DD.Cvalues[5,7]:=Cmul(Cmul(Cmul(rrdu,CDiv(C_ONE,auxdodu1)),auxdodu1DW),Trigup3);
         DD.Cvalues[5,6]:=Cmul(Cmul(Cmul(nndu,auxdodu2),auxdodu2DW),Trigup3);
         DD.Cvalues[4,7]:=Cmul(Cmul(Cmul(nndu,CDiv(C_ONE,auxdodu2)),auxdodu2DW),Trigup3);
         {bloc matrice -+}
         DD.Cvalues[6,4]:=Cmul(Cmul(Cmul(rrud,auxdoud1),auxdoud1DW),Trigup2);
         DD.Cvalues[7,5]:=Cmul(Cmul(Cmul(rrud,CDiv(C_ONE,auxdoud1)),auxdoud1DW),Trigup2);
         DD.Cvalues[7,4]:=Cmul(Cmul(Cmul(nnud,auxdoud2),auxdoud2DW),Trigup2);
         DD.Cvalues[6,5]:=Cmul(Cmul(Cmul(nnud,CDiv(C_ONE,auxdoud2)),auxdoud2DW),Trigup2);
         {bloc matrice --}
         DD.Cvalues[6,6]:=Cmul(Cmul(Cmul(rruu,auxdouu1),auxdouu1DW),Trigup1);
         DD.Cvalues[7,7]:=Cmul(Cmul(Cmul(rruu,CDiv(C_ONE,auxdouu1)),auxdouu1DW),Trigup1);
         DD.Cvalues[7,6]:=Cmul(Cmul(Cmul(nnuu,auxdouu2),auxdouu2DW),Trigup1);
         DD.Cvalues[6,7]:=Cmul(Cmul(Cmul(nnuu,CDiv(C_ONE,auxdouu2)),auxdouu2DW),Trigup1);
        { DD.Cvalues[0,4]:=0;
         DD.Cvalues[0,5]:=0;  }
         DD.scale(0.5);  { ne marche pas correctement }
         Maux.Mul(M,DD);
         M.copy(Maux);
         qquu1:=qquu0;
         qqdu1:=qqdu0;
         qqud1:=qqud0;
         qqdd1:=qqdd0;

      end;
      aux11:=M.Cvalues[0,0];aux13:=M.Cvalues[0,2];
      aux21:=M.Cvalues[1,0];aux23:=M.Cvalues[1,2];
      aux31:=M.Cvalues[2,0];aux33:=M.Cvalues[2,2];
      aux41:=M.Cvalues[3,0];aux43:=M.Cvalues[3,2];
      aux55:=M.Cvalues[4,4];aux57:=M.Cvalues[4,6];
      aux65:=M.Cvalues[5,4];aux67:=M.Cvalues[5,6];
      aux75:=M.Cvalues[6,4];aux77:=M.Cvalues[6,6];
      aux85:=M.Cvalues[7,4];aux87:=M.Cvalues[7,6];
      auxupD:=Csub(Cmul(aux11,aux33),Cmul(aux13,aux31));
      auxdoD:=Csub(Cmul(aux55,aux77),Cmul(aux57,aux75));
      Ruu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux21,aux33),Cmul(aux23,aux31)),auxupD)));
      Rdu.values[i]:=Sqr(Cabs(CDiv(Csub(Cmul(aux41,aux33),Cmul(aux43,aux31)),auxupD)));
      Rdd.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux65,aux77),Cmul(aux67,aux75)),auxdoD)));
      Rud.values[i]:=Sqr(Cabs(CDiv(Csub(Cmul(aux85,aux77),Cmul(aux87,aux75)),auxdoD)));

   end;
   freeit(DD,M,Maux);
   freeit(quu,qdu,qud,qdd);
   freeIt(Equu,Eqdd,Eim);
   freeIt(Equd,Eqdu);
end;

{ version 1.3    F. Ott                                       }
{ procédure de calcul de la réflectivité (neut. polarisés)    }
{ on ne tient pas compte des k parallleles                    }
{prise en compte de la rugosite                               }
{ INPUT: longueur d'onde en nm                                }
{        angles, vecteur des angles d'incidence (en  °)       }
{        z et V, vecteur des épaisseurs et des potentiels des }
{        différentes couches                                  }
{ OUTPUT: un vecteur contenant les réflectivités correspondant aux }
{         différents angles d'incidence                       }
{ REM: aucun paramètre implicite!!                            }
{      dans cette version, les matrices auxilaires D,P,M et Maux       }
{      sont créées à l'initialisation (gain de 10% seulement)          }
{ * V10 to V11 : simplification du calcul des variables qquu1 et qqdd1                         }
{                les noms de ces variables sont vaseux qquu1 et qqdd1 seraient mieux!          }
{ * V11 to V12 : lambda n'est plus passé comme paramètre        }
{                = prmemière version utilisable avec la V1.4    }
{ * V12 to V13 : optimisation de la vitesse de calcul           }
{                actuellement, l'utilisation de variables complexes auxilaires ralentit le calcul : Cplx(-sig*sig/2,0), Cplx(0.5,0) }
{                ca doit être du à la librairie  Cplx   }
{                * D.scale(0.5) gain de 0.5s seulement! }
{                * sig2=-sig*sig/2 -> perte de 10% en vitesse !! }
{               POUR LE MOMENT PAS D'EXPLICATION A CE COMPORTEMENT BIZARRE      }
{       m'activation de l'optimisation n'a aucun effet                          }
{ * V13 to V14 : * calcul avec des matrices 4*4                                 }
{                  correction d'une fuite mémoire sur Vaux                      }
{                  = gain de 5% seulement                                       }
{                * précalcul des termes trogonométrique                         }
{                  = gain de 17%                                                }
procedure CalculReflecSF_HF(var angles:Tvec;Qz:Tvec;var z,Vu,Vd,Vi,sigma,BeffMag,BeffTheta,BeffPhi:Tvec;var Ruu,Rdu,Rud,Rdd:Tvec);
var i,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    zz:double;
    E0,q0,q02:double;
    qquu0,qqdu0,qqud0,qqdd0,qquu1,qqdu1,qqud1,qqdd1,zz0:TCplx;
    quu,qdu,qud,qdd:Tvec;            { vecteur d'onde perpendiculaire dans chaque couche }
    Equu,Eqdd,Equd,Eqdu,Eim:Tvec;    { Energie cinétique perp. dans chaque couche }
    Vaux:Tvec;
    ppuu,ppud,ppdu,ppdd,mmuu,mmud,mmdu,mmdd:Tcplx; {coefficient p et m (Gibaud p.101)}
    rruu,rrud,rrdu,rrdd,nnuu,nnud,nndu,nndd:Tcplx;
    auxupuu1,auxupuu2,auxupud1,auxupud2,auxupdu1,auxupdu2,auxupdd1,auxupdd2:Tcplx; {coefficient de propagation}
    auxdouu1,auxdouu2,auxdoud1,auxdoud2,auxdodu1,auxdodu2,auxdodd1,auxdodd2:Tcplx;
    sig:double;
    auxupuu1DW,auxupuu2DW,auxupud1DW,auxupud2DW:Tcplx; {coefficient de rugosité}
    auxupdu1DW,auxupdu2DW,auxupdd1DW,auxupdd2DW:Tcplx;
    auxdouu1DW,auxdouu2DW,auxdoud1Dw,auxdoud2Dw:Tcplx;
    auxdodu1DW,auxdodu2DW,auxdodd1DW,auxdodd2DW:Tcplx;
    DDu,DDd,Mu,Md,MauxU,Mauxd:TMtx;
    auxupD,aux11,aux13,aux21,aux23,aux31,aux33,aux41,aux43:TCplx;
    auxdoD,aux55,aux57,aux65,aux67,aux75,aux77,aux85,aux87:TCplx;
    tt1,tt2,phi1,phi2:TCplx;
    trigup1,trigup2,trigup3,trigup4:TVec;
{    sig2:TCplx; }
begin
   createIt(Equu,Eqdd,Eim);
   createIt(Equd,Eqdu);
   createIt(Vaux);
   createIt(DDu,DDd);DDu.size(4,4,true);DDd.size(4,4,true);
   createIt(Mu,Md);Mu.size(4,4,true);Md.size(4,4,true);
   createIt(MauxU,MauxD);MauxU.size(4,4,true);MauxD.size(4,4,true);
   createIt(trigup1,trigup2,trigup3,trigup4);
   nbrLayers := z.length;
   nbrPoints := Qz.length;
   Vaux.Size(nbrLayers,false);
   Equu.Size(nbrLayers,false);Equd.Size(nbrLayers,false);
   Eqdu.Size(nbrLayers,false);Eqdd.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0);
   createIt(quu,qdu,qud,qdd);
   quu.Size(nbrLayers,true);
   qud.Size(nbrLayers,true);
   qdu.Size(nbrLayers,true);
   qdd.Size(nbrLayers,true);
   trigup1.Size(nbrLayers,true);trigup2.Size(nbrLayers,true);trigup3.Size(nbrLayers,true);trigup4.Size(nbrLayers,true);

   { précalcul des termes "trigonométriques" }
   for n:=nbrLayers-2 downto 0 do   (* boucle du vide vers le substrat *)
   begin
         tt1:= Cplx(BeffTheta[n+1]/2,0);
         tt2:= Cplx(BeffTheta[n]/2,0);
         phi1:= Cplx(BeffPhi[n+1],0);
         phi2:= Cplx(BeffPhi[n],0);
         trigup1.Cvalues[n]:= Cadd(Cmul(cos(tt2),cos(tt1)),Cmul(Exp(Cmul(C_I,Csub(phi2,phi1))),Cmul(sin(tt1),sin(tt2))));
         trigup4.Cvalues[n]:= Cadd(Cmul(cos(tt2),cos(tt1)),Cmul(Exp(Cmul(C_I,Csub(phi1,phi2))),Cmul(sin(tt1),sin(tt2))));
         trigup2.Cvalues[n]:= Csub(Cmul(Exp(Cmul(Neg(C_I),phi1)),Cmul(cos(tt2),sin(tt1))),Cmul(Exp(Cmul(Neg(C_I),phi2)),Cmul(cos(tt1),sin(tt2))));
         trigup3.Cvalues[n]:= Csub(Cmul(Exp(Cmul(C_I,phi2)),Cmul(cos(tt1),sin(tt2))),Cmul(Exp(Cmul(C_I,phi1)),Cmul(cos(tt2),sin(tt1))));
   end;

   Ruu.Size(Qz);Rdu.Size(Qz);Rud.Size(Qz);Rdd.Size(Qz);
   for i:=0 to nbrPoints-1 do
   begin
      DDu.eye(4,4,true);DDd.eye(4,4,true);
      Mu.eye(4,4,true);Md.eye(4,4,true); { initialisation de la récursion }
      q0 := Qz.values[i]/2;   { on a supprimé *DEG }
      q02:=q0*q0;

      Vaux.setVal(Cplx(mun*BeffMag[nbrlayers-1]*Cos(angles.values[i])*Cos(angles.values[i]),0)); { a quoi correspond ce potentiel ?? }

      {calcul des q++,q+-,....}
      Equu.setVal(hb2over2mn);Eqdd.setVal(hb2over2mn);
      Equd.setVal(hb2over2mn);Eqdu.setVal(hb2over2mn);
      Equu.scale(q02);Eqdd.scale(q02);
      Equd.scale(q02);Eqdu.scale(q02);

      Eim.setval(0);Eim.sub(Vi); { NEW : to take into account the absorption }

      Equu.sub(Vu);Eqdd.sub(Vd);
      Equd.sub(Vu);Eqdu.sub(Vd);
      Equu.sub(Vaux);Eqdd.add(Vaux);
      Equd.add(Vaux);Eqdu.sub(Vaux);
      quu.RealToCplx(Equu,Eim);quu.scale(1/hb2over2mn);quu.sqrt; { q peut-etre complexe, racine }
      qdu.RealToCplx(Eqdu,Eim);qdu.scale(1/hb2over2mn);qdu.sqrt; { q peut-etre complexe, racine }
      qdd.RealToCplx(Eqdd,Eim);qdd.scale(1/hb2over2mn);qdd.sqrt;
      qud.RealToCplx(Equd,Eim);qud.scale(1/hb2over2mn);qud.sqrt;

      { initialisation de variables pour la récursion }
      qquu1:=quu.Cvalues[nbrLayers-1]; qqdu1:=qdd.Cvalues[nbrLayers-1]; { ces lignes avaient été mises en commentaires }
      qqud1:=quu.Cvalues[nbrLayers-1]; qqdd1:=qdd.Cvalues[nbrLayers-1];
      zz:=0;

      for n:=nbrLayers-2 downto 0 do   (* boucle du vide vers le substrat *)
      begin
         qquu0:=quu.Cvalues[n];
         qqdu0:=qdu.Cvalues[n];
         qqud0:=qud.Cvalues[n];
         qqdd0:=qdd.Cvalues[n];

         sig:=sigma.values[n]; { rugosité de l'interface j+1/j }
{         sig2:=Cplx(-sig*sig/2,0); }
         zz:=zz+z.values[n+1]; { on fixe le séro à la surface de la couche }
         zz0:=Cplx(zz,0);      { passage en complexe pour faire les multiplications }

         {détermination des coefficients de rugosité}
         auxupuu1DW :=exp(Cmul(Csqr(Csub(qquu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupuu2DW :=exp(Cmul(Csqr(Cadd(qquu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupud1DW :=exp(Cmul(Csqr(Csub(qqdu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupud2DW :=exp(Cmul(Csqr(Cadd(qqdu0,qquu1)),Cplx(-sig*sig/2,0)));
         auxupdu1DW :=exp(Cmul(Csqr(Csub(qquu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdu2DW :=exp(Cmul(Csqr(Cadd(qquu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdd1DW :=exp(Cmul(Csqr(Csub(qqdu0,qqdu1)),Cplx(-sig*sig/2,0)));
         auxupdd2DW :=exp(Cmul(Csqr(Cadd(qqdu0,qqdu1)),Cplx(-sig*sig/2,0)));
         { détermination de la matrice de transfert }
         { coefficient p ++} ppuu:=Cdiv(Cadd(qquu1,qquu0),qquu1);
         { coefficient p +-} ppud:=Cdiv(Cadd(qquu1,qqdu0),qquu1);
         { coefficient p -+} ppdu:=Cdiv(Cadd(qqdu1,qquu0),qqdu1);
         { coefficient p --} ppdd:=Cdiv(Cadd(qqdu1,qqdu0),qqdu1);
         { coefficient m ++} mmuu:=Cdiv(Csub(qquu1,qquu0),qquu1);
         { coefficient m +-} mmud:=Cdiv(Csub(qquu1,qqdu0),qquu1);
         { coefficient m -+} mmdu:=Cdiv(Csub(qqdu1,qquu0),qqdu1);
         { coefficient m --} mmdd:=Cdiv(Csub(qqdu1,qqdu0),qqdu1);
         { paramètres de "propagation"}
         auxupuu1 := Exp(Cmul(C_I,Cmul(Csub(qquu0,qquu1),zz0)));
         auxupuu2 := Exp(Cmul(C_I,Cmul(Cadd(qquu0,qquu1),zz0)));
         auxupud1 := Exp(Cmul(C_I,Cmul(Csub(qqdu0,qquu1),zz0)));
         auxupud2 := Exp(Cmul(C_I,Cmul(Cadd(qqdu0,qquu1),zz0)));
         auxupdu1 := Exp(Cmul(C_I,Cmul(Csub(qquu0,qqdu1),zz0)));
         auxupdu2 := Exp(Cmul(C_I,Cmul(Cadd(qquu0,qqdu1),zz0)));
         auxupdd1 := Exp(Cmul(C_I,Cmul(Csub(qqdu0,qqdu1),zz0)));
         auxupdd2 := Exp(Cmul(C_I,Cmul(Cadd(qqdu0,qqdu1),zz0)));

         {bloc matrice ++}
         DDu.Cvalues[0,0]:=Cmul(Cmul(Cmul(ppuu,auxupuu1),auxupuu1DW),trigup1.Cvalues[n]);
         DDu.Cvalues[1,1]:=Cmul(Cmul(Cmul(ppuu,CDiv(C_ONE,auxupuu1)),auxupuu1DW),trigup1.Cvalues[n]);
         DDu.Cvalues[1,0]:=Cmul(Cmul(Cmul(mmuu,auxupuu2),auxupuu2DW),trigup1.Cvalues[n]);
         DDu.Cvalues[0,1]:=Cmul(Cmul(Cmul(mmuu,CDiv(C_ONE,auxupuu2)),auxupuu2DW),trigup1.Cvalues[n]);
         {bloc matrice +-}
         DDu.Cvalues[0,2]:=Cmul(Cmul(Cmul(ppud,auxupud1),auxupud1DW),trigup2.Cvalues[n]);
         DDu.Cvalues[1,3]:=Cmul(Cmul(Cmul(ppud,CDiv(C_ONE,auxupud1)),auxupud1DW),trigup2.Cvalues[n]);
         DDu.Cvalues[1,2]:=Cmul(Cmul(Cmul(mmud,auxupud2),auxupud2DW),trigup2.Cvalues[n]);
         DDu.Cvalues[0,3]:=Cmul(Cmul(Cmul(mmud,CDiv(C_ONE,auxupud2)),auxupud2DW),trigup2.Cvalues[n]);
         {bloc matrice -+}
         DDu.Cvalues[2,0]:=Cmul(Cmul(Cmul(ppdu,auxupdu1),auxupdu1DW),trigup3.Cvalues[n]);
         DDu.Cvalues[3,1]:=Cmul(Cmul(Cmul(ppdu,CDiv(C_ONE,auxupdu1)),auxupdu1DW),trigup3.Cvalues[n]);
         DDu.Cvalues[3,0]:=Cmul(Cmul(Cmul(mmdu,auxupdu2),auxupdu2DW),trigup3.Cvalues[n]);
         DDu.Cvalues[2,1]:=Cmul(Cmul(Cmul(mmdu,CDiv(C_ONE,auxupdu2)),auxupdu2DW),trigup3.Cvalues[n]);
         {bloc matrice --}
         DDu.Cvalues[2,2]:=Cmul(Cmul(Cmul(ppdd,auxupdd1),auxupdd1DW),trigup4.Cvalues[n]);
         DDu.Cvalues[3,3]:=Cmul(Cmul(Cmul(ppdd,CDiv(C_ONE,auxupdd1)),auxupdd1DW),trigup4.Cvalues[n]);
         DDu.Cvalues[3,2]:=Cmul(Cmul(Cmul(mmdd,auxupdd2),auxupdd2DW),trigup4.Cvalues[n]);
         DDu.Cvalues[2,3]:=Cmul(Cmul(Cmul(mmdd,CDiv(C_ONE,auxupdd2)),auxupdd2DW),trigup4.Cvalues[n]);
    {            on remplace pp par rr et mm par nn.....  }
         auxdouu1DW :=Exp(Cmul(Csqr(Csub(qqud0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdouu2DW :=Exp(Cmul(Csqr(Cadd(qqud0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdoud1DW :=Exp(Cmul(Csqr(Csub(qqdd0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdoud2DW :=Exp(Cmul(Csqr(Cadd(qqdd0,qqud1)),Cplx(-sig*sig/2,0)));
         auxdodu1DW :=Exp(Cmul(Csqr(Csub(qqud0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodu2DW :=Exp(Cmul(Csqr(Cadd(qqud0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodd1DW :=Exp(Cmul(Csqr(Csub(qqdd0,qqdd1)),Cplx(-sig*sig/2,0)));
         auxdodd2DW :=Exp(Cmul(Csqr(Cadd(qqdd0,qqdd1)),Cplx(-sig*sig/2,0)));

         { coefficient p ++} rruu:=Cdiv(Cadd(qqud1,qqud0),qqud1);
         { coefficient p +-} rrud:=Cdiv(Cadd(qqud1,qqdd0),qqud1);
         { coefficient p -+} rrdu:=Cdiv(Cadd(qqdd1,qqud0),qqdd1);
         { coefficient p --} rrdd:=Cdiv(Cadd(qqdd1,qqdd0),qqdd1);
         { coefficient m ++} nnuu:=Cdiv(Csub(qqud1,qqud0),qqud1);
         { coefficient m +-} nnud:=Cdiv(Csub(qqud1,qqdd0),qqud1);
         { coefficient m -+} nndu:=Cdiv(Csub(qqdd1,qqud0),qqdd1);
         { coefficient m --} nndd:=Cdiv(Csub(qqdd1,qqdd0),qqdd1);
         { initialisation of the D matrix}
         auxdouu1 := Exp(Cmul(C_I,Cmul(Csub(qqud0,qqud1),zz0)));
         auxdouu2 := Exp(Cmul(C_I,Cmul(Cadd(qqud0,qqud1),zz0)));
         auxdoud1 := Exp(Cmul(C_I,Cmul(Csub(qqdd0,qqud1),zz0)));
         auxdoud2 := Exp(Cmul(C_I,Cmul(Cadd(qqdd0,qqud1),zz0)));
         auxdodu1 := Exp(Cmul(C_I,Cmul(Csub(qqud0,qqdd1),zz0)));
         auxdodu2 := Exp(Cmul(C_I,Cmul(Cadd(qqud0,qqdd1),zz0)));
         auxdodd1 := Exp(Cmul(C_I,Cmul(Csub(qqdd0,qqdd1),zz0)));
         auxdodd2 := Exp(Cmul(C_I,Cmul(Cadd(qqdd0,qqdd1),zz0)));

         {bloc matrice ++}
         DDd.Cvalues[0,0]:=Cmul(Cmul(Cmul(rrdd,auxdodd1),auxdodd1DW),Trigup4.Cvalues[n]);
         DDd.Cvalues[1,1]:=Cmul(Cmul(Cmul(rrdd,CDiv(C_ONE,auxdodd1)),auxdodd1DW),Trigup4.Cvalues[n]);
         DDd.Cvalues[1,0]:=Cmul(Cmul(Cmul(nndd,auxdodd2),auxdodd2DW),Trigup4.Cvalues[n]);
         DDd.Cvalues[0,1]:=Cmul(Cmul(Cmul(nndd,CDiv(C_ONE,auxdodd2)),auxdodd2DW),Trigup4.Cvalues[n]);
         {bloc matrice +-}
         DDd.Cvalues[0,2]:=Cmul(Cmul(Cmul(rrdu,auxdodu1),auxdodu1DW),Trigup3.Cvalues[n]);
         DDd.Cvalues[1,3]:=Cmul(Cmul(Cmul(rrdu,CDiv(C_ONE,auxdodu1)),auxdodu1DW),Trigup3.Cvalues[n]);
         DDd.Cvalues[1,2]:=Cmul(Cmul(Cmul(nndu,auxdodu2),auxdodu2DW),Trigup3.Cvalues[n]);
         DDd.Cvalues[0,3]:=Cmul(Cmul(Cmul(nndu,CDiv(C_ONE,auxdodu2)),auxdodu2DW),Trigup3.Cvalues[n]);
         {bloc matrice -+}
         DDd.Cvalues[2,0]:=Cmul(Cmul(Cmul(rrud,auxdoud1),auxdoud1DW),Trigup2.Cvalues[n]);
         DDd.Cvalues[3,1]:=Cmul(Cmul(Cmul(rrud,CDiv(C_ONE,auxdoud1)),auxdoud1DW),Trigup2.Cvalues[n]);
         DDd.Cvalues[3,0]:=Cmul(Cmul(Cmul(nnud,auxdoud2),auxdoud2DW),Trigup2.Cvalues[n]);
         DDd.Cvalues[2,1]:=Cmul(Cmul(Cmul(nnud,CDiv(C_ONE,auxdoud2)),auxdoud2DW),Trigup2.Cvalues[n]);
         {bloc matrice --}
         DDd.Cvalues[2,2]:=Cmul(Cmul(Cmul(rruu,auxdouu1),auxdouu1DW),Trigup1.Cvalues[n]);
         DDd.Cvalues[3,3]:=Cmul(Cmul(Cmul(rruu,CDiv(C_ONE,auxdouu1)),auxdouu1DW),Trigup1.Cvalues[n]);
         DDd.Cvalues[3,2]:=Cmul(Cmul(Cmul(nnuu,auxdouu2),auxdouu2DW),Trigup1.Cvalues[n]);
         DDd.Cvalues[2,3]:=Cmul(Cmul(Cmul(nnuu,CDiv(C_ONE,auxdouu2)),auxdouu2DW),Trigup1.Cvalues[n]);

         DDu.scale(0.5);DDd.scale(0.5);  { ne marche pas correctement }
         MauxU.Mul(Mu,DDu);MauxD.Mul(Md,DDd);
         Mu.copy(MauxU);Md.copy(MauxD);
         qquu1:=qquu0;
         qqdu1:=qqdu0;
         qqud1:=qqud0;
         qqdd1:=qqdd0;

      end;
      aux11:=Mu.Cvalues[0,0];aux13:=Mu.Cvalues[0,2];
      aux21:=Mu.Cvalues[1,0];aux23:=Mu.Cvalues[1,2];
      aux31:=Mu.Cvalues[2,0];aux33:=Mu.Cvalues[2,2];
      aux41:=Mu.Cvalues[3,0];aux43:=Mu.Cvalues[3,2];
      aux55:=Md.Cvalues[0,0];aux57:=Md.Cvalues[0,2];
      aux65:=Md.Cvalues[1,0];aux67:=Md.Cvalues[1,2];
      aux75:=Md.Cvalues[2,0];aux77:=Md.Cvalues[2,2];
      aux85:=Md.Cvalues[3,0];aux87:=Md.Cvalues[3,2];
      auxupD:=Csub(Cmul(aux11,aux33),Cmul(aux13,aux31));
      auxdoD:=Csub(Cmul(aux55,aux77),Cmul(aux57,aux75));
      Ruu.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux21,aux33),Cmul(aux23,aux31)),auxupD)));
      Rdu.values[i]:=Sqr(Cabs(CDiv(Csub(Cmul(aux41,aux33),Cmul(aux43,aux31)),auxupD)));
      Rdd.values[i]:=Sqr(Cabs(Cdiv(Csub(Cmul(aux65,aux77),Cmul(aux67,aux75)),auxdoD)));
      Rud.values[i]:=Sqr(Cabs(CDiv(Csub(Cmul(aux85,aux77),Cmul(aux87,aux75)),auxdoD)));

   end;
   freeit(DDu,Mu,MauxU);freeit(DDd,Md,MauxD);
   freeit(quu,qdu,qud,qdd);
   freeIt(Equu,Eqdd,Eim);
   freeIt(Equd,Eqdu);
   freeIt(Vaux);
   freeIt(trigup1,trigup2,trigup3,trigup4);
end;


{ calcul de la réflectivité en tenant compte de la resolution             }
{ neutrons NON POLARISES                                                  }
{ INPUT : le spectre de resolution en Qz et sa distribution               }
{         nota : on supose que la distribution est la meme qqsoit Qz      }
procedure CalculateReflectivityNP(var DataCalc:TDataReflectivity;var z,V,Vi,sigma:Tvec;var R,T:Tvec);
var i,j:integer;
    Raux,Taux:Tvec;
    QzSPectrumAux,anglesAux:Tvec;
    QzFastMin,dtheta,l:double;
    qBot,Rbot,Tbot,qTop,Rtop,Ttop,q,QzCrt,sig,QzInt,dq,TotalIntensity,Raux1,Taux1:double;
    IdxBot,IdxTop,dQNbr:integer;
begin
    CreateIt(QzSpectrumAux,anglesAux);
    CreateIt(Raux,Taux);
    Raux.size(DataCalc.angles);Raux.setval(0);
    Taux.size(DataCalc.angles);Taux.setval(0);
    { on a différent type de convolution possibles }
    if DataCalc.convolutionCalculation = 'QzGaussian' { calcul dans l'espace Qz }
    then for i:=0 to DataCalc.QzIntensity.length-1 do
         begin
           CalculReflecNP(DataCalc.QzSpectrum[i],z,V,Vi,sigma,R,T);
           R.scale(DataCalc.QzIntensity[i]);T.scale(DataCalc.QzIntensity[i]);
           Raux.add(R);Taux.add(T);
         end
    else if DataCalc.convolutionCalculation = 'full' then
         { calcul dans le cas ou on veux faire une convolution complète }
         { ce calcul est valable pour une mesure à longueur d'onde fixe uniquement }
         for i:=0 to DataCalc.thetaSpectrum.length-1 do
         begin
           dTheta:=dataCalc.thetaSpectrum.values[i];
           anglesAux.copy(dataCalc.Angles);anglesAux.Offset(dTheta);
           for j:=0 to DataCalc.lambdaSpectrum.length-1 do
           begin
             l:=DataCalc.lambdaSpectrum.values[j];
             { il faut créer un vecteur Qz à partir des valeurs dTheta et l }
             QzSpectrumAux.Sin(AnglesAux);
             QzSpectrumAux.Scale(4*3.1416);QzSpectrumAux.Scale(1/l);
             CalculReflecNP(QzSpectrumAux,z,V,Vi,sigma,R,T);
             R.scale(DataCalc.lambdaIntensity.values[j]);T.scale(DataCalc.lambdaIntensity.values[j]);
             R.scale(DataCalc.thetaIntensity.values[i]);T.scale(DataCalc.thetaIntensity.values[i]);
             Raux.add(R);Taux.add(T);
           end;
         end
         else begin { cas de la convolution fast }
              QzFastMin:=DataCalc.QzFast.min;
              R.size(DataCalc.QzFast);Raux.setval(0);R.setval(0);
              T.size(DataCalc.QzFast);Taux.setval(0);T.setval(0);
              CalculReflecNP(DataCalc.QzFast,z,V,Vi,sigma,R,T); { calcul du spectre suréchantillonné }
              { reconstruction du spectre a partir de ces valeurs suréchantillonnées }
              for i:=0 to dataCalc.Qz.length-1 do
                begin
                  QzCrt:=dataCalc.Qz.values[i];sig:=dataCalc.QzSig.values[i];
                  dQNbr:=trunc(dataCalc.dQzNbr/2);
                  if dQNbr>0 then dq:=DataCalc.QzMult*sig*sqrt(2*ln(2))/dQNbr else dq:=0;
                  { generation de la gaussienne ayant cette largeur }
                  TotalIntensity:=0;Raux1:=0;Taux1:=0;
                  for j:=-dqNbr to dqNbr do
                  begin { il y a un problème quand qz est négatif }
                      q:=QzCrt+j*dq;{ NOTA abs est un bricolo pour que q soit toujours positif ... }
                      QzInt:=1/(sig*sqrt(2*pi))*exp(-(sqr(j*dq)/(2*sig*sig)));
                      TotalIntensity:=TotalIntensity+QzInt;
                      idxBot:=DataCalc.IdxInQzFastBot[round(abs(q-QzFastMin)*1e-6)];idxTop:=idxBot+1; { NEW offset de la valeur QzFast min ! }
                      qBot:=DataCalc.QzFast.values[IdxBot];qTop:=DataCalc.QzFast.values[IdxTop];
                      Rbot:=R.values[IdxBot];RTop:=R.values[IdxTop];
                      Tbot:=T.values[IdxBot];TTop:=T.values[IdxTop];
                      Raux1:=Raux1+(Rbot+(Rtop-Rbot)*(q-qBot)/(qTop-qBot))*QzInt;
                      Taux1:=Taux1+(Tbot+(Ttop-Tbot)*(q-qBot)/(qTop-qBot))*QzInt;
                  end;
                  Raux1:=Raux1/TotalIntensity;Taux1:=Taux1/TotalIntensity;
                  Raux.values[i]:=Raux1;Taux.values[i]:=Taux1; { pas nécessaire en principe si la gaussienne est correctement calculée !! }
                end;
              end;

    R.copy(Raux);T.copy(Taux);
    FreeIt(Raux,Taux);
    FreeIt(QzSpectrumAux,anglesAux);
end;

{ calcul de la réflectivité en tenant compte des différents cas possibles }
{ neutrons POLARISES                                                  }
{ ACHTUNG: entree implicite "lambda"                                  }
procedure CalculateReflectivityP(var DataCalc:TDataReflectivity;var z,Vup,Vdo,Vi,sigma:Tvec;var Rup,Rdo,Tup,Tdo:Tvec);
begin
     CalculateReflectivityNP(DataCalc,z,Vup,Vi,sigma,Rup,Tup);
     CalculateReflectivityNP(DataCalc,z,Vdo,Vi,sigma,Rdo,Tdo)
end;

{ ============================================================================= }
{ calcul de la réflectivité en SF en tenant compte des différents cas possibles }
{ neutrons POLARISES avec calcul du SF                                          }
procedure CalculateReflectivitySF(var DataCalc:TDataReflectivity;var z,Vup,Vdo,Vi,phi:Tvec;var Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu:Tvec);
var i,j:integer;
    Raux,Taux:Tvec;
    RuuAux,TuuAux,RddAux,TddAux,RudAux,TudAux,RduAux,TduAux:Tvec;
    QzSPectrumAux,anglesAux:Tvec;
    dtheta,l:double;
    QzFastMin,QzCrt,sig,q,dq,TotalIntensity,QzInt,qBot,qTop:double;
    dQNbr,idxBot,idxTop:integer;
    Ruuaux1,Rddaux1,Rudaux1,Rduaux1:double;
    RuuBot,RddBot,RudBot,RduBot:double;
    RuuTop,RddTop,RudTop,RduTop:double;
begin
    CreateIt(QzSpectrumAux,anglesAux);
    createIt(RuuAux,RddAux,RudAux,RduAux);createIt(TuuAux,TddAux,TudAux,TduAux);
    RuuAux.size(DataCalc.angles);RuuAux.setval(0);TuuAux.size(DataCalc.angles);TuuAux.setval(0);
    RddAux.size(DataCalc.angles);RddAux.setval(0);TddAux.size(DataCalc.angles);TddAux.setval(0);
    RudAux.size(DataCalc.angles);RudAux.setval(0);TudAux.size(DataCalc.angles);TudAux.setval(0);
    RduAux.size(DataCalc.angles);RduAux.setval(0);TduAux.size(DataCalc.angles);TduAux.setval(0);

    { on a différents types de convolution possibles }
    if DataCalc.convolutionCalculation = 'QzGaussian' { calcul dans l'espace Qz }
    then for i:=0 to DataCalc.QzIntensity.length-1 do
         begin
         CalculReflecSF(DataCalc.QzSpectrum[i],z,Vup,Vdo,Vi,phi,Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu);
         Ruu.scale(DataCalc.QzIntensity[i]);Tuu.scale(DataCalc.QzIntensity[i]);
         Rdd.scale(DataCalc.QzIntensity[i]);Tdd.scale(DataCalc.QzIntensity[i]);
         Rud.scale(DataCalc.QzIntensity[i]);Tud.scale(DataCalc.QzIntensity[i]);
         Rdu.scale(DataCalc.QzIntensity[i]);Tdu.scale(DataCalc.QzIntensity[i]);
         RuuAux.add(Ruu);TuuAux.add(Tuu);RddAux.add(Rdd);TddAux.add(Tdd);
         RudAux.add(Rud);TudAux.add(Tud);RduAux.add(Rdu);TduAux.add(Tdu);
         end
    else if DataCalc.convolutionCalculation = 'full' then { calcul dans le cas ou on veux faire une convolution complète }
         { ce calcul est valable pour une mesure à longueur d'onde fixe }
         for i:=0 to DataCalc.thetaSpectrum.length-1 do
         begin
           dTheta:=dataCalc.thetaSpectrum.values[i];
           anglesAux.copy(dataCalc.Angles);anglesAux.Offset(dTheta);
           for j:=0 to DataCalc.lambdaSpectrum.length-1 do
           begin
             l:=DataCalc.lambdaSpectrum.values[j];
             { il faut créer un vecteur Qz à partir des valeurs dTheta et l }
             QzSpectrumAux.Sin(AnglesAux);
             QzSpectrumAux.Scale(4*3.1416);QzSpectrumAux.Scale(1/l);
             CalculReflecSF(QzSpectrumAux,z,Vup,Vdo,Vi,phi,Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu);
             Ruu.scale(DataCalc.lambdaIntensity.values[j]);Tuu.scale(DataCalc.lambdaIntensity.values[j]);
             Rdd.scale(DataCalc.lambdaIntensity.values[j]);Tdd.scale(DataCalc.lambdaIntensity.values[j]);
             Rud.scale(DataCalc.lambdaIntensity.values[j]);Tud.scale(DataCalc.lambdaIntensity.values[j]);
             Rdu.scale(DataCalc.lambdaIntensity.values[j]);Tdu.scale(DataCalc.lambdaIntensity.values[j]);
             Ruu.scale(DataCalc.thetaIntensity.values[i]);Tuu.scale(DataCalc.thetaIntensity.values[i]);
             Rdd.scale(DataCalc.thetaIntensity.values[i]);Tdd.scale(DataCalc.thetaIntensity.values[i]);
             Rud.scale(DataCalc.thetaIntensity.values[i]);Tud.scale(DataCalc.thetaIntensity.values[i]);
             Rdu.scale(DataCalc.thetaIntensity.values[i]);Tdu.scale(DataCalc.thetaIntensity.values[i]);
             RuuAux.add(Ruu);TuuAux.add(Tuu);RddAux.add(Rdd);TddAux.add(Tdd);
             RudAux.add(Rud);TudAux.add(Tud);RduAux.add(Rdu);TduAux.add(Tdu);
           end;
         end
         else begin { cas de la convolution fast }
              QzFastMin:=DataCalc.QzFast.min;
              Ruu.size(DataCalc.QzFast);Ruu.setval(0);RuuAux.setval(0);
              Rdd.size(DataCalc.QzFast);Rdd.setval(0);RddAux.setval(0);
              Rud.size(DataCalc.QzFast);Rud.setval(0);RudAux.setval(0);
              Rdu.size(DataCalc.QzFast);Rdu.setval(0);RduAux.setval(0);

              CalculReflecSF(DataCalc.QzFast,z,Vup,Vdo,Vi,phi,Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu); { calcul du spectre suréchantillonné }
              { reconstruction du spectre a partir de ces valeurs suréchantillonnées }
              for i:=0 to dataCalc.Qz.length-1 do
                begin
                  QzCrt:=dataCalc.Qz.values[i];sig:=dataCalc.QzSig.values[i];
                  dQNbr:=trunc(dataCalc.dQzNbr/2);
                  if dQNbr>0 then dq:=DataCalc.QzMult*sig*sqrt(2*ln(2))/dQNbr else dq:=0;
                  { generation de la gaussienne ayant cette largeur }
                  TotalIntensity:=0;
                  Ruuaux1:=0;Rddaux1:=0;Rudaux1:=0;Rduaux1:=0;
                  for j:=-dqNbr to dqNbr do
                  begin { il y a un problème quand qz est négatif }
                      q:=QzCrt+j*dq;{ NOTA abs est un bricolo pour que q soit toujours positif ... }
                      QzInt:=1/(sig*sqrt(2*pi))*exp(-(sqr(j*dq)/(2*sig*sig)));
                      TotalIntensity:=TotalIntensity+QzInt;
                      idxBot:=DataCalc.IdxInQzFastBot[round(abs(q-QzFastMin)*1e-6)];idxTop:=idxBot+1; { NEW offset de la valeur QzFast min ! }
                      qBot:=DataCalc.QzFast.values[IdxBot];qTop:=DataCalc.QzFast.values[IdxTop];
                      RuuBot:=Ruu.values[IdxBot];RuuTop:=Ruu.values[IdxTop];RuuAux1:=RuuAux1+(Ruubot+(RuuTop-RuuBot)*(q-qBot)/(qTop-qBot))*QzInt;
                      RddBot:=Rdd.values[IdxBot];RddTop:=Rdd.values[IdxTop];RddAux1:=RddAux1+(Rddbot+(RddTop-RddBot)*(q-qBot)/(qTop-qBot))*QzInt;
                      RudBot:=Rud.values[IdxBot];RudTop:=Rud.values[IdxTop];RudAux1:=RudAux1+(Rudbot+(RudTop-RudBot)*(q-qBot)/(qTop-qBot))*QzInt;
                      RduBot:=Rdu.values[IdxBot];RduTop:=Rdu.values[IdxTop];RduAux1:=RduAux1+(Rdubot+(RduTop-RduBot)*(q-qBot)/(qTop-qBot))*QzInt;
                  end;
                  Ruuaux1:=Ruuaux1/TotalIntensity;Ruuaux.values[i]:=Ruuaux1; { pas nécessaire en principe si la gaussienne est correctement calculée !! }
                  Rddaux1:=Rddaux1/TotalIntensity;Rddaux.values[i]:=Rddaux1; { pas nécessaire en principe si la gaussienne est correctement calculée !! }
                  Rudaux1:=Rudaux1/TotalIntensity;Rudaux.values[i]:=Rudaux1; { pas nécessaire en principe si la gaussienne est correctement calculée !! }
                  Rduaux1:=Rduaux1/TotalIntensity;Rduaux.values[i]:=Rduaux1; { pas nécessaire en principe si la gaussienne est correctement calculée !! }
                end;
              end;

   Ruu.copy(RuuAux);Tuu.copy(TuuAux);Rdd.copy(RddAux);Tdd.copy(TddAux);
   Rud.copy(RudAux);Tud.copy(TudAux);Rdu.copy(RduAux);Tdu.copy(TduAux);
   freeIt(RuuAux,RddAux,RudAux,RduAux);freeIt(TuuAux,TddAux,TudAux,TduAux);
   FreeIt(QzSpectrumAux,anglesAux);
end;

{ calcul de la réflectivité en polarisé high field (P.Humbert)en tenant compte des différents cas possibles }
{ ACHTUNG: entree implicite "lambda"                                  }
procedure CalculateReflectivitySF_HF(var DataCalc:TDataReflectivity;var z,Vup,Vdo,Vi,BeffMag,BeffTheta,BeffPhi,sigma:Tvec;var Ruu,Rdd,Rud,Rdu,Tuu,Tdd,Tud,Tdu:Tvec);
var i,j:integer;
    Raux,Taux:Tvec;
    RuuAux,TuuAux,RddAux,TddAux,RudAux,TudAux,RduAux,TduAux:Tvec;
    QzSPectrumAux,anglesAux:Tvec;
    dtheta,l:double;
    QzFastMin,QzCrt,sig,q,dq,TotalIntensity,QzInt,qBot,qTop:double;
    dQNbr,idxBot,idxTop:integer;
    Ruuaux1,Rddaux1,Rudaux1,Rduaux1:double;
    RuuBot,RddBot,RudBot,RduBot:double;
    RuuTop,RddTop,RudTop,RduTop:double;
begin
    CreateIt(QzSpectrumAux,anglesAux);
    createIt(RuuAux,RddAux,RudAux,RduAux);createIt(TuuAux,TddAux,TudAux,TduAux);
    RuuAux.size(DataCalc.angles);RuuAux.setval(0);TuuAux.size(DataCalc.angles);TuuAux.setval(0);
    RddAux.size(DataCalc.angles);RddAux.setval(0);TddAux.size(DataCalc.angles);TddAux.setval(0);
    RudAux.size(DataCalc.angles);RudAux.setval(0);TudAux.size(DataCalc.angles);TudAux.setval(0);
    RduAux.size(DataCalc.angles);RduAux.setval(0);TduAux.size(DataCalc.angles);TduAux.setval(0);

    { on a différents types de convolution possibles }
    if DataCalc.convolutionCalculation = 'QzGaussian' { calcul dans l'espace Qz }
    then for i:=0 to DataCalc.QzIntensity.length-1 do
         begin
         { NOTA: pour le moment, on passe le vecteur des angles d'incidences }
         {       ce vecteur est utilisé à une occasion dans la procédure CalculReflecSF_HF_V12 }
         {       cependant, cela n'influe pas sur le calcul de la convolution                   }
         CalculReflecSF_HF(DataCalc.Angles,DataCalc.QzSpectrum[i],z,Vup,Vdo,Vi,sigma,BeffMag,BeffTheta,BeffPhi,Ruu,Rdu,Rud,Rdd);
         Ruu.scale(DataCalc.QzIntensity[i]);Tuu.scale(DataCalc.QzIntensity[i]);
         Rdd.scale(DataCalc.QzIntensity[i]);Tdd.scale(DataCalc.QzIntensity[i]);
         Rud.scale(DataCalc.QzIntensity[i]);Tud.scale(DataCalc.QzIntensity[i]);
         Rdu.scale(DataCalc.QzIntensity[i]);Tdu.scale(DataCalc.QzIntensity[i]);
         RuuAux.add(Ruu);TuuAux.add(Tuu);RddAux.add(Rdd);TddAux.add(Tdd);
         RudAux.add(Rud);TudAux.add(Tud);RduAux.add(Rdu);TduAux.add(Tdu);
         end
    else if DataCalc.convolutionCalculation = 'full' then { calcul dans le cas ou on veux faire une convolution complète }
         { ce calcul est valable pour une mesure à longueur d'onde fixe }
         for i:=0 to DataCalc.thetaSpectrum.length-1 do
         begin
           dTheta:=dataCalc.thetaSpectrum.values[i];
           anglesAux.copy(dataCalc.Angles);anglesAux.Offset(dTheta);
           for j:=0 to DataCalc.lambdaSpectrum.length-1 do
           begin
             l:=DataCalc.lambdaSpectrum.values[j];
             { il faut créer un vecteur Qz à partir des valeurs dTheta et l }
             QzSpectrumAux.Sin(AnglesAux);
             QzSpectrumAux.Scale(4*3.1416);QzSpectrumAux.Scale(1/l);
             {CalculReflecSF_HF_V11(l,AnglesAux,z,Vup,Vdo,sigma,BeffMag,BeffTheta,BeffPhi,Ruu,Rdu,Rud,Rdd); }
             CalculReflecSF_HF(AnglesAux,QzSpectrumAux,z,Vup,Vdo,Vi,sigma,BeffMag,BeffTheta,BeffPhi,Ruu,Rdu,Rud,Rdd);
             Ruu.scale(DataCalc.lambdaIntensity.values[j]);Tuu.scale(DataCalc.lambdaIntensity.values[j]);
             Rdd.scale(DataCalc.lambdaIntensity.values[j]);Tdd.scale(DataCalc.lambdaIntensity.values[j]);
             Rud.scale(DataCalc.lambdaIntensity.values[j]);Tud.scale(DataCalc.lambdaIntensity.values[j]);
             Rdu.scale(DataCalc.lambdaIntensity.values[j]);Tdu.scale(DataCalc.lambdaIntensity.values[j]);
             Ruu.scale(DataCalc.thetaIntensity.values[i]);Tuu.scale(DataCalc.thetaIntensity.values[i]);
             Rdd.scale(DataCalc.thetaIntensity.values[i]);Tdd.scale(DataCalc.thetaIntensity.values[i]);
             Rud.scale(DataCalc.thetaIntensity.values[i]);Tud.scale(DataCalc.thetaIntensity.values[i]);
             Rdu.scale(DataCalc.thetaIntensity.values[i]);Tdu.scale(DataCalc.thetaIntensity.values[i]);
             RuuAux.add(Ruu);TuuAux.add(Tuu);RddAux.add(Rdd);TddAux.add(Tdd);
             RudAux.add(Rud);TudAux.add(Tud);RduAux.add(Rdu);TduAux.add(Tdu);
           end;
         end
         else begin { cas de la convolution fast }
              QzFastMin:=DataCalc.QzFast.min;
              Ruu.size(DataCalc.QzFast);Ruu.setval(0);RuuAux.setval(0);
              Rdd.size(DataCalc.QzFast);Rdd.setval(0);RddAux.setval(0);
              Rud.size(DataCalc.QzFast);Rud.setval(0);RudAux.setval(0);
              Rdu.size(DataCalc.QzFast);Rdu.setval(0);RduAux.setval(0);
              { NOTA: pour le moment, on passe le vecteur des angles d'incidences }
              {       ce vecteur est utilisé à une occasion dans la procédure CalculReflecSF_HF_V12 }
              {       cependant, cela n'influe pas sur le calcul de la convolution                   }
              { IMPRECISION : DataCalc n'est pas compatible avec QzFast !!! }
              CalculReflecSF_HF(DataCalc.Angles,DataCalc.QzFast,z,Vup,Vdo,Vi,sigma,BeffMag,BeffTheta,BeffPhi,Ruu,Rdu,Rud,Rdd);
              { reconstruction du spectre a partir de ces valeurs suréchantillonnées }
              for i:=0 to dataCalc.Qz.length-1 do
                begin
                  QzCrt:=dataCalc.Qz.values[i];sig:=dataCalc.QzSig.values[i];
                  dQNbr:=trunc(dataCalc.dQzNbr/2);
                  if dQNbr>0 then dq:=DataCalc.QzMult*sig*sqrt(2*ln(2))/dQNbr else dq:=0;
                  { generation de la gaussienne ayant cette largeur }
                  TotalIntensity:=0;
                  Ruuaux1:=0;Rddaux1:=0;Rudaux1:=0;Rduaux1:=0;
                  for j:=-dqNbr to dqNbr do
                  begin { il y a un problème quand qz est négatif }
                      q:=QzCrt+j*dq;{ NOTA abs est un bricolo pour que q soit toujours positif ... }
                      QzInt:=1/(sig*sqrt(2*pi))*exp(-(sqr(j*dq)/(2*sig*sig)));
                      TotalIntensity:=TotalIntensity+QzInt;
                      idxBot:=DataCalc.IdxInQzFastBot[round(abs(q-QzFastMin)*1e-6)];idxTop:=idxBot+1; { NEW offset de la valeur QzFast min ! }
                      qBot:=DataCalc.QzFast.values[IdxBot];qTop:=DataCalc.QzFast.values[IdxTop];
                      RuuBot:=Ruu.values[IdxBot];RuuTop:=Ruu.values[IdxTop];RuuAux1:=RuuAux1+(Ruubot+(RuuTop-RuuBot)*(q-qBot)/(qTop-qBot))*QzInt;
                      RddBot:=Rdd.values[IdxBot];RddTop:=Rdd.values[IdxTop];RddAux1:=RddAux1+(Rddbot+(RddTop-RddBot)*(q-qBot)/(qTop-qBot))*QzInt;
                      RudBot:=Rud.values[IdxBot];RudTop:=Rud.values[IdxTop];RudAux1:=RudAux1+(Rudbot+(RudTop-RudBot)*(q-qBot)/(qTop-qBot))*QzInt;
                      RduBot:=Rdu.values[IdxBot];RduTop:=Rdu.values[IdxTop];RduAux1:=RduAux1+(Rdubot+(RduTop-RduBot)*(q-qBot)/(qTop-qBot))*QzInt;
                  end;
                  Ruuaux1:=Ruuaux1/TotalIntensity;Ruuaux.values[i]:=Ruuaux1; { pas nécessaire en principe si la gaussienne est correctement calculée !! }
                  Rddaux1:=Rddaux1/TotalIntensity;Rddaux.values[i]:=Rddaux1; { pas nécessaire en principe si la gaussienne est correctement calculée !! }
                  Rudaux1:=Rudaux1/TotalIntensity;Rudaux.values[i]:=Rudaux1; { pas nécessaire en principe si la gaussienne est correctement calculée !! }
                  Rduaux1:=Rduaux1/TotalIntensity;Rduaux.values[i]:=Rduaux1; { pas nécessaire en principe si la gaussienne est correctement calculée !! }
                end;
              end;

   Ruu.copy(RuuAux);Tuu.copy(TuuAux);Rdd.copy(RddAux);Tdd.copy(TddAux);
   Rud.copy(RudAux);Tud.copy(TudAux);Rdu.copy(RduAux);Tdu.copy(TduAux);
   freeIt(RuuAux,RddAux,RudAux,RduAux);freeIt(TuuAux,TddAux,TudAux,TduAux);
   FreeIt(QzSpectrumAux,anglesAux);
end;

{ procedure générale de calcul de la réflectivité }
procedure CalculateReflectivity(var MyExpPar:TExperimentalParameters;var MyML:TMultiLayer;var MyDC:TDataReflectivity);
begin
   if MyExpPar.polarisationType='non polarised'
   then
      if MyML.MLsystem
      then begin
           MyML.GenerateMLParameters(MyExpPar);
           CalculateReflectivityNP(MyDC,MyML.zML,MyML.VML,MyML.ViML,MyML.sigmaML,MyDC.R,MyDC.T);
           end
      else CalculateReflectivityNP(MyDC,MyML.z,MyML.V,MyML.Vi,MyML.sigma,MyDC.R,MyDC.T);
   if MyExpPar.polarisationType='polarised'
   then
      if MyML.MLsystem
      then begin
           MyML.GenerateMLParameters(MyExpPar);
           CalculateReflectivityP(MyDC,MyML.zML,MyML.VupML,MyML.VdoML,MyML.ViML,MyML.sigmaML,MyDC.Rup,MyDC.Rdo,MyDC.Tup,MyDC.Tdo);
           end
      else CalculateReflectivityP(MyDC,MyML.z,MyML.Vup,MyML.Vdo,MyML.Vi,MyML.sigma,MyDC.Rup,MyDC.Rdo,MyDC.Tup,MyDC.Tdo);
   if MyExpPar.polarisationType='with spin-flip'
   then
      if MyML.MLsystem
      then begin
           MyML.GenerateMLParameters(MyExpPar);
           CalculateReflectivitySF(MyDC,MyML.zML,MyML.VupML,MyML.VdoML,MyML.ViML,MyML.phiML,
                MyDC.Ruu,MyDC.Rdd,MyDC.Rud,MyDC.Rdu,
                MyDC.Tuu,MyDC.Tdd,MyDC.Tud,MyDC.Tdu);
           end
      else CalculateReflectivitySF(MyDC,MyML.z,MyML.Vup,MyML.Vdo,MyML.Vi,MyML.phi,
                MyDC.Ruu,MyDC.Rdd,MyDC.Rud,MyDC.Rdu,
                MyDC.Tuu,MyDC.Tdd,MyDC.Tud,MyDC.Tdu);
   if MyExpPar.polarisationType='with SF high field'
   then
      if MyML.MLsystem
      then begin
           MyML.GenerateMLParameters(MyExpPar);
           CalculateReflectivitySF_HF(MyDC,MyML.zML,MyML.VupML,MyML.VdoML,MyML.ViML,
                MyML.BeffMagML,MyML.BeffThetaML,MyML.BeffPhiML,MyML.sigmaML,
                MyDC.Ruu,MyDC.Rdd,MyDC.Rud,MyDC.Rdu,
                MyDC.Tuu,MyDC.Tdd,MyDC.Tud,MyDC.Tdu);
           end
      else CalculateReflectivitySF_HF(MyDC,MyML.z,MyML.Vup,MyML.Vdo,MyML.Vi,
                MyML.BeffMag,MyML.BeffTheta,MyML.BeffPhi,MyML.sigma,
                MyDC.Ruu,MyDC.Rdd,MyDC.Rud,MyDC.Rdu,
                MyDC.Tuu,MyDC.Tdd,MyDC.Tud,MyDC.Tdu);

end;

{ version 2.7                                                                   }
{ procédure de calcul de la réflectivité (neut. non polarisés)                  }
{ NEW : au lieu d'utilser l'angle d'incidence,, on va rentrer un vecteur de vecteurs de diffusion }
{       c'est en vue des modifications pour le calcul des convolutions expérimentales             }
{ Cette procédure est fonctionnelle quelle que soit la situation NML ou ML      }
{ INPUT: Qz est le vecteur de diffusion 4pi/lambda * sin(theta)                 }
{        z et V, vecteur des épaisseurs et des potentiels des  différentes couches }
{ OUTPUT: un vecteur contenant les vecteurs d'ondes dans les différentes couches}
{         deux vecteurs AR et AT donnant les amplitudes dans les différentes ocuches    }
{         couche i : phi = AT[i] * exp(i*kz[i]*z) + AR[i] * exp(-i*kz[i]*z)         }
procedure CalculReflecNPamplitudes(Qz:double;var z,V,Vi,sigma:Tvec;var kz,AR,AT:Tvec);
var i,n,nn:integer;
    nbrLayers,nbrML:integer; { nombre de couches dans le système }
    nbrPoints:integer; { nombre d'angles d'incidence       }
    zz:double;          { norme du vecteur d'onde incident  }
    q0:double;         { composante perpendiculaire du vecteur d'onde pour les différents angles d'incidence }
    qq0,qq1,zz0,aux,auxC1,auxC2:TCplx;
    rp,rm,exp1,exp2,pp,mm,M1,M2:Tcplx;
    aux1,aux2,aux1DW,aux2DW:TCplx;
    sig:double;
    aux1r,aux2r:double;
    Eq,Eim,AD,ADaux:Tvec;
    DD,M,Maux,MM3:TMtx;
begin
   nbrLayers := z.length;
   createIt(DD,M,Maux,MM3);
   createIt(Eq,Eim);
   DD.size(2,2,true);M.size(2,2,true);Maux.size(2,2,true);MM3.size(2,2,true);
   createit(AD,ADaux);AD.size(2,true);ADaux.size(2,true);
   AR.size(nbrlayers,true);AT.size(nbrlayers,true);
   kz.Size(nbrLayers,true); { les composantes de q peuvent etre complexes }
   Eq.Size(nbrLayers,false);
   Eim.Size(nbrLayers,false);Eim.setVal(0); { a quoi ca sert ???? }

   begin
      M.eye(2,2,true); { initialisation de la récursion }
      q0 := Qz/2; { Qz est le vecteur de diffusion; q0 est la projection de k0 }
      Eq.setVal(hb2over2mn);Eq.scale(q0*q0);Eq.sub(V);
      Eim.setval(0);Eim.sub(Vi); { *** NEW *** prise en compte de la partie complexe du potentiel *** }
      kz.RealToCplx(Eq,Eim);kz.scale(1/hb2over2mn);kz.sqrt;kz.scale(-1); { q peut-etre complexe, racine }
      { on applique le facteur -1 pour des raisons de convention de  l'axe Oz }
      { initialisation de variables pour la récursion }
      qq1:=kz.Cvalues[nbrLayers-1];
      zz:=0;
      n:=nbrLayers-2;
      for n:=nbrLayers-2 downto 0 do   (* boucle du vide vers le substrat *)
      begin
         sig:=sigma.values[n]; { rugosité de l'interface j+1/j }
         qq0:=kz.Cvalues[n];
         zz:=z.values[n]; { epaisseur de la couche considérée }
         zz0:=Cplx(zz,0);      { passage en complexe pour faire les multiplications }
         { détermination de la matrice de transfert }
         { coefficients p et m }
         pp:=Cdiv(Cadd(qq1,qq0),qq1);mm:=Cdiv(Csub(qq1,qq0),qq1);
         { TEMPORAIRE TANT QUE SCALE N'EST PAS REGLE }
         pp:=CMul(Cplx(0.5,0),pp);mm:=CMul(Cplx(0.5,0),mm);
         { calcul des coefficients de rugosité }
         aux1DW :=exp(Cmul(Csqr(Csub(qq0,qq1)),Cplx(-sig*sig/2,0)));
         aux2DW :=exp(Cmul(Csqr(Cadd(qq0,qq1)),Cplx(-sig*sig/2,0)));
         { initialisation of the D matrix}
         aux1 := Exp(Cmul(C_I,Cmul(qq0,zz0)));        { coefficients de la matrice de propagation }
         aux2 := CDiv(C_ONE,aux1);                    { le signe est ok parce que qq0 est negatif ! }
         DD.Cvalues[0,0]:=Cmul(Cmul(pp,aux1),aux1DW);
         DD.Cvalues[0,1]:=Cmul(Cmul(mm,aux2),aux2DW);
         DD.Cvalues[1,0]:=Cmul(Cmul(mm,aux1),aux2DW);
         DD.Cvalues[1,1]:=Cmul(Cmul(pp,aux2),aux1DW);
{         DD.scale(Cplx(0.5,0)); } { CA MARCHE PAS !!!! }
         Maux.Mul(M,DD);M.copy(Maux);
         qq1:=qq0;
      end;
   end;
   { initialisation de la récursion }
   M1:=M.Cvalues[1,0];M2:=M.Cvalues[0,0];
   AR.Cvalues[nbrlayers-1]:=Cdiv(M1,M2);
   AR.Cvalues[0]:=C_ZERO;
   AT.Cvalues[nbrlayers-1]:=C_ONE;      { faiseau incident }
   AT.Cvalues[0]:=Cdiv(C_one,M2);

{ LES VALEURS DE DAR ET DAT SONT CALCULEES pour la partie "basse" de la couche  }
{ N.B. : cela suppose une interface virtuelle pour le substrat                  }
{calcul du vecteur des amplitudes a chaque interface, [0] pour R et [1] pour T  }
{initialisation de la récursion}
   AD.Cvalues[0]:=AT.Cvalues[0]; { les indices 0 et 1 correspondent au vecteur [T,R] }
   AD.Cvalues[1]:=AR.Cvalues[0];
{recursion pour calculer les amplitudes a chaque interface}


   for i:=0 to nbrlayers-2 do {boucle du substrat vers le vide} { ON PEUT ECONOMISER LE DERNIER CALCUL; SERT DE VERIFICATION }
    begin
    {calculation of Fresnel coefficients for each layer}
    {Fresnel reflected and transmitted}
    pp:=Cdiv(Cadd(kz.Cvalues[i+1],kz.Cvalues[i]),Cmul(Cplx(2,0),kz.Cvalues[i+1]));
    mm:=Cdiv(Csub(kz.Cvalues[i+1],kz.Cvalues[i]),Cmul(Cplx(2,0),kz.Cvalues[i+1]));
    {we introduce the roughness by r1 and r2}
    sig:=sigma.values[i];
    rp:=Cmul(pp,exp(Cmul(Csqr(Csub(kz.Cvalues[i+1],kz.Cvalues[i])),Cplx(-sig*sig/2,0))));
    rm:=Cmul(mm,exp(Cmul(Csqr(Cadd(kz.Cvalues[i+1],kz.Cvalues[i])),Cplx(-sig*sig/2,0))));
    {elements of the propagation matrix}
    exp1:=exp(Cmul(C_I,Cmul(kz.Cvalues[i],Cplx(z.values[i],0))));
    exp2:=Cdiv(C_one,exp1);
    {on passe par les matrices pour obtenir DAR[1]}
    MM3.Cvalues[0,0]:=Cmul(rp,exp1);
    MM3.Cvalues[0,1]:=Cmul(rm,exp2);
    MM3.Cvalues[1,0]:=Cmul(rm,exp1);
    MM3.Cvalues[1,1]:=Cmul(rp,exp2);
    {calcul des amplitudes a l'interface}
    ADaux.Cvalues[0]:=Cadd(Cmul(MM3.Cvalues[0,0],AD.Cvalues[0]),Cmul(MM3.Cvalues[0,1],AD.Cvalues[1]));
    ADaux.Cvalues[1]:=Cadd(Cmul(MM3.Cvalues[1,0],AD.Cvalues[0]),Cmul(MM3.Cvalues[1,1],AD.Cvalues[1]));

    AT.Cvalues[i+1]:=ADaux.Cvalues[0];
    AR.Cvalues[i+1]:=ADaux.Cvalues[1];
    AD.copy(ADaux);
    end;

freeIt(DD,M,Maux,MM3);freeIt(Eq,Eim);freeIt(AD,ADaux);

end;


end.


