{*************************************************************
*       Calculate Function Gamma with a complex argument     *
* ---------------------------------------------------------- *
* EXPLANATION:                                               *
* Purpose: This program computes the gamma function G(z)     *
*          or Ln[G(z)] for a complex argument using          *
*          subroutine CGAMA                                  *
* Input :  x  --- Real part of z                             *
*          y  --- Imaginary part of z                        *
*          KF --- Function code                              *
*          KF=0 for Ln[G(z)]                                 *
*          KF=1 for G(z)                                     *
* Output:  GR --- Real part of Ln[G(z)] or G(z)              *
*          GI --- Imaginary part of Ln[G(z)] or G(z)         *
* Examples:                                                  *
*    x         y           Re[G(z)]           Im[G(z)]       *
*  --------------------------------------------------------  *
*   2.50      5.00     .2267360319E-01    -.1172284404E-01   *
*   5.00     10.00     .1327696517E-01     .3639011746E-02   *
*   2.50     -5.00     .2267360319E-01     .1172284404E-01   *
*   5.00    -10.00     .1327696517E-01    -.3639011746E-02   *
*                                                            *
*    x         y          Re[LnG(z)]         Im[LnG(z)]      *
*  --------------------------------------------------------  *
*   2.50      5.00    -.3668103262E+01     .5806009801E+01   *
*   5.00     10.00    -.4285507444E+01     .1911707090E+02   *
*   2.50     -5.00    -.3668103262E+01    -.5806009801E+01   *
*   5.00    -10.00    -.4285507444E+01    -.1911707090E+02   *
* ---------------------------------------------------------- *
* SAMPLE RUNS:                                               *
*                                                            *
* Please enter KF, x and y: 1 2.5 5.0                        *
*                                                            *
*    x         y           Re[G(z)]           Im[G(z)]       *
*  --------------------------------------------------------  *
*   2.50      5.00       0.0226736032      -0.0117228440     *
*                                                            *
* Please enter KF, zx and zy: 0 2.5 5.0                      *
*                                                            *
*    x         y          Re[LnG(z)]         Im[LnG(z)]      *
*  --------------------------------------------------------  *
*   2.50      5.00      -3.6681032624      5.8060098006      *
*                                                            *
* ---------------------------------------------------------- *
* REFERENCE:                                                 *
*    "Fortran Routines for Computation of Special Functions, *
*     jin.ece.uiuc.edu/routines/routines.html".              *
*                                                            *
*                       Pascal Release By J-P Moreau, Paris. *
*                                (www.jpmoreau.fr)           *
*************************************************************}
  PROGRAM MCGAMA;
    Uses WinCrt;

    Var X,X0,Y,GR,GI: Double;
        KF: Integer;

    {auxiliary functions}
    Function SINH(x:double):double;
    Var expx: double;
    Begin
      expx := exp(x);
      SINH := 0.5*(expx-1.0/expx)
    End;

    Function COSH(x:double):double;
    Var expx: double;
    Begin
      expx := exp(x);
      COSH := 0.5*(expx+1.0/expx)
    End;

    Function Power(x:Double; n:Integer): Double;
    {Calculate x power n}
    var i,m : integer;
        result :double;
    begin
      result := 1.0;
      if n=0 then
      begin
        Power:=result;
        exit;
      end;
      m:=  n;
      if n<0 then m:=-n;
      for i:=1 to m do result :=x*result;
      Power :=result;
      if n<0 then Power:=1.0/result;
    end;


  Procedure CGAMA(X,Y:Double; KF:Integer; Var GR,GI:Double);
{ ===========================================================
        Purpose: Compute the gamma function G(z) or Ln[G(z)]
                 for a complex argument
        Input :  x  --- Real part of z
                 y  --- Imaginary part of z
                 KF --- Function code
                        KF=0 for Ln[G(z)]
                        KF=1 for G(z)
        Output:  GR --- Real part of Ln[G(z)] or G(z)
                 GI --- Imaginary part of Ln[G(z)] or G(z)
  =========================================================== }
  Label Return;
  Var
      A: Array[1..10] of Double;
      G0,GR1,GI1,SR,SI,T,TH,TH1,TH2,X1,Y1,Z1,Z2: Double;
      J,K,NA:Integer;
  Begin
    {Initialize table A}
    A[1]:= 8.333333333333333E-02; A[2]:= -2.777777777777778E-03;
    A[3]:= 7.936507936507937E-04; A[4]:= -5.952380952380952E-04;
    A[5]:= 8.417508417508418E-04; A[6]:= -1.917526917526918E-03;
    A[7]:= 6.410256410256410E-03; A[8]:= -2.955065359477124E-02;
    A[9]:= 1.796443723688307E-01; A[10]:=-1.39243221690590;

    IF (Y = 0.0) AND (X = INT(X)) AND (X <= 0.0) THEN
    begin
       GR:=1E+300;  {arbitrary big number}
       GI:=0.0;
       goto RETURN
    end
    ELSE IF X < 0.0 THEN
    begin
       X1:=X;
       Y1:=Y;
       X:=-X;
       Y:=-Y
    end;
    X0:=X;
    IF X <= 7.0 THEN
    begin
      NA:=Round(7.0-X);
      X0:=X+NA
    end;
    Z1:=SQRT(X0*X0+Y*Y);
    TH:=ArcTan(Y/X0);
    GR:=(X0-0.5)*Ln(Z1)-TH*Y-X0+0.5*Ln(2.0*PI);
    GI:=TH*(X0-0.5)+Y*Ln(Z1)-Y;
    For K:=1 to 10 do
    begin
      T:=Power(Z1,(1-2*K));
      GR:=GR+A[K]*T*COS((2.0*K-1.0)*TH);
      GI:=GI-A[K]*T*SIN((2.0*K-1.0)*TH)
    end;
    IF X <= 7.0 THEN
    begin
      GR1:=0.0;
      GI1:=0.0;
      For J:=0 to NA-1 do
      begin
        GR1:=GR1+0.5*Ln(Sqr(X+J)+Y*Y);
        GI1:=GI1+ArcTan(Y/(X+J))
      end;
      GR:=GR-GR1;
      GI:=GI-GI1
    end;
    IF X1 < 0.0 THEN
    begin
      Z1:=SQRT(X*X+Y*Y);
      TH1:=ArcTan(Y/X);
      SR:=-SIN(PI*X)*COSH(PI*Y);
      SI:=-COS(PI*X)*SINH(PI*Y);
      Z2:=SQRT(SR*SR+SI*SI);
      TH2:=ArcTan(SI/SR);
      IF SR < 0.0 Then TH2:=PI+TH2;
      GR:=Ln(PI/(Z1*Z2))-GR;
      GI:=-TH1-TH2-GI;
      X:=X1;
      Y:=Y1
    end;
    IF KF = 1 THEN
    begin
      G0:=EXP(GR);
      GR:=G0*COS(GI);
      GI:=G0*SIN(GI)
    end;
  Return: End;


  {main program}
  BEGIN
    WRITELN;
    WRITE('  Please enter KF, x and y: ');
    READLN(KF, X, Y);
    WRITELN;
    IF KF = 1 THEN
      WRITELN('        x         y           Re[G(z)]           Im[G(z)]')
    ELSE
      WRITELN('        x         y          Re[LnG(z)]          Im[LnG(z)]');
    WRITELN('    ---------------------------------------------------------');

    CGAMA(X,Y,KF,GR,GI);

    WRITELN(' ',X:10:2,Y:10:2,GR:19:10,GI:19:10);
    WRITELN;

    ReadKey;
    DoneWinCrt

  END.

{end of file mcgama.pas}