{-----------------------------------------------------------------------------
 Unit Name: Gaunt_recursive
 Author:    kaetzel
 Purpose:   Provides the recursive calculation of Gaunt coefficients because
            direct calculation is not always exact and very time consuming
            (days-->weeks)
-----------------------------------------------------------------------------}

unit Gaunt_recursive;

interface

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

        {                       Public Functions                        }

function GauntRecursiveSingleValue(m,n,mu,nu,psoll:integer):extended;
procedure AllValuesRecursive(Nl,n:ShortInt;Filenamed:String);


implementation

type
 EFileExistsError=class(Exception);


{-----------------------------------------------------------------------------
  Procedure: Alpha_p
  Purpose:   Aid Function for the recursive calculation of Gaunt Coefficients
  Comments:
-----------------------------------------------------------------------------}

function Alpha_p(n,nu,p:integer):real;
begin
  Result:=(p*p-Sqr(n-nu))*(p*p-Sqr(n+nu+1))/(4*p*p-1);
end;


{-----------------------------------------------------------------------------
  Procedure: Ap
  Purpose:   Aid function for the recursive calculation of Gaunt coefficients
  Comments:
-----------------------------------------------------------------------------}

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


{-----------------------------------------------------------------------------
  Procedure: GammaFunction
  Purpose:   Calculates the integral t^(x-1) exp(-t) dt from zero to infinity
  Comments:  Based on Matlab-Code, Version 6.5
             This function was used to calculate the Pochhammer-Symbol for
             non-integer values in another kind of recursion formulae
-----------------------------------------------------------------------------}

function GammaFunction(x:extended):extended;
const
  p:array[1..8] of extended = (-1.71618513886549492533811e+0,
                                2.47656508055759199108314e+1,
                               -3.79804256470945635097577e+2,
                                6.29331155312818442661052e+2,
                                8.66966202790413211295064e+2,
                               -3.14512729688483675254357e+4,
                               -3.61444134186911729807069e+4,
                                6.64561438202405440627855e+4);
  q:array[1..8] of extended = (-3.08402300119738975254353e+1,
                                3.15350626979604161529144e+2,
                               -1.01515636749021914166146e+3,
                               -3.10777167157231109440444e+3,
                                2.25381184209801510330112e+4,
                                4.75584627752788110767815e+3,
                               -1.34659959864969306392456e+5,
                               -1.15132259675553483497211e+5);
  c:array[1..7] of extended = (-1.910444077728e-3,
                                8.4171387781295e-4,
                               -5.952379913043012e-4,
                                7.93650793500350248e-4,
                               -2.777777777777681622553e-3,
                                8.333333333333333331554247e-2,
                                5.7083835261e-3);
  spi=0.9189385332046727417803297;
var
  y1,i,j,xn:integer;
  IsNegative,IsSmallerOne,IsSmallerTwelve:bool;
  y,x1:extended;
  res,factor,z,xnum,xden,sum:extended;
begin
   IsNegative:=False;
   IsSmallerOne:=False;
   IsSmallerTwelve:=False;
   xn:=0;
   x1:=0;
   factor:=0;
   res:=0;
   Result:=0;
   if x<=0 then
     begin
       IsNegative:=True;
       y:=-x;
       y1:=trunc(y);
       res:=y-y1;
       factor:=-Pi/(sin(Pi*res))*(1-2*(y1 mod 2));
       x:=y+1;
     end;
   if x<1 then
     begin
       IsSmallerOne:=True;
       x1:=x;
       x:=x1+1;
     end;
   if x<12 then
     begin
       IsSmallerTwelve:=True;
       xn:=trunc(x)-1;
       x:=x-xn;
     end;
   if IsSmallerTwelve then
     begin
       z:=x-1;
       xnum:=0*z;
       xden:=xnum+1;
       for i:=1 to 8 do
         begin
           xnum:=(xnum+p[i])*z;
           xden:=(xden*z+q[i]);
         end;
       res:=xnum/xden+1;
     end;
   if IsSmallerOne then
     res:=res/x1;
   for j:=1 to xn do
     begin
       res:=res*x;
       x:=x+1;
     end;
   if x>=12 then
     begin
       y:=x;
       sum:=c[7];
       for i:=1 to 6 do
         sum:=sum/(y*y)+c[i];
       sum:=sum/y-y+spi;
       sum:=sum+(y-0.5)*ln(y);
       res:=exp(sum);
     end;
   if isnegative then
     res:=factor/res;
   Result:=res;
