// Gauss-Legendre Quadrature
// Adapted from "Applied Numerical Methods," Carnahan, Luther, Wilkes,
// John Wiley & Sons, 1969, pp. 101-112.
{modifié par TDD afin de pouvoir faire l'integration de fonctions a plusieurs parametres(ici 9)}

UNIT GaussLegendre;

INTERFACE

  TYPE
  TGLFunction = FUNCTION(x: double;par1,par2,par3,par4,par5,par6,par7,par8,par9,par10,par11,par12,par13:double):  DOUBLE;
  TGLFunction0 = FUNCTION(x: double;par1,par2,par3,par4,par5,par6:double):  DOUBLE;

  Procedure GaussLegendreQuadrature(CONST A, B:  DOUBLE;
                                    CONST N:  INTEGER;
                                    CONST F:  TGLFunction;
                                    CONST par1,par2,par3,par4,par5,par6,par7,par8,par9,par10,par11,par12,par13:double;
                                    var result:double);
  Procedure GaussLegendreQuadrature0(CONST A, B:  DOUBLE;
                                    CONST N:  INTEGER;
                                    CONST F:  TGLFunction0;
                                    CONST par1,par2,par3,par4,par5,par6:double;
                                    var result:double);
  { intégration de foncionts oscillantes }
  Procedure GaussLegendreOsc(CONST A, B, C:  DOUBLE;
                                    CONST N:  INTEGER;
                                    CONST F:  TGLFunction0;
                                    CONST par1,par2,par3,par4,par5,par6:double;
                                    var result:double);

  function FCorr(x,sigvali,sigvalj,xsivali,xsivalj,xsizval,hvali,hvalj,zi,zj:double):double; { TGLFunction0 }
  function integrandSQ(x,sigval2,xsival,hval,qx,qz2,p6:double):double;         { TGLFunction0 }
  function integrandSQre(x,sigval2,xsival,hval,qx,qz2re,qz2im:double):double;  { TGLFunction0 }
  function integrandSQim(x,sigval2,xsival,hval,qx,qz2re,qz2im:double):double;  { TGLFunction0 }
  function integrandSQ_IPC(x,sigval2,xsival,hval,qx,qz2,zphase,xsiz,p8,p9,p10,p11,p12,p13:double):double;   { TGLFunction }
  function integrandSQre_IPC(x,sigval2,xsival,hval,qx,qz2re,qz2im,zphase,xsiz,p9,p10,p11,p12,p13:double):double; { TGLFunction }
  function integrandSQim_IPC(x,sigval2,xsival,hval,qx,qz2re,qz2im,zphase,xsiz,p9,p10,p11,p12,p13:double):double; { TGLFunction }
  function Fre(x,sigvali,sigvalj,xsivali,xsivalj,xsizval,hvali,hvalj,zi,zj,qx,qq2r,qq2i,res:double):double; { TGLFunction }
  function Fim(x,sigvali,sigvalj,xsivali,xsivalj,xsizval,hvali,hvalj,zi,zj,qx,qq2r,qq2i,res:double):double; { TGLFunction }
  { Percus - Yevick }
  function integrandPercusYevick(x,Q,sigma,etha,alpha,beta,gamma:double):double;  { TGLFunction0 }


IMPLEMENTATION

  USES    SysUtils,Math  {Exception};

  TYPE
    EQuadratureError = CLASS(Exception);

  CONST
    PointCount:  ARRAY[1..7] OF INTEGER = (2, 3, 4, 5, 6, 10, 15);
    PointKey  :  ARRAY[1..8] OF INTEGER = (1, 2, 4, 6, 9, 12, 17, 25);
    Deg=pi/180;

{n-Point Formula}
    Root:  ARRAY[1..24] OF DOUBLE =
 {2}    (0.577350269189626,

 {3}     0.0,               0.774596669241483,

 {4}     0.339981043584856, 0.861136311594053,

 {5}     0.0,               0.538469310105683, 0.906179845938664,

 {6}     0.238619186083197, 0.661209386466265, 0.932469514203152,

{10}     0.148874338981631, 0.433395394129247, 0.679409568299024,
         0.865063366688985, 0.973906528517172,

{15}     0.0,               0.201194093997435, 0.394151347077563,
         0.570972172608539, 0.724417731360170, 0.848206583410427,
         0.937273392400706, 0.987992518020485);

    Weight:  ARRAY[1..24] OF DOUBLE =
 {2}    (1.0,

 {3}     0.888888888888889, 0.555555555555556,

 {4}     0.652145154862546, 0.347854845137454,

 {5}     0.568888888888889, 0.478628670499366, 0.236926885056189,

 {6}     0.467913934572691, 0.360761573048139, 0.171324492379170,

{10}     0.295524224714753, 0.269266719309996, 0.219086362515982,
         0.149451349150581, 0.066671344308688,

{15}     0.202578241925561, 0.198431485327111, 0.186161000115562,
         0.166269205816994, 0.139570677926154, 0.107159220467172,
         0.070366047488108, 0.030753241996117);

{ ========================================== }
{ fonctions de correlation in-plane pures }
function integrandSQ(x,sigval2,xsival,hval,qx,qz2,p6:double):double;
var fcor:double;
begin
fcor:=sigval2*exp(-Power(x/xsival,2*hval));
result:=cos(x*qx)*(exp(qz2*FCor)-1);
end;
function integrandSQre(x,sigval2,xsival,hval,qx,qz2re,qz2im:double):double;
var fcor:double;
begin
fcor:=sigval2*exp(-Power(x/xsival,2*hval));
result:=cos(x*qx)*(cos(qz2im*Fcor)*exp(qz2re*FCor)-1);
end;
function integrandSQim(x,sigval2,xsival,hval,qx,qz2re,qz2im:double):double;
var fcor:double;
begin
fcor:=sigval2*exp(-Power(x/xsival,2*hval));
result:=cos(x*qx)*(sin(qz2im*Fcor)*exp(qz2re*FCor)); { il n' y a pas de -1 !!! }
end;
{ ========================================== }
{ fonctions de correlation in-plane pures }
function integrandSQ_IPC(x,sigval2,xsival,hval,qx,qz2,zphase,xsiz,p8,p9,p10,p11,p12,p13:double):double;
var fcor:double;
begin
fcor:=sigval2*exp(-Power(x/xsival,2*hval))/exp(Abs(zphase)/xsiz);
result:=cos(x*qx)*(exp(qz2*FCor)-1);
end;
function integrandSQre_IPC(x,sigval2,xsival,hval,qx,qz2re,qz2im,zphase,xsiz,p9,p10,p11,p12,p13:double):double;
var fcor:double;
begin
fcor:=sigval2*exp(-Power(x/xsival,2*hval))/exp(Abs(zphase)/xsiz);
result:=cos(x*qx)*(cos(qz2im*Fcor)*exp(qz2re*FCor)-1);
end;
function integrandSQim_IPC(x,sigval2,xsival,hval,qx,qz2re,qz2im,zphase,xsiz,p9,p10,p11,p12,p13:double):double;
var fcor:double;
begin
fcor:=sigval2*exp(-Power(x/xsival,2*hval))/exp(Abs(zphase)/xsiz);
result:=cos(x*qx)*(sin(qz2im*Fcor)*exp(qz2re*FCor)); { il n' y a pas de -1 !!! }
end;
{ ========================================== }
{ ========================================== }
{ fonctions de correlation AVEC CORRELATIONs ENTRE INTERFACES  }
function FCorr(x,sigvali,sigvalj,xsivali,xsivalj,xsizval,hvali,hvalj,zi,zj:double):double;
var gi,hi,gj,hj:double;
begin
gi:=x/xsivali;
hi:=2*hvali;
gj:=x/xsivalj;
hj:=2*hvalj;
{result:=sqr(sigval)/(exp(Power(gi,hi)*exp(Abs(zi-zj)/xsizval)));}
result:=sigvali*sigvalj/exp(Power(gi,hi))/exp(Abs(zi-zj)/xsizval);
end;

function Fre(x,sigvali,sigvalj,xsivali,xsivalj,xsizval,hvali,hvalj,zi,zj,qx,qq2r,qq2i,res:double):double;
var fcor,expfcor:double;
begin
fcor:=FCorr(x,sigvali,sigvalj,xsivali,xsivalj,xsizval,hvali,hvalj,zi,zj);
expfcor:=exp(qq2r*Fcor);
result:=cos(x*qx)*(cos(qq2i*FCor)*expFCor-1);
end;

function Fim(x,sigvali,sigvalj,xsivali,xsivalj,xsizval,hvali,hvalj,zi,zj,qx,qq2r,qq2i,res:double):double;
var fcor,expFcor:double;
begin
fcor:=FCorr(x,sigvali,sigvalj,xsivali,xsivalj,xsizval,hvali,hvalj,zi,zj);
expfcor:=exp(qq2r*Fcor);
result:=cos(x*qx)*sin(qq2i*FCor)*expFCor;
end;

function integrandPercusYevick(x,Q,sigma,etha,alpha,beta,gamma:double):double;  { TGLFunction0 }
begin
  result:=x*x*sin(x*Q*sigma)/(x*Q*sigma)*(alpha+beta*x+gamma*x*x*x);
end;

{ ========================================== }

Procedure GaussLegendreQuadrature(CONST A, B:  DOUBLE;
                                  CONST N:  INTEGER;
                                  CONST F:  TGLFunction;
                                  CONST par1,par2,par3,par4,par5,par6,par7,par8,par9,par10,par11,par12,par13:double;
                                  var result:double);

    VAR
      C     :  DOUBLE;
      D     :  DOUBLE;
      found :  BOOLEAN;
      i     :  INTEGER;
      j     :  INTEGER;
      jFirst:  INTEGER;
      jLast :  INTEGER;
      Sum   :  DOUBLE;

  BEGIN
    i := 1;
    REPEAT
      found := (N = PointCount[i]);
      IF   NOT found
      THEN INC(i)
    UNTIL found OR (i > High(PointCount));

    IF   NOT found
    THEN RAISE EQuadratureError.Create('GaussLegendre:  ' +
               IntToStr(N) + '-point formula not implemented.');

    jFirst := PointKey[i];
    jLast  := PointKey[i+1]-1;
    C := (B - A) / 2.0;
    D := (B + A) / 2.0;

    Sum := 0;
    FOR j := jFirst TO jLast DO
    BEGIN
      IF   Root[j] = 0.0
      THEN
      begin
      Sum:= Sum + Weight[j]*F(D,par1,par2,par3,par4,par5,par6,par7,par8,par9,par10,par11,par12,par13)
      end
      ELSE
      begin
      Sum:= Sum + Weight[j]*(F( Root[j]*C + D,par1,par2,par3,par4,par5,par6,par7,par8,par9,par10,par11,par12,par13)+F(-Root[j]*C + D,par1,par2,par3,par4,par5,par6,par7,par8,par9,par10,par11,par12,par13))
      end;
    END;
  RESULT := C * SUM
  END {GaussLegendreQuadrature};

{ version allégée de Gauss-Legendre          }
{ on limite le nombre de parametres a passer }
Procedure GaussLegendreQuadrature0(CONST A, B:  DOUBLE;
                                  CONST N:  INTEGER;
                                  CONST F:  TGLFunction0;
                                  CONST par1,par2,par3,par4,par5,par6:double;
                                  var result:double);
    VAR
      C     :  DOUBLE;
      D     :  DOUBLE;
      found :  BOOLEAN;
      i     :  INTEGER;
      j     :  INTEGER;
      jFirst:  INTEGER;
      jLast :  INTEGER;
      Sum   :  DOUBLE;

  BEGIN
    i := 1;
    REPEAT
      found := (N = PointCount[i]);
      IF   NOT found
      THEN INC(i)
    UNTIL found OR (i > High(PointCount));

    IF   NOT found
    THEN RAISE EQuadratureError.Create('GaussLegendre:  ' +
               IntToStr(N) + '-point formula not implemented.');

    jFirst := PointKey[i];
    jLast  := PointKey[i+1]-1;
    C := (B - A) / 2.0;
    D := (B + A) / 2.0;

    Sum := 0;
    FOR j := jFirst TO jLast DO
    BEGIN
      IF   Root[j] = 0.0
      THEN
      begin
      Sum:= Sum + Weight[j]*F(D,par1,par2,par3,par4,par5,par6)
      end
      ELSE
      begin
      Sum:= Sum + Weight[j]*(F( Root[j]*C + D,par1,par2,par3,par4,par5,par6)+F(-Root[j]*C + D,par1,par2,par3,par4,par5,par6))
      end;
    END;
  RESULT := C * SUM
  END {GaussLegendreQuadrature};

  { intégration de foncionts oscillantes }
  { l'utilisation de GL standard même d'ordre élevé ne donne pas de bons résultats }
  { INPUT : on rajoute en entrée, le facteur A de la fonction oscillante sin ( C x ) par exemple }
  { on fait des intégrations par tranches de demi-pétriodes d'oscillations = Pi/C  }
  Procedure GaussLegendreOsc(CONST A, B, C:  DOUBLE;
                                    CONST N:  INTEGER;
                                    CONST F:  TGLFunction0;
                                    CONST par1,par2,par3,par4,par5,par6:double;
                                    var result:double);
  var int,intervalle,A1,B1:double; { resultat intermédiare de l'intégration }
  begin
    result:=0;
    intervalle:=pi/C;
    A1:=A;
    while (A1+intervalle)<B do
    begin B1:=A1+intervalle;
          GaussLegendreQuadrature0(A1,B1,N,F,par1,par2,par3,par4,par5,par6,int);
          result:=result+int;
          A1:=B1; end;
    B1:=B;
    GaussLegendreQuadrature0(A1,B1,N,F,par1,par2,par3,par4,par5,par6,int);
    result:=result+int;
  end;


END.
