{-----------------------------------------------------------------------------
 Unit Name: MieParam
 Author:    kaetzel
 Purpose:   Contains the calculation of the Mie coefficients an and bn and
            Calculation of the Spherical Bessel and Hankel Functions
 Comments:  Originally written by Crawley (1994)
            Has been adapted for the calculation for many spheres
            Bugs have been fixed and Comments added
-----------------------------------------------------------------------------}

unit MieParam;

interface

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

        {                     public procedures                       }

procedure Miecalc;
procedure SpherBes (n:integer; z:float; var FK,SK:float);
function SpherHankFK(n:integer;z:float):TComplexCart;

        //all global variables are declared in main.pas

implementation

uses Main, Functions;

const
        //Integer bessel function terms:
  BfhMaxTerms = 200;
        //highest order of an,bn required
  nmax=50;

type
  ArrEx = array [0..nmax+200] of extended;
  Yvec=array [0..nmax+100] of extended;
  Yvecptr=^Yvec;
  ArrPtr=^ArrEx;

var
  n1A_r,n1B_r,n1B_i,n3A_r,n3A_i:Yvecptr;
  BfhMaxOrder, BfhMinOrder:integer;
  JaA,JfA,PsiA,chiA:ArrPtr;


{-----------------------------------------------------------------------------
  Procedure: BesselFractional
  Purpose:   Calculation of Bessel functions of fractional order
  Comments:  according to Crawley (1994) and Lenz (1976) by using continued
             fractions.
             BfhMaxOrder is the maximum order of a Bessel function that is
             calculated, because it calculates with both forward and backward
             recursion, this number must be higher than the order wanted to
             provide sufficient precision.
-----------------------------------------------------------------------------}

procedure BesselFractional (z:extended; BfhMaxOrder,BfhMaxTerms:integer;
          var Ja,Jf:ArrPtr);
var
        //Maximum order of the fractional Bessel function
  FractMaxOrder:extended;
  i: integer;
  a,b,Jj,Jjm: ArrPtr;

begin
  new(a);new(b);new(Jj);new(Jjm);

  {simplified expressions for low order bessel functions : }

  {Ja^[0] = J (.5,z) } //Start values
  Ja^[0]:=Sqrt(2/pi/z)*sin(z); //Start values

  {Jf^[0] = Y (.5,z) } //here a bug was fixed, the sign has to be altered
                        //in earlier versions of this unit
  Jf^[0]:=-Sqrt(2/pi/z)*cos(z);

  FractMaxOrder:=BfhMaxOrder+0.5;

  for i:=1 to BfhMaxTerms do
    begin
      a^[i]:=2*(FractMaxOrder+i)/z;
      if odd(i) then a^[i]:=-a^[i];
    end;

  b^[BfhMaxTerms]:=a^[BfhMaxTerms];
  for i:=BfhMaxTerms downto 2 do
    b^[i-1]:=a^[i-1]+1/b^[i];

  {Jj^=ratio of BesselFunctions}
  Jj^[BfhMaxOrder]:=2*FractMaxOrder/z + 1/b^[1];
  for i:=BfhMaxOrder downto 1 do
    Jj^[i-1]:=2*(i-0.5)/z-1/Jj^[i];

  Jjm^[0]:=Jj^[0];
  for i:=1 to BfhMaxOrder do
    Jjm^[i]:=2*(-i+0.5)/z-1/Jjm^[i-1];

  for i:=1 to BfhMaxOrder do
    begin
      Ja^[i]:=Ja^[i-1]/Jj^[i];
      Jf^[i]:=Jf^[i-1]*Jjm^[i];
    end;
  dispose(a);dispose(b);dispose(Jj);dispose(Jjm); //Memory release
end;
        //The values have been tested with tables in Abramowitz/Stegun
        //Handbook of Mathematical functions


{-----------------------------------------------------------------------------
  Procedure: RicattiBesselCalculation
  Purpose:   Calculation of Ricatti-Bessel functions for the size parameter
  Comments:
-----------------------------------------------------------------------------}

procedure RicattiBesselCalculation (sizep:extended);
var
  i:integer;
begin
  for i := 0 to BfhMaxOrder do
  begin
        //JaA is the fractional Bessel function
    psiA^[i]:=sqrt(pi*sizep/2.0)*JaA^[i];
        //JfA is the fractional Neumann function
    chiA^[i]:=-1*sqrt(pi*sizep/2.0)*Jfa^[i];
        //JfA is calculated with alternating sign, this is changed here
    if odd (i) then chiA^[i]:=- chiA^[i];
  end;
end;


{-----------------------------------------------------------------------------
  Procedure: calcn3A
  Purpose:   Calculates eta3_n which is related to the ratio of the derivative
             of the zeta function to the zeta function
  Comments:  see Kerker (1969)
-----------------------------------------------------------------------------}

