      SUBROUTINE CEXPLI (MO, A, W)
C     ******************************************************************
C     FORTRAN SUBROUTINE FOR COMPLEX EXPONENTIAL INTEGRAL
C     ******************************************************************
C     MO = MODE OF OPERATION (0 OR 1)
C     A  = ARGUMENT (COMPLEX NUMBER)
C     W  = FUNCTION (COMPLEX NUMBER)
C     -------------------
      COMPLEX A, W
      DIMENSION AZ(2), ZL(2)
      DIMENSION TS(2), SR(2), SM(2), TM(2), QF(2), CD(18), CE(18)
      DATA CD(1) /0.00000000000000E00/,  CD(2) /3.11105957086528E-2/,
     1     CD(3) /1.03661260539112E-1/,  CD(4) /2.16532335244554E-1/,
     2     CD(5) /3.69931427960192E-1/,  CD(6) /5.66766259990589E-1/,
     3     CD(7) /8.14042066324748E-1/,  CD(8) /1.12384247540813E00/,
     4     CD(9) /1.51400478148512E00/,  CD(10)/2.00886795032284E00/,
     5     CD(11)/2.64052411823592E00/,  CD(12)/3.45098449933392E00/,
     6     CD(13)/4.49583360763202E00/,  CD(14)/5.85058263409822E00/,
     7     CD(15)/7.62273501463380E00/,  CD(16)/9.97814501584578E00/,
     8     CD(17)/1.32122064896408E+1/,  CD(18)/1.80322948376021E+1/
      DATA CE(1) /8.50156516121093E-3/,  CE(2) /5.05037465849058E-2/,
     1     CE(3) /8.36817368956407E-2/,  CE(4) /1.07047582417607E-1/,
     2     CE(5) /1.20424719029462E-1/,  CE(6) /1.25096631582229E-1/,
     3     CE(7) /1.22314435224685E-1/,  CE(8) /1.12621417553907E-1/,
     4     CE(9) /9.63419407392582E-2/,  CE(10)/7.47398422757511E-2/,
     5     CE(11)/5.08596135953441E-2/,  CE(12)/2.90822706773628E-2/,
     6     CE(13)/1.32201640530101E-2/,  CE(14)/4.43802939829067E-3/,
     7     CE(15)/9.92612478987576E-4/,  CE(16)/1.26579795112011E-4/,
     8     CE(17)/7.02150908253350E-6/,  CE(18)/9.10281532564632E-8/
C     -------------------
      AZ(1)=REAL(A)
      AZ(2)=AIMAG(A)
  001 ZS=AZ(1)*AZ(1)+AZ(2)*AZ(2)
  002 IF(ZS-1.0)010,010,003
  003 IF(ZS-1600.0)004,005,005
  004 IF(-AZ(1)+0.064*AZ(2)*AZ(2))010,010,008
  005 SM(1)=0.0
      SM(2)=0.0
      QF(1)=+AZ(1)/ZS
      QF(2)=-AZ(2)/ZS
      TM(1)=QF(1)
      TM(2)=QF(2)
      PM=0.0
      GO TO 007
  006 PM=PM+1.0
      TS(1)=TM(1)*QF(1)-TM(2)*QF(2)
      TS(2)=TM(1)*QF(2)+TM(2)*QF(1)
      TM(1)=PM*TS(1)
      TM(2)=PM*TS(2)
      IF(ABS(SM(1))+ABS(TM(1)).NE.ABS(SM(1)))GO TO 007
      IF(ABS(SM(2))+ABS(TM(2)).EQ.ABS(SM(2)))GO TO 014
  007 SM(1)=SM(1)+TM(1)
      SM(2)=SM(2)+TM(2)
      IF(PM.LT.40.0)GO TO 006
      GO TO 014
  008 SM(1)=0.0
      SM(2)=0.0
      DO 009 I=1,18
      TS(1)=AZ(1)-CD(I)
      TS(2)=AZ(2)
      SS=TS(1)*TS(1)+TS(2)*TS(2)
      TM(1)=+CE(I)*TS(1)/SS
      TM(2)=-CE(I)*TS(2)/SS
      SM(1)=SM(1)+TM(1)
      SM(2)=SM(2)+TM(2)
  009 CONTINUE
      GO TO 014
  010 ZL(1)=0.5*ALOG(ZS)
      ZL(2)=ATAN2(-AZ(2),-AZ(1))
      SR(1)=0.577215664901533+ZL(1)+AZ(1)
      SR(2)=ZL(2)+AZ(2)
      SM(1)=0.0
      SM(2)=0.0
      TM(1)=AZ(1)
      TM(2)=AZ(2)
      PM=1.0
  011 PM=PM+1.0
      TS(1)=TM(1)*AZ(1)-TM(2)*AZ(2)
      TS(2)=TM(1)*AZ(2)+TM(2)*AZ(1)
      TM(1)=TS(1)/PM
      TM(2)=TS(2)/PM
      TS(1)=TM(1)/PM
      TS(2)=TM(2)/PM
      IF(ABS(SM(1))+ABS(TS(1)).NE.ABS(SM(1)))GO TO 012
      IF(ABS(SM(2))+ABS(TS(2)).EQ.ABS(SM(2)))GO TO 013
  012 SM(1)=SM(1)+TS(1)
      SM(2)=SM(2)+TS(2)
      GO TO 011
  013 SM(1)=SR(1)+SM(1)
      SM(2)=SR(2)+SM(2)
      GO TO 015
  014 IF(MO.NE.0)GO TO 017
      QM=EXP(AZ(1))
      QF(1)=QM*COS(AZ(2))
      QF(2)=QM*SIN(AZ(2))
      GO TO 016
  015 IF(MO.EQ.0)GO TO 017
      QM=EXP(-AZ(1))
      QF(1)=QM*COS(-AZ(2))
      QF(2)=QM*SIN(-AZ(2))
  016 TS(1)=QF(1)*SM(1)-QF(2)*SM(2)
      TS(2)=QF(1)*SM(2)+QF(2)*SM(1)
      SM(1)=TS(1)
      SM(2)=TS(2)
  017 W=CMPLX(SM(1),SM(2))
      RETURN
      END
