MINICERN/mathlib/gen/c/adigam.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:13  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #if defined(CERLIB_OBSOLETE)
  10 #include "gen/pilot.h"
  11 #if defined(CERNLIB_DOUBLE)
  12       FUNCTION ADIGAM(RX)
  13 #include "gen/def64.inc"
  14      +      DDIGAM,DX,X,C(6),ZERO,ONE,HALF,C0,PI,A,V,H,R,P
  15       CHARACTER NAME*6
  16
  17       DATA ZERO /0D0/, ONE /1D0/, HALF /0.5D0/, C0 /15D0/
  18
  19       DATA PI /3.14159 26535 89793 24D0/
  20       DATA C(1) / 8.33333 33333 33333 33D-2/
  21       DATA C(2) /-8.33333 33333 33333 33D-3/
  22       DATA C(3) / 3.96825 39682 53968 25D-3/
  23       DATA C(4) /-4.16666 66666 66666 67D-3/
  24       DATA C(5) / 7.57575 75757 57575 76D-3/
  25       DATA C(6) /-2.10927 96092 79609 28D-2/
  26
  27       X=RX
  28       NAME='ADIGAM'
  29       ADIGAM = ZERO
  30     9 A=ABS(X)
  31       IF(-A .EQ. INT(X)) THEN
  32        WRITE(6,100) NAME,X
  33        RETURN
  34       END IF
  35       V=A
  36       H=ZERO
  37       IF(A .LT. C0) THEN
  38        H=ONE/V
  39        DO 1 I = 1,14-INT(A)
  40        V=V+ONE
  41     1  H=H+ONE/V
  42        V=V+ONE
  43       END IF
  44       R=ONE/V**2
  45       P=R*C(1)
  46       DO 2 I = 6,1,-1
  47     2 P=R*(C(I)+P)
  48       H=LOG(V)-HALF/V-P-H
  49       IF(X .LT. ZERO) H=H+ONE/A+PI/TAN(PI*A)
  50       IF(NAME .EQ. 'DDIGAM') THEN
  51       GO TO 10
  52       ELSE
  53        ADIGAM=SNGL(H+(H-DBLE(SNGL(H))))
  54       END IF
  55       RETURN
  56
  57       ENTRY DDIGAM(DX)
  58       X=DX
  59       NAME='DDIGAM'
  60       DDIGAM = ZERO
  61       GO TO 9
  62   10  DDIGAM = H
  63
  64
  65   100 FORMAT(7X,'***** CERN C317 ',A6,' ... ARGUMENT EQUALS ',
  66      1          'NON-POSITIVE INTEGER =',E10.1)
  67       END
  68 #endif
  69 #if defined(CERNLIB_SINGLE)
  70       FUNCTION ADIGAM(X)
  71       DIMENSION B(6)
  72       DATA PI /3.14159 26535 89793/
  73       DATA B /+8.33333 33333 33333E-2, -8.33333 33333 33333E-3,
  74      1        +3.96825 39682 53968E-3, -4.16666 66666 66667E-3,
  75      2        +7.57575 75757 57576E-3, -2.10927 96092 79609E-2/
  76       A=ABS(X)
  77       IF(-A .EQ. AINT(X)) GO TO 4
  78       V=A
  79       H=0.
  80       IF(A .GE. 15.0) GO TO 3
  81       N=14-INT(A)
  82       H=1.0/V
  83       IF(N .EQ. 0) GO TO 2
  84       DO 1 I = 1,N
  85       V=V+1.0
  86     1 H=H+1.0/V
  87     2 V=V+1.0
  88     3 R=1.0/V**2
  89       ADIGAM=LOG(V)-0.5/V-R*(B(1)+R*(B(2)+R*(B(3)+R*(B(4)+R*(B(5)+R*
  90      1 (B(6)+R*B(1)))))))-H
  91       IF(X .GE. 0.0) RETURN
  92       H=PI*A
  93       ADIGAM=ADIGAM+1.0/A+PI*COS(H)/SIN(H)
  94       RETURN
  95     4 WRITE(6,100)X
  96       ADIGAM=0.
  97       RETURN
  98   100 FORMAT(' ADIGAM ... ARGUMENT IS NON-POSITIVE INTEGER =',F20.2)
  99       END
 100 #endif
 101 #endif