]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 |