#### /UAFWLIS2.4A/XING/YVFK.F

FORTRAN Legacy | 465 lines | 236 code | 132 blank | 97 comment | 9 complexity | d7c0cab968ef0bd8496d2f59e3da4f84 MD5 | raw file
1. COMPLEX*16 FUNCTION YVFK(x_m,y_m,z_m,alpha_mn,
2. & l_m,b_m,x_n,y_n,z_n,l_n,b_n,k,d,w,EPSILON_
3. & ,STORE_KRHO,NUM_KRHO,T_L,T_U,
4. & RHO_MAX, curr_flag_m,
5. & l_n2)
6. ********************************************************************************
7. * THIS FUNCTION IS TO EVALUATE YY_0(ALPHA_MN=0) USING EXACT REACTION FORM
8. * Y MEANS RECTANGULAR BASED HALF CELL
9. ********************************************************************************
10. DOUBLE PRECISION, INTENT(IN) :: X_M,X_N,Y_M,Y_N,z_m,z_n,l_m,l_n
11. DOUBLE PRECISION, INTENT(IN) :: l_n2, ALPHA_MN
12. DOUBLE PRECISION, INTENT(IN) :: b_m,b_n,d,W,RHO_MAX
13. COMPLEX*16, INTENT(IN) :: k,EPSILON_
14. INTEGER,INTENT(IN) :: curr_flag_m, T_L, T_U
15. INTEGER, INTENT(INOUT) :: NUM_KRHO
16. INTEGER :: NUM_GREEN_ARRAY
17. PARAMETER (NUM_GREEN_ARRAY=20000)
18. COMPLEX*16, INTENT(INOUT) :: STORE_KRHO(0:NUM_GREEN_ARRAY)
19. ************* LOCAL VARIABLE **********************
20. INTEGER P, P_STEP,P_START,P_MAX
21. INTEGER E_L, E_U
22. INTEGER M_B,N_B, M_T, M_E
23. DOUBLE PRECISION THETA_MN, COS_ALPHA_MN,SIN_ALPHA_MN
24. DOUBLE PRECISION R_MN_ARRAY(-1:1,-1:1,-1:1,-1:1)
25. DOUBLE PRECISION THETA_ARRAY(-1:1,-1:1,-1:1,-1:1)
26. COMPLEX*16 :: A1MN(-1:1,-1:1,-1:1,-1:1),
27. & A2MN(-1:1,-1:1,-1:1,-1:1)
28. COMPLEX*16 :: A1, A2
29. COMPLEX*16 :: IM, KRHO
30. INTEGER MNUMA
31. PARAMETER (MNUMA=12)
32. INTEGER NUMA,PCKEXP(MNUMA)
33. COMPLEX*16 A(MNUMA),OMEGA_PLUS(MNUMA),OMEGA_MINUS(MNUMA),Z,
34. & H01,H02,H11,H12
35. DOUBLE PRECISION PHZ
36. DOUBLE PRECISION RHO_MN,ZMI,ZNI,RHO_MIN
37. ******************************************************************************
38. INTEGER SD
39. INTEGER :: SD_ILHI_HV
40. DOUBLE PRECISION :: RELERR_ILHI_HV,RELERR_RES_HV
41. COMMON /SIGDIG/ SD
42. COMMON/RELERR_HV/RELERR_ILHI_HV,RELERR_RES_HV
43. COMMON/SD_HV/SD_ILHI_HV
44. ************************************
45. double precision :: RELERR
46. complex*16 :: K_2, term_sum
47. integer :: ECOUNT, COUNT, CHOICE
48. DOUBLE PRECISION :: KRHO_MAX, K_Z
49. double precision :: R_MN_ARRAY_VL(-1:1,-1:1,-1:1,-1:1)
50. double precision :: R_MN_ARRAY_VR(-1:1,-1:1,-1:1,-1:1)
51. double precision :: THETA_ARRAY_VL(-1:1,-1:1,-1:1,-1:1)
52. double precision :: THETA_ARRAY_VR(-1:1,-1:1,-1:1,-1:1)
53. COMPLEX*16 :: A1MN_VL(-1:1,-1:1,-1:1,-1:1),
54. & A2MN_VL(-1:1,-1:1,-1:1,-1:1)
55. COMPLEX*16 :: A1MN_VR(-1:1,-1:1,-1:1,-1:1),
56. & A2MN_VR(-1:1,-1:1,-1:1,-1:1)
57. LOGICAL :: AMN_NE1_VL(-1:1,-1:1,-1:1,-1:1),
58. & AMN_NE2_VL(-1:1,-1:1,-1:1,-1:1)
59. LOGICAL :: AMN_NE1_VR(-1:1,-1:1,-1:1,-1:1),
60. & AMN_NE2_VR(-1:1,-1:1,-1:1,-1:1)
61. DOUBLE PRECISION :: A1MN_SQP1_VL(-1:1,-1:1,-1:1,-1:1),
62. & A2MN_SQP1_VL(-1:1,-1:1,-1:1,-1:1)
63. DOUBLE PRECISION :: A1MN_SQP1_VR(-1:1,-1:1,-1:1,-1:1),
64. & A2MN_SQP1_VR(-1:1,-1:1,-1:1,-1:1)
65. ******************************************************************************
66. COMPLEX*16 YVYK
67. EXTERNAL YVYK
68. EXTERNAL OMEGA,CHANKL, CHANKLH
69. EXTERNAL RHO_THETA1, SGN
70. EXTERNAL FIND_KRHO
71. ****************************************
72. IM = (0.D0, 1.D0)
73. RELERR=5.D0*10.D0**(-SD-1)
74. K_2=K*K
75. P_MAX=20000
76. COS_ALPHA_MN=COS(alpha_mn)
77. SIN_ALPHA_MN=SIN(alpha_mn)
78. IF (ABS((z_m-z_n)/z_m)<1.D-4
79. & .and. T_L==-1 .and. T_U==1 ) THEN
80. WRITE(*,*) 'A HOR CELL IS AT CENTER'
81. WRITE(*,*) 'A FULL VER CELL IS AT CENTER'
82. WRITE(*,*) 'WHAT TO DO NOW?'
84. STOP
85. ENDIF
86. IF (ABS((z_n-d/2.D0)/z_n)<1.D-4) THEN
87. * CHECK TO SEE IF THE HOR. (EXPANSION) CURRENT IS IN THE
88. * CENTER OF THE WAVEGUIDE. THEN EVEN MODES ARE NOT COMPUTED
89. P_STEP=2
90. P_START=1
91. ELSEIF ( ABS((z_m-d/2.D0)/z_m)<1.D-4
92. & .and. T_L==-1 .and. T_U==1) THEN
93. * CHECK TO SEE IF THE FULL VER. (TEST) CURRENT IS IN THE
94. * CENTER OF THE WAVEGUIDE. THEN ODD MODES ARE NOT COMPUTED
95. P_STEP=2
96. P_START=2
97. ELSE
98. * NEED ALL MODES
99. P_STEP=1
100. P_START=1
101. ENDIF
102. ** COMPUTE SOME COMMON USED VARIABLES
103. **** CALCULATE AN ARRAY OF THE REQUIRED DISTANCES AND ANGLES.
104. *--------------------case VL combination--------------
105. CALL RHO_THETA1(x_m,y_m,x_n,y_n,l_m,l_n,b_m,b_n,
106. & alpha_mn,R_MN_ARRAY,THETA_ARRAY,RHO_MIN
107. & ,-1,1,-1,1,-1,1,1,1,COS_ALPHA_MN,SIN_ALPHA_MN)
108. A1MN = -IM*sin(THETA_ARRAY)
109. A2MN = -IM*cos(THETA_ARRAY)
110. *** Compute and Store the repeatedly used Parameters
111. R_MN_ARRAY_VL = R_MN_ARRAY
112. THETA_ARRAY_VL = THETA_ARRAY
113. A1MN_VL = A1MN
114. A2MN_VL = A2MN
115. DO m_t=-1,1,1
116. DO m_e=-1,1,1
117. DO m_b=-1,1,1
118. DO n_b=-1,1,1
119. A1=A1MN(M_T,M_E,M_B,N_B)
120. A2=A2MN(M_T,M_E,M_B,N_B)
121. A1MN_SQP1_VL(M_T,M_E,M_B,N_B)=DBLE(a1*a1+1.D0)
122. A2MN_SQP1_VL(M_T,M_E,M_B,N_B)=DBLE(a2*a2+1.D0)
123. IF (ABS(a1)/=1.D0) THEN
124. AMN_NE1_VL(M_T,M_E,M_B,N_B)=.TRUE.
125. ELSE
126. AMN_NE1_VL(M_T,M_E,M_B,N_B)=.FALSE.
127. END IF
128. IF (ABS(a2)/=1.D0) THEN
129. AMN_NE2_VL(M_T,M_E,M_B,N_B)=.TRUE.
130. ELSE
131. AMN_NE2_VL(M_T,M_E,M_B,N_B)=.FALSE.
132. END IF
133. END DO
134. END DO
135. END DO
136. END DO
137. *--------------------case VL combination ends--------------
138. *--------------------case VR combination--------------
139. CALL RHO_THETA1(x_m,y_m,x_n,y_n,l_m,l_n2,b_m,b_n,
140. & alpha_mn,R_MN_ARRAY,THETA_ARRAY,RHO_MIN
141. & ,-1,1,-1,1,-1,1,1,1,COS_ALPHA_MN,SIN_ALPHA_MN)
142. A1MN = -IM*sin(THETA_ARRAY)
143. A2MN = -IM*cos(THETA_ARRAY)
144. *** Compute and Store the repeatedly used Parameters
145. R_MN_ARRAY_VR = R_MN_ARRAY
146. THETA_ARRAY_VR = THETA_ARRAY
147. A1MN_VR = A1MN
148. A2MN_VR = A2MN
149. DO m_t=-1,1,1
150. DO m_e=-1,1,1
151. DO m_b=-1,1,1
152. DO n_b=-1,1,1
153. A1=A1MN(M_T,M_E,M_B,N_B)
154. A2=A2MN(M_T,M_E,M_B,N_B)
155. A1MN_SQP1_VR(M_T,M_E,M_B,N_B)=DBLE(a1*a1+1.D0)
156. A2MN_SQP1_VR(M_T,M_E,M_B,N_B)=DBLE(a2*a2+1.D0)
157. IF (ABS(a1)/=1.D0) THEN
158. AMN_NE1_VR(M_T,M_E,M_B,N_B)=.TRUE.
159. ELSE
160. AMN_NE1_VR(M_T,M_E,M_B,N_B)=.FALSE.
161. END IF
162. IF (ABS(a2)/=1.D0) THEN
163. AMN_NE2_VR(M_T,M_E,M_B,N_B)=.TRUE.
164. ELSE
165. AMN_NE2_VR(M_T,M_E,M_B,N_B)=.FALSE.
166. END IF
167. END DO
168. END DO
169. END DO
170. END DO
171. NUMA=4
172. *------------------case VR combination ends--------------------
173. *-------------------------------------------------------------
174. E_L = -1
175. E_U = 1
176. cx IF (RHO_MIN .GT. -LOG(1.D-5)/KRHO_MAX/1.5D0) THEN
177. IF (RHO_MIN .GT. RHO_MAX/1.5D0) then
178. * USE approximation
179. YVFK =0.D0
180. P=1
181. COUNT=-1
182. ECOUNT=0
183. DO
184. TERM_sum =0.D0
185. if (p>P_max) exit
186. END DO
187. *--------------------------------------
188. *======================================
189. * NEED USE EXACT FORM OF REACTION
190. ELSE
191. p_max=20000
192. YVFK =0.D0
193. *------------------------------------------------------------
194. *============================================================
195. * ------------begin pole loop for exponential term summation ---------
196. P=1
197. COUNT=-1
198. ECOUNT=0
199. CHOICE =1
200. DO
201. CALL FIND_KRHO(KRHO,K,K_Z,P,NUM_KRHO,P_STEP,STORE_KRHO
202. & ,D, NUM_GREEN_ARRAY)
203. TERM_sum =0.D0
204. if (p>P_max) exit
205. if(1) then
206. * ----------- 1. VL combination begins----------------
207. TERM_sum = TERM_sum+
208. & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
209. & ,z_n,l_n,b_n,k,d,w,EPSILON_,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_L
210. & ,RHO_MIN, RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
211. & R_MN_ARRAY_VL, THETA_ARRAY_VL, A1MN_VL, A2MN_VL,
212. & AMN_NE1_VL, AMN_NE2_VL, A1MN_SQP1_VL, A2MN_SQP1_VL )
213. end if
214. *-----------------1. VL combination ends------------
215. if (1) then
216. *-----------------2. VR combination begins----------
217. TERM_sum = TERM_sum+
218. & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
219. & ,z_n,l_n2,b_n,k,d,w,EPSILON_,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_U
220. & ,RHO_MIN,RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
221. & R_MN_ARRAY_VR, THETA_ARRAY_VR, A1MN_VR, A2MN_VR,
222. & AMN_NE1_VR, AMN_NE2_VR, A1MN_SQP1_VR, A2MN_SQP1_VR )
223. *-----------------2. VR combination ends------------
224. end if
225. YVFK = YVFK + TERM_sum
226. IF (( ABS(TERM_sum/YVFK) .LT. RELERR_RES_HV) .or.
227. & (TERM_sum==0 .and. YVFK==0) ) THEN
228. IF (p .EQ. COUNT+p_step) THEN
229. * SINCE THE ODD MODES YIELD ZERO RESULTS WHEN EITHER TEST OR EXPANSION
230. * FUNCTION IS AT THE MIDDLE OF THE PARALLEL PLATES. WE NEED TO MAKE SURE
231. * THAT THE RELATIVE ERROR BOUND IS SATISFIED AT LEAST 3 TIMES
232. * LATER COMMENT: NOTE THAT THIS IS POSSIBLY TAKEN CARE OF BY ONLY SUMMING
233. * THE EVEN RESIDUES FOR TEST AND EXPANSION FUNCTIONS SITTING AT THE CENTER
234. * OF THE PARALLEL PLATE WAVEGUIDE
235. IF(ECOUNT==1) EXIT
236. ECOUNT=ECOUNT+1
237. ELSE
238. ECOUNT=0
239. ENDIF
240. COUNT=p
241. END IF
242. * GO BACK TO CALCULATE NEW RESIDUE VALUE FOR EXPONENTIALS
243. p=p+p_step
244. END DO
245. c print *, "p=", p, YVFK
246. *-------------------------------END of Exponential term summation-------------
247. *
248. *
249. *------------------------------BEGIN pole loop for ILHI term summation-------------------
250. P=1
251. COUNT=-1
252. ECOUNT=0
253. CHOICE =2
254. DO
255. CALL FIND_KRHO(KRHO,K,K_Z,P,NUM_KRHO,P_STEP,STORE_KRHO
256. & ,D, NUM_GREEN_ARRAY)
257. TERM_sum =0.D0
258. if (p>P_max) exit
259. if(1) then
260. * ----------- 1. VL combination begins----------------
261. TERM_sum = TERM_sum+
262. & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
263. & ,z_n,l_n,b_n,k,d,w,EPSILON_,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_L
264. & ,RHO_MIN,RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
265. & R_MN_ARRAY_VL, THETA_ARRAY_VL, A1MN_VL, A2MN_VL,
266. & AMN_NE1_VL, AMN_NE2_VL, A1MN_SQP1_VL, A2MN_SQP1_VL )
267. end if
268. *-----------------1. VL combination ends------------
269. if (1) then
270. *-----------------2. VR combination begins----------
271. TERM_sum = TERM_sum+
272. & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
273. & ,z_n,l_n2,b_n,k,d,w,EPSILON_,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_U
274. & ,RHO_MIN,RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
275. & R_MN_ARRAY_VR, THETA_ARRAY_VR, A1MN_VR, A2MN_VR,
276. & AMN_NE1_VR, AMN_NE2_VR, A1MN_SQP1_VR, A2MN_SQP1_VR )
277. *-----------------2. VR combination ends------------
278. end if
279. YVFK = YVFK + TERM_sum
280. IF (( ABS(TERM_sum/YVFK) .LT. RELERR_RES_HV) .or.
281. & (TERM_sum==0 .and. YVFK==0) ) THEN
282. IF (p .EQ. COUNT+p_step) THEN
283. * SINCE THE ODD MODES YIELD ZERO RESULTS WHEN EITHER TEST OR EXPANSION
284. * FUNCTION IS AT THE MIDDLE OF THE PARALLEL PLATES. WE NEED TO MAKE SURE
285. * THAT THE RELATIVE ERROR BOUND IS SATISFIED AT LEAST 3 TIMES
286. * LATER COMMENT: NOTE THAT THIS IS POSSIBLY TAKEN CARE OF BY ONLY SUMMING
287. * THE EVEN RESIDUES FOR TEST AND EXPANSION FUNCTIONS SITTING AT THE CENTER
288. * OF THE PARALLEL PLATE WAVEGUIDE
289. IF(ECOUNT==1) EXIT
290. ECOUNT=ECOUNT+1
291. ELSE
292. ECOUNT=0
293. ENDIF
294. COUNT=p
295. END IF
296. * GO BACK TO CALCULATE NEW RESIDUE VALUE FOR EXPONENTIALS
297. p=p+p_step
298. END DO
299. c print *, "p=", p, YVFK
300. *-------------------------------END pole loop for ILHI term summation-------------
301. *------------------------------BEGIN Brach point calculation-------------------
302. KRHO = 0
303. k_z = k
304. TERM_sum =0.D0
305. CHOICE =3
306. * ----------- 1. VL combination----------------
307. if (1) then
308. TERM_sum = TERM_sum+
309. & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
310. & ,z_n,l_n,b_n,k,d,w,EPSILON_,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_L
311. & ,RHO_MIN,RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
312. & R_MN_ARRAY_VL, THETA_ARRAY_VL, A1MN_VL, A2MN_VL,
313. & AMN_NE1_VL, AMN_NE2_VL, A1MN_SQP1_VL, A2MN_SQP1_VL )
314. *-----------------1. VL combination ends------------
315. end if
316. if (1) then
317. *-----------------2. VR combination begins----------
318. TERM_sum = TERM_sum+
319. & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
320. & ,z_n,l_n2,b_n,k,d,w,EPSILON_,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_U
321. & ,RHO_MIN,RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
322. & R_MN_ARRAY_VR, THETA_ARRAY_VR, A1MN_VR, A2MN_VR,
323. & AMN_NE1_VR, AMN_NE2_VR, A1MN_SQP1_VR, A2MN_SQP1_VR )
324. *-----------------2. VR combination ends------------
325. end if
326. YVFK = YVFK + TERM_sum
327. c print *, TERM_sum
328. *-------------------------------END of Branch Point summation
329. *
330. END IF
331. * above endif is the judge for exact/approximation
332. c print *, YVFK
333. END FUNCTION YVFK