{ ================================================================================================ }
{ procedures de calcul du diffus }
{ ================================================================================================ }
{ 27.02.02 : FO : creations de 2 nouvelles variables globales : xlimit, xsiz    }
{                 ces variables sont mises a jour dans RefreshParameters        }

unit CalculDiffus;

interface

uses Qdialogs,constants, dataReflec, ndspl, Math387, MathAux, MtxVec, GaussLegendre;

procedure CalculAmpSpeculaire(lambda,theta:double;var z,V,Vi,sigma,xsi,exph:Tvec;var kz,DAR,DAT:Tvec);
procedure CalculNP_Born_NIPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var Born:double);
procedure CalculNP_Born_SA_NIPC(typeCorr:integer;lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var Born:double);
procedure CalculNP_Born_IPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var Born:double);
procedure CalculNP_Born_SA_IPC(typeCorr:integer;lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var Born:double);
procedure CalculNP_DWBA_II_NIPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var IDWBA:double);
procedure CalculNP_DWBA_II_IPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var IDWBA:double);
procedure CalculNP_DWBA_III_NIPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var IDWBA:double);
procedure CalculNP_DWBA_III_SA_NIPC(typeCorr:integer;lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var IDWBA:double);
procedure CalculNP_DWBA_III_IPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var IDWBA:double);

var
   gaussOrder:integer;

implementation

{ ============================================================================= }
{ Calculation of Diffuse Amplitudes for the Reflection and Transmission         }
{ INPUT : theta en °                                                            }
{         lambda en nm                                                          }
{ IMPLICITE : zdiff, Vdiff, sigmadiff, xsidiff, hdiff:Tvec convention Stettner  }
{ OUTPUT: kz, vecteurs d'ondes dans chaque couche (indexation Stettner          }
{         signe de kz ?                                                         }
{        DAR et DAT, coefficients de Fresnel complexes a toutes les interfaces  }
{ MAJ 20.08.01 : sig:=sigma.values[i]; au lieur de sig:=sigma.values[i-1];      }
{                dans le calcul des DAR et DAT intermediaires, correction des indices aussi }
{                on ne calculait pas les bons DAT et DAR                        }
{                Plus inversion du couple (DAR;DAT) en (DAT;DAR) voir formalisme }
{ ============================================================================= }
procedure CalculAmpSpeculaire(lambda,theta:double;var z,V,Vi,sigma,xsi,exph:Tvec;var kz,DAR,DAT:Tvec);
var
   i,j,nbrlayers:integer;
   k0z,sig:double;
   M1,M2:Tcplx;
   kzre,kzim:Tvec;
   MI1,MM1,MM2,MI2,MM3:TMtx;
   pp,mm,rp,rm,exp1,exp2:Tcplx;
   AD,ADaux:Tvec; { vecteur contenant les amplitudes complxes a une interface; variable intermediaire }

begin
{Vdiff est un reel, mais on veut kz complexe, on cree donc ces 2 vecteurs (partie reelle et imaginaire) }
createit(kzre,kzim);
createit(MI1,MM1);createit(MM2);createit(MI2,MM3);
createit(AD,ADaux);AD.size(2,true);ADaux.size(2,true);
nbrlayers:=z.length;  { valeur de "N" }
kzre.Size(nbrLayers,false);kzim.Size(nbrLayers,false);kz.size(nbrlayers,true);
DAR.size(nbrlayers,true);DAT.size(nbrlayers,true);

k0z:=2*Pi/(lambda)*1e9*sin(theta*Deg);
kzre.setval(hb2over2mn);kzre.scale(k0z*k0z);kzre.sub(V);
kzim.setval(0);kzim.sub(Vi);
kz.RealToCplx(kzre,kzim);kz.scale(1/hb2over2mn);kz.sqrt;kz.scale(-1);

MI1.size(2,2,true);MM1.size(2,2,true);MM2.size(2,2,true);
MI1.eye(2,2,true); { matrice identite(eye) pour initialisation de la récursion }
for i:=nbrlayers-2 downto 0 do { boucle du vide vers le substrat }
    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]}
    MM1.Cvalues[0,0]:=Cmul(rp,exp1);
    MM1.Cvalues[0,1]:=Cmul(rm,exp2);
    MM1.Cvalues[1,0]:=Cmul(rm,exp1);
    MM1.Cvalues[1,1]:=Cmul(rp,exp2);
    MM2.Mul(MI1,MM1);MI1.copy(MM2);
    end;

{ initialisation de la récursion }
M1:=MI1.Cvalues[1,0];M2:=MI1.Cvalues[0,0];
DAR.Cvalues[nbrlayers-1]:=Cdiv(M1,M2);
DAR.Cvalues[0]:=C_ZERO;
DAT.Cvalues[nbrlayers-1]:=C_ONE;      { faiseau incident }
DAT.Cvalues[0]:=Cdiv(C_one,M2);
{ LES VALEURS DE DAR ET DAT SONT CALCULEES pour la partie "basse" de la couche }
{calcul du vecteur des amplitudes a chaque interface, [0] pour R et [1] pour T }
{initialisation de la récursion}
AD.Cvalues[0]:=DAT.Cvalues[0];
AD.Cvalues[1]:=DAR.Cvalues[0];
{recursion pour calculer les amplitudes a chaque interface}
MM3.size(2,2,true);

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]));
    DAT.Cvalues[i+1]:=ADaux.Cvalues[0];
    DAR.Cvalues[i+1]:=ADaux.Cvalues[1];
    AD.copy(ADaux);
    end;

freeit(kzre,kzim);freeit(MI1,MM1,MM2);freeit(AD,ADaux);freeit(MI2,MM3);
end;

{ version 1 : pas de corrélations entre les couches                    }
{ NIPC no interplane correlations                                      }
{ calcul dans l'approximation de Born                                  }
{ INPUT : lambda (nm); thetai, thetaf (deg)                            }
{ OUTPUT : Born, l'intensité diffusée (complexe)  dsigma/domega        }
{          on laisse les corrections pour la procedure appelante       }
{          pour le moement on laisse aussi tomber le facteur G         }
{ REM : il faut thetain <> thetaout sinon singularité                  }
{ on applique la formule de stettner 2.70 telle quelle (avec k0^4)     }
procedure CalculNP_Born_NIPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var Born:double);
var i,j,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    k0,deltatheta:double;
    k,kr,ki,kz,DAR,DAT:Tvec;
    qx,qz,qz2,sqrk,qqr,qsig,SFRe,SFR:double;
    intensite:double;
    sigvali,sigvali2,xsivali,hvali:double;
    nlimit:integer;
    xint,xi,xf,xlimit:double;
Begin
createit(k,kr,ki);
nbrLayers := z.length;
k.size(nbrlayers,true);kr.size(nbrlayers,false);ki.size(nbrlayers,false);

qx:=2*pi/(lambda)*1e9*(cos(thetai*Deg)-cos(thetaf*Deg)); { signe etait faux, si on prend (Ox) suivant le faisceau incident }
qz:=2*Pi/(lambda)*1e9*(sin(thetai*Deg)+sin(thetaf*Deg));qz2:=qz*qz;
{ verifier si c'est juste }
k0:=2*Pi/(lambda)*1e9;
kr.setval(hb2over2mn);kr.scale(k0*k0);kr.sub(V);
ki.setval(0);ki.sub(Vi);
k.realtocplx(kr,ki);k.scale(1/hb2over2mn);k.sqrt; { pourquoi pas multiplie par -1?? }

{ on se balance du signe, c'est uniquement le module qui nous intéresse }

deltatheta:=abs(thetai-thetaf); { en radians ! }
intensite:=0; { initialisation de la somme de diffus }
if deltaTheta<0.001 { spéculaire }
then begin { calcul du spéculaire }
     createit(kz,DAR,DAT);
     kz.size(nbrlayers,true);DAR.size(nbrlayers,true);DAT.size(nbrlayers,true);
     CalculAmpSpeculaire(lambda,thetai,z,V,Vi,sigma,xsi,exph,kz,DAR,DAT);
     Born:=sqr(Cabs(DAR.Cvalues[nbrlayers-1]));
     FreeIt(kz,DAR,DAT);
     end
else begin
     for i:=1 to nbrLayers-1 do { sommes croisées de la première couche (excepte substrat) au vide }
      begin
          sigvali:=sigma.values[i-1];xsivali:=xsi.values[i-1];hvali:=exph.values[i-1];
          sigvali2:=sigvali*sigvali;
          {prefacteur au facteur de structure}
          {sqrk: difference des carres des k; remplace K0^4 (n-n)(n-n)*}
          sqrk:=Sqr(Cabs(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])))); { indices ok }
          qsig:=exp(-qz*qz*sigvali2);
          qqr:=4*pi/qz/qz;
          { calcul de S(Q)  facteur de structure }
              xlimit:=xsivali*Power(-ln(critereIntegral),1/(2*hvali));
              SFRE:=0;
              xint:=abs(2*pi/qx); { on va integrer par périodes pour avoir plus de precision }
              if xint<(xlimit) then begin xi:=0;xf:=xint;nlimit:=floor(xlimit/xint)+1; end
                               else begin xi:=0;xf:=xlimit;nlimit:=1;end;
              for n:=1 to nlimit do begin
                                    GaussLegendreQuadrature0(xi,xf,gaussOrder,integrandSQ,sigvali2,xsivali,hvali,qx,qz2,0,SFR);
                                    SFRe:=SFRe+SFR;
                                    xi:=n*xint;xf:=(n+1)*xint;
                                    end;
          intensite:=intensite+sqrk*qsig*qqr*SFRe;
      end; { fin des boucle i }
      Born:=intensite/16/pi2;
      end;
   FreeIt(k,kr,ki);
end;

