      SUBROUTINE CGAMMA (MO, Z, W)
C     ******************************************************************
C     EVALUATION OF THE COMPLEX GAMMA AND LOGGAMMA FUNCTIONS
C     ******************************************************************
C     MO IS AN INTEGER, Z A COMPLEX ARGUMENT, AND W A COMPLEX VARIABLE.
C     WHEN CGAMMA IS CALLED, W IS ASSIGNED THE VALUE GAMMA(Z) IF MO = 0
C     AND THE VALUE LOG(GAMMA(Z)) IF MO IS NOT ZERO.
C     -------------------
      COMPLEX Z, W
      REAL COEFF(18)
C     -------------------
C     DPI = 2.0 * PI
C     ALPI = LOG(PI)
C     HL2P = 0.5 * LOG(DPI)
C     -------------------
      DATA PI/3.14159265358979324/
      DATA DPI/6.28318530717958648/
      DATA ALPI/1.14472988584940017/
      DATA HL2P/.918938533204672742/
C     -------------------
C     COEFFICIENTS OF A RATIONAL APPROXIMATION FOR 1.0 - EXP(U)
C     WHEN -0.1054 .LE. U .LE. 0.0
C     -------------------
      DATA F0/840.07385296052619/,F1/20.001230821894200/
      DATA D0/1680.1477059210524/,D1/180.01477047052042/
C     -------------------
C     COEFFICIENTS FOR THE ASYMPTOTIC EXPANSION
C     -------------------
      DATA COEFF(1) /+8.333333333333333E-2/,
     1     COEFF(2) /-2.777777777777777E-3/,
     2     COEFF(3) /+7.936507936507937E-4/,
     3     COEFF(4) /-5.952380952380952E-4/,
     4     COEFF(5) /+8.417508417508418E-4/,
     5     COEFF(6) /-1.917526917526918E-3/,
     6     COEFF(7) /+6.410256410256410E-3/,
     7     COEFF(8) /-2.955065359477124E-2/,
     8     COEFF(9) /+1.796443723688306E-1/,
     9     COEFF(10)/-1.392432216905901E00/,
     1     COEFF(11)/+1.340286404416839E+1/,
     2     COEFF(12)/-1.568482846260020E+2/,
     3     COEFF(13)/+2.193103333333333E+3/,
     4     COEFF(14)/-3.610877125372499E+4/,
     5     COEFF(15)/+6.914722688513131E+5/,
     6     COEFF(16)/-1.523822153940742E+7/,
     7     COEFF(17)/+3.829007513914141E+8/,
     8     COEFF(18)/-1.088226603578439E10/
C     -------------------
      X = REAL(Z)
      Y = AIMAG(Z)
      IF (X .GE. 0.0) GO TO 50
C     ------------------------------------------------------------------
C     CASE WHEN THE REAL PART OF Z IS NEGATIVE
C     ------------------------------------------------------------------
      Y = ABS(Y)
      T = -PI*Y
      ET = EXP(T)
      E2T = ET*ET
C
C     SET  A1 = (1 + E2T)/2  AND  A2 = (1 - E2T)/2
C
      A1 = 0.5*(1.0 + E2T)
      IF (T .GE. -.0527) GO TO 10
      A2 = 0.5*(1.0 - E2T)
      GO TO 20
   10 U = T + T
      U2 = U*U
      C = -U*(F0 + F1*U2)
      A2 = C/((U2 + D1)*U2 + D0 + C)
C
C     COMPUTE SIN(PI*X) AND COS(PI*X)
C
   20 K = -X
      U = X + K
      K = MOD(K,2)
      IF (U .GT. -0.5) GO TO 21
      U = U + 1.0
      K = K + 1
   21 U = PI*U
      SN = SIN(U)
      CN = COS(U)
      IF (K .NE. 1) GO TO 30
      SN = -SN
      CN = -CN
C
C     SET  H1 + H2*I  TO  PI/SIN(PI*Z)  OR  LOG(PI/SIN(PI*Z))
C
   30 A1 = SN*A1
      A2 = CN*A2
      A = A1*A1 + A2*A2
      IF (MO .NE. 0) GO TO 40
C
      H1 = A1/A
      H2 = -A2/A
      C = PI*ET
      H1 = C*H1
      H2 = C*H2
      GO TO 41
C
   40 H1 = (ALPI + T) - 0.5*ALOG(A)
      H2 = -ATAN2(A2,A1)
   41 IF (AIMAG(Z) .LT. 0.0) GO TO 42
      X = 1.0 - X
      Y = -Y
      GO TO 50
   42 H2 = -H2
      X = 1.0 - X