procedure calcn3A(zr:extended);
{G:=P+iC}
var
  denom:extended;
  i:integer;
begin
  for i:=1 to BfhMaxOrder-1 do
    begin
      denom:=sqr(psiA^[i])+sqr(chiA^[i]);
      n3A_r^[i]:=(psiA^[i]*psiA^[i-1]+chiA^[i]*chiA^[i-1])/denom-i/zr;
      n3A_i^[i]:=(psiA^[i]*chiA^[i-1]-psiA^[i-1]*chiA^[i])/denom;
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: calcn1A
  Purpose:   Calculates eta1_n which is related to the ratio of the derivative
             of the psi function to the psi function
  Comments:  see Kerker (1969)
-----------------------------------------------------------------------------}

procedure calcn1A(zr:extended);
var
  i:integer;
begin
  for i:=1 to BfhMaxOrder-1 do
    n1A_r^[i]:=psiA^[i-1]/psiA^[i]-i/zr;
end;


{-----------------------------------------------------------------------------
  Procedure: SimplifyComplexFraction
  Purpose:     c
	    ------    --->   u + iv
	    a + ib
  Comments:
-----------------------------------------------------------------------------}

procedure SimplifyComplexFraction(a,b,c:extended;var u,v:extended);
var
  den:extended;
begin
  den:=sqr(a)+sqr(b);
  u:=c*a/den;
  v:=c*(-b)/den;
end;


{-----------------------------------------------------------------------------
  Procedure: calcn1B
  Purpose:   calculation of eta1_n with complex Argument
  Comments:
-----------------------------------------------------------------------------}

procedure calcn1B(zr,zi:extended);
var
  FractMaxOrder,br_,bi_,tempr,tempi:extended;
  i:integer;
  ar_,ai_,Jr,Ji:ArrPtr;
begin
  new(ar_);new(ai_);new(Jr);new(Ji);
  FractMaxorder := BfhMaxOrder+0.5;
  for i:=1 to BfhMaxTerms do
    begin
      SimplifyComplexFraction(zr,zi,2.0*(FractMaxOrder+i),ar_^[i],ai_^[i]);
      if odd(i) then
        begin
	  ar_^[i]:=-ar_^[i];
	  ai_^[i]:=-ai_^[i];
        end;
    end;
  br_:=ar_^[BfhMaxTerms];
  bi_:=ai_^[BfhMaxTerms];
  for i:=BfhMaxTerms downto 2 do
    begin
      SimplifyComplexFraction(br_,bi_,1.0,tempr,tempi);
      br_:=ar_^[i-1]+tempr;
      bi_:=ai_^[i-1]+tempi;
    end;
  SimplifyComplexFraction(br_,bi_,1.0,tempr,tempi);
  {Jx = ratio of BesselFunctions}
  SimplifyComplexFraction(zr,zi,2.0*FractMaxOrder,Jr^[BfhMaxOrder],Ji^[BfhMaxOrder]);
  Jr^[BfhMaxOrder]:=Jr^[BfhMaxOrder]+tempr;
  Ji^[BfhMaxOrder]:=Ji^[BfhMaxOrder]+tempi;
  for i:=BfhMaxOrder downto 1 do
    begin
      SimplifyComplexFraction(Jr^[i],Ji^[i],1.0,tempr,tempi);
      SimplifyComplexFraction(zr,zi,2.0*(i - 0.5),Jr^[i-1],Ji^[i-1]);
      Jr^[i-1]:=Jr^[i-1]-tempr;
      Ji^[i-1]:=Ji^[i-1]-tempi;
    end;
  for i:=0 to BfhMaxOrder do
    begin
      SimplifyComplexFraction(zr,zi,i,tempr,tempi);
      n1B_r^[i] := Jr^[i] - tempr;
      n1B_i^[i] := Ji^[i] - tempi;
    end;
  dispose(ar_);dispose(ai_);dispose(Jr);dispose(Ji);
end;


{-----------------------------------------------------------------------------
  Procedure: calc_anbn
  Purpose:   Calculates the Mie parameters for a given sizeparameter and
             refractive index
  Comments:  //according to Kerker (1969) equations 3.6.25 and 3.6.26
-----------------------------------------------------------------------------}

procedure calc_anbn (c:integer; mr,mi:extended;
  		var anr,ani,bnr,bni:extended);
var
  denomA,denomB,aX,aY,aU,aV,bX,bY,bU,bV,ar,ai,br,bi:extended;