end;

{-----------------------------------------------------------------------------
  Procedure: Pochhammer
  Purpose:   Calculates Pochhammers Symbol for integer values
  Comments:  Refers to Abramowitz/Stegun Handbook of Mathematical Functions
             Original Definition is: (z)n=Gamma(z+n)/Gamma(z)
-----------------------------------------------------------------------------}

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


{-----------------------------------------------------------------------------
  Procedure: GauntRecursiveSingleValue
  Purpose:   Calculates a specified Gaunt Coefficient by backward recursion
  Comments:  The stability of this algorithm has only been tested for all
             possible values for n,nu<=20, actually most values are still
             correct until n=45 but especially values with Ap(p+4)=0 are
             unstable.
-----------------------------------------------------------------------------}

function GauntRecursiveSingleValue(m,n,mu,nu,psoll:integer):extended;
var
  qsoll,q,qmax,p,p1,p2,pmin,n4:integer;
  c0,c1,c2,c3:extended;
  a:array of extended;
begin
  qsoll:=(n+nu-psoll)div 2;
        //All values for q which are not integers do not exist
  if ((n+nu-psoll) mod 2)<>0 then
    begin
      Result:=0;
      exit;
    end;
  qmax:=min(min(n,nu),((n+nu-(abs(m+mu)))div 2));
        //All values bigger than qmax are exactly zero
  if qsoll>qmax then
    begin
      Result:=0;
      exit;
    end;
  n4:=n+nu-m-mu;
  SetLength(a,qmax+1);

        //Backward Recursion

        //Starting value
  a[0]:=Pochhammer(n+1,n)*Pochhammer(nu+1,nu)*Factorials(n+nu-m-mu)/
        (Pochhammer(n+nu+1,n+nu)*Factorials(n-m)*Factorials(nu-mu));

  if qsoll>=1 then
        //Second Starting value
    a[1]:=a[0]*(2*n+2*nu-3)/2*(1-(2*n+2*nu-1)/(n4*(n4-1))*((m-n)*(m-n+1)/
             (2*n-1)+(mu-nu)*(mu-nu+1)/(2*nu-1)));

        //Recursion
  for q:=2 to qsoll do
    begin
      p:=n+nu-2*q;
      p1:=p-m-mu;
      p2:=p+m+mu;
        //Two term recurrence relation
      if (((mu=m) and (nu=n)) or ((m=0) and (mu=0))) then
        a[q]:=(p+1)*(p2+2)*Alpha_p(n,nu,p+2)/((p+2)*(p1+1)*Alpha_p(n,nu,p+1))*
              a[q-1]
      else
        begin
                //Three-term recurrence relation
          if Ap(m,n,mu,nu,p+4)<>0 then
            begin
              c0:=(p+2)*(p+3)*(p1+1)*(p1+2)*Ap(m,n,mu,nu,p+4)*Alpha_p(n,nu,p+1);
              c1:=Ap(m,n,mu,nu,p+2)*Ap(m,n,mu,nu,p+3)*Ap(m,n,mu,nu,p+4)+
                  (p+1)*(p+3)*(p1+2)*(p2+2)*Ap(m,n,mu,nu,p+4)*Alpha_p(n,nu,p+2)+
                  (p+2)*(p+4)*(p1+3)*(p2+3)*Ap(m,n,mu,nu,p+2)*Alpha_p(n,nu,p+3);
              c2:=-(p+2)*(p+3)*(p2+3)*(p2+4)*Ap(m,n,mu,nu,p+2)*Alpha_p(n,nu,p+4);
              a[q]:=c1/c0*a[q-1]+c2/c0*a[q-2];
            end
          else
            begin
                //Four term recurrence relation
              if q=2 then
                        //Third starting value
                  a[q]:=a[0]*(2*n+2*nu-1)*(2*n+2*nu-7)/4 * ( (2*n+2*nu-3)/(n4*(n4-1))
                        * ( (2*n+2*nu-5)/(2*(n4-2)*(n4-3)) * ( (m-n)*(m-n+1)*
                        (m-n+2)*(m-n+3)/((2*n-1)*(2*n-3)) + 2*(m-n)*(m-n+1)*(mu-nu)*
                        (mu-nu+1)/((2*n-1)*(2*nu-1)) + (mu-nu)*(mu-nu+1)*(mu-nu+2)*
                        (mu-nu+3)/((2*nu-1)*(2*nu-3)) ) - (m-n)*(m-n+1)/(2*n-1) -
                        (mu-nu)*(mu-nu+1)/(2*nu-1) ) + 0.5 )
              else
                begin
                  c0:=(p+2)*(p+3)*(p+5)*(p1+1)*(p1+2)*(p1+4)*Ap(m,n,mu,nu,p+6)*
                      Alpha_p(n,nu,p+1);
                  c1:=(p+5)*(p1+4)*Ap(m,n,mu,nu,p+6)*(Ap(m,n,mu,nu,p+2)*
                      Ap(m,n,mu,nu,p+3)+(p+1)*(p+3)*(p1+2)*(p2+2)*Alpha_p(n,nu,p+2));
                  c2:=(p+2)*(p2+3)*Ap(m,n,mu,nu,p+2)*(Ap(m,n,mu,nu,p+5)*
                      Ap(m,n,mu,nu,p+6)+(p+4)*(p+6)*(p1+5)*(p2+5)*Alpha_p(n,nu,p+5));
                  c3:=-1*(p+2)*(p+4)*(p+5)*(p2+3)*(p2+5)*(p2+6)*Ap(m,n,mu,nu,p+2)*
                      Alpha_p(n,nu,p+6);
                  a[q]:=c1/c0*a[q-1]+c2/c0*a[q-2]+c3/c0*a[q-3];
                end;
            end;
        end;
    end;
  Result:=a[qsoll];