C     ------------------------------------------------------------------
C     CASE WHEN THE REAL PART OF Z IS NONNEGATIVE
C     ------------------------------------------------------------------
   50 W1 = 0.0
      W2 = 0.0
      N = 0
      T = X
      Y2 = Y*Y
      A = T*T + Y2
      IF (A .GE. 36.0) GO TO 80
   51 N = N + 1
      T = T + 1.0
      A = T*T + Y2
      IF (A .LT. 36.0) GO TO 51
C
C     LET S1 + S2*I BE THE PRODUCT OF THE TERMS (Z+J)/(Z+N)
C
      U1 = (X*T + Y2)/A
      U2 = Y/A
      S1 = U1
      S2 = N*U2
      IF (N .LT. 2) GO TO 70
      U = T/A
      NM1 = N - 1
      DO 60 J = 1,NM1
      V1 = U1 + J*U
      V2 = (N - J)*U2
      C = S1*V1 - S2*V2
      D = S1*V2 + S2*V1
      S1 = C
   60 S2 = D
C
C     SET  W1 + W2*I = LOG(S1 + S2*I)  WHEN MO IS NONZERO
C
   70 S = S1*S1 + S2*S2
      IF (MO .EQ. 0) GO TO 80
      W1 = 0.5 * ALOG(S)
      W2 = ATAN2(S2,S1)
C
C     SET  V1 + V2*I = (Z - 0.5) * LOG(Z + N) - Z
C
   80 T1 = 0.5 * ALOG(A) - 1.0
      T2 = ATAN2(Y,T)
      U = X - 0.5
      V1 = (U*T1 - 0.5) - Y*T2
      V2 = U*T2 + Y*T1
C
C     LET A1 + A2*I BE THE ASYMPTOTIC SUM
C
      U1 = T/A
      U2 = -Y/A
      Q1 = U1*U1 - U2*U2
      Q2 = 2.0*U1*U2
      A1 = 0.0
      A2 = 0.0
      SM1 = 0.0
      SM2 = 0.0
      DO 91 J = 1,18
      T1 = COEFF(J) * U1
      T2 = COEFF(J) * U2
      SM1 = SM1 + T1
      SM2 = SM2 + T2
      IF (A1 .NE. SM1) GO TO 90
      IF (A2 .EQ. SM2) GO TO 100
   90 A1 = SM1
      A2 = SM2
      T1 = U1*Q1 - U2*Q2
      T2 = U1*Q2 + U2*Q1
      U1 = T1
   91 U2 = T2
C     ------------------------------------------------------------------
C     GATHERING TOGETHER THE RESULTS
C     ------------------------------------------------------------------
  100 W1 = (((A1 + HL2P) - W1) + V1) - N
      W2 = (A2 - W2) + V2
      IF (REAL(Z) .LT. 0.0) GO TO 120
      IF (MO .NE. 0) GO TO 110
C
C     CASE WHEN THE REAL PART OF Z IS NONNEGATIVE AND MO = 0
C
      A = EXP(W1)
      W1 = A * COS(W2)
      W2 = A * SIN(W2)
      IF (N .EQ. 0) GO TO 140
      C = (S1*W1 + S2*W2)/S
      D = (S1*W2 - S2*W1)/S
      W1 = C
      W2 = D
      GO TO 140
C
C     CASE WHEN THE REAL PART OF Z IS NONNEGATIVE AND MO IS NONZERO
C     HERE THE ANGLE W2 IS REDUCED TO THE INTERVAL -PI.LT.W2.LE.PI
C
  110 IF (W2 .GT. PI) GO TO 111
      K = 0.5 - W2/DPI
      AK = K
      W2 = W2 + DPI*AK
      GO TO 140
  111 R = W2/DPI - 0.5
      K = R
      AK = K
      IF (R .GT. AK) GO TO 112
      W2 = PI
      GO TO 140
  112 W2 = W2 - DPI*(AK + 1.0)
      GO TO 140
C
C     CASE WHEN THE REAL PART OF Z IS NEGATIVE AND MO IS NONZERO
C
  120 IF (MO .EQ. 0) GO TO 130
      W1 = H1 - W1
      W2 = H2 - W2
      GO TO 110
C
C
C     CASE WHEN THE REAL PART OF Z IS NEGATIVE AND MO = 0
C
  130 A = EXP(-W1)
      T1 = A * COS(-W2)
      T2 = A * SIN(-W2)
      W1 = H1*T1 - H2*T2
      W2 = H1*T2 + H2*T1
      IF (N .EQ. 0) GO TO 140
      C = W1*S1 - W2*S2
      D = W1*S2 + W2*S1
      W1 = C
      W2 = D
C
C     TERMINATION
C
  140 W = CMPLX(W1,W2)
      RETURN
      END