begin
        //real part of the denominator of the Mie equation
  aU:=psiA^[c]*(n1B_r^[c]-mr*n3A_r^[c]+mi*n3A_i^[c])
      +chiA^[c]*(-n1B_i^[c]+n3A_r^[c]*mi+mr*n3A_i^[c]);
        //imaginary part of the denominator
  aV:=-psiA^[c]*(-n1B_i^[c]+n3A_r^[c]*mi+mr*n3A_i^[c])
      +chiA^[c]*(n1B_r^[c]-mr*n3A_r^[c]+mi*n3A_i^[c]);
        //real part of the numerator of 3.6.25
  aX:=psiA^[c]*(n1B_r^[c]-mr*n1A_r^[c]);
        //imaginary part of the numerator
  aY:=psiA^[c]*(n1B_i^[c]-mi*n1A_r^[c]);
  denomA:=sqr(aU)+sqr(aV);
  ar:=(aX*aU+aY*aV)/denomA;
  ai:=(aY*aU-aX*aV)/denomA;

        //same for equation 3.6.26
  bU:=psiA^[c]*(n3A_r^[c]-mr*n1B_r^[c]+mi*n1B_i^[c])
      +chiA^[c]*(-n3A_i^[c]+n1B_r^[c]*mi+mr*n1B_i^[c]);
  bV:=-psiA^[c]*(-n3A_i^[c]+n1B_r^[c]*mi+mr*n1B_i^[c])
      +chiA^[c]*(n3A_r^[c]-mr*n1B_r^[c]+mi*n1B_i^[c]);
  bX:=psiA^[c]*(n1A_r^[c]-mr*n1B_r^[c]+mi*n1B_i^[c]);
  bY:=-psiA^[c]*(mi*n1B_r^[c]+mr*n1B_i^[c]);
  denomB:=sqr(bU)+sqr(bV);
  br:=(bX*bU+bY*bV)/denomB;
  bi:=(bY*bU-bX*bV)/denomB;

  anr:=ar;
  ani:=ai;
  bnr:=br;
  bni:=bi;
end;


{-----------------------------------------------------------------------------
  Procedure: Miecalc
  Purpose:   Calculates all necessary Mie parameters for each sphere of the
             agglomerate using the Wiscombe criterion and stores the values in a
             three-dimensional array
  Comments:
-----------------------------------------------------------------------------}

procedure Miecalc;
var
  betar,betai : extended;
  i,j,k:integer;
begin
  SetLength(MieParameter,Nl,4,50);
  new(JaA);new(JfA);new(PsiA);
  new(chiA);
        {new(n1A_r)};//new(n1B_r);new(n1B_i);
  new(n3A_r);new(n3A_i);
  for i:=0 to Nl-1 do
    begin
      new(n1A_r);new(n1B_r);new(n1B_i);
        //new(n3A_r);new(n3A_i);
        {order of Bessel Functions must be higher than alpha}
      BfhMaxOrder:=trunc(WiscombeNumber[i])+100;
      BfhMinOrder:=trunc(WiscombeNumber[i])+50;
      BesselFractional (x[i],BfhMaxOrder,BfhMaxTerms,JaA,JfA);
      RicattiBesselCalculation(x[i]);
      betar:= RelRefrIndex[i,0]*x[i];
      betai:= RelRefrIndex[i,1]*x[i];
      calcn1A(x[i]);
      calcn1B(betar,betai);
      calcn3A(x[i]);
      for j:=0 to 3 do
        begin
          for k:=0 to 49 do
            MieParameter[i,j,k]:=0;
        end;
      for j:=1 to WiscombeNumber[i] do
        begin
          calc_anbn (j, RelRefrIndex[i,0], RelRefrIndex[i,1],MieParameter[i,0,j],MieParameter[i,1,j],
          MieParameter[i,2,j],MieParameter[i,3,j]);
        end;
      dispose(n1A_r);//dispose(n1B_r);dispose(n1B_i);
        //dispose(n3A_r);dispose(n3A_i);
    end;
  dispose(JaA);dispose(JfA);
  dispose(chiA);
  dispose(n3A_r);dispose(n3A_i);
  {dispose(n1A_r);}dispose(n1B_r);dispose(n1B_i);
end;


{-----------------------------------------------------------------------------
  Procedure: SpherBes
  Purpose:   Calculates the spherical Bessel functions of the first and
             second kind
  Comments:
-----------------------------------------------------------------------------}

procedure SpherBes (n:integer; z:float; var FK,SK:float);
begin
  new(JaA);new(Jfa);
  BesselFractional(z,n+1,BfhMaxTerms,JaA,Jfa);
  FK:=Sqrt(PiOn2/z)*JaA^[n];
  SK:=Sqrt(PiOn2/z)*Jfa^[n];
        //change of sign according to alternating series for JfA
  if odd(n) then SK:=-SK;
  dispose(JaA);dispose(Jfa);
end;


{-----------------------------------------------------------------------------
  Procedure: SpherHankFK
  Purpose:   Calculates the spherical Hankel function of the first kind
  Comments:
-----------------------------------------------------------------------------}

function SpherHankFK(n:integer;z:float):TComplexCart;
var
  Real,Imaginary:float;
begin
  SpherBes(n,z,Real,Imaginary);
  Result.Re:=Real;
  Result.Im:=Imaginary;
end;

end.
