MINICERN/mathlib/gen/v/rvnspc.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:58  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       FUNCTION RVNSPC(R,RHO,D)
  12
  13       DIMENSION AZ(4),BZ(4),PZ(4),QZ(4)
  14
  15       PARAMETER (Z1 = 1)
  16       PARAMETER (C1 = 4*Z1/3, C2 = 2*Z1/3, C3 = 4*Z1/9, C4 = Z1/3)
  17       PARAMETER (PI = 3.14159 265D0, CL = 1.38629 436D0)
  18       PARAMETER (SF = 4*PI/3, SFH = 2*PI/3, C0 = 2*PI/3-8*Z1/9)
  19       PARAMETER (PIH = PI/2, Z10 = 10, HF = Z1/2, CA = Z10**(-4))
  20
  21       DATA AZ
  22      1/1.45133 8556D-2, 3.74253 9571D-2, 3.58998 0090D-2,
  23      2 9.66633 8350D-2/
  24       DATA BZ
  25      1/4.41839 8230D-3, 3.32852 1016D-2, 6.88029 5505D-2,
  26      2 1.24985 9468D-1/
  27       DATA PZ
  28      1/1.73631 4854D-2, 4.75740 4429D-2, 6.26076 1942D-2,
  29      2 4.43251 5145D-1/
  30       DATA QZ
  31      1/5.26378 9328D-3, 4.06946 8414D-2, 9.20010 9374D-2,
  32      2 2.49983 6641D-1/
  33
  34       RS=ABS(R)
  35       RC=ABS(RHO)
  36       DA=ABS(D)
  37       DR=RS-RC
  38       RS2=RS**2
  39       RS3=RS*RS2
  40       IF(RC .EQ. 0 .OR. RS .EQ. 0 .OR. DA .GE. RS+RC) THEN
  41        V=0
  42       ELSEIF(DR .LE. DA .AND. DA .LE. -DR) THEN
  43        V=SF*RS3
  44       ELSEIF(DA .EQ. RC .AND. RS .EQ. 2*DA) THEN
  45        V=C0*RS3
  46       ELSEIF(DA .EQ. 0) THEN
  47        V=RS3
  48        IF(RS .GT. RC) V=V-SQRT(RS2-RC**2)**3
  49        V=SF*V
  50       ELSE
  51        BP=DA+RC
  52        BM=DA-RC
  53        BP2=BP**2
  54        A=MAX(RS2,BP2)
  55        B=MIN(RS2,BP2)
  56        C=BM**2
  57        AB=A-B
  58        AC=A-C
  59        BC=B-C
  60        S=BP*BM
  61        IF(DA .NE. DR) THEN
  62         Y=AB/AC
  63         YL=LOG(Y)
  64         PA=AZ(1)
  65         PB=BZ(1)
  66         DO 1 I = 2,4
  67         PA=PA*Y+AZ(I)
  68     1   PB=PB*Y+BZ(I)
  69         HK=CL+PA*Y-YL*(HF+PB*Y)
  70         PA=PZ(1)
  71         PB=QZ(1)
  72         DO 2 I = 2,4
  73         PA=PA*Y+PZ(I)
  74     2   PB=PB*Y+QZ(I)
  75         HE=1+(PA-YL*PB)*Y
  76         IF(C .NE. 0) THEN
  77          YKP=SQRT(AB/AC)
  78          EE=YKP
  79          AM0=1
  80          PP=B/C
  81          IF(PP .GT. 0) THEN
  82           CC=1
  83           PP=SQRT(PP)
  84           DD=1/PP
  85          ELSE
  86           GG=1-PP
  87           FF=YKP**2-PP
  88           PP=SQRT(FF/GG)
  89           DD=-BC/(AC*GG*PP)
  90           CC=0
  91          ENDIF
  92     3    FF=CC
  93          CC=DD/PP+CC
  94          GG=EE/PP
  95          DD=2*(FF*GG+DD)
  96          PP=GG+PP
  97          GG=AM0
  98          AM0=YKP+AM0
  99          IF(ABS(GG-YKP) .GT. CA*GG) THEN
 100           YKP=2*SQRT(EE)
 101           EE=YKP*AM0
 102           GO TO 3
 103          ENDIF
 104          H3=PIH*(CC*AM0+DD)/(AM0*(AM0+PP))
 105         ENDIF
 106        ENDIF
 107        IF(DA .LT. DR) THEN
 108         IF(C .EQ. 0) THEN
 109          V=SFH*RS3+C3*SQRT(A)*(AB*HK-2*(A+AB)*HE)
 110         ELSE
 111          V=C1*(H3*A**2*S/C-HK*(A*S-C4*AB*AC)-HE*AC*(S+C2*(AB+AC)))/
 112      1     SQRT(AC)
 113          IF(RC .GT. DA) V=V+SF*RS3
 114         ENDIF
 115        ELSEIF(DA .EQ. DR) THEN
 116         V=C1*(RS3*ATAN2(2*SQRT(DA*RC),BM)-SQRT(AC)*(S+C2*AC))
 117        ELSE
 118         IF(C .EQ. 0) THEN
 119          V=SFH*RS3+C3*(AB*(B-2*AB)*HK+2*A*(AB-B)*HE)/SQRT(A)
 120         ELSE
 121          V=C1*(H3*B**2*S/C+HK*(S*(AB-B)+C4*AB*(BC-2*AB))
 122      1    -HE*AC*(S-C2*(AB-BC)))/SQRT(AC)
 123          IF(RC .GT. DA) V=V+SF*RS3
 124         ENDIF
 125        ENDIF
 126       ENDIF
 127       RVNSPC=V
 128       RETURN
 129       END
 130 #endif