MINICERN/mathlib/gen/c/relikc64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:11  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       FUNCTION DELIKC(AK)
  12 C
  13 #include "gen/imp64.inc"
  14 C
  15       CHARACTER*(*) NAMEK,NAMEE
  16       PARAMETER(NAMEK='RELIKC/DELIKC',NAMEE='RELIEC/DELIEC')
  17 #endif
  18 #if !defined(CERNLIB_DOUBLE)
  19       FUNCTION RELIKC(AK)
  20 C
  21       CHARACTER*(*) NAMEK,NAMEE
  22       PARAMETER(NAMEK='RELIKC',NAMEE='RELIEC')
  23 #endif
  24       CHARACTER*80 ERRTXT
  25 C
  26 C     Based on
  27 C      W.J. Cody, Chebyshev approximations for the complete elliptic
  28 C      integrals K and E, Math. Comp. l9 (1965) 105-112
  29 C
  30       DIMENSION A(8),B(8),P(8),Q(8)
  31
  32       PARAMETER (Z1 = 1, HF = Z1/2, C = 1.38629 43611 19890 6D0)
  33       DATA A
  34      1/6.49984 43329 39018 0D-4, 6.69055 09906 89793 6D-3,
  35      2 1.38556 01247 15656 0D-2, 1.12089 18554 64409 2D-2,
  36      3 9.65875 79861 75311 3D-3, 1.49789 88178 70462 9D-2,
  37      4 3.08855 73486 75269 4D-2, 9.65735 90797 58901 8D-2/
  38       DATA B
  39      1/1.50491 81783 60188 3D-4, 3.18313 09927 86288 6D-3,
  40      2 1.41053 80776 15804 8D-2, 2.71898 61116 78825 0D-2,
  41      3 3.70683 98934 15542 2D-2, 4.88180 58565 40395 2D-2,
  42      4 7.03124 26464 62736 1D-2, 1.24999 99994 11792 3D-1/
  43       DATA P
  44      1/7.09809 64089 98722 9D-4, 7.33561 64974 29036 5D-3,
  45      2 1.53771 02528 55201 9D-2, 1.30341 46073 73143 2D-2,
  46      3 1.25105 92410 84464 4D-2, 2.18762 20647 18619 8D-2,
  47      4 5.68056 57874 69535 8D-2, 4.43147 18112 15580 6D-1/
  48       DATA Q
  49      1/1.64272 10797 04802 5D-4, 3.48386 79435 89649 2D-3,
  50      2 1.55251 29948 04072 1D-2, 3.03027 47728 41284 8D-2,
  51      3 4.23828 07456 94790 0D-2, 5.85828 39536 55902 4D-2,
  52      4 9.37499 20249 68011 3D-2, 2.49999 99993 61762 2D-1/
  53
  54 #if defined(CERNLIB_DOUBLE)
  55       ENTRY DELLIK(AK)
  56 #endif
  57 #if !defined(CERNLIB_DOUBLE)
  58       ENTRY ELLICK(AK)
  59 #endif
  60
  61       U=AK**2
  62       IF(U .LT. 1) THEN
  63        Y=1-U
  64        PA=A(1)
  65        PB=B(1)
  66        DO 1 I = 2,8
  67        PA=PA*Y+A(I)
  68     1  PB=PB*Y+B(I)
  69        H=C+PA*Y-LOG(Y)*(HF+PB*Y)
  70       ELSE
  71        H=0
  72        WRITE(ERRTXT,105) AK
  73        CALL MTLPRT(NAMEK,'C347.5',ERRTXT)
  74       ENDIF
  75       GO TO 9
  76
  77 #if defined(CERNLIB_DOUBLE)
  78       ENTRY DELIEC(AK)
  79       ENTRY DELLIE(AK)
  80 #endif
  81 #if !defined(CERNLIB_DOUBLE)
  82       ENTRY RELIEC(AK)
  83       ENTRY ELLICE(AK)
  84 #endif
  85
  86       U=AK**2
  87       IF(U .LT. 1) THEN
  88        Y=1-U
  89        PA=P(1)
  90        PB=Q(1)
  91        DO 2 I = 2,8
  92        PA=PA*Y+P(I)
  93     2  PB=PB*Y+Q(I)
  94        H=1+(PA-LOG(Y)*PB)*Y
  95       ELSE IF(U .EQ. 1) THEN
  96        H=1
  97       ELSE
  98        H=0
  99        WRITE(ERRTXT,105) AK
 100        CALL MTLPRT(NAMEE,'C347.6',ERRTXT)
 101       ENDIF
 102 #if defined(CERNLIB_DOUBLE)
 103     9 DELIKC=H
 104 #endif
 105 #if !defined(CERNLIB_DOUBLE)
 106     9 RELIKC=H
 107 #endif
 108       RETURN
 109   105 FORMAT('ILLEGAL  AK = ',1P,D15.8)
 110       END