{-----------------------------------------------------------------------------
 Unit Name: Functions
 Author:    kaetzel
 Purpose:   In this Unit all the mathematical functions that are needed for the
            calculation of the scattering parameters are defined
-----------------------------------------------------------------------------}

unit Functions;

interface

uses Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
     Dialogs, StdCtrls,JvDSADialogs,JclMath,Math,JclUnitConv,JclBase,
     Gaunt;

type
        //Type of extended in which the Gaunt coefficients are saved
  Gaunt=packed record
    GC:extended;
  end;
  TGaunt=file of Gaunt;
        //Integer Type in which the position numbers are saved in the
        //PosMatrix.dat files
  pos=packed record
    position_value:integer;
  end;
  TPosition=file of pos;
        //Error Class for the check whether the calculated
        //Gaunt coefficients are sufficient
  EExceedingError=class(Exception);
  EFileNotFoundError=class(Exception);


        {Public Functions}

function  Overlap:bool;
function  CenterOfGravity (Coordinate:array of float):float;
procedure PointOfOriginTransformation (xs,ys,zs:real);
procedure CartesianToSpherical;
procedure SizeParameter(Diameter:array of real);
procedure WiscombeCriterion(x:array of extended);
procedure RelativeRefractiveIndex(NpReal,NpImag:array of real;Nmed:real);
function  Absorption:Boolean;
procedure DistanceMatrix;
procedure Matrix_Construction;
procedure CalcVector_amnl_bmnl;
function  ExtinctionCoefficient:extended;
function  ScatteringCoefficient:extended;
function  AsymmetryParameter(Csca:extended):extended;
procedure GetRotationAngles(const i:integer; const NoO:integer;
                            out theta_x:float; out theta_y:float);
procedure RotateAgglomerate (X0coord,Y0coord,Z0coord:array of float;
                           theta_x:float; theta_y:float);
function  AverageValue (values:array of extended):extended;

        {all global variables are declared in main.pas}

implementation

uses Main, NumSpheres, MieParam, ComplexNumbers;


{-----------------------------------------------------------------------------
  Procedure: Distance
  Purpose:   Calculates the shortest distance between two points in a 3D-space
  Comments:
-----------------------------------------------------------------------------}

function Distance(x1,y1,z1,x2,y2,z2:real):real;
begin
 Result:=Sqrt(Sqr(x2-x1)+Sqr(y2-y1)+Sqr(z2-z1));
end;


{-----------------------------------------------------------------------------
  Procedure: Overlap
  Purpose:   Controls whether the spheres of a defined agglomerate do intersect
  Comments:  The shortest distance between two spheres is calculated and then
             compared to the sum of the radii of the two spheres
-----------------------------------------------------------------------------}

function Overlap:bool;
var
  i,j:integer;
  dist:real;
begin
  try
    for i:=0 to (Nl-2) do
      begin
        for j:=i+1 to (Nl-1) do
          begin
            dist:=Distance(Xcoord[i],Ycoord[i],Zcoord[i],
                           Xcoord[j],Ycoord[j],Zcoord[j]);
            if ((dia[i]/2+dia[j]/2)-dist)>1e-6 then abort;
          end;
      end;
    Result:=False;
  except
    Result:=True;
  end;
end;


{-----------------------------------------------------------------------------
  Procedure: CenterOfGravity
  Purpose:   This function calculates the geometrical centre of gravity for a
             given agglomerate to ensure that the program works with least
             distances
  Comments:  Saves calculation time for whole process
-----------------------------------------------------------------------------}

function CenterOfGravity (Coordinate:array of float):float;
var
  i:integer;
  sum_dia,sum_weight:real;
begin
  sum_dia:=0;
  sum_weight:=0;
  for i:=0 to Nl-1 do
    sum_dia:=sum_dia+dia[i];
  for i:=0 to Nl-1 do
    sum_weight:=sum_weight+dia[i]*Coordinate[i];
  Result:=sum_weight/sum_dia;
end;


{-----------------------------------------------------------------------------
  Procedure: PointOfOriginTransformation
  Purpose:   This procedure moves the origin of the coordinate system which was
             defined in the agglomerate configuration to the geometrical centre
             of gravity
  Comments:
-----------------------------------------------------------------------------}

procedure PointOfOriginTransformation (xs,ys,zs:real);
var
  i:integer;
begin
  for i:=0 to Nl-1 do
    begin
      Xcoord[i]:=Xcoord[i]-xs;
      Ycoord[i]:=Ycoord[i]-ys;
      Zcoord[i]:=Zcoord[i]-zs;
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: CartesianToSpherical
  Purpose:   Calculates the spherical coordinates of the spheres of the
             agglomerate in the centered coordinate system
  Comments:  Uses embedded function of JclUnitConv
-----------------------------------------------------------------------------}

procedure CartesianToSpherical;
var
  i:integer;
begin
  for i:=0 to Nl-1 do
    begin
      CartesianToSpheric(Xcoord[i],Ycoord[i],Zcoord[i],Rcoord[i],Phicoord[i],
                         Thetacoord[i]);
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: SizeParameter
  Purpose:   Calculates the dimensionless size parameter of each sphere as a
             preparation for the calculation of the Mie coefficients
  Comments:
-----------------------------------------------------------------------------}

procedure SizeParameter(Diameter:array of real);
var
  i:integer;
begin
  for i:=0 to Nl-1 do
    x[i]:=Pi*dia[i]*1000/(WLength/RefrIndMed);
        //Factor 1000 because Wavelength is given in nm and dia in micron
end;


