{-----------------------------------------------------------------------------
 Unit Name: Gaunt
 Author:    kaetzel
 Purpose:   Calculation of the Gaunt coefficients for the calculation of
            multisphere scattering problems
-----------------------------------------------------------------------------}

unit Gaunt;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  Math, StdCtrls, Menus, ExtCtrls, ComCtrls, JvComCtrls, JvxSlider,JvDSADialogs,
  JclMath, JvWaitingProgress;

        {                   Public Functions                            }

function AssociatedLegendre(n,m:integer;x:extended):extended;
function Factorials (z:integer):extended;

type
  TForm1 = class(TForm)
    StartButton:TButton;
    CloseButton: TButton;
    Headline:TLabel;
    EnterSingle: TLabel;
    Panel1:TPanel;
    Panel2:TPanel;
    Panel3: TPanel;
    Value_n:TEdit;
    Value_m:TEdit;
    Value_nu:TEdit;
    Value_mu:TEdit;
    Value_p: TEdit;
    Label1:Tlabel;
    Label2:Tlabel;
    Label3:Tlabel;
    Label4:Tlabel;
    Label5:Tlabel;
    Label6:Tlabel;
    Label7:Tlabel;
    Label8:Tlabel;
    Label9:Tlabel;
    Label10:Tlabel;
    GauntCoefficient: TEdit;
    Label11: TLabel;
    Label12: TLabel;
    Label13: TLabel;
    Label14: TLabel;
    Label15: TLabel;
    Label16: TLabel;
    Label17: TLabel;
    Label18: TLabel;
    Label19: TLabel;
    Label20: TLabel;
    Label21: TLabel;
    Label23: TLabel;
    Label24: TLabel;
    WiscombeNumber:TEdit;
    Filename: TEdit;
    Label25: TLabel;
    Calculating: TLabel;
    Finished: TLabel;
    JvxSlider1: TJvxSlider;
    ResetButton: TButton;
    Edit1: TEdit;
    Label26: TLabel;
    Label27: TLabel;
    procedure StartButtonClick(Sender: TObject);
    procedure CloseButtonClick(Sender: TObject);
    procedure ResetButtonClick(Sender: TObject);
    procedure FormActivate(Sender: TObject);
    procedure JvxSlider1Changed(Sender: TObject);

  private
    { Déclarations privées }
  public
    { Déclarations publiques }
  end;

var
  Form1: TForm1;
  ButtonClick,ExitTrue:bool;

implementation

uses Gaunt_recursive;

{$R *.DFM}


{-----------------------------------------------------------------------------
  Procedure: Factorials
  Purpose:   Calculation of the Factorial of an integer
  Comments:  according to the model of Lanczos
             Lanczos, C.J. ; SIAM. Numer. Anal. Ser. B1, 86 (1964)
-----------------------------------------------------------------------------}

function Factorials (z:integer):extended;
const
  c0:real = 0.1E+01;
  c1:double = 0.16427423239836267E+05;
  c2:double = -0.48589401600331902E+05;
  c3:double = 0.55557391003815523E+05;
  c4:double = -0.30964901015912058E+05;
  c5:double = 0.87287202992571788E+04;
  c6:double = -0.11714474574532352E+04;
  c7:double = 0.63103078123601037E+02;
  c8:double = -0.93060589791758878E+00;
  c9:double = 0.13919002438227877E-02;
  c10:double = -0.45006835613027859E-08;
  c11:double = 0.13069587914063262E-09;
begin
  Result:=(Exp(Ln(2*Pi)/2-(z+10+0.5)+(z+0.5)*Ln(z+10+0.5)+Ln(c0+c1/(z+1)+
  c2/(z+2)+c3/(z+3)+c4/(z+4)+c5/(z+5)+c6/(z+6)+c7/(z+7)+c8/(z+8)+c9/(z+9)+
  c10/(z+10)+c11/(z+11))));
end;


{-----------------------------------------------------------------------------
  Procedure: AssociatedLegendre
  Purpose:   Calculation of the associated Legendre Polynomial of order n and
             degree m
  Comments:  According to "Numerical Recipes in Pascal"
             Stable backward recursion
-----------------------------------------------------------------------------}

function AssociatedLegendre (n,m: integer; x: extended): extended;
var
   IsNegative:bool;
   fact,pll,pmm,pmmp1,somx2: real;
   i,ll: integer;
begin
  IsNegative:=False;
  if (m < 0) then
    begin
      IsNegative:=True;
      m:=-m;
    end;
  pmm:=1.0;
  if (m>0) then
    begin
      somx2:=sqrt((1.0-x)*(1.0+x));
      fact:=1.0;
      for i:=1 to m do
        begin
          pmm:=-pmm*fact*somx2;
          fact:=fact+2.0
        end;
    end;
  if (n=m) then
    begin
     Result:=pmm;
    end
  else
    begin
      pmmp1:=x*(2*m+1)*pmm;
      if (n=m+1) then
        begin
          Result:=pmmp1;
        end
      else
        begin
          for ll:=m+2 to n do
            begin
              pll:=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
              pmm:=pmmp1;
              pmmp1:=pll;
            end;
          Result:=pll;
        end;
    end;
  Result:=Result*Power(-1,m);
  if IsNegative then
    Result:=Result*Power(-1,m)*Factorials(n-m)/Factorials(n+m);
