MINICERN/mathlib/gen/c/dclaus64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.2  1996/04/02 16:23:45  mclareni
   6 * More precise dclaus64 (C326), test added and C344 removed from TESTALL
   7 *
   8 * Revision 1.1.1.1  1996/04/01 15:02:03  mclareni
   9 * Mathlib gen
  10 *
  11 *
  12 #include "gen/pilot.h"
  13 #if defined(CERNLIB_DOUBLE)
  14       FUNCTION DCLAUS(X)
  15       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  16 #endif
  17 #if !defined(CERNLIB_DOUBLE)
  18       FUNCTION RCLAUS(X)
  19 #include "gen/imp64.inc"
  20 #endif
  21
  22       DIMENSION A(0:8),B(0:13)
  23
  24       PARAMETER (R1 = 1, HF =R1/2)
  25       PARAMETER (PI = 3.14159 26535 89793 24D0)
  26       PARAMETER (PI2 = 2*PI, PIH = PI/2, RPIH = 2/PI)
  27
  28       DATA A( 0) / 0.02795 28319 73575 6613D0/
  29       DATA A( 1) / 0.00017 63088 74389 8116D0/
  30       DATA A( 2) / 0.00000 12662 74146 1157D0/
  31       DATA A( 3) / 0.00000 00117 17181 8134D0/
  32       DATA A( 4) / 0.00000 00001 23006 4129D0/
  33       DATA A( 5) / 0.00000 00000 01395 2729D0/
  34       DATA A( 6) / 0.00000 00000 00016 6908D0/
  35       DATA A( 7) / 0.00000 00000 00000 2076D0/
  36       DATA A( 8) / 0.00000 00000 00000 0027D0/
  37
  38       DATA B( 0) / 0.63909 70888 57265 341D0/
  39       DATA B( 1) /-0.05498 05693 01851 716D0/
  40       DATA B( 2) /-0.00096 12619 45950 606D0/
  41       DATA B( 3) /-0.00003 20546 86822 550D0/
  42       DATA B( 4) /-0.00000 13294 61695 426D0/
  43       DATA B( 5) /-0.00000 00620 93601 824D0/
  44       DATA B( 6) /-0.00000 00031 29600 656D0/
  45       DATA B( 7) /-0.00000 00001 66351 954D0/
  46       DATA B( 8) /-0.00000 00000 09196 527D0/
  47       DATA B( 9) /-0.00000 00000 00524 004D0/
  48       DATA B(10) /-0.00000 00000 00030 580D0/
  49       DATA B(11) /-0.00000 00000 00001 820D0/
  50       DATA B(12) /-0.00000 00000 00000 110D0/
  51       DATA B(13) /-0.00000 00000 00000 007D0/
  52
  53       V=MOD(ABS(X),PI2)
  54       S=SIGN(R1,X)
  55       IF(V .GT. PI) THEN
  56        V=PI2-V
  57        S=-S
  58       ENDIF
  59       IF(V .EQ. 0 .OR. V .EQ. PI) THEN
  60        H=0
  61       ELSEIF(V .LT. PIH) THEN
  62        U=RPIH*V
  63        H=2*U**2-1
  64        ALFA=H+H
  65        B1=0
  66        B2=0
  67        DO 1 I = 8,0,-1
  68        B0=A(I)+ALFA*B1-B2
  69        B2=B1
  70     1  B1=B0
  71        H=V*(1-LOG(V)+HF*V**2*(B0-H*B2))
  72       ELSE
  73        U=RPIH*V-2
  74        H=2*U**2-1
  75        ALFA=H+H
  76        B1=0
  77        B2=0
  78        DO 2 I = 13,0,-1
  79        B0=B(I)+ALFA*B1-B2
  80        B2=B1
  81     2  B1=B0
  82        H=(PI-V)*(B0-H*B2)
  83       ENDIF
  84       DCLAUS=S*H
  85       RETURN
  86       END