{-----------------------------------------------------------------------------
  Procedure: WiscombeCriterion
  Purpose:   Calculates the stop criterion of the series expansion in Mie
             calculations and also for the multisphere solution
  Comments:  This parameter may not exceed 20 because the number of calculated
             Gaunt coefficients will not suffice for further computations
-----------------------------------------------------------------------------}

procedure WiscombeCriterion(x:array of extended);
var
  i:integer;
begin
  for i:=0 to Nl-1 do
    begin
      WiscombeNumber[i]:=Ceil(x[i]+4*Power(x[i],(1/3))+2);
      if WiscombeNumber[i]>20 then
        Raise EExceedingError.Create('Warning! Wiscombe Parameter of Sphere '+
              IntToStr(i+1)+' exceeds 20! For this case no Gaunt coefficients'+
              ' have been calculated! Please redefine agglomerate or calculate'+
              ' the required coefficients with MSS_DividedCalc.exe!');
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: RelativeRefractiveIndex
  Purpose:   Computes the relative refractive index for each sphere as it is
             needed for the calculation of the Mie parameters
  Comments:
-----------------------------------------------------------------------------}

procedure RelativeRefractiveIndex(NpReal,NpImag:array of real;Nmed:real);
var
  i:integer;
begin
  for i:=0 to Nl-1 do
    begin
      RelRefrIndex[i,0]:=NpReal[i]/Nmed;
      RelRefrIndex[i,1]:=NpImag[i]/Nmed;
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: Absorption
  Purpose:   Checks whether there is any absorption in the agglomerate
  Comments:
-----------------------------------------------------------------------------}

function Absorption:Boolean;
var
  i:integer;
begin
  try
    for i:=0 to Nl-1 do
      begin
        if RefrIndexImag[i]<>0 then abort;
      end;
    Result:=False;
  except
    Result:=True;
  end;
end;


{-----------------------------------------------------------------------------
  Procedure: PositionRead
  Purpose:   Reads the required Position of the Gaunt coefficients in their
             files.
  Comments:  The files have to be in the same path as the executing program.
-----------------------------------------------------------------------------}

function PositionRead(n,m,nu,mu,p:integer;var Filename:TPosition):integer;
var
  Data:pos;
  loi,DataPointer,sum,i:integer;
        //loi=length of interesting area
  aidv,lowerbound,upperbound,pmax:integer;
begin
        //fixing the length of equal m
  loi:=FileSize(Filename) div (2*n+1);
        //this gives zero for m=-n
  aidv:=n+m;
        //first dataset with value m
  lowerbound:=loi*aidv;
        //the value for nu is not directly calculable so one has to sum up all
        //the values which are not of interest in the actual purpose
  sum:=0;
  if nu=1 then
    sum:=0
  else
    begin
      for i:=1 to (nu-1) do
        sum:=sum+i;
    end;
        //first dataset with value nu
  lowerbound:=lowerbound+(n+1)*(2*sum+nu-1);
        //last dataset with value nu
  upperbound:=lowerbound+(n+1)*(2*nu+1)-1;
        //fixing length of equal mu
  loi:=(upperbound+1-lowerbound) div (2*nu+1);
        //gives zero for mu=-nu
  aidv:=nu+mu;
        //first dataset with value mu
  lowerbound:=lowerbound+loi*aidv;
        //highest possible value of p
  pmax:=n+nu;
        //Pointer to the required Data in the file
  DataPointer:=lowerbound+(pmax-p)div 2;
        //set the filepointer to the position
  seek(Filename,DataPointer);
  read(Filename,Data);
  Result:=Data.position_value;
end;


{-----------------------------------------------------------------------------
  Procedure: GauntCoeffRead
  Purpose:   Reads the Gaunt Coefficients in the already calculated files
  Comments:  Unfortunately the order of these is not recursive to derive
             directly the position in the file, therefore another file which
             contains the position numbers must be read (PosMatrix*.dat)
-----------------------------------------------------------------------------}

function GauntCoeffRead(n,m,nu,mu,p:integer;var FilenameGaunt:TGaunt;
                        var FilenamePos:TPosition):extended;
var
  Datapointer:integer;
  Data:Gaunt;
begin
  if p<abs(m+mu) then
      Result:=0
                //This value contains several Gaunt Coefficients which are
                //called from TranslationCoeff_B but are not included in
                //the files, furthermore this way is faster
  else
    begin
                //calls the position read
      Datapointer:=PositionRead(n,m,nu,mu,p,FilenamePos);
                //seeks the position in the file
      seek(FilenameGaunt,Datapointer);
      read(FilenameGaunt,Data);
      Result:=Data.GC;
      if Result>1E100 then
        ShowMessage('WARNING! Gaunt Coefficient a('+IntToStr(m)+','+IntToStr(n)+
                    ','+IntToStr(mu)+','+IntToStr(nu)+','+IntToStr(p)+
                    ') is >1E100!!'+
                    ' Please Verify with MSS_Dividedcalc.exe!',dckActiveForm);
      if ((Result>0) and (Result<1E-100)) then
        ShowMessage('WARNING! Gaunt Coefficient a('+IntToStr(m)+','+IntToStr(n)+
                    ','+IntToStr(mu)+','+IntToStr(nu)+','+IntToStr(p)+
                    ') is <1E-100!!'+
                    ' Please Verify with MSS_Dividedcalc.exe!',dckActiveForm);
                                //These two messages have been included to check
                                //all calculated files for computation errors
                                //if this message occurs, normally the file has
                                //to be recalculated
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: Pochhammer
  Purpose:   Calculates Pochhammer's Symbol as it is defined in the Abramowitz
  Comments:
