PageRenderTime 50ms CodeModel.GetById 16ms RepoModel.GetById 0ms app.codeStats 0ms

/UAFWLIS2.4A/XING/YVFK.F

http://ua-fwlis.googlecode.com/
FORTRAN Legacy | 465 lines | 236 code | 132 blank | 97 comment | 9 complexity | d7c0cab968ef0bd8496d2f59e3da4f84 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. 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?'
  83. WRITE(*,*) 'PLEASE CHECK YVYK.FOR'
  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