      SUBROUTINE ELLPCI (ARC, PK, QK, FF, FE)
C     ******************************************************************
C     FORTRAN SUBROUTINE FOR ELLIPTIC INTEGRALS OF FIRST AND SECOND KIND
C     ******************************************************************
C     ARC = ARGUMENT = PHI              (ABS(PHI).LE.PI/2.0)
C     PK = MODULUS = K                  (ABS(PK).LE.1.0)
C     QK = COMODULUS = SQRT (1 - K*K)   (ABS(QK).LE.1.0)
C     FF = ELLIPTIC INTEGRAL OF FIRST KIND = F(PHI, K)
C     FE = ELLIPTIC INTEGRAL OF SECOND KIND = E(PHI, K)
C
  001 PS=PK*PK
      QS=QK*QK
      SN=SIN(ARC)
      CN=COS(ARC)
      SS=SN*SN
      PX=ABS(PK*SN)
      QX=ABS(PK*CN)
      IF(PX.GE.7.61594155955765E-01)GO TO 006
C
C     SERIES EXPANSION FOR ABS(K*SIN(PHI)).LE.TANH(1.0)
C
  002 RA=ARC
      RF=1.0
      S1=0.0
      S2=0.0
      PN=1.0
      QN=2.0
      TR=ARC*SS
      TS=SN*CN
      GO TO 004
  003 PN=QN+1.0
      QN=PN+1.0
      TR=SS*TR
      TS=SS*TS
  004 R0=S1
      RA=(PN*RA-TS)/QN
      RF=PS*RF/QN
      S2=S2+RF*RA
      RF=PN*RF
      S1=S1+RF*RA
      IF(ABS(TR).LT.ABS(RA))GO TO 005
      IF(ABS(S1).GT.ABS(R0))GO TO 003
  005 FF=ARC+S1
      FE=ARC-S2
      RETURN
C
C     SERIES EXPANSION FOR ABS(K*SIN(PHI)).GT.TANH(1.0)
C
  006 DS=QK*QK+QX*QX
      DR=SQRT(DS)
      SI=1.0
      SJ=1.0
      SK=0.0
      RM=0.0
      RN=0.0
      S1=0.0
      S2=0.0
      S3=0.0
      S4=0.0
      TD=QX*DR
      DN=2.0
      GO TO 008
  007 SI=RI
      SJ=RJ
      SK=RK
      DN=DN+2.0
      TD=DS*TD
  008 PN=(DN-1.0)/DN
      QN=(DN+1.0)/(DN+2.0)
      RI=PN*SI
      RJ=PN*PN*QS*SJ
      RK=SK+2.0/(DN*(DN-1.0))
      R0=TD/DN
      RM=QN*QN*QS*(RM-R0*RI)
      RN=PN*QN*QS*(RN-R0*SI)
      D1=RJ
      D2=QN*RJ
      D3=RM-RJ*RK
      D4=RN-PN*QS*SJ*RK+QS*SJ/(DN*DN)
      R0=S3
      S1=S1+D1
      S2=S2+D2
      S3=S3+D3
      S4=S4+D4
      IF(S3.LT.R0)GO TO 007
  009 TR=DR+QX
      RQ=TR/4.0
      PL=-ALOG(RQ)
      TS=1.0+PX
      SQ=TS/2.0
      QL=+ALOG(SQ)
      T1=(1.0+S1)*PL+QX/DR*QL
      T2=(0.5+S2)*QS*PL+1.0-QX*DR/TS
      FF=T1+S3
      FE=T2+S4
      IF(ARC.GE.0.0)RETURN
      FF=-FF
      FE=-FE
      RETURN
      END