end;


{-----------------------------------------------------------------------------
  Procedure: AllValuesRecursive
  Purpose:   Calculates all possible Gaunt Coefficients with known number of
             Wiscombe summations
  Comments:  Backward recursion stability only tested for values Nl<20, but
             seems to work also until 25, 29 is known to be unstable
-----------------------------------------------------------------------------}

procedure AllValuesRecursive(Nl,n:ShortInt;Filenamed:String);
type DataSet=packed record
  GC1:extended;
end;
var
  Gaunts2:file of DataSet;
  data:DataSet;
  m2,n2,mu2,nu2,p,p1,p2,q,qmax,qmax_check,pmin,n4:SmallInt;
  GC:extended;
  c0,c1,c2,c3,actrl,actrl2,realeps:extended;
  a:array of extended;
begin
  try
    AssignFile(Gaunts2,Filenamed);
    if FileExists(Filenamed) then
      Raise EFileExistsError.CreateFmt('File already exists! %s'+
                                       'We will not do the same work again! %s',
                                       [#10#13#10#13,#10#13] )
    else
      begin
        ReWrite(Gaunts2);
        n2:=n;
          for m2:=-n2 to n2 do
            begin
              for nu2:=1 to Nl do
                begin
                  for mu2:=-nu2 to nu2 do
                    begin
                      Application.ProcessMessages;
                      if ButtonClick then
                        begin
                          CloseFile(Gaunts2);
                          DeleteFile(Filenamed);
                          ExitTrue:=True;
                          exit;
                        end;
                        //Two different calculations for qmax have been
                        //implemented to assure the same file structure as
                        //with the integration algorithm
                      qmax:=Min(n2,Min(nu2,Round((n2+nu2-abs(m2+mu2))/2)));
                      SetLength(a,qmax+1);
                      n4:=n2+nu2-mu2-m2;
                      qmax_check:=Min(n2,Min(nu2,(n2+nu2-abs(m2+mu2)) div 2));
                      if qmax_check<qmax then
                        begin
                          a[qmax]:=0;
                          qmax:=qmax_check;
                        end;
                      pmin:=n2+nu2-2*qmax;

                        //Backward Recursion

                        //Starting value
                      a[0]:=Pochhammer(n2+1,n2)*Pochhammer(nu2+1,nu2)*
                            Factorials(n2+nu2-m2-mu2)/(Pochhammer(n2+nu2+1,n2+nu2)*
                            Factorials(n2-m2)*Factorials(nu2-mu2));
                      if qmax>=1 then
                                //Second starting value
                        a[1]:=a[0]*(2*n2+2*nu2-3)/2*(1-(2*n2+2*nu2-1)/(n4*(n4-1))*
                                   ((m2-n)*(m2-n2+1)/
                                   (2*n-1)+(mu2-nu2)*(mu2-nu2+1)/(2*nu2-1)));
                        //Recursion
                      for q:=2 to qmax do
                        begin
                          p:=n2+nu2-2*q;
                          p1:=p-m2-mu2;
                          p2:=p+m2+mu2;
                          if (((mu2=m2) and (nu2=n2)) or ((m2=0) and (mu2=0))) then
                            a[q]:=(p+1)*(p2+2)*Alpha_p(n2,nu2,p+2)/
                                  ((p+2)*(p1+1)*Alpha_p(n2,nu2,p+1))*a[q-1]
                          else
                            begin
                              if Ap(m2,n2,mu2,nu2,p+4)<>0 then
                                begin
                                  c0:=(p+2)*(p+3)*(p1+1)*(p1+2)*
                                      Ap(m2,n2,mu2,nu2,p+4)*Alpha_p(n2,nu2,p+1);
                                  c1:=Ap(m2,n2,mu2,nu2,p+2)*Ap(m2,n2,mu2,nu2,p+3)*
                                      Ap(m2,n2,mu2,nu2,p+4)+(p+1)*(p+3)*(p1+2)*
                                      (p2+2)*Ap(m2,n2,mu2,nu2,p+4)*
                                      Alpha_p(n2,nu2,p+2)+(p+2)*(p+4)*(p1+3)*
                                      (p2+3)*Ap(m2,n2,mu2,nu2,p+2)*
                                      Alpha_p(n2,nu2,p+3);
                                  c2:=-(p+2)*(p+3)*(p2+3)*(p2+4)*
                                      Ap(m2,n2,mu2,nu2,p+2)*Alpha_p(n2,nu2,p+4);
                                  a[q]:=c1/c0*a[q-1]+c2/c0*a[q-2];
                                end
                              else
                                begin
                                  if q=2 then
                                    a[q]:=a[0]*(2*n2+2*nu2-1)*(2*n2+2*nu2-7)/4 *
                                          ( (2*n2+2*nu2-3)/(n4*(n4-1))
                                          * ( (2*n2+2*nu2-5)/(2*(n4-2)*(n4-3)) *
                                          ( (m2-n2)*(m2-n2+1)*
                                          (m2-n2+2)*(m2-n2+3)/((2*n2-1)*(2*n2-3)) +
                                          2*(m2-n2)*(m2-n2+1)*(mu2-nu2)*
                                          (mu2-nu2+1)/((2*n2-1)*(2*nu2-1)) +
                                          (mu2-nu2)*(mu2-nu2+1)*(mu2-nu2+2)*
                                          (mu2-nu2+3)/((2*nu2-1)*(2*nu2-3)) ) -
                                          (m2-n2)*(m2-n2+1)/(2*n2-1) -
                                          (mu2-nu2)*(mu2-nu2+1)/(2*nu2-1) ) + 0.5 )
                                  else
                                    begin
                                      c0:=(p+2)*(p+3)*(p+5)*(p1+1)*(p1+2)*(p1+4)*
                                          Ap(m2,n2,mu2,nu2,p+6)*Alpha_p(n2,nu2,p+1);
                                      c1:=(p+5)*(p1+4)*Ap(m2,n2,mu2,nu2,p+6)*
                                          (Ap(m2,n2,mu2,nu2,p+2)*
                                          Ap(m2,n2,mu2,nu2,p+3)+
                                          (p+1)*(p+3)*(p1+2)*(p2+2)*
                                          Alpha_p(n2,nu2,p+2));
                                      c2:=(p+2)*(p2+3)*Ap(m2,n2,mu2,nu2,p+2)*
                                          (Ap(m2,n2,mu2,nu2,p+5)*
                                          Ap(m2,n2,mu2,nu2,p+6)+
                                          (p+4)*(p+6)*(p1+5)*(p2+5)*
                                          Alpha_p(n2,nu2,p+5));
                                      c3:=-1*(p+2)*(p+4)*(p+5)*(p2+3)*(p2+5)*(p2+6)*
                                          Ap(m2,n2,mu2,nu2,p+2)*Alpha_p(n2,nu2,p+6);
                                      a[q]:=c1/c0*a[q-1]+c2/c0*a[q-2]+c3/c0*a[q-3];
                                    end;
                                end;
                            end;
                        end;


                        //Last value of the recursion can be calculated
                        //analytically, this is used to check the stability
                        //of the backward recursion algorithm but can also be
                        //further developed to a combined forward and backward
                        //recursion scheme

                        if pmin=(n2-nu2) then
                          actrl:=Power(-1,mu2)*Pochhammer(nu2+1,nu2)*
                                 Factorials(n2+m2)*Factorials(2*pmin+1)/
                                 (Pochhammer(n2+1,n2+1)*Factorials(nu2-mu2)*
                                 Factorials(n2-nu2)*Factorials(pmin+m2+mu2))
                        else
                          begin
                            if pmin=(nu2-n2) then
                              actrl:=Power(-1,m2)*Pochhammer(n2+1,n2)*
                                     Factorials(nu2+mu2)*Factorials(2*pmin+1)/
                                     (Pochhammer(nu2+1,nu2+1)*Factorials(n2-m2)*
                                     Factorials(nu2-n2)*Factorials(pmin+m2+mu2))
                            else
                              begin
                                if pmin=m2+mu2 then
                                  actrl:=Power(-1,n2+m2-qmax)*(2*pmin+1)*
                                         Pochhammer(qmax+1,qmax)*
                                         Factorials(n2+nu2-qmax)*
                                         Factorials(n2+m2)*Factorials(nu2+mu2)/
                                         (Factorials(n2-qmax)*Factorials(nu2-qmax)*
                                         Factorials(n2-m2)*Factorials(nu2-mu2)*
                                         Factorials(n2+nu2+pmin+1))
                                else
                                  begin
                                    if pmin=-m2-mu2 then
                                      actrl:=Power(-1,nu2+mu2-qmax)*(2*pmin+1)*
                                             Pochhammer(qmax+1,qmax)*
                                             Factorials(n2+nu2-qmax)*
                                             Factorials(pmin-m2-mu2)/
                                             (Factorials(n2-qmax)*
                                             Factorials(nu2-qmax)*
                                             Factorials(n2+nu2+pmin+1))
                                    else
                                      begin
                                        if pmin=m2+mu2+1 then
                                          actrl:=Power(-1,n2+m2-qmax)*
                                                 Ap(m2,n2,mu2,nu2,pmin)*(2*pmin+1)*
                                                 Pochhammer(qmax+1,qmax)/
                                                 ((pmin-1)*Factorials(n2+nu2+pmin+1))*
                                                 Factorials(n2+nu2-qmax)*
                                                 Factorials(n2+m2)*
                                                 Factorials(nu2+mu2)/
                                                 (Factorials(n2-qmax)*
                                                 Factorials(nu2-qmax)*
                                                 Factorials(n2-m2)*
                                                 Factorials(nu2-mu2))
                                        else
                                          actrl:=Power(-1,nu2+mu2-qmax)*
                                                 Ap(m2,n2,mu2,nu2,pmin)*
                                                 (2*pmin+1)*Pochhammer(qmax+1,qmax)/
                                                 ((pmin-1)*Factorials(n2+nu2+pmin+1))*
                                                 Factorials(n2+nu2-qmax)*
                                                 Factorials(pmin-m2-mu2)/
                                                 (Factorials(n2-qmax)*
                                                 Factorials(nu2-qmax));
                                      end;
                                  end;
                              end;
                          end;
                        //Calculation of the deviation from the real value at
                        //the end of the recursion
                      if actrl=0 then
                        realeps:=a[qmax]
                      else
                        realeps:=abs((a[qmax]-actrl)/actrl);
                      a[qmax]:=actrl;
                      if realeps>1E-5 then abort;
                        //Writes the whole array to the file
                      for q:=0 to High(a) do
                        begin
                          data.GC1:=a[q];
                          Write(Gaunts2,(data));
                        end;
                    end;
                end;
            end;
        CloseFile(Gaunts2);
      end;
  except
    on EFileExistsError do
      ShowMessage(Format('File already exists! %s'+
                         'We will not do the same work again! %s',
                         [#10#13#10#13,#10#13]),dckActiveForm )
    else
      begin
        ShowMessage(Format('WARNING! %s'+'Backward Recursion for Gaunt Coefficients'+
                           ' failed (eps>10e-5)! %s'+'Process aborted! %s',
                           [#10#13#10#13,#10#13,#10#13] ),dckActiveForm);
        CloseFile(Gaunts2);
        DeleteFile(Filenamed);
        Form1.Calculating.Hide;
        Form1.Label27.Show;
      end;
  end;
end;

end.
