1 : C CYLWRK SOURCE PV 11/03/07 21:16:15 6885
2 : SUBROUTINE CYLWRK(RF,ZF,KFLAG,U,UR,UZ,URR,URZ,UZZ,ITYPE,
3 : $NPN,EDICON,ICOO,ITRAVA)
4 : C
5 : C EVALUATE POTENTIAL AND ITS DERIVATIVES FOR CYLINDER COORDINATES
6 : C-----------------------------------------------------------------------
7 : C
8 : IMPLICIT INTEGER(I-N)
9 : IMPLICIT REAL*8(A-H,O-Z)
10 : C POTENTIAL VALUES
11 : SEGMENT EDICON
12 : INTEGER KSTRT, KSTEP, NMIR, IS
13 : REAL*8 CROT, SROT, SYMFCT
14 : LOGICAL LREAL, LIMAG
15 : ENDSEGMENT
16 : SEGMENT ICOO
17 : REAL*8 X(MV),Y(MV),P(MV),WNODE(MV)
18 : INTEGER LISVO(MV)
19 : ENDSEGMENT
20 : SEGMENT ITRAVA
21 : REAL*8 KENN(M42,2),SIGMA(M42),DELRHO(M42),C(M50,M50)
22 : REAL*8 AK(M50),UM(M50),RM(M50)
23 : INTEGER IL(M50)
24 : ENDSEGMENT
25 : C COMMON // NEU, NPN, NPU
26 : C
27 : C CONSTANTS DESCRIBING SYMMETRY CONDITIONS
28 : C
29 : C
30 : C SOURCE TERMS
31 : C
32 : C TABLES AND WORKING STORE FOR POTENTIAL AND FIELD EDIT
33 : C
34 : C-----------------------------------------------------------------------
35 : C
36 : M50=RM(/1)
37 : NPU = MIN(NPN,19)
38 : KEND = KSTRT + KSTEP * (NPU - IS - 2)
39 : IFLAG = 0
40 : J = 0
41 : C
42 : DO 90 JJ = 1, NPN
43 : C
44 : IF(WNODE(JJ) .EQ. 0.0) GO TO 90
45 : RN = X(JJ)
46 : ZN = Y(JJ)
47 : UN = P(JJ)
48 : C
49 : C SET UP HARMONIC FUNCTIONS FOR LOCATION (RN,ZN)
50 : C NOTE. U(SUBSCRIPT M) IS STORED IN UM(M+3)
51 : C R(SUBSCRIPT M) IS STORED IN RM(M+3)
52 : C THE DO LOOP INDEX K IS (M+1)
53 : C
54 : UM(3) = 1.0
55 : IF(IS .NE. 0) UM(3) = RN
56 : RM(3) = FLOAT(IS+IS) * ZN
57 : C
58 : DO 40 K = 2, KEND
59 : UM(K+2) = (RN*RM(K+1) + ZN*UM(K+1))
60 : RM(K+2) = (ZN*RM(K+1) - RN*UM(K+1)) * FLOAT(K+IS+IS-1) / FLOAT(K)
61 : 40 CONTINUE
62 : C
63 : C SET UP COEFFICIENT MATRIX FOR FIT
64 : C
65 : J = J + 1
66 : I = IS
67 : C(J,1) = WNODE(JJ)
68 : C
69 : IF(KEND .LT. KSTRT) GO TO 80
70 : C
71 : DO 70 K = KSTRT, KEND, KSTEP
72 : I = I + 1
73 : C(J,I) = WNODE(JJ) * UM(3) * UM(K+2)
74 : 70 CONTINUE
75 : C
76 : 80 C(J,I+1) = WNODE(JJ) * UN
77 : C
78 : 90 CONTINUE
79 : C WRITE(6,FMT= ' ('' I ET J MAX '',/,(2I6))') I,J
80 : C WRITE(6,FMT= ' ('' KSTRT KSTEP KEND ITYPE '',/,(4I5))') KSTRT,
81 : C $KSTEP,KEND,ITYPE
82 : C
83 : C-----------------------------------------------------------------------
84 : C PERFORM LEAST-SQUARES FIT
85 : C
85 : CALL GLSQ(C,AK,IL,J,I,ALFA,1.0D-10,1.0D-10,M50)
86 : C
87 : C-----------------------------------------------------------------------
88 : C COMPUTE POTENTIAL AND DERIVATIVES
89 : C
90 : UM(1) = 0.0
91 : UM(2) = 0.0
92 : UM(3) = 1.0
93 : IF(IS .NE. 0) UM(3) = RF
94 : RM(1) = 0.0
95 : RM(2) = FLOAT(IS)
96 : RM(3) = 2.0 * RM(2) * ZF
97 : C
98 : DO 110 K = 2, KEND
99 : UM(K+2) = (RF*RM(K+1) + ZF*UM(K+1))
100 : RM(K+2) = (ZF*RM(K+1) - RF*UM(K+1)) * FLOAT(K+IS+IS-1) / FLOAT(K)
101 : 110 CONTINUE
102 : C
103 : IFLAG = KFLAG
104 : U = AK(1) * FLOAT(IS)
105 : UR = 0.0
106 : UZ = 0.0
107 : URZ = 0.0
108 : UZZ = 0.0
109 : I = IS
110 : C
111 : IF(KEND .LT. KSTRT) GO TO 150
112 : C
113 : DO 120 K = KSTRT, KEND, KSTEP
114 : I = I + 1
115 : U = U + AK(I) * UM(3) * UM(K+2)
116 : XK = AK(I) * FLOAT(K+IS+IS-1)
117 : UR = UR + XK * RM(K+1)
118 : UZ = UZ + XK * UM(K+1)
119 : XK = XK * FLOAT(K+IS+IS-2)
120 : URZ = URZ + XK * RM(K)
121 : UZZ = UZZ + XK * UM(K)
122 : 120 CONTINUE
123 : C
124 : 150 KFLAG = IFLAG
125 : URR = - UZZ
126 : RETURN
127 : C***********************************************************************
128 : END
129 :
130 :
© Cast3M 2003 - All rights reserved.
Disclaimer