end;


{-----------------------------------------------------------------------------
  Procedure: GauntCo
  Purpose:   Direct calculation of the Gaunt coefficient a(m,n,mu,nu,p)
  Comments:  Integral calculation of three associated Legendre Polynomials
             with Simpson's rule
             Because of Singularities of the Associated Legendre Polynomial at
             the integration boundaries, this calculation is not always exact.
             That is why the purpose of this procedure is only to check the
             recursive calculation.
-----------------------------------------------------------------------------}

function GauntCo(m,n,mu,nu,p:integer):extended;
const
  upper=1;
  lower=-1;
  delta:real=0.001;
var
  i:integer;
  x:real;
  f,flower,odd_sum,even_sum,integral:extended;
begin
  if p<abs(m+mu) then
    Result:=0
  else
    begin
      flower:=AssociatedLegendre(n,m,-1)*AssociatedLegendre(nu,mu,-1)*
      AssociatedLegendre(p,(m+mu),-1);
      odd_sum:=0;
      even_sum:=0;
      for i:=1 to (round((upper-lower)/delta)) do
      begin
        x:=lower+i*delta;
        f:=AssociatedLegendre(n,m,x)*AssociatedLegendre(nu,mu,x)*
        AssociatedLegendre(p,(m+mu),x);
        if odd(i) then
          odd_sum:=odd_sum+f
        else
          even_sum:=even_sum+f;
      end;
      integral:=delta/3*(flower+4*odd_sum+2*even_sum);
      Result:=(2*p+1)/2*Factorials(p-m-mu)/(Factorials(p+m+mu))*integral;
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: Delay
  Purpose:   Stops the program for a defined time
  Comments:
-----------------------------------------------------------------------------}

procedure Delay(msec: LongInt);
var start, stop: LongInt;
begin
  start := GetTickCount;
  repeat
    stop := GetTickCount;
    Application.ProcessMessages;
  until (stop-start)>=msec;
end;


{-----------------------------------------------------------------------------
  Procedure: TForm1.StartButtonClick
  Purpose:   Performs the calculation of the Gaunt Coefficient specified or all
             possible values, depends on the position of the slider
  Comments:
-----------------------------------------------------------------------------}

procedure TForm1.StartButtonClick(Sender: TObject);
var
  n,m,mu,nu,p,Nl:integer;
  Filenamed:string[22];