{ version 1 : pas de corrélations entre les couches                    }
{ NIPC no interplane correlations                                      }
{ calcul dans l'approximation de Born et simple Fourier transform      }
{ de la fonction de corrélation                                        }
{ dérivé de la procedure ci-dessus dans laquelle l'intégration est analytique }
{ INPUT : lambda (nm); thetai, thetaf (deg)                            }
{         typeCorr, type de la fonction de corrélation                 }
{ OUTPUT : Born, l'intensité diffusée (complexe)  dsigma/domega        }
{          on laisse les corrections pour la procedure appelante       }
{          pour le moement on laisse aussi tomber le facteur G         }
{ REM : il faut thetain <> thetaout sinon singularité                  }
{ on applique la formule de stettner 2.70 telle quelle (avec k0^4)     }
procedure CalculNP_Born_SA_NIPC(typeCorr:integer;lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var Born:double);
var i,j,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    k0,deltatheta:double;
    k,kr,ki,kz,DAR,DAT:Tvec;
    qx,qz,qz2,sqrk,qqr,qsig,SFRe,SFR:double;
    intensite:double;
    sigvali,sigvali2,xsivali,hvali:double;
    nlimit:integer;
    xint,xi,xf,xlimit:double;
Begin
createit(k,kr,ki);
nbrLayers := z.length;
k.size(nbrlayers,true);kr.size(nbrlayers,false);ki.size(nbrlayers,false);

qx:=2*pi/(lambda)*1e9*(cos(thetai*Deg)-cos(thetaf*Deg)); { signe etait faux, si on prend (Ox) suivant le faisceau incident }
qz:=2*Pi/(lambda)*1e9*(sin(thetai*Deg)+sin(thetaf*Deg));qz2:=qz*qz;
{ verifier si c'est juste }
k0:=2*Pi/(lambda)*1e9;
kr.setval(hb2over2mn);kr.scale(k0*k0);kr.sub(V);
ki.setval(0);ki.sub(Vi);
k.realtocplx(kr,ki);k.scale(1/hb2over2mn);k.sqrt; { pourquoi pas multiplie par -1?? }
{ on se balance du signe, c'est uniquement le module qui nous intéresse }

deltatheta:=abs(thetai-thetaf); { en radians ! }
intensite:=0; { initialisation de la somme de diffus }
if deltaTheta<0.001 { spéculaire }
then begin { calcul du spéculaire }
     createit(kz,DAR,DAT);
     kz.size(nbrlayers,true);DAR.size(nbrlayers,true);DAT.size(nbrlayers,true);
     CalculAmpSpeculaire(lambda,thetai,z,V,Vi,sigma,xsi,exph,kz,DAR,DAT);
     Born:=sqr(Cabs(DAR.Cvalues[nbrlayers-1]));
     FreeIt(kz,DAR,DAT);
     end
else begin
     for i:=1 to nbrLayers-1 do { sommes croisées de la première couche (excepte substrat) au vide }
      begin
          sigvali:=sigma.values[i-1];xsivali:=xsi.values[i-1];hvali:=exph.values[i-1];
          sigvali2:=sigvali*sigvali;
          {prefacteur au facteur de structure}
          {sqrk: difference des carres des k; remplace K0^4 (n-n)(n-n)*}
          sqrk:=Sqr(Cabs(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])))); { indices ok }
          qsig:=exp(-qz*qz*sigvali2);
          qqr:=2*pi; { LE FACTEUR 2 DISPARAIT : IL PASSE DANS LA TF }
          { calcul de S(Q)  facteur de structure }
          case typeCorr of
          0:SFRe:=3.14;
          1:SFRe:=sqrt(pi)*xsivali*sigvali2*exp(-qx*qx/4*xsivali*xsivali); { cas gaussian }
          2:SFRe:=2*sigvali2*xsivali/(xsivali*xsivali*qx*qx-1);
          3:SFRe:=3.14;
          end;
          intensite:=intensite+sqrk*qsig*qqr*SFRe;
      end; { fin des boucle i }
      Born:=intensite/16/pi2;
      end;
   FreeIt(k,kr,ki);
end;

{ version 1 : corrélations entre les couches                    }
{ LIMITATION : on suppose que les fonctions de correlations sont les memes dans chaque couche !!! }
{ ca simplifie le calcul des fonctions de cross correlations :             }
{ il suffit d'appliquer le facteur de damping exp(-|zi-zj|/xsi)            }
{ NIPC no interplane correlations                                      }
{ calcul dans l'approximation de Born                                  }
{ INPUT : lambda (nm); thetai, thetaf (deg)                            }
{ OUTPUT : Born, l'intensité diffusée (complexe)  dsigma/domega        }
{          on laisse les corrections pour la procedure appelante       }
{          pour le moement on laisse aussi tomber le facteur G         }
{ REM : il faut thetain <> thetaout sinon singularité                  }
{ on applique la formule de stettner 2.70 telle quelle (avec k0^4)     }
{ legere optimisation possible si on calcule explicitement les termes Cii }
procedure CalculNP_Born_IPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var Born:double);
var i,j,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    k0,deltatheta:double;
    k,kr,ki,kz,DAR,DAT:Tvec;
    qx,qz,qz2,qqr,qsig,SFRe,SFR:double;
    sqrkR:double;
    sqrk,qphase:Tcplx;
    intensite:double;
    zzi,zzj,sigvali,sigvali2,xsivali,hvali,sigvalj,sigvalj2,xsivalj,hvalj:double;
    nlimit:integer;
    xint,xi,xf,xlimit:double;
    zz:array[0..nbrLayersMax] of double; { tableau des cotes des interfaces }
    xsiIPC:double;
Begin
xsiIPC:=xsiZ; { passage en SI }
createit(k,kr,ki);
nbrLayers := z.length;
k.size(nbrlayers,true);kr.size(nbrlayers,false);ki.size(nbrlayers,false);

{ creation d'un tableau des profondeurs }
zz[nbrLayers-1]:=0; { tableua des cotes des interfaces } { la surface est mise a la cote 0 }
for i:=nbrLayers-2 downto 0 do zz[i]:=zz[i+1]-z.values[i];

qx:=2*pi/(lambda)*1e9*(cos(thetai*Deg)-cos(thetaf*Deg)); { signe etait faux, si on prend (Ox) suivant le faisceau incident }
qz:=2*Pi/(lambda)*1e9*(sin(thetai*Deg)+sin(thetaf*Deg));qz2:=qz*qz;
{ verifier si c'est juste }
k0:=2*Pi/(lambda)*1e9;
kr.setval(hb2over2mn);kr.scale(k0*k0);kr.sub(V);
ki.setval(0);ki.sub(Vi);
k.realtocplx(kr,ki);k.scale(1/hb2over2mn);k.sqrt; { pourquoi pas multiplie par -1?? }
{ on se balance du signe, c'est uniquement le module qui nous intéresse }

deltatheta:=abs(thetai-thetaf); { en radians ! }
intensite:=0; { initialisation de la somme de diffus }
if deltaTheta<0.001 { spéculaire }
then begin { calcul du spéculaire }
     createit(kz,DAR,DAT);
     kz.size(nbrlayers,true);DAR.size(nbrlayers,true);DAT.size(nbrlayers,true);
     CalculAmpSpeculaire(lambda,thetai,z,V,Vi,sigma,xsi,exph,kz,DAR,DAT);
     Born:=sqr(Cabs(DAR.Cvalues[nbrlayers-1]));
     FreeIt(kz,DAR,DAT);
     end
else begin
     for i:=1 to nbrLayers-1 do { sommes croisées des termes i=j }
     begin
          zzi:=zz[i]; { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          sigvali:=sigma.values[i-1];sigvali2:=sigvali*sigvali; { indices OK }
          xsivali:=xsi.values[i-1]; { indices OK }
          hvali:=exph.values[i-1]; { indices OK }
          {prefacteur au facteur de structure}
          {sqrk: difference des carres des k; remplace K0^4 (n-n)(n-n)*}
          sqrkR:=sqr(Cabs(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])))); { indices ok }
          { facteur de phase complexe }
          qsig:=exp(-qz*qz*sigvali2);
          qqr:=4*pi/qz/qz;
          { calcul de S(Q)  facteur de structure }
              xlimit:=xsivali*Power(-ln(critereIntegral),1/(2*hvali));
              SFRe:=0;
              xint:=abs(2*pi/qx); { on va integrer par périodes pour avoir plus de precision }
              if xint<(xlimit) then begin xi:=0;xf:=xint;nlimit:=floor(xlimit/xint)+1; end
                               else begin xi:=0;xf:=xlimit;nlimit:=1;end;
              for n:=1 to nlimit do begin
                                    GaussLegendreQuadrature0(xi,xf,gaussOrder,integrandSQ,sigvali2,xsivali,hvali,qx,qz2,0,SFR);
                                    SFRe:=SFRe+SFR;
                                    xi:=n*xint;xf:=(n+1)*xint;
                                    end;
          intensite:=intensite+sqrkR*qsig*qqr*SFRe;
     end;
     for i:=1 to nbrLayers-1 do { sommes croisées de la première couche (excepte substrat) au vide }
       for j:=1 to i-1 do       { mais avec i<>j }
       begin
          zzi:=zz[i];zzj:=zz[j]; { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          sigvali:=sigma.values[i-1];sigvali2:=sigvali*sigvali;sigvalj:=sigma.values[j-1];sigvalj2:=sigvalj*sigvalj; { indices OK }
          xsivali:=xsi.values[i-1];xsivalj:=xsi.values[j-1]; { indices OK }
          hvali:=exph.values[i-1];hvalj:=exph.values[j-1]; { indices OK }

          {prefacteur au facteur de structure}
          {sqrk: difference des carres des k; remplace K0^4 (n-n)(n-n)*}
          sqrk:=Cmul(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])),Conj(Csub(Csqr(k.Cvalues[j]),Csqr(k.Cvalues[j-1])))); { indices ok }
          { facteur de phase complexe }
          qphase:=exp(Cplx(0,-qz*(zz[i]-zz[j])));
          qsig:=exp(-0.5*qz*qz*(sigvali2+sigvalj2));
          qqr:=4*pi/qz/qz;
          { calcul de S(Q)  facteur de structure }
              xlimit:=xsivali*Power(-ln(critereIntegral),1/(2*hvali));
              SFRe:=0;
              xint:=abs(2*pi/qx); { on va integrer par périodes pour avoir plus de precision }
              if xint<(xlimit) then begin xi:=0;xf:=xint;nlimit:=floor(xlimit/xint)+1; end
                               else begin xi:=0;xf:=xlimit;nlimit:=1;end;
              for n:=1 to nlimit do begin
                                    GaussLegendreQuadrature(xi,xf,gaussOrder,integrandSQ_IPC,sigvali2,xsivali,hvali,qx,qz2,zz[i]-zz[j],xsiIPC,0,0,0,0,0,0,SFR);
                                    SFRe:=SFRe+SFR;
                                    xi:=n*xint;xf:=(n+1)*xint;
                                    end;
          intensite:=intensite+2*RealPart(Cmul(sqrk,qphase))*qsig*qqr*SFRe;
      end; { fin des boucle i }
      Born:=intensite/16/pi2;
      end;
   FreeIt(k,kr,ki);
end;

{ version 1 : corrélations entre les couches                    }
{ calcul dans l'approximation de Born et simple Fourier transform      }
{ de la fonction de corrélation                                        }
{ dérivé de la procedure ci-dessus dans laquelle l'intégration est analytique }
{ LIMITATION : on suppose que les fonctions de correlations sont les memes dans chaque couche !!! }
{ ca simplifie le calcul des fonctions de cross correlations :             }
{ il suffit d'appliquer le facteur de damping exp(-|zi-zj|/xsi)            }
{ NIPC no interplane correlations                                      }
{ calcul dans l'approximation de Born                                  }
{ INPUT : lambda (nm); thetai, thetaf (deg)                            }
{         typeCorr, type de la fonction de corrélation                 }
{ OUTPUT : Born, l'intensité diffusée (complexe)  dsigma/domega        }
{          on laisse les corrections pour la procedure appelante       }
{          pour le moement on laisse aussi tomber le facteur G         }
{ REM : il faut thetain <> thetaout sinon singularité                  }
{ on applique la formule de stettner 2.70 telle quelle (avec k0^4)     }
{ legere optimisation possible si on calcule explicitement les termes Cii }
procedure CalculNP_Born_SA_IPC(typeCorr:integer;lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var Born:double);
var i,j,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    k0,deltatheta:double;
    k,kr,ki,kz,DAR,DAT:Tvec;
    qx,qz,qz2,qqr,qsig,SFRe,SFR:double;
    sqrkR:double;
    sqrk,qphase:Tcplx;
    intensite:double;
    zzi,zzj,sigvali,sigvali2,xsivali,hvali,sigvalj,sigvalj2,xsivalj,hvalj:double;
    nlimit:integer;
    xint,xi,xf,xlimit:double;
    zz:array[0..nbrLayersMax] of double; { tableau des cotes des interfaces }
    xsiIPC:double;
Begin
xsiIPC:=xsiz; { passage en SI }
createit(k,kr,ki);
nbrLayers := z.length;
k.size(nbrlayers,true);kr.size(nbrlayers,false);ki.size(nbrlayers,false);

{ creation d'un tableau des profondeurs }
zz[nbrLayers-1]:=0; { tableua des cotes des interfaces } { la surface est mise a la cote 0 }
for i:=nbrLayers-2 downto 0 do zz[i]:=zz[i+1]-z.values[i];

qx:=2*pi/(lambda)*1e9*(cos(thetai*Deg)-cos(thetaf*Deg)); { signe etait faux, si on prend (Ox) suivant le faisceau incident }
qz:=2*Pi/(lambda)*1e9*(sin(thetai*Deg)+sin(thetaf*Deg));qz2:=qz*qz;
{ verifier si c'est juste }
k0:=2*Pi/(lambda)*1e9;
kr.setval(hb2over2mn);kr.scale(k0*k0);kr.sub(V);
ki.setval(0);ki.sub(Vi);
k.realtocplx(kr,ki);k.scale(1/hb2over2mn);k.sqrt; { pourquoi pas multiplie par -1?? }
{ on se balance du signe, c'est uniquement le module qui nous intéresse }

deltatheta:=abs(thetai-thetaf); { en radians ! }
intensite:=0; { initialisation de la somme de diffus }
if deltaTheta<0.001 { spéculaire }
then begin { calcul du spéculaire }
     createit(kz,DAR,DAT);
     kz.size(nbrlayers,true);DAR.size(nbrlayers,true);DAT.size(nbrlayers,true);
     CalculAmpSpeculaire(lambda,thetai,z,V,Vi,sigma,xsi,exph,kz,DAR,DAT);
     Born:=sqr(Cabs(DAR.Cvalues[nbrlayers-1]));
     FreeIt(kz,DAR,DAT);
     end
else begin
     for i:=1 to nbrLayers-1 do { sommes croisées des termes i=j }
     begin
          zzi:=zz[i]; { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          sigvali:=sigma.values[i-1];sigvali2:=sigvali*sigvali; { indices OK }
          xsivali:=xsi.values[i-1]; { indices OK }
          hvali:=exph.values[i-1]; { indices OK }
          {prefacteur au facteur de structure}
          {sqrk: difference des carres des k; remplace K0^4 (n-n)(n-n)*}
          sqrkR:=sqr(Cabs(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])))); { indices ok }
          { facteur de phase complexe }
          qsig:=exp(-qz*qz*sigvali2);
          qqr:=2*pi;
          { calcul de S(Q)  facteur de structure }
          case typeCorr of
          0:SFRe:=3.14;
          1:SFRe:=sqrt(pi)*xsivali*sigvali2*exp(-qx*qx/4*xsivali*xsivali); { cas gaussian }
          2:SFRe:=2*sigvali2*xsivali/(xsivali*xsivali*qx*qx-1);
          3:SFRe:=3.14;
          end;
          intensite:=intensite+sqrkR*qsig*qqr*SFRe;
     end;
     for i:=1 to nbrLayers-1 do { sommes croisées de la première couche (excepte substrat) au vide }
       for j:=1 to i-1 do       { mais avec i<>j }
       begin
          zzi:=zz[i];zzj:=zz[j]; { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          sigvali:=sigma.values[i-1];sigvali2:=sigvali*sigvali;sigvalj:=sigma.values[j-1];sigvalj2:=sigvalj*sigvalj; { indices OK }
          xsivali:=xsi.values[i-1];xsivalj:=xsi.values[j-1]; { indices OK }
          hvali:=exph.values[i-1];hvalj:=exph.values[j-1]; { indices OK }

          {prefacteur au facteur de structure}
          {sqrk: difference des carres des k; remplace K0^4 (n-n)(n-n)*}
          sqrk:=Cmul(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])),Conj(Csub(Csqr(k.Cvalues[j]),Csqr(k.Cvalues[j-1])))); { indices ok }
          { facteur de phase complexe }
          qphase:=exp(Cplx(0,-qz*(zz[i]-zz[j])));
          qsig:=exp(-0.5*qz*qz*(sigvali2+sigvalj2));
          qqr:=2*pi;
          { calcul de S(Q)  facteur de structure }
          case typeCorr of
          0:SFRe:=3.14;
          1:SFRe:=sqrt(pi)*xsivali*sigvali2*exp(-qx*qx/4*xsivali*xsivali); { cas gaussian }
          2:SFRe:=2*sigvali2*xsivali/(xsivali*xsivali*qx*qx-1);
          3:SFRe:=3.14;
          end;
          intensite:=intensite+2*RealPart(Cmul(sqrk,qphase))*qsig*qqr*SFRe;
      end; { fin des boucle i }
      Born:=intensite/16/pi2;
      end;
   FreeIt(k,kr,ki);
end;

{ =============================================================================}
{ calcul de l'intensité diffusée en DWBA                                       }
{ ref Stettner p.8, formalisme II                                              }
{ Cette version se limite a NIPC                                               }
{ INPUT:  lambdadiff (en nm)                                                   }
{         theta in et theta final en °                                         }
{         il faut thetain <> thetaf                                            }
{ OUTPUT: IDWBA l'intensité diffusée (partie réelle ??)                        }
{ N.B.: explication des variables (a comparer avec Stettner p.25)              }
{ MAJ 20.08.01 : inversion du signe de qx                                      }
{                correctiopns d'indices multiples                              }
{ on suit le formalisme II de Stettner, cela permet de ne pas avoir a distingeur entre le cas dioptre et pas dioptre }
{ les DAR et les DAT sont definis correctement (bottom de la couche)           }
{ QUID des signes des k ..., on peut reprendre Stettner, il a les memes conventions de (Oz)    }
{ on prend le calcul poussé au bout pour acceler le process de calcul (voir classeur }
{ on gagne un facteur 4 en ne calculant pas systématiquement les 16 inegrales complexes }
{ DOUTE SUR LES SIGNES DES Z !!   cotes ou epaisseurs des couches ??           }
{ =============================================================================}
procedure CalculNP_DWBA_II_NIPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var IDWBA:double);
var i,j,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    zzi,zzj:TCplx;
    DARi,DATi,DARf,DATf:Tvec;
    k,kr,ki,kzi,kzf:Tvec;
    Gmi,Gmbi,qmi:array [0..3] of array [0..nbrlayersmax] of Tcplx;
    k0,qx,qq2r,qq2i,SFRe,SFIm,integ,integR,integI:double;
    qsig,qqr,qz2,sqrk,intensite,sommeMN:double;
    qsigC,qz2C,qqrC:Tcplx;
    xlimit:double;
    sigvali,sigvali2,xsivali,hvali,sigvalj,xsivalj,hvalj,zi,zj:double;
    ti,tf,a1,deltatheta:double;
    nlimit:integer;
    xtot,xint,xi,xf,SFR,SFI:double;
    dthetai,dthetaf,res:double;
    zz:array[0..nbrLayersMax] of double; { tableau des cotes des interfaces }

{ la perte de temps induite par cet appel est non mesurable within 5% }
function GLintegration(F:TGLFunction0;qx,qz2R,qz2I:double):double;
var n:integer;
begin
    SFRE:=0;xint:=abs(2*pi/qx);
    if xint<(xlimit) then begin xi:=0;xf:=xint;nlimit:=floor(xlimit/xint)+1; end
                     else begin xi:=0;xf:=xlimit;nlimit:=1; end;
    for n:=1 to nlimit do begin GaussLegendreQuadrature0(xi,xf,gaussOrder,F,sigvali2,xsivali,hvali,qx,qz2R,qz2I,SFR);
                                SFRe:=SFRe+SFR;xi:=n*xint;xf:=(n+1)*xint; end;
    result:=SFRe;
end;

begin
createit(k,kr,ki,kzi);createit(kzf);createit(DARi,DARf,DATi,DATf);
nbrLayers := z.length;
k.size(nbrlayers,true);kr.size(nbrlayers,false);ki.size(nbrlayers,false);kzi.size(nbrlayers,true);kzf.size(nbrlayers,true);
DARi.size(nbrlayers,true);DATi.size(nbrlayers,true);DARf.size(nbrlayers,true);DATf.size(nbrlayers,true);

{ calcul des vecteurs d'ondes dans chaque couche -> kzin kzfinal }
{ avec nos conventions kzi est negatif et kzf est negatif (au-dessus de l'angle critique) }
{ calcul des coefficients de Fresnel complexes a chaque interface du système, dans la partie basse de la couche }
CalculAmpSpeculaire(lambda,thetai,z,V,Vi,sigma,xsi,exph,kzi,DARi,DATi);CalculAmpSpeculaire(lambda,thetaf,z,V,Vi,sigma,xsi,exph,kzf,DARf,DATf);
qx:=2*pi/(lambda)*1e9*(cos(thetai*Deg)-cos(thetaf*Deg)); { signe etait faux, si on prend (Ox) suivant le faisceau incident }
{calcul des vecteurs k pour toutes les couches}
k0:=2*Pi/(lambda)*1e9;
kr.setval(hb2over2mn);kr.scale(k0*k0);kr.sub(V);
ki.setval(0);ki.sub(Vi);
k.realtocplx(kr,ki);k.scale(1/hb2over2mn);k.sqrt; { pourquoi pas multiplie par -1 ??? }

deltatheta:=abs(thetai-thetaf); { en radians ! }
{ creation d'un tableau des profondeurs }
zz[nbrLayers-1]:=0; { tableua des cotes des interfaces }
for i:=nbrLayers-2 downto 0 do zz[i]:=zz[i+1]-z.values[i];

intensite:=0; { initilaisation de la somme de diffus }
if deltatheta<0.001     { limite en radians ?! }
then IDWBA:=sqr(Cabs(DARi.Cvalues[nbrlayers-1])) { on met lm'intensité spéculaire si on est trop proche du spéculaire , sinon point singulier }
else begin { sommes croisées sur les différentes couches de 0 à N }
       for i:=1 to nbrLayers-1 do { sommes croisées de la première couche (excepte substrat) au vide }
          begin
          qmi[0,i]:=Cadd(kzi.Cvalues[i-1],kzf.Cvalues[i-1]);qmi[1,i]:=Csub(kzi.Cvalues[i-1],kzf.Cvalues[i-1]);qmi[2,i]:=neg(qmi[1,i]);qmi[3,i]:=neg(qmi[0,i]);
          { facteur de pahse pour les coefficients G ~ }
          { on definit le zero au niveau de la surface de l'échantillon : interface vide/première couche <=> onde incidente a une phase nulle au niveau de la surface }
          { "profondeur" dans l'echantillon }
          zzi:=Cplx(zz[i],0);  { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          Gmi[0,i]:=Cmul(DATi.Cvalues[i-1],DATf.Cvalues[i-1]); Gmbi[0,i]:=Cmul(Gmi[0,i],exp(Cdiv(Cmul(qmi[0,i],zzi),C_I)));
          Gmi[1,i]:=Cmul(DATi.Cvalues[i-1],DARf.Cvalues[i-1]); Gmbi[1,i]:=Cmul(Gmi[1,i],exp(Cdiv(Cmul(qmi[1,i],zzi),C_I)));
          Gmi[2,i]:=Cmul(DARi.Cvalues[i-1],DATf.Cvalues[i-1]); Gmbi[2,i]:=Cmul(Gmi[2,i],exp(Cdiv(Cmul(qmi[2,i],zzi),C_I)));
          Gmi[3,i]:=Cmul(DARi.Cvalues[i-1],DARf.Cvalues[i-1]); Gmbi[3,i]:=Cmul(Gmi[3,i],exp(Cdiv(Cmul(qmi[3,i],zzi),C_I)));
          sigvali:=sigma.values[i-1];xsivali:=xsi.values[i-1];hvali:=exph.values[i-1];
          sigvali2:=sigvali*sigvali;

          sqrk:=Sqr(Cabs(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])))); { valeur reelle ! }
          SommeMN:=0;
          xlimit:=xsivali*Power(-ln(critereIntegral),1/(2*hvali)); { EVITER L'APPEL A UN CONTROEL GRAPHIQUE .value !!! }
          { sommation sur les indices m et n                              }
          { on fait la somme de manière explicite pour booster le calcul  }
          { Q_00, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[0,i].re)-sqr(qmi[0,i].im)));
          qz2:=sqr(Cabs(qmi[0,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2;
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*(Sqr(Cabs(Gmbi[0,i]))+Sqr(Cabs(Gmbi[3,i])));
          { Q_11, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[1,i].re)-sqr(qmi[1,i].im)));
          qz2:=sqr(Cabs(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2;
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*(Sqr(Cabs(Gmbi[1,i]))+Sqr(Cabs(Gmbi[2,i])));
          { Q_03, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[0,i].re)-sqr(qmi[0,i].im)));
          qz2:=-sqr(Cabs(qmi[0,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2; { attention jonglage avec les signes, voir calculs papier }
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*2*RealPart(Cmul(Gmbi[0,i],conj(Gmbi[3,i])));
          { Q_12, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[1,i].re)-sqr(qmi[1,i].im)));
          qz2:=-sqr(Cabs(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2;
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*2*RealPart(Cmul(Gmbi[1,i],conj(Gmbi[2,i])));
          { Q_01, une intégrale complexe a calculer }
          qsigC:=exp(Cmul(Cplx(-0.5*sigvali*sigvali,0),Cadd(Csqr(qmi[0,i]),Csqr(conj(qmi[1,i])))));
          qz2C:=Cmul(qmi[0,i],conj(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqrC:=Cdiv(C_one,qz2C);
          integR:=GLintegration(integrandSQre,qx,qz2C.re,qz2C.im);
          integI:=GLintegration(integrandSQim,qx,qz2C.re,qz2C.im);
          sommeMN:=sommeMN+4*pi*2*
              RealPart(Cmul(Cadd(Cmul(Gmbi[0,i],conj(Gmbi[1,i])),Cmul(Gmbi[3,i],conj(Gmbi[2,i]))),
                            Cmul(Cmul(qsigC,qqrC),Cplx(integR,integI))));
          { Q_02, une intégrale complexe a calculer }
          qsigC:=exp(Cmul(Cplx(-0.5*sigvali*sigvali,0),Cadd(Csqr(qmi[0,i]),Csqr(conj(qmi[2,i])))));
          qz2C:=Cmul(qmi[0,i],conj(qmi[2,i])); { valeur de q pour faire le calcul de l'integrale }
          qqrC:=Cdiv(C_one,qz2C);
          integR:=GLintegration(integrandSQre,qx,qz2C.re,qz2C.im);
          integI:=GLintegration(integrandSQim,qx,qz2C.re,qz2C.im);
          sommeMN:=sommeMN+4*pi*2*
              RealPart(Cmul(Cadd(Cmul(Gmbi[0,i],conj(Gmbi[2,i])),Cmul(Gmbi[3,i],conj(Gmbi[1,i]))),
                            Cmul(Cmul(qsigC,qqrC),Cplx(integR,integI))));
          intensite:=intensite+sqrk*SommeMN;
     end; { fin des boucle i }
     IDWBA:=intensite/16/pi2;
   end;
   FreeIt(k,kr,ki);FreeIt(kzi,kzf);FreeIt(DARi,DARf,DATi,DATf);
end;


{ =============================================================================}
{ calcul de l'intensité diffusée en DWBA                                       }
{ ref Stettner p.8, formalisme II                                              }
{ version IPC, on tient compte de correlations interplanaires                  }
{ INPUT:  lambdadiff (en nm)                                                   }
{         theta in et theta final en °                                         }
{         il faut thetain <> thetaf                                            }
{ OUTPUT: IDWBA l'intensité diffusée (partie réelle ??)                        }
{ N.B.: explication des variables (a comparer avec Stettner p.25)              }
{ MAJ 20.08.01 : inversion du signe de qx                                      }
{                correctiopns d'indices multiples                              }
{ on suit le formalisme II de Stettner, cela permet de ne pas avoir a distingeur entre le cas dioptre et pas dioptre }
{ les DAR et les DAT sont definis correctement (bottom de la couche)           }
{ QUID des signes des k ..., on peut reprendre Stettner, il a les memes conventions de (Oz)    }
{ =============================================================================}
procedure CalculNP_DWBA_II_IPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var IDWBA:double);

var i,j,f,mindx,nindx:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    zzi,zzj:TCplx;
    DARi,DATi,DARf,DATf:Tvec;
    k,kr,ki,kzi,kzf:Tvec;
    Gmi,Gmbi,qmi:array [0..3] of array [0..nbrlayersmax] of Tcplx;
    k0,qx,qq2r,qq2i,SFRe,SFIm:double;
    qz,sqrk,Gmb2,qq2,qqr,qsig,aux1,sommeMN,intensite:Tcplx;
    sommeTest:Tcplx; { ?? }
    xlimit:double;
    sigvali,xsivali,hvali,sigvalj,xsivalj,hvalj,xsizval,zi,zj:double;
    ti,tf,a1,deltatheta:double;
    nlimit:integer;
    xtot,xint,xi,xf,SFR,SFI:double;
    dthetai,dthetaf,res:double;
    zz:array[0..nbrLayersMax] of double; { tableua des cotes des interfaces }


begin
createit(k,kr,ki,kzi);createit(kzf);
createit(DARi,DARf,DATi,DATf);
nbrLayers := z.length;
k.size(nbrlayers,true);kr.size(nbrlayers,false);ki.size(nbrlayers,false);
kzi.size(nbrlayers,true);kzf.size(nbrlayers,true);
DARi.size(nbrlayers,true);DATi.size(nbrlayers,true);
DARf.size(nbrlayers,true);DATf.size(nbrlayers,true);

{ calcul des vecteurs d'ondes dans chaque couche -> kzin kzfinal }
{ avec nos conventions kzi est negatif et kzf est negatif (au-dessus de l'angle critique) }
{ calcul des coefficients de Fresnel complexes a chaque interface du système, dans la partie basse de la couche }
CalculAmpSpeculaire(lambda,thetai,z,V,Vi,sigma,xsi,exph,kzi,DARi,DATi);
CalculAmpSpeculaire(lambda,thetaf,z,V,Vi,sigma,xsi,exph,kzf,DARf,DATf);
qx:=2*pi/(lambda)*1e9*(cos(thetai*Deg)-cos(thetaf*Deg)); { signe etait faux, si on prend (Ox) suivant le faisceau incident }
{calcul des vecteurs k pour toutes les couches}
k0:=2*Pi/(lambda)*1e9;
kr.setval(hb2over2mn);kr.scale(k0*k0);kr.sub(V);
ki.setval(0);ki.sub(Vi);
k.realtocplx(kr,ki);k.scale(1/hb2over2mn);k.sqrt; { pourquoi pas multiplie par -1 ??? }

deltatheta:=abs(thetai-thetaf); { en radians ! }
{ creation d'un tableau des profondeurs }
zz[nbrLayers-1]:=0; { tableau des cotes des interfaces}
for i:=nbrLayers-2 downto 0 do zz[i]:=zz[i+1]-z.values[i];

intensite:=c_zero; { initilaisation de la somme de diffus }

xsizval:=xsiz;

if deltatheta<0.001     { limite en radians ?! }
then IDWBA:=sqr(Cabs(DARi.Cvalues[nbrlayers-1])) { on met lm'intensité spéculaire si on est trop proche du spéculaire , sinon point singulier }
else begin { sommes croisées sur les différentes couches de 0 à N }
      { for i:=layerindexi.value to layerindexj.value do
          for j:=layerindexi.value to layerindexj.value do }
       for i:=1 to nbrLayers-1 do { sommes croisées de la première couche (excepte substrat) au vide }
          for j:=1 to nbrLayers-1 do { tres inefficace car on recalcule tout plusieurs fois }
          begin
          qmi[0,i]:=Cadd(kzi.Cvalues[i-1],kzf.Cvalues[i-1]);qmi[1,i]:=Csub(kzi.Cvalues[i-1],kzf.Cvalues[i-1]);qmi[2,i]:=neg(qmi[1,i]);qmi[3,i]:=neg(qmi[0,i]);

          Gmi[0,i]:=Cmul(DATi.Cvalues[i-1],DATf.Cvalues[i-1]);
          Gmi[1,i]:=Cmul(DATi.Cvalues[i-1],DARf.Cvalues[i-1]);
          Gmi[2,i]:=Cmul(DARi.Cvalues[i-1],DATf.Cvalues[i-1]);
          Gmi[3,i]:=Cmul(DARi.Cvalues[i-1],DARf.Cvalues[i-1]);

{ Il faut mettre un facteur de phase }
{ on definit le zero au niveau de la surface de l'échantillon : interface vide/première couche <=> onde incidente a une phase nulle au niveau de la surface }
{ "profondeur" dans l'echantillon }
          zzi:=Cplx(zz[i],0); zzj:=Cplx(zz[j],0); { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          Gmbi[0,i]:=Cmul(Gmi[0,i],exp(Cdiv(Cmul(qmi[0,i],zzi),C_I)));
          Gmbi[1,i]:=Cmul(Gmi[1,i],exp(Cdiv(Cmul(qmi[1,i],zzi),C_I)));
          Gmbi[2,i]:=Cmul(Gmi[2,i],exp(Cdiv(Cmul(qmi[2,i],zzi),C_I)));
          Gmbi[3,i]:=Cmul(Gmi[3,i],exp(Cdiv(Cmul(qmi[3,i],zzi),C_I)));

          sigvali:=sigma.values[i-1];sigvalj:=sigma.values[j-1]; { indices OK }
          xsivali:=xsi.values[i-1];xsivalj:=xsi.values[j-1]; { indices OK }
          hvali:=exph.values[i-1];hvalj:=exph.values[j-1]; { indices OK }

          SommeMN:=C_zero; { sommes partielles sur les indices m et n }
          {sqrk: difference des carres des k; remplace K0^4 (n-n)(n-n)*}
          sqrk:=Cmul(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])),Conj(Csub(Csqr(k.Cvalues[j]),Csqr(k.Cvalues[j-1])))); { indices ok }
{          for mindx:=mindex.value to nindex.value do
              for nindx:=mindex.value to nindex.value do }
          for mindx:=0 to 3 do
              for nindx:=0 to 3 do
              begin
              { fot tests purposes only }
              if (mindx<0) or (mindx>3) then begin MessageDlg('problème avec midx', mtInformation,[mbOk], 0);halt;halt end;
              if (nindx<0) or (nindx>3) then begin MessageDlg('problème avec nidx', mtInformation,[mbOk], 0);halt;halt end;
              {Gmj multiplie par Gnj"tild"}
              Gmb2:=Cmul(Gmbi[mindx,i],Conj(Gmbi[nindx,j]));
              qsig:=Cmul(Cplx(-0.5,0),Cadd(Csqr(Cmul(qmi[mindx,i],Cplx(sigvali,0))),Csqr(Cmul(conj(qmi[nindx,j]),Cplx(sigvalj,0)))));
              {prefacteur au facteur de structure}
              aux1:=Cmul(Gmb2,exp(qsig));
              {calcul de facteur de structure}
              {qq2: qmj multiplie par conjugue de qnj}
              qq2:=Cmul(qmi[mindx,i],conj(qmi[nindx,j]));qq2r:=qq2.re;qq2i:=qq2.im;
              {prefacteur de l'integrale du facteur de structure}
              qqr:=Cdiv(Cplx(4*pi,0),qq2);
              zi:=z.values[i];zj:=z.values[j]; { for the vertical correlations }
              { calcul du facter de structure }
              xlimit:=xsivali*Power(-ln(critereintegral),1/(2*hvali)); { POURQUOI NE PRENDRE Q4UN CRITERE SUR XSIVALI ??? }
              SFRE:=0;SFIm:=0;
              xint:=abs(2*pi/qx); { on va integrer par périodes pour avoir plus de precision }
              if xint<(xlimit)
              then begin xi:=0;xf:=xint;nlimit:=floor(xlimit/xint)+1; end
              else begin xi:=0;xf:=xlimit;nlimit:=1;end;

              for f:=1 to nlimit do
              begin
                   GaussLegendreQuadrature(xi,xf,gaussOrder,Fre,sigvali,sigvalj,xsivali,xsivalj,xsizval,hvali,hvalj,zi,zj,qx,qq2r,qq2i,res,SFR);
                   GaussLegendreQuadrature(xi,xf,gaussOrder,Fim,sigvali,sigvalj,xsivali,xsivalj,xsizval,hvali,hvalj,zi,zj,qx,qq2r,qq2i,res,SFI);
                   SFRe:=SFRe+SFR;SFIm:=SFIm+SFI;
                   xi:=f*xint;xf:=(f+1)*xint;
              end;
              SommeMN:=Cadd(SommeMN,Cmul(aux1,Cmul(qqr,Cplx(SFRe,SFIm))));
          end; { fin des boucles m,n }
          intensite:=Cadd(intensite,Cmul(sqrk,SommeMN))
     end; { fin des boucle i,j }
     IDWBA:=Cabs(intensite)/16/pi/pi;
     end;
FreeIt(k,kr,ki);FreeIt(kzi,kzf);FreeIt(DARi,DARf,DATi,DATf);
end;

{ =============================================================================}
{ calcul de l'intensité diffusée en DWBA                                       }
{ ref Stettner p.8, formalisme III                                             }
{ Cette version se limite a NIPC                                               }
{ INPUT:  lambdadiff (en nm)                                                   }
{         theta in et theta final en °                                         }
{         il faut thetain <> thetaf                                            }
{ OUTPUT: IDWBA l'intensité diffusée (partie réelle ??)                        }
{ N.B.: explication des variables (a comparer avec Stettner p.25)              }
{ MAJ 20.08.01 : inversion du signe de qx                                      }
{                correctiopns d'indices multiples                              }
{ on suit le formalisme II de Stettner, cela permet de ne pas avoir a distingeur entre le cas dioptre et pas dioptre }
{ les DAR et les DAT sont definis correctement (bottom de la couche)           }
{ QUID des signes des k ..., on peut reprendre Stettner, il a les memes conventions de (Oz)    }
{ on prend le calcul poussé au bout pour acceler le process de calcul (voir classeur }
{ on gagne un facteur 4 en ne calculant pas systématiquement les 16 inegrales complexes }
{ DOUTE SUR LES SIGNES DES Z !!   cotes ou epaisseurs des couches ??           }
{ =============================================================================}
{ les tableaux qmi sont surdefinis                                             }
procedure CalculNP_DWBA_III_NIPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var IDWBA:double);
var i,j,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    zzi,zzj:TCplx;
    DARi,DATi,DARf,DATf:Tvec;
    k,kr,ki,kzi,kzf:Tvec;
    Gmi,Gmbi,qmi:array [0..3] of array [0..nbrlayersmax] of Tcplx;
    k0,qx,qq2r,qq2i,SFRe,SFIm,integ,integR,integI:double;
    qsig,qqr,qz2,sqrk,intensite,sommeMN:double;
    qsigC,qz2C,qqrC:Tcplx;
    xlimit:double;
    sigvali,sigvali2,xsivali,hvali,sigvalj,xsivalj,hvalj,zi,zj:double;
    ti,tf,a1,deltatheta:double;
    nlimit:integer;
    xtot,xint,xi,xf,SFR,SFI:double;
    dthetai,dthetaf,res:double;
    zz:array[0..nbrLayersMax] of double; { tableau des cotes des interfaces }

{ la perte de temps induite par cet appel est non mesurable within 5% }
function GLintegration(F:TGLFunction0;qx,qz2R,qz2I:double):double;
var n:integer;
begin
    SFRE:=0;xint:=abs(2*pi/qx);
    if xint<(xlimit) then begin xi:=0;xf:=xint;nlimit:=floor(xlimit/xint)+1; end
                     else begin xi:=0;xf:=xlimit;nlimit:=1; end;
    for n:=1 to nlimit do begin GaussLegendreQuadrature0(xi,xf,gaussOrder,F,sigvali2,xsivali,hvali,qx,qz2R,qz2I,SFR);
                                SFRe:=SFRe+SFR;xi:=n*xint;xf:=(n+1)*xint; end;
    result:=SFRe;
end;

begin
createit(k,kr,ki,kzi);createit(kzf);createit(DARi,DARf,DATi,DATf);
nbrLayers := z.length;
k.size(nbrlayers,true);kr.size(nbrlayers,false);ki.size(nbrlayers,false);kzi.size(nbrlayers,true);kzf.size(nbrlayers,true);
DARi.size(nbrlayers,true);DATi.size(nbrlayers,true);DARf.size(nbrlayers,true);DATf.size(nbrlayers,true);

{ calcul des vecteurs d'ondes dans chaque couche -> kzin kzfinal }
{ avec nos conventions kzi est negatif et kzf est negatif (au-dessus de l'angle critique) }
{ calcul des coefficients de Fresnel complexes a chaque interface du système, dans la partie basse de la couche }
CalculAmpSpeculaire(lambda,thetai,z,V,Vi,sigma,xsi,exph,kzi,DARi,DATi);CalculAmpSpeculaire(lambda,thetaf,z,V,Vi,sigma,xsi,exph,kzf,DARf,DATf);
qx:=2*pi/(lambda)*1e9*(cos(thetai*Deg)-cos(thetaf*Deg)); { signe etait faux, si on prend (Ox) suivant le faisceau incident }
{calcul des vecteurs k pour toutes les couches}
k0:=2*Pi/(lambda)*1e9;
kr.setval(hb2over2mn);kr.scale(k0*k0);kr.sub(V);
ki.setval(0);ki.sub(Vi);
k.realtocplx(kr,ki);k.scale(1/hb2over2mn);k.sqrt; { pourquoi pas multiplie par -1 ??? }

deltatheta:=abs(thetai-thetaf); { en radians ! }
{ creation d'un tableau des profondeurs }
zz[nbrLayers-1]:=0; { tableua des cotes des interfaces }
for i:=nbrLayers-2 downto 0 do zz[i]:=zz[i+1]-z.values[i];

intensite:=0; { initilaisation de la somme de diffus }
if deltatheta<0.001     { limite en radians ?! }
then IDWBA:=sqr(Cabs(DARi.Cvalues[nbrlayers-1])) { on met lm'intensité spéculaire si on est trop proche du spéculaire , sinon point singulier }
else begin { sommes croisées sur les différentes couches de 0 à N }
       for i:=1 to nbrLayers-1 do { sommes croisées de la première couche (excepte substrat) au vide }
          begin
          qmi[0,i]:=Cadd(kzi.Cvalues[i],kzf.Cvalues[i]);qmi[1,i]:=Csub(kzi.Cvalues[i],kzf.Cvalues[i]);qmi[2,i]:=neg(qmi[1,i]);qmi[3,i]:=neg(qmi[0,i]);
          { facteur de pahse pour les coefficients G ~ }
          { on definit le zero au niveau de la surface de l'échantillon : interface vide/première couche <=> onde incidente a une phase nulle au niveau de la surface }
          { "profondeur" dans l'echantillon }
          zzi:=Cplx(zz[i],0);  { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          Gmi[0,i]:=Cmul(DATi.Cvalues[i],DATf.Cvalues[i]); Gmbi[0,i]:=Cmul(Gmi[0,i],exp(Cdiv(Cmul(qmi[0,i],zzi),C_I)));
          Gmi[1,i]:=Cmul(DATi.Cvalues[i],DARf.Cvalues[i]); Gmbi[1,i]:=Cmul(Gmi[1,i],exp(Cdiv(Cmul(qmi[1,i],zzi),C_I)));
          Gmi[2,i]:=Cmul(DARi.Cvalues[i],DATf.Cvalues[i]); Gmbi[2,i]:=Cmul(Gmi[2,i],exp(Cdiv(Cmul(qmi[2,i],zzi),C_I)));
          Gmi[3,i]:=Cmul(DARi.Cvalues[i],DARf.Cvalues[i]); Gmbi[3,i]:=Cmul(Gmi[3,i],exp(Cdiv(Cmul(qmi[3,i],zzi),C_I)));
          sigvali:=sigma.values[i-1];xsivali:=xsi.values[i-1];hvali:=exph.values[i-1];
          sigvali2:=sigvali*sigvali;

          sqrk:=Sqr(Cabs(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])))); { valeur reelle ! }
          SommeMN:=0;
          xlimit:=xsivali*Power(-ln(critereIntegral),1/(2*hvali)); { EVITER L'APPEL A UN CONTROEL GRAPHIQUE .value !!! }
          { sommation sur les indices m et n                              }
          { on fait la somme de manière explicite pour booster le calcul  }
          { Q_00, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[0,i].re)-sqr(qmi[0,i].im)));
          qz2:=sqr(Cabs(qmi[0,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2;
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*(Sqr(Cabs(Gmbi[0,i]))+Sqr(Cabs(Gmbi[3,i])));
          { Q_11, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[1,i].re)-sqr(qmi[1,i].im)));
          qz2:=sqr(Cabs(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2;
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*(Sqr(Cabs(Gmbi[1,i]))+Sqr(Cabs(Gmbi[2,i])));
          { Q_03, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[0,i].re)-sqr(qmi[0,i].im)));
          qz2:=-sqr(Cabs(qmi[0,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2; { attention jonglage avec les signes, voir calculs papier }
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*2*RealPart(Cmul(Gmbi[0,i],conj(Gmbi[3,i])));
          { Q_12, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[1,i].re)-sqr(qmi[1,i].im)));
          qz2:=-sqr(Cabs(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2;
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*2*RealPart(Cmul(Gmbi[1,i],conj(Gmbi[2,i])));
          { Q_01, une intégrale complexe a calculer }
          qsigC:=exp(Cmul(Cplx(-0.5*sigvali*sigvali,0),Cadd(Csqr(qmi[0,i]),Csqr(conj(qmi[1,i])))));
          qz2C:=Cmul(qmi[0,i],conj(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqrC:=Cdiv(C_one,qz2C);
          integR:=GLintegration(integrandSQre,qx,qz2C.re,qz2C.im);
          integI:=GLintegration(integrandSQim,qx,qz2C.re,qz2C.im);
          sommeMN:=sommeMN+4*pi*2*
              RealPart(Cmul(Cadd(Cmul(Gmbi[0,i],conj(Gmbi[1,i])),Cmul(Gmbi[3,i],conj(Gmbi[2,i]))),
                            Cmul(Cmul(qsigC,qqrC),Cplx(integR,integI))));
          { Q_02, une intégrale complexe a calculer }
          qsigC:=exp(Cmul(Cplx(-0.5*sigvali*sigvali,0),Cadd(Csqr(qmi[0,i]),Csqr(conj(qmi[2,i])))));
          qz2C:=Cmul(qmi[0,i],conj(qmi[2,i])); { valeur de q pour faire le calcul de l'integrale }
          qqrC:=Cdiv(C_one,qz2C);
          integR:=GLintegration(integrandSQre,qx,qz2C.re,qz2C.im);
          integI:=GLintegration(integrandSQim,qx,qz2C.re,qz2C.im);
          sommeMN:=sommeMN+4*pi*2*
              RealPart(Cmul(Cadd(Cmul(Gmbi[0,i],conj(Gmbi[2,i])),Cmul(Gmbi[3,i],conj(Gmbi[1,i]))),
                            Cmul(Cmul(qsigC,qqrC),Cplx(integR,integI))));
          intensite:=intensite+sqrk*SommeMN;
     end; { fin des boucle i }
     IDWBA:=intensite/16/pi2;
   end;
   FreeIt(k,kr,ki);FreeIt(kzi,kzf);FreeIt(DARi,DARf,DATi,DATf);
end;

{ =============================================================================}
{ calcul de l'intensité diffusée en DWBA                                       }
{ ref Stettner p.8, formalisme III                                             }
{ Cette version se limite a NIPC                                               }
{ cas particulier d el'approximation petits angles : calcu analytiqe de la FT  }
{ INPUT:  lambdadiff (en nm)                                                   }
{         theta in et theta final en °                                         }
{         il faut thetain <> thetaf                                            }
{ OUTPUT: IDWBA l'intensité diffusée (partie réelle ??)                        }
{ N.B.: explication des variables (a comparer avec Stettner p.25)              }
{ MAJ 20.08.01 : inversion du signe de qx                                      }
{                correctiopns d'indices multiples                              }
{ on suit le formalisme II de Stettner, cela permet de ne pas avoir a distingeur entre le cas dioptre et pas dioptre }
{ les DAR et les DAT sont definis correctement (bottom de la couche)           }
{ QUID des signes des k ..., on peut reprendre Stettner, il a les memes conventions de (Oz)    }
{ on prend le calcul poussé au bout pour acceler le process de calcul (voir classeur }
{ on gagne un facteur 4 en ne calculant pas systématiquement les 16 inegrales complexes }
{ DOUTE SUR LES SIGNES DES Z !!   cotes ou epaisseurs des couches ??           }
{ =============================================================================}
{ les tableaux qmi sont surdefinis                                             }
procedure CalculNP_DWBA_III_SA_NIPC(typeCorr:integer;lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var IDWBA:double);
var i,j,n:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    zzi,zzj:TCplx;
    DARi,DATi,DARf,DATf:Tvec;
    k,kr,ki,kzi,kzf:Tvec;
    Gmi,Gmbi,qmi:array [0..3] of array [0..nbrlayersmax] of Tcplx;
    k0,qx,qq2r,qq2i,SFRe,SFIm,integ,integR,integI:double;
    qsig,qqr,qz2,sqrk,intensite,sommeMN:double;
    qsigC,qz2C,qqrC:Tcplx;
    xlimit:double;
    sigvali,sigvali2,xsivali,hvali,sigvalj,xsivalj,hvalj,zi,zj:double;
    ti,tf,a1,deltatheta:double;
    nlimit:integer;
    xtot,xint,xi,xf,SFR,SFI:double;
    dthetai,dthetaf,res:double;
    zz:array[0..nbrLayersMax] of double; { tableau des cotes des interfaces }

begin
createit(k,kr,ki,kzi);createit(kzf);createit(DARi,DARf,DATi,DATf);
nbrLayers := z.length;
k.size(nbrlayers,true);kr.size(nbrlayers,false);ki.size(nbrlayers,false);kzi.size(nbrlayers,true);kzf.size(nbrlayers,true);
DARi.size(nbrlayers,true);DATi.size(nbrlayers,true);DARf.size(nbrlayers,true);DATf.size(nbrlayers,true);

{ calcul des vecteurs d'ondes dans chaque couche -> kzin kzfinal }
{ avec nos conventions kzi est negatif et kzf est negatif (au-dessus de l'angle critique) }
{ calcul des coefficients de Fresnel complexes a chaque interface du système, dans la partie basse de la couche }
CalculAmpSpeculaire(lambda,thetai,z,V,Vi,sigma,xsi,exph,kzi,DARi,DATi);CalculAmpSpeculaire(lambda,thetaf,z,V,Vi,sigma,xsi,exph,kzf,DARf,DATf);
qx:=2*pi/(lambda)*1e9*(cos(thetai*Deg)-cos(thetaf*Deg)); { signe etait faux, si on prend (Ox) suivant le faisceau incident }
{calcul des vecteurs k pour toutes les couches}
k0:=2*Pi/(lambda)*1e9;
kr.setval(hb2over2mn);kr.scale(k0*k0);kr.sub(V);
ki.setval(0);ki.sub(Vi);
k.realtocplx(kr,ki);k.scale(1/hb2over2mn);k.sqrt; { pourquoi pas multiplie par -1 ??? }

deltatheta:=abs(thetai-thetaf); { en radians ! }
{ creation d'un tableau des profondeurs }
zz[nbrLayers-1]:=0; { tableua des cotes des interfaces }
for i:=nbrLayers-2 downto 0 do zz[i]:=zz[i+1]-z.values[i];

intensite:=0; { initilaisation de la somme de diffus }
if deltatheta<0.001     { limite en radians ?! }
then IDWBA:=sqr(Cabs(DARi.Cvalues[nbrlayers-1])) { on met lm'intensité spéculaire si on est trop proche du spéculaire , sinon point singulier }
else begin { sommes croisées sur les différentes couches de 0 à N }
       for i:=1 to nbrLayers-1 do { sommes croisées de la première couche (excepte substrat) au vide }
          begin
          qmi[0,i]:=Cadd(kzi.Cvalues[i],kzf.Cvalues[i]);qmi[1,i]:=Csub(kzi.Cvalues[i],kzf.Cvalues[i]);qmi[2,i]:=neg(qmi[1,i]);qmi[3,i]:=neg(qmi[0,i]);
          { facteur de pahse pour les coefficients G ~ }
          { on definit le zero au niveau de la surface de l'échantillon : interface vide/première couche <=> onde incidente a une phase nulle au niveau de la surface }
          { "profondeur" dans l'echantillon }
          zzi:=Cplx(zz[i],0);  { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          Gmi[0,i]:=Cmul(DATi.Cvalues[i],DATf.Cvalues[i]); Gmbi[0,i]:=Cmul(Gmi[0,i],exp(Cdiv(Cmul(qmi[0,i],zzi),C_I)));
          Gmi[1,i]:=Cmul(DATi.Cvalues[i],DARf.Cvalues[i]); Gmbi[1,i]:=Cmul(Gmi[1,i],exp(Cdiv(Cmul(qmi[1,i],zzi),C_I)));
          Gmi[2,i]:=Cmul(DARi.Cvalues[i],DATf.Cvalues[i]); Gmbi[2,i]:=Cmul(Gmi[2,i],exp(Cdiv(Cmul(qmi[2,i],zzi),C_I)));
          Gmi[3,i]:=Cmul(DARi.Cvalues[i],DARf.Cvalues[i]); Gmbi[3,i]:=Cmul(Gmi[3,i],exp(Cdiv(Cmul(qmi[3,i],zzi),C_I)));
          sigvali:=sigma.values[i-1];xsivali:=xsi.values[i-1];hvali:=exph.values[i-1];
          sigvali2:=sigvali*sigvali;

          sqrk:=Sqr(Cabs(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])))); { valeur reelle ! }
          SommeMN:=0;
          { sommation sur les indices m et n                              }
          { on fait la somme de manière explicite pour booster le calcul  }
          { Q_00, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[0,i].re)-sqr(qmi[0,i].im)));
          qz2:=sqr(Cabs(qmi[0,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=2*pi; { LE FACTEUR 2 DISPARAIT : IL PASSE DANS LA TF }
          case typeCorr of
          0:SFRe:=3.14;
          1:SFRe:=sqrt(pi)*xsivali*sigvali2*exp(-qx*qx/4*xsivali*xsivali); { cas gaussian }
          2:SFRe:=2*sigvali2*xsivali/(xsivali*xsivali*qx*qx-1);
          3:SFRe:=3.14;
          end;
          integ:=SFRe;
          sommeMN:=sommeMN+qqr*qsig*integ*(Sqr(Cabs(Gmbi[0,i]))+Sqr(Cabs(Gmbi[3,i])));
          { Q_11, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[1,i].re)-sqr(qmi[1,i].im)));
          qz2:=sqr(Cabs(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=2*pi; { LE FACTEUR 2 DISPARAIT : IL PASSE DANS LA TF }
          case typeCorr of
          0:SFRe:=3.14;
          1:SFRe:=sqrt(pi)*xsivali*sigvali2*exp(-qx*qx/4*xsivali*xsivali); { cas gaussian }
          2:SFRe:=2*sigvali2*xsivali/(xsivali*xsivali*qx*qx-1);
          3:SFRe:=3.14;
          end;
          integ:=SFRe;
          sommeMN:=sommeMN+qqr*qsig*integ*(Sqr(Cabs(Gmbi[1,i]))+Sqr(Cabs(Gmbi[2,i])));
          { Q_03, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[0,i].re)-sqr(qmi[0,i].im)));
          qz2:=-sqr(Cabs(qmi[0,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2; { attention jonglage avec les signes, voir calculs papier }
{          integ:=GLintegration(integrandSQ,qx,qz2,0); }
          sommeMN:=sommeMN+qqr*qsig*integ*2*RealPart(Cmul(Gmbi[0,i],conj(Gmbi[3,i])));
          { Q_12, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[1,i].re)-sqr(qmi[1,i].im)));
          qz2:=-sqr(Cabs(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2;
{          integ:=GLintegration(integrandSQ,qx,qz2,0); }
          sommeMN:=sommeMN+qqr*qsig*integ*2*RealPart(Cmul(Gmbi[1,i],conj(Gmbi[2,i])));
          { Q_01, une intégrale complexe a calculer }
          qsigC:=exp(Cmul(Cplx(-0.5*sigvali*sigvali,0),Cadd(Csqr(qmi[0,i]),Csqr(conj(qmi[1,i])))));
          qz2C:=Cmul(qmi[0,i],conj(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqrC:=Cdiv(C_one,qz2C);
{          integR:=GLintegration(integrandSQre,qx,qz2C.re,qz2C.im); }
{          integI:=GLintegration(integrandSQim,qx,qz2C.re,qz2C.im); }
          sommeMN:=sommeMN+4*pi*2*
              RealPart(Cmul(Cadd(Cmul(Gmbi[0,i],conj(Gmbi[1,i])),Cmul(Gmbi[3,i],conj(Gmbi[2,i]))),
                            Cmul(Cmul(qsigC,qqrC),Cplx(integR,integI))));
          { Q_02, une intégrale complexe a calculer }
          qsigC:=exp(Cmul(Cplx(-0.5*sigvali*sigvali,0),Cadd(Csqr(qmi[0,i]),Csqr(conj(qmi[2,i])))));
          qz2C:=Cmul(qmi[0,i],conj(qmi[2,i])); { valeur de q pour faire le calcul de l'integrale }
          qqrC:=Cdiv(C_one,qz2C);
{          integR:=GLintegration(integrandSQre,qx,qz2C.re,qz2C.im);
          integI:=GLintegration(integrandSQim,qx,qz2C.re,qz2C.im); }
          sommeMN:=sommeMN+4*pi*2*
              RealPart(Cmul(Cadd(Cmul(Gmbi[0,i],conj(Gmbi[2,i])),Cmul(Gmbi[3,i],conj(Gmbi[1,i]))),
                            Cmul(Cmul(qsigC,qqrC),Cplx(integR,integI))));
          intensite:=intensite+sqrk*SommeMN;
     end; { fin des boucle i }
     IDWBA:=intensite/16/pi2;
   end;
   FreeIt(k,kr,ki);FreeIt(kzi,kzf);FreeIt(DARi,DARf,DATi,DATf);
end;


{ =============================================================================}
{ calcul de l'intensité diffusée en DWBA                                       }
{ ref Stettner p.8, formalisme III                                             }
{ Cette version se limite a NIPC                                               }
{ INPUT:  lambdadiff (en nm)                                                   }
{         theta in et theta final en °                                         }
{         il faut thetain <> thetaf                                            }
{ OUTPUT: IDWBA l'intensité diffusée (partie réelle ??)                        }
{ N.B.: explication des variables (a comparer avec Stettner p.25)              }
{ MAJ 20.08.01 : inversion du signe de qx                                      }
{                correctiopns d'indices multiples                              }
{ on suit le formalisme II de Stettner, cela permet de ne pas avoir a distingeur entre le cas dioptre et pas dioptre }
{ les DAR et les DAT sont definis correctement (bottom de la couche)           }
{ QUID des signes des k ..., on peut reprendre Stettner, il a les memes conventions de (Oz)    }
{ on prend le calcul poussé au bout pour acceler le process de calcul (voir classeur }
{ on gagne un facteur 4 en ne calculant pas systématiquement les 16 inegrales complexes }
{ DOUTE SUR LES SIGNES DES Z !!   cotes ou epaisseurs des couches ??           }
{ =============================================================================}
{ remplir les tableaux une seule fois                                          }
procedure CalculNP_DWBA_III_IPC(lambda,thetai,thetaf:double;var z,V,Vi,sigma,xsi,exph:Tvec;var IDWBA:double);
var i,j,mm,nn,midx,nidx:integer;
    nbrLayers:integer; { nombre de couches dans le système }
    zzi,zzj:TCplx;
    DARi,DATi,DARf,DATf:Tvec;
    k,kr,ki,kzi,kzf:Tvec;
    Gmi,Gmbi,qmi:array [0..3] of array [0..nbrlayersmax] of Tcplx;
    k0,qx,qq2r,qq2i,SFRe,SFIm,integ,integR,integI:double;
    qsig,qqr,qz2,sqrk,intensite,sommeMN:double;
    qsigC,qz2C,qqrC,nij:Tcplx;
    xlimit:double;
    sigvali,sigvali2,xsivali,hvali,sigvalj,sigvalj2,xsivalj,hvalj,zi,zj:double;
    ti,tf,a1,deltatheta:double;
    nlimit:integer;
    xtot,xint,xi,xf,SFR,SFI:double;
    dthetai,dthetaf,res:double;
    zz:array[0..nbrLayersMax] of double; { tableau des cotes des interfaces }
    xsiIPC:double;

{ la perte de temps induite par cet appel est non mesurable within 5% }
function GLintegration(F:TGLFunction0;qx,qz2R,qz2I:double):double;
var n:integer;
begin
    SFRE:=0;xint:=abs(2*pi/qx);
    if xint<(xlimit) then begin xi:=0;xf:=xint;nlimit:=floor(xlimit/xint)+1; end
                     else begin xi:=0;xf:=xlimit;nlimit:=1; end;
    for n:=1 to nlimit do begin GaussLegendreQuadrature0(xi,xf,gaussOrder,F,sigvali2,xsivali,hvali,qx,qz2R,qz2I,SFR);
                                SFRe:=SFRe+SFR;xi:=n*xint;xf:=(n+1)*xint; end;
    result:=SFRe;
end;

function GLintegration_IPC(F:TGLFunction;qx,qz2R,qz2I,zphase,xsiz:double):double;
var n:integer;
begin
    SFRE:=0;xint:=abs(2*pi/qx);
    if xint<(xlimit) then begin xi:=0;xf:=xint;nlimit:=floor(xlimit/xint)+1; end
                     else begin xi:=0;xf:=xlimit;nlimit:=1; end;
    for n:=1 to nlimit do begin GaussLegendreQuadrature(xi,xf,gaussOrder,F,sigvali2,xsivali,hvali,qx,qz2R,qz2I,zphase,xsiz,0,0,0,0,0,SFR);
                                SFRe:=SFRe+SFR;xi:=n*xint;xf:=(n+1)*xint; end;
    result:=SFRe;
end;

begin
xsiIPC:=xsiz;              { en nm }
createit(k,kr,ki,kzi);createit(kzf);createit(DARi,DARf,DATi,DATf);
nbrLayers := z.length;
k.size(nbrlayers,true);kr.size(nbrlayers,false);ki.size(nbrlayers,false);kzi.size(nbrlayers,true);kzf.size(nbrlayers,true);
DARi.size(nbrlayers,true);DATi.size(nbrlayers,true);DARf.size(nbrlayers,true);DATf.size(nbrlayers,true);

{ calcul des vecteurs d'ondes dans chaque couche -> kzin kzfinal }
{ avec nos conventions kzi est negatif et kzf est negatif (au-dessus de l'angle critique) }
{ calcul des coefficients de Fresnel complexes a chaque interface du système, dans la partie basse de la couche }
CalculAmpSpeculaire(lambda,thetai,z,V,Vi,sigma,xsi,exph,kzi,DARi,DATi);CalculAmpSpeculaire(lambda,thetaf,z,V,Vi,sigma,xsi,exph,kzf,DARf,DATf);
qx:=2*pi/(lambda)*1e9*(cos(thetai*Deg)-cos(thetaf*Deg)); { signe etait faux, si on prend (Ox) suivant le faisceau incident }
{calcul des vecteurs k pour toutes les couches}
k0:=2*Pi/(lambda)*1e9;
kr.setval(hb2over2mn);kr.scale(k0*k0);kr.sub(V);
ki.setval(0);ki.sub(Vi);
k.realtocplx(kr,ki);k.scale(1/hb2over2mn);k.sqrt; { pourquoi pas multiplie par -1 ??? }

deltatheta:=abs(thetai-thetaf); { en radians ! }
{ creation d'un tableau des profondeurs }
zz[nbrLayers-1]:=0; { tableau des cotes des interfaces }
for i:=nbrLayers-2 downto 0 do zz[i]:=zz[i+1]-z.values[i];

intensite:=0; { initilaisation de la somme de diffus }
if deltatheta<0.001     { limite en radians ?! }
then IDWBA:=sqr(Cabs(DARi.Cvalues[nbrlayers-1])) { on met lm'intensité spéculaire si on est trop proche du spéculaire , sinon point singulier }
else begin { sommes croisées sur les différentes couches de 0 à N }
     { creation des tableaux des Gm, Gmbi et des qmi }
     for i:=1 to nbrLayers-1 do
     begin
          zzi:=Cplx(zz[i],0);  { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          qmi[0,i]:=Cadd(kzi.Cvalues[i],kzf.Cvalues[i]);qmi[1,i]:=Csub(kzi.Cvalues[i],kzf.Cvalues[i]);
          qmi[2,i]:=neg(qmi[1,i]);qmi[3,i]:=neg(qmi[0,i]);
          Gmi[0,i]:=Cmul(DATi.Cvalues[i],DATf.Cvalues[i]); Gmbi[0,i]:=Cmul(Gmi[0,i],exp(Cdiv(Cmul(qmi[0,i],zzi),C_I)));
          Gmi[1,i]:=Cmul(DATi.Cvalues[i],DARf.Cvalues[i]); Gmbi[1,i]:=Cmul(Gmi[1,i],exp(Cdiv(Cmul(qmi[1,i],zzi),C_I)));
          Gmi[2,i]:=Cmul(DARi.Cvalues[i],DATf.Cvalues[i]); Gmbi[2,i]:=Cmul(Gmi[2,i],exp(Cdiv(Cmul(qmi[2,i],zzi),C_I)));
          Gmi[3,i]:=Cmul(DARi.Cvalues[i],DARf.Cvalues[i]); Gmbi[3,i]:=Cmul(Gmi[3,i],exp(Cdiv(Cmul(qmi[3,i],zzi),C_I)));
     end;
     for i:=1 to nbrLayers-1 do { termes avec i = j  } { on conserve ce qui a été ecrit dans DWBA III NIPC }
     begin
          { facteur de pahse pour les coefficients G ~ }
          { on definit le zero au niveau de la surface de l'échantillon : interface vide/première couche <=> onde incidente a une phase nulle au niveau de la surface }
          { "profondeur" dans l'echantillon }
          zzi:=Cplx(zz[i],0);  { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          sigvali:=sigma.values[i-1];xsivali:=xsi.values[i-1];hvali:=exph.values[i-1];
          sigvali2:=sigvali*sigvali;

          sqrk:=Sqr(Cabs(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])))); { valeur reelle ! }
          SommeMN:=0;
          xlimit:=xsivali*Power(-ln(critereIntegral),1/(2*hvali)); { EVITER L'APPEL A UN CONTROEL GRAPHIQUE .value !!! }
          { sommation sur les indices m et n                              }
          { on fait la somme de manière explicite pour booster le calcul  }
          { Q_00, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[0,i].re)-sqr(qmi[0,i].im)));
          qz2:=sqr(Cabs(qmi[0,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2;
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*(Sqr(Cabs(Gmbi[0,i]))+Sqr(Cabs(Gmbi[3,i])));
          { Q_11, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[1,i].re)-sqr(qmi[1,i].im)));
          qz2:=sqr(Cabs(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2;
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*(Sqr(Cabs(Gmbi[1,i]))+Sqr(Cabs(Gmbi[2,i])));
          { Q_03, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[0,i].re)-sqr(qmi[0,i].im)));
          qz2:=-sqr(Cabs(qmi[0,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2; { attention jonglage avec les signes, voir calculs papier }
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*2*RealPart(Cmul(Gmbi[0,i],conj(Gmbi[3,i])));
          { Q_12, une intégrale reelle a calculer }
          qsig:=exp(-sigvali*sigvali*(sqr(qmi[1,i].re)-sqr(qmi[1,i].im)));
          qz2:=-sqr(Cabs(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqr:=4*pi/qz2;
          integ:=GLintegration(integrandSQ,qx,qz2,0);
          sommeMN:=sommeMN+qqr*qsig*integ*2*RealPart(Cmul(Gmbi[1,i],conj(Gmbi[2,i])));
          { Q_01, une intégrale complexe a calculer }
          qsigC:=exp(Cmul(Cplx(-0.5*sigvali*sigvali,0),Cadd(Csqr(qmi[0,i]),Csqr(conj(qmi[1,i])))));
          qz2C:=Cmul(qmi[0,i],conj(qmi[1,i])); { valeur de q pour faire le calcul de l'integrale }
          qqrC:=Cdiv(C_one,qz2C);
          integR:=GLintegration(integrandSQre,qx,qz2C.re,qz2C.im);
          integI:=GLintegration(integrandSQim,qx,qz2C.re,qz2C.im);
          sommeMN:=sommeMN+4*pi*2*
              RealPart(Cmul(Cadd(Cmul(Gmbi[0,i],conj(Gmbi[1,i])),Cmul(Gmbi[3,i],conj(Gmbi[2,i]))),
                            Cmul(Cmul(qsigC,qqrC),Cplx(integR,integI))));
          { Q_02, une intégrale complexe a calculer }
          qsigC:=exp(Cmul(Cplx(-0.5*sigvali*sigvali,0),Cadd(Csqr(qmi[0,i]),Csqr(conj(qmi[2,i])))));
          qz2C:=Cmul(qmi[0,i],conj(qmi[2,i])); { valeur de q pour faire le calcul de l'integrale }
          qqrC:=Cdiv(C_one,qz2C);
          integR:=GLintegration(integrandSQre,qx,qz2C.re,qz2C.im);
          integI:=GLintegration(integrandSQim,qx,qz2C.re,qz2C.im);
          sommeMN:=sommeMN+4*pi*2*
              RealPart(Cmul(Cadd(Cmul(Gmbi[0,i],conj(Gmbi[2,i])),Cmul(Gmbi[3,i],conj(Gmbi[1,i]))),
                            Cmul(Cmul(qsigC,qqrC),Cplx(integR,integI))));
          intensite:=intensite+sqrk*SommeMN;
     end; { fin des boucle i }
     for i:=1 to nbrLayers-1 do { termes avec i <> j  }
       for j:=1 to i-1 do
          begin
          SommeMN:=0;
          { "profondeur" dans l'echantillon }
          zzi:=Cplx(zz[i],0);  { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          sigvali:=sigma.values[i-1];xsivali:=xsi.values[i-1];hvali:=exph.values[i-1];
          sigvali2:=sigvali*sigvali;
          zzj:=Cplx(zz[j],0);  { c'est un facteur de phase, pas une épaisseur!, on utilise le tableau de cotes zz }
          sigvalj:=sigma.values[j-1];xsivalj:=xsi.values[j-1];hvalj:=exph.values[j-1];
          sigvalj2:=sigvalj*sigvalj;

          nij:=Cmul(Csub(Csqr(k.Cvalues[i]),Csqr(k.Cvalues[i-1])),Csub(Csqr(k.Cvalues[j]),Csqr(k.Cvalues[j-1])));
          xlimit:=xsivali*Power(-ln(critereIntegral),1/(2*hvali));
          for mm:=0 to 3 do { somme sur les couples m = n }
          begin
             midx:=mm;
             qsigC:=exp(Cadd(Cmul(Cplx(sigvali/sqrt(2),0),qmi[midx,i]),Cmul(Cplx(sigvalj/sqrt(2),0),conj(qmi[midx,j]))));
             qz2C:=Cmul(qmi[midx,i],conj(qmi[midx,j])); { valeur de q pour faire le calcul de l'integrale }
             qqrC:=Cdiv(C_one,qz2C);
             integR:=GLintegration_IPC(integrandSQre_IPC,qx,qz2C.re,qz2C.im,zz[i]-zz[j],xsiIPC);
             integI:=GLintegration_IPC(integrandSQim_IPC,qx,qz2C.re,qz2C.im,zz[i]-zz[j],xsiIPC);
             sommeMN:=sommeMN+4*pi*2*
               RealPart(Cmul(Cmul(Cmul(Cmul(nij,Cmul(Gmbi[midx,i],conj(Gmbi[midx,j]))),qsigC),qqrC),Cplx(integR,integI)));
          end;
          for mm:=0 to 3 do { corriger le signe de exp (q q sig sig) }
             for nn:=mm+1 to 3 do
             begin
                midx:=mm;nidx:=nn;
                qsigC:=exp(Cadd(Cmul(Cplx(sigvali/sqrt(2),0),qmi[midx,i]),Cmul(Cplx(sigvalj/sqrt(2),0),conj(qmi[nidx,j]))));
                qz2C:=Cmul(qmi[midx,i],conj(qmi[nidx,j])); { valeur de q pour faire le calcul de l'integrale }
                qqrC:=Cdiv(C_one,qz2C);
                integR:=GLintegration_IPC(integrandSQre_IPC,qx,qz2C.re,qz2C.im,zz[i]-zz[j],xsiIPC);
                integI:=GLintegration_IPC(integrandSQim_IPC,qx,qz2C.re,qz2C.im,zz[i]-zz[j],xsiIPC);
                sommeMN:=sommeMN+4*pi*2*
                RealPart(Cmul(Cmul(Cmul(Cmul(nij,Cmul(Gmbi[midx,i],conj(Gmbi[nidx,j]))),qsigC),qqrC),Cplx(integR,integI)));

                qsigC:=exp(Cadd(Cmul(Cplx(sigvali/sqrt(2),0),qmi[nidx,i]),Cmul(Cplx(sigvalj/sqrt(2),0),conj(qmi[midx,j]))));
                qz2C:=Cmul(qmi[nidx,i],conj(qmi[midx,j])); { valeur de q pour faire le calcul de l'integrale }
                qqrC:=Cdiv(C_one,qz2C);
                integR:=GLintegration_IPC(integrandSQre_IPC,qx,qz2C.re,qz2C.im,zz[i]-zz[j],xsiIPC);
                integI:=GLintegration_IPC(integrandSQim_IPC,qx,qz2C.re,qz2C.im,zz[i]-zz[j],xsiIPC);
                sommeMN:=sommeMN+4*pi*2*
                RealPart(Cmul(Cmul(Cmul(Cmul(nij,Cmul(Gmbi[nidx,i],conj(Gmbi[midx,j]))),qsigC),qqrC),Cplx(integR,integI)));
             end;
         intensite:=intensite+SommeMN; { on cumule les sommes MN pour un couple (i,j) donné) }
     end; { fin des boucles i,j }
     IDWBA:=intensite/16/pi2;
   end;
   FreeIt(k,kr,ki);FreeIt(kzi,kzf);FreeIt(DARi,DARf,DATi,DATf);
end;

end.