-----------------------------------------------------------------------------}

function Pochhammer (z,n:integer):extended;
begin
  Result:=Factorials(z+n-1)/Factorials(z-1);
end;


{-----------------------------------------------------------------------------
  Procedure: DistanceMatrix
  Purpose:   Calculates the Distance Matrix which contains the spherical
             coordinates of each sphere when the centre of the coordinate
             system is the centre of the j-th sphere
  Comments:
-----------------------------------------------------------------------------}

procedure DistanceMatrix;
var
  xjl,yjl,zjl:float;
  l,j:integer;
begin
  for j:=1 to Nl do
    begin
      for l:=1 to Nl do
        begin
          xjl:=Xcoord[l-1]-Xcoord[j-1];
          yjl:=Ycoord[l-1]-Ycoord[j-1];
          zjl:=Zcoord[l-1]-Zcoord[j-1];
          CartesianToSpheric(xjl,yjl,zjl,
                             djl[j-1,l-1],phijl[j-1,l-1],thetajl[j-1,l-1]);
        end;
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: VectorTranslationCoeff_A
  Purpose:   Calculates the vector translation coefficient A for Vector
             Spherical Harmonics
  Comments:  It is already defined as a complex function to avoid splitting the
             whole calculations in Real and Imaginary parts.
             Refers to Xu, Yu-Lin J.Comp.Phys. 139, 137-165 (1998)
             Xu, Yu-Lin www.astro.ufl.edu/~xu, description.pdf Eq. 15a
             and Xu, Yu-Lin Appl. Opt. 34, 4573-4588 (1995)
-----------------------------------------------------------------------------}

function VectorTranslationCoeff_A (m,n,mu,nu,l,j:integer; tilde:bool;
                                   var Filename:TGaunt;
                                   var Fileposition:TPosition):TComplexCart;
var
  C0,sumfactor,c:extended;
  p,q,qmax:integer;
  i:TComplexCart;
  FirstPart:TComplexPolar;
  FirstPartCart:TComplexCart;
  sum:TComplexCart;
  aq,FK,SK:float;
begin
  if (tilde and (j=l)) then
    begin
      if ((m=mu) and (n=nu)) then
        Result:=RealOne
      else
        Result:=ComplexZeroCart;
      exit;
    end;
  C0:=0.5*Sqrt(((2*n+1)*(2*nu+1)*Factorials(n+m)*Factorials(nu-mu))/
          (n*(n+1)*nu*(nu+1)*Factorials(n-m)*Factorials(nu+mu)));
  FirstPart.r:=C0;
  FirstPart.arg:=(mu-m)*phijl[j-1,l-1];
  FirstPart:=ComplexRealMul(FirstPart,Power(-1,m));
  qmax:=Min(n,Min(nu,Trunc((n+nu-abs(mu-m))/2)));
  sum:=ComplexZeroCart;
  for q:=0 to qmax do
    begin
      p:=n+nu-2*q;
      i:=ImaginaryIntegerPower(p);
      aq:=GauntCoeffRead(n,-m,nu,mu,p,Filename,FilePosition);
      c:=cos(thetajl[j-1,l-1]);
      if abs(c)<0.000000000000001 then
        c:=0;
      if tilde then
        begin
          SpherBes(p,Wavenumber*djl[j-1,l-1],FK,SK);
          sumfactor:=(n*(n+1)+nu*(nu+1)-p*(p+1))*aq*
                      AssociatedLegendre(p,(mu-m),c)*FK;
        end
      else
        begin
          sumfactor:=(n*(n+1)+nu*(nu+1)-p*(p+1))*aq*
                      AssociatedLegendre(p,(mu-m),c);
          i:=ComplexMul(i,(SpherHankFK(p,Wavenumber*djl[j-1,l-1])));
        end;
      i:=ComplexRealMul(i,sumfactor);
      sum:=ComplexAdd(sum,i);
    end;
  FirstPartCart:=ComplexPolarToCart(FirstPart);
  Result:=ComplexMul(FirstPartCart,sum);
end;


{-----------------------------------------------------------------------------
  Procedure: Factor_Ap
  Purpose:   Auxiliary variable for the calculation of the vector translation
             coefficient B
  Comments:
-----------------------------------------------------------------------------}

function Factor_Ap (n,m,nu,mu,p:integer):float;
begin
  Result:=-p*(p-1)*(m+mu)+(m-mu)*(n-nu)*(n+nu+1);
end;


{-----------------------------------------------------------------------------
  Procedure: Factor_Alpha
  Purpose:   Auxiliary variable for VTC B
  Comments:
-----------------------------------------------------------------------------}

function Factor_Alpha (n,nu,p:integer):float;
begin
  Result:=(Power(p,2)-Power(n-nu,2))*(Power(p,2)-Power(n+nu+1,2))/
          (4*Power(p,2)-1);
end;


{-----------------------------------------------------------------------------
  Procedure: VectorTranslationCoeff_B
  Purpose:   Calculates the vector translation coefficient B for Vector
             Spherical Harmonics
  Comments:  It is already defined as a complex function to avoid splitting the
             whole calculations in Real and Imaginary parts.
             Refers to Xu, Yu-Lin J.Comp.Phys. 139, 137-165 (1998)
             Xu, Yu-Lin www.astro.ufl.edu/~xu, description.pdf Eq. 15b
             and Xu, Yu-Lin Appl. Opt. 34, 4573-4588 (1995)
-----------------------------------------------------------------------------}

function VectorTranslationCoeff_B (m,n,mu,nu,l,j:integer; tilde:bool;
                                   var Filename:TGaunt;
                                   var FilePosition:TPosition):TComplexCart;