begin
  ButtonClick:=False;
  ExitTrue:=False;
  Nl:=0;
  Finished.Hide;
  Label27.Hide;
        //calculation of a single Gaunt coefficient
  if JvxSlider1.Value=0 then
        //reading the indizes of the coefficient
    begin
      try
        n:=StrToInt(Value_n.Text);
      except
        ShowMessage('n must be an integer value!',dckActiveForm);
        Value_n.SetFocus;
        exit;
      end;
      if n<1 then
        begin
          ShowMessage('n must be greater than Zero!',dckActiveForm);
          Value_n.SetFocus;
          exit;
        end
      else
        begin
          try
            m:=StrToInt(Value_m.Text);
          except
            ShowMessage('m must be an integer value!',dckActiveForm);
            Value_m.SetFocus;
            exit;
          end;
          if abs(m)>n then
            begin
              ShowMessage('Modulus of m must be smaller than or equal to n!'
                          ,dckActiveForm);
              Value_m.SetFocus;
              exit;
            end
          else
            begin
              try
                nu:=StrToInt(Value_nu.Text);
              except
                ShowMessage('nu must be an integer value!',dckActiveForm);
                Value_nu.SetFocus;
                exit;
              end;
              if nu<1 then
                begin
                  ShowMessage('nu must be greater than Zero!',dckActiveForm);
                  Value_nu.SetFocus;
                  exit;
                end
              else
                begin
                  try
                    mu:=StrToInt(Value_mu.Text);
                  except
                    ShowMessage('mu must be an integer value!',dckActiveForm);
                    Value_mu.SetFocus;
                    exit;
                  end;
                  if abs(mu)>nu then
                    begin
                      ShowMessage('Modulus of mu must be smaller than or equal'+
                                  ' to nu!',dckActiveForm);
                      Value_mu.SetFocus;
                      exit;
                    end
                  else
                    begin
                      try
                        p:=StrToInt(Value_p.Text);
                      except
                        ShowMessage('p must be an integer value!',dckActiveForm);
                        Value_p.SetFocus;
                        exit;
                      end;
                      if ((p<0) or (p>(n+nu))) then
                        begin
                          ShowMessage('p must be greater than or equal to Zero'+
                                      ' and smaller'+
                                      ' than n+nu!',dckActiveForm);
                          Value_p.SetFocus;
                          exit;
                        end;
                    end;
                end;
            end;
        end;
        Calculating.Show;
        Delay(1);
        GauntCoefficient.Text:=FloatToStrF(GauntRecursiveSingleValue(m,n,mu,nu,p)
                                           ,ffExponent,11,3);
        Calculating.Hide;
        Finished.Show;
    end;
        //This is the option of calculating and saving all possible Gaunt
        //coefficients for a known number of Wiscombe summations
    if JvxSlider1.Value=1 then
      begin
        try
          Nl:=StrToInt(WiscombeNumber.Text);
          n:=StrToInt(Edit1.Text);
          if ((n<1) or (n>N)) then abort;
        except
          ShowMessage('Nl and n must be an integer values! 0<n<=N!',dckActiveForm);
          WiscombeNumber.SetFocus;
        end;
        if ((Nl<0)or (Nl>20)) then
          begin
            ShowMessage(Format(' Nl should be in the range of'+
            ' 0<=Nl<=20'+'%s Stability of backward recursion can not be guaranteed! %s',
            [#10#13#10#13,#10#13] ),dckActiveForm);
            WiscombeNumber.SetFocus;
          end
        else
          begin
            Filenamed:='GauntCoeff_'+IntToStr(Nl)+'_'+IntToStr(n)+'.dat';
            Filename.Text:=Filenamed;
            Calculating.Show;
            Delay(1);
            AllValuesRecursive(Nl,n,Filenamed);
            if ExitTrue then
              begin
                Calculating.Hide;
                Label27.Show;
                exit;
              end;
            Calculating.Hide;
            Finished.Show;
          end;
      end;
end;


{-----------------------------------------------------------------------------
  Procedure: TForm1.CloseButtonClick
  Purpose:   Closes the application
  Comments:
-----------------------------------------------------------------------------}

procedure TForm1.CloseButtonClick(Sender: TObject);
begin
  ButtonClick:=True;
  Form1.Close;
end;


{-----------------------------------------------------------------------------
  Procedure: TForm1.ResetButtonClick
  Purpose:   Reinitializes all values in the Edit Components and stops
             calculation
  Comments:
-----------------------------------------------------------------------------}

procedure TForm1.ResetButtonClick(Sender: TObject);
begin
  ButtonClick:=True;
  Calculating.Visible:=False;
  Finished.Visible:=False;
  Value_p.Text:='';
  Value_m.Text:='';
  Value_n.Text:='';
  Value_mu.Text:='';
  Value_nu.Text:='';
  GauntCoefficient.Text:='';
  Filename.Text:='';
  WiscombeNumber.Text:='';
  Edit1.Text:='';
  if JvxSlider1.Value=0 then
    Value_n.SetFocus
  else
    WiscombeNumber.SetFocus;
end;


{-----------------------------------------------------------------------------
  Procedure: TForm1.FormActivate
  Purpose:   Standard Values
  Comments:
-----------------------------------------------------------------------------}

procedure TForm1.FormActivate(Sender: TObject);
begin
  JvxSlider1.Value:=0;
  WiscombeNumber.Enabled:=False;
  Edit1.Enabled:=False;
  Value_n.SetFocus;
end;


{-----------------------------------------------------------------------------
  Procedure: TForm1.JvxSlider1Changed
  Purpose:   Change of standard values if Slider position is changed
  Comments:
-----------------------------------------------------------------------------}

procedure TForm1.JvxSlider1Changed(Sender: TObject);
begin
  if JvxSlider1.Value=1 then
    begin
      WiscombeNumber.Enabled:=True;
      Edit1.Enabled:=True;
      WiscombeNumber.SetFocus;
      Value_n.Text:='';
      Value_n.Enabled:=False;
      Value_m.Text:='';
      Value_m.Enabled:=False;
      Value_nu.Text:='';
      Value_nu.Enabled:=False;
      Value_mu.Text:='';
      Value_mu.Enabled:=False;
      Value_p.Text:='';
      Value_p.Enabled:=False;
      GauntCoefficient.Text:='';
      Finished.Hide;
      Calculating.Hide;
    end
  else
    begin
      WiscombeNumber.Enabled:=False;
      Edit1.Enabled:=False;
      WiscombeNumber.Text:='';
      Filename.Text:='';
      Value_n.Enabled:=True;
      Value_n.SetFocus;
      Value_m.Enabled:=True;
      Value_nu.Enabled:=True;
      Value_mu.Enabled:=True;
      Value_p.Enabled:=True;
      Finished.Hide;
      Calculating.Hide;
    end;
end;

end.
