MINICERN/mathlib/gen/v/rvnspc64.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 DVNSPC(R,RHO,D)
  12 C
  13 #include "gen/imp64.inc"
  14 C
  15       CHARACTER*(*) NAME
  16       PARAMETER(NAME='DVNSPC')
  17 #endif
  18 #if !defined(CERNLIB_DOUBLE)
  19       FUNCTION RVNSPC(R,RHO,D)
  20 C
  21       CHARACTER*(*) NAME
  22       PARAMETER(NAME='RVNSPC')
  23 #endif
  24 C
  25 C     Based on
  26 C     F. LAMARCHE and C. LEROY, Evaluation of the volume of
  27 C     a sphere with a cylinder by elliptic integrals,
  28 C     Computer Phys. Comm. 59 (1990) 359-369
  29 C
  30       PARAMETER (Z1 = 1)
  31       PARAMETER (C1 = 4*Z1/3, C2 = 2*Z1/3, C3 = 4*Z1/9, C4 = Z1/3)
  32       PARAMETER (PI = 3.14159 26535 89793 24D0)
  33       PARAMETER (SF = 4*PI/3, SFH = 2*PI/3, C0 = 2*PI/3-8*Z1/9)
  34
  35       RS=ABS(R)
  36       RC=ABS(RHO)
  37       DA=ABS(D)
  38       DR=RS-RC
  39       RS2=RS**2
  40       RS3=RS2*RS
  41       IF(RC .EQ. 0 .OR. RS .EQ. 0 .OR. DA .GE. RS+RC) THEN
  42        V=0
  43       ELSEIF(DR .LE. DA .AND. DA .LE. -DR) THEN
  44        V=SF*RS3
  45       ELSEIF(DA .EQ. RC .AND. RS .EQ. 2*DA) THEN
  46        V=C0*RS3
  47       ELSEIF(DA .EQ. 0) THEN
  48        V=RS3
  49        IF(RS .GT. RC) V=V-SQRT(RS2-RC**2)**3
  50        V=SF*V
  51       ELSE
  52        BP=DA+RC
  53        BM=DA-RC
  54        BP2=BP**2
  55        A=MAX(RS2,BP2)
  56        B=MIN(RS2,BP2)
  57        C=BM**2
  58        AB=A-B
  59        AC=A-C
  60        BC=B-C
  61        S=BP*BM
  62        XK2=BC/AC
  63        XK=SQRT(XK2)
  64        IF(DA .LT. DR) THEN
  65         IF(C .EQ. 0) THEN
  66 #if defined(CERNLIB_DOUBLE)
  67          V=SFH*RS3+C3*SQRT(A)*(AB*DELIKC(XK)-2*(A+AB)*DELIEC(XK))
  68         ELSE
  69          V=C1*(DELI3C(SQRT(1-XK2),XK2,B/C)*A**2*S/C
  70      1    -DELIKC(XK)*(A*S-C4*AB*AC)
  71      2    -DELIEC(XK)*AC*(S+C2*(AB+AC)))/SQRT(AC)
  72 #endif
  73 #if !defined(CERNLIB_DOUBLE)
  74          V=SFH*RS3+C3*SQRT(A)*(AB*RELIKC(XK)-2*(A+AB)*RELIEC(XK))
  75         ELSE
  76          V=C1*(RELI3C(SQRT(1-XK2),XK2,B/C)*A**2*S/C
  77      1    -RELIKC(XK)*(A*S-C4*AB*AC)
  78      2    -RELIEC(XK)*AC*(S+C2*(AB+AC)))/SQRT(AC)
  79 #endif
  80          IF(RC .GT. DA) V=V+SF*RS3
  81         ENDIF
  82        ELSEIF(DA .EQ. DR) THEN
  83         V=C1*(RS3*ATAN2(2*SQRT(DA*RC),BM)-SQRT(AC)*(S+C2*AC))
  84        ELSE
  85         IF(C .EQ. 0) THEN
  86 #if defined(CERNLIB_DOUBLE)
  87          V=SFH*RS3+C3*(AB*(B-2*AB)*DELIKC(XK)+2*A*(AB-B)*DELIEC(XK))/
  88      1     SQRT(A)
  89         ELSE
  90          V=C1*(DELI3C(SQRT(1-XK2),XK2,B/C)*B**2*S/C
  91      1    +DELIKC(XK)*(S*(AB-B)+C4*AB*(BC-2*AB))
  92      2    -DELIEC(XK)*AC*(S-C2*(AB-BC)))/SQRT(AC)
  93 #endif
  94 #if !defined(CERNLIB_DOUBLE)
  95          V=SFH*RS3+C3*(AB*(B-2*AB)*RELIKC(XK)+2*A*(AB-B)*RELIEC(XK))/
  96      1     SQRT(A)
  97         ELSE
  98          V=C1*(RELI3C(SQRT(1-XK2),XK2,B/C)*B**2*S/C
  99      1    +RELIKC(XK)*(S*(AB-B)+C4*AB*(BC-2*AB))
 100      2    -RELIEC(XK)*AC*(S-C2*(AB-BC)))/SQRT(AC)
 101 #endif
 102          IF(RC .GT. DA) V=V+SF*RS3
 103         ENDIF
 104        ENDIF
 105       ENDIF
 106 #if defined(CERNLIB_DOUBLE)
 107       DVNSPC=V
 108 #endif
 109 #if !defined(CERNLIB_DOUBLE)
 110       RVNSPC=V
 111 #endif
 112       RETURN
 113       END