1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
FUNCTION DKLS(NP,I,L,R,Z)
C-----
C THIS ROUTINE CALCULATES THE SINGLE PRECISION INTEGRALS FOR
C AXISYMMETRIC SOLIDS IN EMG
C
C INPUT
C NP = NUMBER OF POINTS (3 OR 4)
C I,L= THE INTEGRAL DESIRED (I SERIES STARTS WITH -1)
C R = RADIUS ARRAY (NP LONG)
C Z = Z-CORD ARRAY (NP LONG)
C
C OUTPUT
C DKL = DESIRED INTEGRAL
C
C-----
INTEGER NAM(2)
REAL R(3),Z(3)
DATA EPS /.01/, NAM /4HDKLS ,1H /
DATA ZERO, ONE, TWO / 0., 1., 2. /
C
DKLS= ZERO
L1 = L+1
L2 = L + 2
DL1 = L1
K = I+1
C
C . LOOP ON NUMBER OF POINTS...
IF (R(1).LE.ZERO) GO TO 300
DO 200 M = 1,NP
J = M+1
IF (M.EQ.NP) J = 1
RA = R(M)
RB = R(J)
ZA = Z(M)
ZB = Z(J)
DR = RB-RA
DZ = ZB-ZA
C
C . TEST IF RADIUS IS .LE. 0 (DRIVER SHOULD FIND THIS)...
IF (RB.LE.ZERO) GO TO 300
GKL = ZERO
PR = RA+RB
AR = PR / TWO
C
C . CHECK FOR APPROXIMATION, DR/AVE(R)...
IF (ABS(DR/AR) .LT. EPS) GO TO 70
C
A = ZA*DR - RA*DZ
BETA = A/DR
C
C . CHECK FOR BETA .EQ. 0 CASE...
IF ( ABS (BETA / AR ) .GT. EPS ) GO TO 10
C
IF (DZ.EQ.ZERO) GO TO 200
LK = L + K + 1
AR = LK
GKL = (DZ/DR)**L1 * (RA**LK-RB**LK) / (DL1*AR)
GO TO 200
C
C . GENERAL CASE...
10 RAK = RA**K
RBK = RB**K
IF ( K ) 300,20,30
C
C . GENERAL CASE, K.EQ.0, CONSTANT TERM...
20 GKL = ALOG(RA/RB)/DL1
GO TO 40
C
C . GENERAL CASE, CONSTANT TERM...
30 AR = K * L1
GKL = (RAK - RBK) / AR
C
C . GENERAL CASE, SUMMATION...
40 IF (DZ.EQ.ZERO) GO TO 65
LFACT = 1
C . CALCULATE FACTORIAL (L+1)...
DO 50 J = 2,L
50 LFACT = LFACT * J
FACTL = LFACT
JFACT = 1
AJ = ONE
DZJ = ONE
LMJF= LFACT * L1
DO 60 J = 1,L1
JFACT = JFACT * J
C . CALCULATE (L+1-J) FACTORIAL IN LMJF...
LMJF = LMJF / (L2-J)
FACT = FACTL / FLOAT (JFACT*LMJF)
DFACT = K + J
DFACT = FACT / DFACT
AJ = AJ * A
RAK = RAK * RA
RBK = RBK * RB
DZJ = DZJ * DZ
60 GKL = GKL + (DFACT * DZJ * (RAK-RBK)) / AJ
C-----
65 GKL = GKL * BETA**L1
GO TO 200
C
C . APPROXIMATE CODE...
70 CONTINUE
IF (DR.EQ.ZERO) GO TO 200
DZJ = L1 * L2
RBK = ZB**L1
J = K - 1
GKL = -DR * AR**J * RBK / DL1
C
IF (DZ.EQ.ZERO) GO TO 200
GKL = GKL + (((2.*RA+RB)/3.)**J *DR*ABS(ZA**L2 - RBK*ZB))/(DZJ*DZ)
C
200 DKLS= DKLS+ GKL
C-----
C
C . ALL DONE
C
210 CONTINUE
RETURN
C
C . ERROR...
C
300 CALL MESAGE (-7,K,NAM)
GO TO 210
END
|
