MINICERN/mathlib/gen/c/celfun64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.2  1997/12/15 16:18:35  mclareni
   6 * Changes for the Portland Group f77 compiler inside cpp define CERNLIB_QFPGF77
   7 *
   8 * Revision 1.1.1.1  1996/04/01 15:02:01  mclareni
   9 * Mathlib gen
  10 *
  11 *
  12 #include "gen/pilot.h"
  13 #if defined(CERNLIB_DOUBLE)
  14       SUBROUTINE WELFUN(W,AK2,SN,CN,DN)
  15 C
  16 #include "gen/imp64.inc"
  17 C
  18       CHARACTER*(*) NAME
  19       PARAMETER(NAME='CELFUN/WELFUN')
  20 #endif
  21 #if !defined(CERNLIB_DOUBLE)
  22       SUBROUTINE CELFUN(W,AK2,SN,CN,DN)
  23 C
  24       CHARACTER*(*) NAME
  25       PARAMETER(NAME='CELFUN')
  26 #endif
  27 C
  28 C     Jacobian Elliptic Functions SN(W,M),CN(W,M),DN(W,M) for
  29 C     complex argument  w = u+i*v.
  30 C     Iterates for parameter m = k**2 or mc= 1-m if mc < m to obtain
  31 C     fastest convergence and uses finally Jacobi's
  32 C     imaginary transformation to go back to sn, cn, dn of m.
  33 C
  34 #include "gen/defc64.inc"
  35      +  W,SN,CN,DN,I
  36       DIMENSION C(4)
  37
  38 C     MACHINE-DEPENDENT: EPS1=2**-(MB/2), EPS2=2**-(MB+3)
  39 C     Where M = Number of bits in mantissa
  40
  41       PARAMETER (MB = 64)
  42       PARAMETER (Z0 = 0, Z1 = 1, Z2 = 2, HF = Z1/2, QU = Z1/4)
  43 #if !defined(CERNLIB_CRAY)
  44       PARAMETER (I = (0D0,1D0))
  45 #endif
  46 #if defined(CERNLIB_CRAY)
  47       PARAMETER (I = (0E0,1E0))
  48 #endif
  49       PARAMETER (PI = 3.14159 26535 89793 24D0, PIH = PI/2)
  50       PARAMETER (EPS1 = Z2**(-MB/2), EPS2 = Z2**(-(MB+3)))
  51
  52       CHARACTER*80 ERRTXT
  53
  54 #if defined(CERNLIB_QFPGF77)
  55       DATA AM0 /-1D20/
  56
  57       SAVE AM0,C,A,L,BIGK
  58 #endif
  59
  60 #if defined(CERNLIB_QF2C)
  61 #include "gen/gcmpfun.inc"
  62 #endif
  63
  64 #if ! defined(CERNLIB_QFPGF77)
  65       DATA AM0 /-1D20/
  66
  67       SAVE AM0,C,A,L,BIGK
  68 #endif
  69
  70 #if !defined(CERNLIB_QF2C)
  71 #include "gen/gcmpfun.inc"
  72 #endif
  73
  74       AM=ABS(AK2)
  75       IF(AM .GT. 1) THEN
  76        WRITE(ERRTXT,101) AM
  77        CALL MTLPRT(NAME,'C320.1',ERRTXT)
  78        RETURN
  79       ENDIF
  80       IF(AM .LE. HF) THEN
  81        U=W
  82        V=GIMAG(W)
  83        IF(AM .EQ. AM0) GO TO 1
  84        XK2=AM
  85        B=SQRT(1-AM)
  86       ELSE
  87        U=GIMAG(W)
  88        V=-W
  89        IF(AM .EQ. AM0) GO TO 1
  90        XK2=1-AM
  91        B=SQRT(AM)
  92       ENDIF
  93
  94       AM0=AM
  95       A=1
  96       L=4
  97       C(L)=QU*XK2
  98
  99 C     Gaussian arithmetic-geometric mean. Skipped if previous M.
 100
 101     2 IF(C(L) .GE. EPS1) THEN
 102        L=L-1
 103        C(L)=(QU*(A-B))**2
 104        A1=HF*(A+B)
 105        B=SQRT(A*B)
 106        A=A1
 107        GO TO 2
 108       ENDIF
 109       A=HF*(A+B)
 110       BIGK=PIH/A
 111
 112 C     Descending Landen-Gauss Trafo for real U
 113
 114     1 X=SIN(A*U)
 115       IF(V .NE. 0 .OR. X .NE. 0) THEN
 116        IF(X .NE. 0) THEN
 117         X=A/X
 118         DO 3 J = L,4
 119         X1=C(J)/X
 120    3    X=X1+X
 121         SU=1/X
 122         DU=1-2*X1*SU
 123         CU=SIGN(SQRT(ABS(1-SU**2)),BIGK-ABS(U))
 124        ENDIF
 125        IF(V .NE. 0) THEN
 126         Y=A/SINH(A*V)
 127         DO 4 J = L,4
 128         Y1=C(J)/Y
 129         Y=Y-Y1
 130         IF(Y .EQ. 0) Y=EPS2
 131     4   CONTINUE
 132         SV=1/Y
 133         Y1=2*Y1*SV
 134         DV=1+Y1
 135         CV=SIGN(SQRT((SV+Y)*SV),Y1)
 136
 137 C     Evaluation of complex sn, cn, dn from real and imaginary ones and
 138 C     Jacobi's imaginary argument trafo if necessary
 139
 140         IF(U .NE. 0) THEN
 141          SS=-SU*SV
 142          D=1/(XK2*SS**2+1)
 143          DV=DV*D
 144          CU=CU*D
 145          CC=CU*CV
 146          DD=DU*DV
 147          DN=GCMPLX(DD,XK2*SS*CC)
 148          CN=GCMPLX(CC,SS*DD)
 149          SN=GCMPLX(SU*CV*DV,SV*CU*DU)
 150         ELSE
 151          SN=GCMPLX(Z0,SV)
 152          CN=CV
 153          DN=DV
 154         ENDIF
 155        ELSE
 156         SN=SU
 157         CN=CU
 158         DN=DU
 159        ENDIF
 160        IF(AM .GT. HF) THEN
 161         CN=1/CN
 162         DN=DN*CN
 163         SN=I*SN*CN
 164        ENDIF
 165       ELSE
 166        SN=0
 167        CN=1
 168        DN=1
 169       ENDIF
 170       RETURN
 171   101 FORMAT('MODULUS AK2 = ',1P,E15.6,' OUT OF RANGE')
 172       END