PageRenderTime 52ms CodeModel.GetById 20ms RepoModel.GetById 1ms app.codeStats 0ms

/UAFWLIS2.3/Backup/VimFWLIS/YVFK.F

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