var
  C0,sumfactor,c:extended;
  p,q,qmax,p1,p2:integer;
  i:TComplexCart;
  FirstPart:TComplexPolar;
  FirstPartCart:TComplexCart;
  sum:TComplexCart;
  bq,FK,SK:float;
begin
  if (tilde and (j=l)) then
    begin
      Result:=ComplexZeroCart;
      exit;
    end;
  C0:=0.5*Sqrt(((2*n+1)*(2*nu+1)*Factorials(n+m)*Factorials(nu-mu))/
          (n*(n+1)*nu*(nu+1)*Factorials(n-m)*Factorials(nu+mu)));
  FirstPart.r:=C0;
  FirstPart.arg:=(mu-m)*phijl[j-1,l-1];
  FirstPart:=ComplexRealMul(FirstPart,Power(-1,m));
  qmax:=Min(n,Min(nu,Trunc((n+nu+1-abs(mu-m))/2)));
  sum:=ComplexZeroCart;
  for q:=1 to qmax do
    begin
      p:=n+nu-2*q;
      if ((Factor_Ap(n,m,nu,mu,p+2)=0)and(Factor_Ap(n,m,nu,mu,p+4)=0)) then
        begin
          Result:=ComplexZeroCart;
          exit;
        end;
      p1:=p+m-mu;
      p2:=p-m+mu;
      if Factor_Ap(n,m,nu,mu,p+2)=0 then
        begin
          bq:=(2*p+3)/(Factor_Ap(n,m,nu,mu,p+4)*(p+3)*(p1+2))*
              ( (Factor_Ap(n,m,nu,mu,p+3)*Factor_Ap(n,m,nu,mu,p+4)+(p+2)*
              (p+4)*(p1+3)*(p2+3)*Factor_Alpha(n,nu,p+3))*
              GauntCoeffRead(n,-m,nu,mu,n+nu-2*(q-1),Filename,Fileposition)
              -(p+2)*(p+3)*(p2+3)*(p2+4)*Factor_Alpha(n,nu,p+4)*
              GauntCoeffRead(n,-m,nu,mu,n+nu-2*(q-2),Filename,FilePosition) );
        end
      else
        begin
          bq:=(2*p+3)/Factor_Ap(n,m,nu,mu,p+2)*((p+2)*(p1+1)*
              Factor_Alpha(n,nu,p+1)*GauntCoeffRead(n,-m,nu,mu,n+nu-2*q,
              Filename,Fileposition)-(p+1)*(p2+2)*Factor_Alpha(n,nu,p+2)*
              GauntCoeffRead(n,-m,nu,mu,n+nu-2*(q-1),Filename,
              FilePosition));
        end;
      i:=ImaginaryIntegerPower(p+1);
      c:=cos(thetajl[j-1,l-1]);
      if abs(c)<0.000000000000001 then
        c:=0;
      if tilde then
        begin
          SpherBes(p+1,Wavenumber*djl[j-1,l-1],FK,SK);
          sumfactor:=bq*AssociatedLegendre((p+1),(mu-m),c)*FK;
        end
      else
        begin
          sumfactor:=bq*AssociatedLegendre(p+1,mu-m,c);
          i:=ComplexMul(i,(SpherHankFK(p+1,Wavenumber*djl[j-1,l-1])));
        end;
      i:=ComplexRealMul(i,sumfactor);
      sum:=ComplexAdd(sum,i);
    end;
  FirstPartCart:=ComplexPolarToCart(FirstPart);
  Result:=ComplexMul(FirstPartCart,sum);
end;


{-----------------------------------------------------------------------------
  Procedure: PrimExpansionCoeff_p
  Purpose:   Calculates the primary expansion coefficient of the incident wave p
  Comments:  Refers to Xu, Yu-Lin www.astro.ufl.edu/~Xu, Eq. 8 and 10
-----------------------------------------------------------------------------}

function PrimExpansionCoeff_p (m,n:integer; beta,zj:float):TComplexCart;
var
  ResultPolar,aidv:TComplexPolar;
begin
  if abs(m)<>1 then
    Result:=ComplexZeroCart
  else
    begin
      if m=1 then
        begin
          ResultPolar.r:=Sqrt(2*n+1)/2;
          ResultPolar.arg:=-beta;
        end
      else
        begin
          ResultPolar.r:=Sqrt(2*n+1)/2;
          ResultPolar.arg:=beta;
          ResultPolar:=ComplexRealMul(ResultPolar,-1);
        end;
      aidv.r:=1;
      aidv.arg:=Wavenumber*zj;
      ResultPolar:=ComplexMul(ResultPolar,aidv);
      Result:=ComplexPolarToCart(ResultPolar);
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: PrimExpansionCoeff_q
  Purpose:   Calculates the primary expansion coefficient of the incident wave q
  Comments:  Refers to Xu, Yu-Lin www.astro.ufl.edu/~Xu Eq. 8 and 10
-----------------------------------------------------------------------------}

function PrimExpansionCoeff_q (m,n:integer; beta,zj:float):TComplexCart;
var
  ResultPolar,aidv:TComplexPolar;
begin
  if abs(m)<>1 then
    Result:=ComplexZeroCart
  else
    begin
      if m=1 then
        begin
          ResultPolar.r:=Sqrt(2*n+1)/2;
          ResultPolar.arg:=-beta;
        end
      else
        begin
          ResultPolar.r:=Sqrt(2*n+1)/2;
          ResultPolar.arg:=beta;
        end;
      aidv.r:=1;
      aidv.arg:=Wavenumber*zj;
      ResultPolar:=ComplexMul(ResultPolar,aidv);
      Result:=ComplexPolarToCart(ResultPolar);
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: Unknowns
  Purpose:   Calculates the number of unknowns in the linear equation system
             for the interactive scattering coefficients a(mn) and b(mn)
             of each sphere
  Comments:
-----------------------------------------------------------------------------}

function Unknowns:integer;
var
  sum,i:integer;
begin
  sum:=0;
  for i:=1 to Nl do
    sum:=sum+2*WiscombeNumber[i-1]*(WiscombeNumber[i-1]+2);
  Result:=sum;
end;


{-----------------------------------------------------------------------------
  Procedure: Matrix_Construction
  Purpose:   Builds the matrix of vector translation coefficients for the
             solution process for the linear equation system Ax=b
             for partial scattering coefficients
  Comments:  ATTENTION!!! THIS COMPUTATION IS MOST TIME AND MEMORY CONSUMING IN
             THE WHOLE PROGRAM BECAUSE IT REQUIRES THE CALCULATION OF ALL
             VECTOR TRANSLATION COEFFICIENTS A AND B WHICH INTEGRATES
             ALL GAUNT COEFFICIENTS
-----------------------------------------------------------------------------}

procedure Matrix_Construction;
var
  Le,h,i,j,l,n,m,nu,mu:integer;
  beta:extended;
  GauntCoeff:array[1..20] of TGaunt;
  PosMatrix:array[1..20] of TPosition;
begin

        {                    File Initialization                       }

  for i:=1 to 20 do
    begin
      if not(FileExists('GauntCoeff_20_'+IntToStr(i)+'.dat')) then
        begin
          MainWindow.Label35.Visible:=False;
          Raise EFileNotFoundError.Create('Error in Matrix_Construction! File with'+
              ' Gaunt Coefficients not found (GauntCoeff_20_'+IntToStr(i)+'.dat)'+
              '. Files must be in same directory as program!');
        end;
                //Reading the files with the Gaunt Coefficients
      AssignFile(GauntCoeff[i],'GauntCoeff_20_'+IntToStr(i)+'.dat');
                //Reading the files with the position numbers of the Gaunt
                //Coeff's in their files because the previous algorithm is not
                //recursive to directly obtain the position of the pointer
      AssignFile(PosMatrix[i],'PosMatrix_20_'+IntToStr(i)+'.dat');
                //open the files
      ReSet(GauntCoeff[i]);
      ReSet(PosMatrix[i]);
    end;
                //Polarization angle Pi/2
  //beta:=PiOn2;
  beta:=0;
  Le:=Unknowns;
                //Length of the vector is equal to number of unknowns
  try
    SetLength(Vector_pq,Le);
    SetLength(Vector_amn_bmn,Le);
    SetLength(Matrix_AB,Le,Le);
  except
    ShowMessage(Format(' Fatal Error!'+
                       '%s Physical Memory not sufficient to calculate problem!'+
                       ' %s ',[#10#13,#10#13] ),dckActiveForm);
  end;

        { Calculation of the resulting vector with primary expansion   }
        { coefficients p(mn) and q(mn)                                 }

  i:=-1;
  for j:=1 to Nl do
    begin
      for n:=1 to WiscombeNumber[j-1] do
        begin
          for m:=-n to n do
            begin
              i:=i+1;
              Vector_pq[i]:=PrimExpansionCoeff_p(m,n,beta,Zcoord[j-1]);
              i:=i+1;
              Vector_pq[i]:=PrimExpansionCoeff_q(m,n,beta,Zcoord[j-1]);
            end;
        end;
    end;

        { Calculation of the coefficient matrix by taking into account }
        { that m,n,l are constant in the rows of the Matrix and        }
        { mu,nu,j are constant in the columns, the values change every }
        { second element because of a change between a(mn) and b(mn)in }
        { every element, therefore the matrix has to be built of small }
        { 2x2 matrixes where m,n,mu,nu,l,j are constant                }

  h:=-2;
  for l:=1 to Nl do
    begin
      for n:=1 to WiscombeNumber[l-1] do
        begin
          for m:=-n to n do
            begin
              h:=h+2;
              i:=0;
              for j:=1 to Nl do
                begin
                  for nu:=1 to WiscombeNumber[j-1] do
                    begin
                      for mu:=-nu to nu do
                        begin
                          if i=h then
                                //This 2x2 matrix contains the diagonal elements
                                //which are the reciprocal of the complex
                                //Mie coefficients
                            begin
                              Matrix_AB[h,i].Re:=MieParameter[j-1,0,n]/
                                              (Sqr(MieParameter[j-1,0,n])
                                              +Sqr(MieParameter[j-1,1,n]));
                              Matrix_AB[h,i].Im:=-MieParameter[j-1,1,n]/
                                              (Sqr(MieParameter[j-1,0,n])
                                              +Sqr(MieParameter[j-1,1,n]));
                              Matrix_AB[h+1,i+1].Re:=MieParameter[j-1,2,n]/
                                                  (Sqr(MieParameter[j-1,2,n])
                                                  +Sqr(MieParameter[j-1,3,n]));
                              Matrix_AB[h+1,i+1].Im:=-MieParameter[j-1,3,n]/
                                                  (Sqr(MieParameter[j-1,2,n])
                                                  +Sqr(MieParameter[j-1,3,n]));
                              Matrix_AB[h,i+1]:=ComplexZeroCart;
                              Matrix_AB[h+1,i]:=ComplexZeroCart;
                              i:=i+2;
                            end
                          else
                            begin
                              if j=l then
                                        //This 2x2 Matrix contains the self-
                                        //interacting coefficients between one
                                        //sphere which are exactly zero,
                                        //since all elements of the matrix
                                        //are complex this is a complex zero
                                begin
                                  Matrix_AB[h,i]:=ComplexZeroCart;
                                  Matrix_AB[h+1,i+1]:=ComplexZeroCart;
                                  Matrix_AB[h+1,i]:=ComplexZeroCart;
                                  Matrix_AB[h,i+1]:=ComplexZeroCart;
                                  i:=i+2;
                                end
                              else
                                        //This 2x2 matrix contains the vector
                                        //translation coefficients A and B in
                                        //the following manner
                                        //   A(m,n,mu,nu,l,j) | B(m,n,mu,nu,l,j)
                                        //   B(m,n,mu,nu,l,j) | A(m,n,mu,nu,l,j)
                                        //Here the handling of 2x2 matrices
                                        //saves about half the calculation time
                                begin
                                  Matrix_AB[h,i]:=VectorTranslationCoeff_A(m,n,
                                                  mu,nu,l,j,false,
                                                  GauntCoeff[n],PosMatrix[n]);
                                  Matrix_AB[h,i+1]:=VectorTranslationCoeff_B(m,
                                                  n,mu,nu,l,j,false,
                                                  GauntCoeff[n],PosMatrix[n]);
                                  Matrix_AB[h+1,i+1]:=Matrix_AB[h,i];
                                  Matrix_AB[h+1,i]:=Matrix_AB[h,i+1];
                                  i:=i+2;
                                end;
                            end;
                        end;
                    end;
                end;
            end;
        end;
    end;
  for i:=1 to 20 do
    begin
      CloseFile(GauntCoeff[i]);
      CloseFile(PosMatrix[i]);
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: CalcVector_amnl_bmnl
  Purpose:   Calculates the vector of interactive scattering coefficients amn(l)
             and bmn(l) which are assigned to the vector translation coeff's
             ~A and ~B with the Bessel function of the first kind as generating
             function
  Comments:  COMPUTATION IS ALSO TIME CONSUMING BUT LESS MEMORY CONSUMING AS
             THE ELEMENTS OF THE VECTOR CAN BE CALCULATED DIRECTLY
-----------------------------------------------------------------------------}

procedure CalcVector_amnl_bmnl;
var
  h,i,j,l,m,n,mu,nu:integer;
  sum1,sum2,VTC_A_Tilde,VTC_B_Tilde:TComplexCart;
  GauntCoeff:array[1..20] of TGaunt;
  PosMatrix:array[1..20] of TPosition;
begin
  SetLength(Vector_amnl_bmnl,length(Vector_amn_bmn));
  for i:=1 to 20 do
    begin
      AssignFile(GauntCoeff[i],'GauntCoeff_20_'+IntToStr(i)+'.dat');
                //Reading the files with the Gaunt Coefficients
      AssignFile(PosMatrix[i],'PosMatrix_20_'+IntToStr(i)+'.dat');
                //Reading the files with the position numbers of the Gaunt
                //Coeff's in their files
      ReSet(GauntCoeff[i]);
      ReSet(PosMatrix[i]);
    end;
  h:=0;
  for l:=1 to Nl do
    begin
      for n:=1 to WiscombeNumber[l-1] do
        begin
          for m:=-n to n do
            begin
              sum1:=ComplexZeroCart;
              sum2:=ComplexZeroCart;
              i:=0;
              for j:=1 to Nl do
                begin
                  for nu:=1 to WiscombeNumber[j-1] do
                    begin
                      for mu:=-nu to nu do
                        begin

                          VTC_A_Tilde:=VectorTranslationCoeff_A(m,n,mu,nu,l,j,
                                       true,GauntCoeff[n],PosMatrix[n]);
                          VTC_B_Tilde:=VectorTranslationCoeff_B(m,n,mu,nu,l,j,
                                       true,GauntCoeff[n],PosMatrix[n]);
                          sum1:=ComplexAdd(sum1,ComplexAdd(ComplexMul(
                                VTC_A_Tilde,Vector_amn_bmn[i]),ComplexMul(
                                VTC_B_Tilde,Vector_amn_bmn[i+1])));
                          sum2:=ComplexAdd(sum2,ComplexAdd(ComplexMul(
                                VTC_A_Tilde,Vector_amn_bmn[i+1]),ComplexMul(
                                VTC_B_Tilde,Vector_amn_bmn[i])));
                          i:=i+2;
                        end;
                    end;
                end;
              Vector_amnl_bmnl[h]:=sum1;
              Vector_amnl_bmnl[h+1]:=sum2;
              h:=h+2;
            end;
        end;
    end;
  for i:=1 to 20 do
    begin
      CloseFile(GauntCoeff[i]);
      CloseFile(PosMatrix[i]);
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: ExtinctionCoefficient
  Purpose:   Calculates the extinction coefficient of one orientation of the
             agglomerate
  Comments:
-----------------------------------------------------------------------------}

function ExtinctionCoefficient:extended;
var
  i:integer;
  sum:extended;
begin
  sum:=0;
  for i:=0 to (high(vector_amn_bmn) div 2) do
        sum:=sum+RealPart(ComplexAdd(ComplexMul(Conjugate(vector_pq[2*i]),
                          vector_amn_bmn[2*i]),ComplexMul(Conjugate(
                          vector_pq[2*i+1]),vector_amn_bmn[2*i+1])));
  Result:=4*Pi/(Power(Wavenumber,2))*sum;
end;


{-----------------------------------------------------------------------------
  Procedure: ScatteringCoefficient
  Purpose:   Calculates the scattering coefficient of one orientation of the
             agglomerate
  Comments:
-----------------------------------------------------------------------------}

function ScatteringCoefficient:extended;
var
  i:integer;
  sum:extended;
begin
  sum:=0;
  for i:=0 to (high(vector_amn_bmn) div 2) do
    sum:=sum+RealPart(ComplexAdd(ComplexMul(Conjugate(vector_amn_bmn[2*i]),
                      vector_amnl_bmnl[2*i]),ComplexMul(
                      Conjugate(vector_amn_bmn[2*i+1]),
                      vector_amnl_bmnl[2*i+1])));
  Result:=4*Pi/(Power(Wavenumber,2))*sum;
end;


{-----------------------------------------------------------------------------
  Procedure: AsymmetryParameter
  Purpose:   Calculates the asymmetry parameter <cos theta> for one orientation
             of the agglomerate
  Comments:
-----------------------------------------------------------------------------}

function AsymmetryParameter(Csca:extended):extended;
var
  i,l,m,n:integer;
  sum,f1,f2,f3:extended;
  amnltilde,bmnltilde:TComplexCart;
begin
  sum:=0;
  i:=0;
  for l:=1 to Nl do
    begin
      for n:=1 to WiscombeNumber[l-1] do
        begin
          for m:=-n to n do
            begin
              if ((n=1)or(abs(m)>(n-1))and not (n=WiscombeNumber[l-1])) then
                begin
                  f1:=m/(n*(n+1));
                  f2:=1/(n+1)*Power((n*(n+2)*(n-m+1)*(n+m+1))/((2*n+1)*(2*n+3)),
                                    0.5);
                  amnltilde:=ComplexAdd(ComplexRealMul(Vector_amnl_bmnl[i+1],f1),
                             ComplexRealMul(Vector_amnl_bmnl[i+2*(2*n+2)],f2));
                  bmnltilde:=ComplexAdd(ComplexRealMul(Vector_amnl_bmnl[i],f1),
                             ComplexRealMul(Vector_amnl_bmnl[i+1+2*(2*n+2)],f2));
                end
              else
                begin
                  if n=WiscombeNumber[l-1] then
                    begin
                      f1:=m/(n*(n+1));
                      if abs(m)>(n-1) then
                        begin
                          amnltilde:=ComplexRealMul(Vector_amnl_bmnl[i+1],f1);
                          bmnltilde:=ComplexRealMul(Vector_amnl_bmnl[i],f1);
                        end
                      else
                        begin
                          f3:=1/n*Power(((n-1)*(n+1)*(n-m)*(n+m))/
                              ((2*n-1)*(2*n+1)),0.5);
                          amnltilde:=ComplexAdd(ComplexRealMul(
                                     Vector_amnl_bmnl[i+1],f1),
                                     ComplexRealMul(
                                     Vector_amnl_bmnl[i-2*(2*(n-1)+2)],f3));
                          bmnltilde:=ComplexAdd(ComplexRealMul(
                                     Vector_amnl_bmnl[i],f1),ComplexRealMul(
                                     Vector_amnl_bmnl[i+1-2*(2*(n-1)+2)],f3));
                        end;
                    end
                  else
                    begin
                      f1:=m/(n*(n+1));
                      f2:=1/(n+1)*Power((n*(n+2)*(n-m+1)*(n+m+1))/
                                       ((2*n+1)*(2*n+3)),0.5);
                      f3:=1/n*Power(((n-1)*(n+1)*(n-m)*(n+m))/
                                   ((2*n-1)*(2*n+1)),0.5);
                      amnltilde:=ComplexAdd(ComplexRealMul(Vector_amnl_bmnl[i+1],
                                 f1),ComplexAdd(ComplexRealMul(Vector_amnl_bmnl
                                 [i+2*(2*n+2)],f2),ComplexRealMul(
                                 Vector_amnl_bmnl[i-2*(2*(n-1)+2)],f3)));
                      bmnltilde:=ComplexAdd(ComplexRealMul(Vector_amnl_bmnl[i],f1),
                                 ComplexAdd(ComplexRealMul(Vector_amnl_bmnl
                                 [i+1+2*(2*n+2)],f2),ComplexRealMul(
                                 Vector_amnl_bmnl[i+1-2*(2*(n-1)+2)],f3)));
                    end;
                end;
              sum:=sum+RealPart(ComplexAdd(ComplexMul(Conjugate(
                                Vector_amn_bmn[i]),amnltilde),ComplexMul(
                                Conjugate(Vector_amn_bmn[i+1]),bmnltilde)));
              i:=i+2;
            end;
        end;
    end;
  Result:=4*Pi/(sqr(Wavenumber)*Csca)*sum;
end;


{-----------------------------------------------------------------------------
  Procedure: GetRotationAngles
  Purpose:   Provides the Euler Angles theta_x and theta_y for the rotation
             of the agglomerate to calculate orientation averaged values
  Comments:  Values depend on the chosen number of different orientations
-----------------------------------------------------------------------------}

procedure GetRotationAngles(const i:integer; const NoO:integer;
                            out theta_x:float; out theta_y:float);
begin
  case NoO of
        //Roughest approximation. The agglomerate is rotated only once by 180°
        //around the x-axis. This is a look from the +z-direction.
    2: begin
         theta_x:=Pi;
         theta_y:=0;
       end;
        //6 orientations. The angles have been calculated from an octaeder where
        //two corners are on the z-axis.
    6: begin
         case i of
           1: begin
                theta_x:=Pi;
                theta_y:=0;
              end;
           2: begin
                theta_x:=PiOn2;
                theta_y:=0;
              end;
           3: begin
                theta_x:=-PiOn2;
                theta_y:=0;
              end;
           4: begin
                theta_x:=0;
                theta_y:=PiOn2;
              end;
           5: begin
                theta_x:=0;
                theta_y:=-PiOn2;
              end;
         end;
       end;
        //8 orientations. The angles have been calculated from the corners of a
        //cube where two corners are on the z-axis.
    8: begin
         case i of
           1: begin
                theta_x:=Pi;
                theta_y:=0;
              end;
           2: begin
                theta_x:=Pi;
                theta_y:=-1.23095941734077;
              end;
           3: begin
                theta_x:=0;
                theta_y:=1.23095941734077;
              end;
           4: begin
                theta_x:=-1.18319964013972;
                theta_y:=-0.490882678289311;
              end;
           5: begin
                theta_x:=1.18319964013972;
                theta_y:=-0.490882678289311;
              end;
           6: begin
                theta_x:=-1.95839301345008;
                theta_y:=0.490882678289311;
              end;
           7: begin
                theta_x:=1.95839301345008;
                theta_x:=0.490882678289311;
              end;
         end;
       end;
        //12 orientations. The angles have been calculated from the corners of a
        //ikosaeder where two corners are on the z-axis.
    12: begin
          case i of
            1:  begin
                  theta_x:=Pi;
                  theta_y:=0;
                end;
            2:  begin
                  theta_x:=1.10714871779409;
                  theta_y:=0;
                end;
            3:  begin
                  theta_x:=0.553574358897046;
                  theta_y:=1.01722196789785;
                end;
            4:  begin
                  theta_x:=-1.01722196789785;
                  theta_y:=0.553574358897045;
                end;
            5:  begin
                  theta_x:=0.553574358897046;
                  theta_y:=-1.01722196789785;
                end;
            6:  begin
                  theta_x:=-1.01722196789785;
                  theta_y:=-0.553574358897045;
                end;
            7:  begin
                  theta_x:=2.12437068569194;
                  theta_y:=0.553574358897045;
                end;
            8:  begin
                  theta_x:=-2.58801829469275;
                  theta_y:=1.01722196789785;
                end;
            9:  begin
                  theta_x:=-2.0344439357957;
                  theta_y:=0;
                end;
            10: begin
                  theta_x:=2.12437068569194;
                  theta_y:=-0.553574358897045;
                end;
            11: begin
                  theta_x:=-2.58801829469275;
                  theta_y:=-1.01722196789785;
                end;
          end;
        end;
        //14 orientations. The angles have been calculated from a equally
        //divided sphere
    14: begin
          case i of
            1:  begin
                  theta_x:=Pi;
                  theta_y:=0;
                end;
            2:  begin
                  theta_x:=PiOn2;
                  theta_y:=0;
                end;
            3:  begin
                  theta_x:=-PiOn2;
                  theta_y:=0;
                end;
            4:  begin
                  theta_x:=0;
                  theta_y:=PiOn2;
                end;
            5:  begin
                  theta_x:=0;
                  theta_y:=-PiOn2;
                end;
            6:  begin
                  theta_x:=PiOn4;
                  theta_y:=0.616;
                end;
            7:  begin
                  theta_x:=PiOn4;
                  theta_y:=-0.616;
                end;
            8:  begin
                  theta_x:=-PiOn4;
                  theta_y:=0.616;
                end;
            9:  begin
                  theta_x:=-PiOn4;
                  theta_y:=-0.616;
                end;
            10: begin
                  theta_x:=0.75*Pi;
                  theta_y:=0.616;
                end;
            11: begin
                  theta_x:=0.75*Pi;
                  theta_y:=-0.616;
                end;
            12: begin
                  theta_x:=-0.75*Pi;
                  theta_y:=0.616;
                end;
            13: begin
                  theta_x:=-0.75*Pi;
                  theta_y:=-0.616;
                end;
          end;
        end;
  end;
end;


{-----------------------------------------------------------------------------
  Procedure: RotateAgglomerate
  Purpose:   Calculates the new coordinates of the agglomerate by Euler Rotation.
             Used for Orientation averaging
  Comments:  Rotation is done first around the x-axis and around the y-axis
             afterwards.
-----------------------------------------------------------------------------}

procedure RotateAgglomerate (X0coord,Y0coord,Z0coord:array of float;
                           theta_x:float; theta_y:float);
var
  XOld,YOld,ZOld:extended;
  i:integer;
begin
  for i:=0 to Nl-1 do
    begin
      XOld:=X0coord[i];
      YOld:=Y0coord[i];
      ZOld:=Z0coord[i];
        //Euler Rotation Matrix
      Xcoord[i]:=XOld*cos(theta_y)+YOld*sin(theta_y)*sin(theta_x)+
                 ZOld*sin(theta_y)*cos(theta_x);
      Ycoord[i]:=YOld*cos(theta_x)-ZOld*sin(theta_x);
      Zcoord[i]:=-XOld*sin(theta_y)+YOld*cos(theta_y)*sin(theta_x)+
                 ZOld*cos(theta_y)*cos(theta_x);
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: AverageValue
  Purpose:   Calculates average value of all elements of an array
  Comments:
-----------------------------------------------------------------------------}

function AverageValue (values:array of extended):extended;
var
  sum:extended;
  i:integer;
begin
  sum:=0;
  for i:=0 to high(values) do
    sum:=sum+values[i];
  Result:=sum/(length(values));
end;

end.
