/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
- COMPLEX*16 FUNCTION YVFK(x_m,y_m,z_m,alpha_mn,
- & l_m,b_m,x_n,y_n,z_n,l_n,b_n,k,d,w,epsilon
- & ,STORE_KRHO,NUM_KRHO,T_L,T_U,
- & RHO_MAX, curr_flag_m,
- & l_n2)
- ********************************************************************************
- * THIS FUNCTION IS TO EVALUATE YY_0(ALPHA_MN=0) USING EXACT REACTION FORM
- * Y MEANS RECTANGULAR BASED HALF CELL
- ********************************************************************************
-
- $DECLARE
-
- DOUBLE PRECISION, INTENT(IN) :: X_M,X_N,Y_M,Y_N,z_m,z_n,l_m,l_n
- DOUBLE PRECISION, INTENT(IN) :: l_n2, ALPHA_MN
- DOUBLE PRECISION, INTENT(IN) :: b_m,b_n,d,W,RHO_MAX
- COMPLEX*16, INTENT(IN) :: k,EPSILON
- INTEGER,INTENT(IN) :: curr_flag_m, T_L, T_U
- INTEGER, INTENT(INOUT) :: NUM_KRHO
-
- INTEGER :: NUM_GREEN_ARRAY
- PARAMETER (NUM_GREEN_ARRAY=20000)
- COMPLEX*16, INTENT(INOUT) :: STORE_KRHO(0:NUM_GREEN_ARRAY)
-
-
-
- ************* LOCAL VARIABLE **********************
- INTEGER P, P_STEP,P_START,P_MAX
- INTEGER E_L, E_U
- INTEGER M_B,N_B, M_T, M_E
- DOUBLE PRECISION THETA_MN, COS_ALPHA_MN,SIN_ALPHA_MN
-
- DOUBLE PRECISION R_MN_ARRAY(-1:1,-1:1,-1:1,-1:1)
- DOUBLE PRECISION THETA_ARRAY(-1:1,-1:1,-1:1,-1:1)
- COMPLEX*16 :: A1MN(-1:1,-1:1,-1:1,-1:1),
- & A2MN(-1:1,-1:1,-1:1,-1:1)
-
- COMPLEX*16 :: A1, A2
- COMPLEX*16 :: IM, KRHO
-
- INTEGER MNUMA
- PARAMETER (MNUMA=12)
- INTEGER NUMA,PCKEXP(MNUMA)
- COMPLEX*16 A(MNUMA),OMEGA_PLUS(MNUMA),OMEGA_MINUS(MNUMA),Z,
- & H01,H02,H11,H12
- DOUBLE PRECISION PHZ
-
- DOUBLE PRECISION RHO_MN,ZMI,ZNI,RHO_MIN
- ******************************************************************************
-
- INTEGER SD
- INTEGER :: SD_ILHI_HV
- DOUBLE PRECISION :: RELERR_ILHI_HV,RELERR_RES_HV
-
- COMMON /SIGDIG/ SD
- COMMON/RELERR_HV/RELERR_ILHI_HV,RELERR_RES_HV
- COMMON/SD_HV/SD_ILHI_HV
-
-
- ************************************
-
- double precision :: RELERR
- complex*16 :: K_2, term_sum
- integer :: ECOUNT, COUNT, CHOICE
- DOUBLE PRECISION :: KRHO_MAX, K_Z
-
-
-
- double precision :: R_MN_ARRAY_VL(-1:1,-1:1,-1:1,-1:1)
- double precision :: R_MN_ARRAY_VR(-1:1,-1:1,-1:1,-1:1)
- double precision :: THETA_ARRAY_VL(-1:1,-1:1,-1:1,-1:1)
- double precision :: THETA_ARRAY_VR(-1:1,-1:1,-1:1,-1:1)
-
-
- COMPLEX*16 :: A1MN_VL(-1:1,-1:1,-1:1,-1:1),
- & A2MN_VL(-1:1,-1:1,-1:1,-1:1)
-
- COMPLEX*16 :: A1MN_VR(-1:1,-1:1,-1:1,-1:1),
- & A2MN_VR(-1:1,-1:1,-1:1,-1:1)
-
- LOGICAL :: AMN_NE1_VL(-1:1,-1:1,-1:1,-1:1),
- & AMN_NE2_VL(-1:1,-1:1,-1:1,-1:1)
-
- LOGICAL :: AMN_NE1_VR(-1:1,-1:1,-1:1,-1:1),
- & AMN_NE2_VR(-1:1,-1:1,-1:1,-1:1)
-
- DOUBLE PRECISION :: A1MN_SQP1_VL(-1:1,-1:1,-1:1,-1:1),
- & A2MN_SQP1_VL(-1:1,-1:1,-1:1,-1:1)
-
- DOUBLE PRECISION :: A1MN_SQP1_VR(-1:1,-1:1,-1:1,-1:1),
- & A2MN_SQP1_VR(-1:1,-1:1,-1:1,-1:1)
-
- ******************************************************************************
- COMPLEX*16 YVYK
- EXTERNAL YVYK
-
- EXTERNAL OMEGA,CHANKL, CHANKLH
- EXTERNAL RHO_THETA1, SGN
- EXTERNAL FIND_KRHO
-
- ****************************************
- IM = (0.D0, 1.D0)
- RELERR=5.D0*10.D0**(-SD-1)
- K_2=K*K
- P_MAX=20000
-
- COS_ALPHA_MN=COS(alpha_mn)
- SIN_ALPHA_MN=SIN(alpha_mn)
-
-
- IF (ABS((z_m-z_n)/z_m)<1.D-4
- & .and. T_L==-1 .and. T_U==1 ) THEN
- WRITE(*,*) 'A HOR CELL IS AT CENTER'
- WRITE(*,*) 'A FULL VER CELL IS AT CENTER'
- WRITE(*,*) 'WHAT TO DO NOW?'
- WRITE(*,*) 'PLEASE CHECK YVYK.FOR'
- STOP
- ENDIF
-
-
- IF (ABS((z_n-d/2.D0)/z_n)<1.D-4) THEN
- * CHECK TO SEE IF THE HOR. (EXPANSION) CURRENT IS IN THE
- * CENTER OF THE WAVEGUIDE. THEN EVEN MODES ARE NOT COMPUTED
- P_STEP=2
- P_START=1
- ELSEIF ( ABS((z_m-d/2.D0)/z_m)<1.D-4
- & .and. T_L==-1 .and. T_U==1) THEN
- * CHECK TO SEE IF THE FULL VER. (TEST) CURRENT IS IN THE
- * CENTER OF THE WAVEGUIDE. THEN ODD MODES ARE NOT COMPUTED
- P_STEP=2
- P_START=2
- ELSE
- * NEED ALL MODES
- P_STEP=1
- P_START=1
- ENDIF
-
-
-
-
-
- ** COMPUTE SOME COMMON USED VARIABLES
- **** CALCULATE AN ARRAY OF THE REQUIRED DISTANCES AND ANGLES.
-
-
- *--------------------case VL combination--------------
-
- CALL RHO_THETA1(x_m,y_m,x_n,y_n,l_m,l_n,b_m,b_n,
- & alpha_mn,R_MN_ARRAY,THETA_ARRAY,RHO_MIN
- & ,-1,1,-1,1,-1,1,1,1,COS_ALPHA_MN,SIN_ALPHA_MN)
-
- A1MN = -IM*sin(THETA_ARRAY)
- A2MN = -IM*cos(THETA_ARRAY)
-
- *** Compute and Store the repeatedly used Parameters
-
- R_MN_ARRAY_VL = R_MN_ARRAY
- THETA_ARRAY_VL = THETA_ARRAY
- A1MN_VL = A1MN
- A2MN_VL = A2MN
-
- DO m_t=-1,1,1
- DO m_e=-1,1,1
- DO m_b=-1,1,1
- DO n_b=-1,1,1
-
- A1=A1MN(M_T,M_E,M_B,N_B)
- A2=A2MN(M_T,M_E,M_B,N_B)
-
- A1MN_SQP1_VL(M_T,M_E,M_B,N_B)=DBLE(a1*a1+1.D0)
- A2MN_SQP1_VL(M_T,M_E,M_B,N_B)=DBLE(a2*a2+1.D0)
-
- IF (ABS(a1)/=1.D0) THEN
- AMN_NE1_VL(M_T,M_E,M_B,N_B)=.TRUE.
- ELSE
- AMN_NE1_VL(M_T,M_E,M_B,N_B)=.FALSE.
- END IF
- IF (ABS(a2)/=1.D0) THEN
- AMN_NE2_VL(M_T,M_E,M_B,N_B)=.TRUE.
- ELSE
- AMN_NE2_VL(M_T,M_E,M_B,N_B)=.FALSE.
- END IF
-
- END DO
- END DO
- END DO
- END DO
-
- *--------------------case VL combination ends--------------
-
- *--------------------case VR combination--------------
-
- CALL RHO_THETA1(x_m,y_m,x_n,y_n,l_m,l_n2,b_m,b_n,
- & alpha_mn,R_MN_ARRAY,THETA_ARRAY,RHO_MIN
- & ,-1,1,-1,1,-1,1,1,1,COS_ALPHA_MN,SIN_ALPHA_MN)
-
- A1MN = -IM*sin(THETA_ARRAY)
- A2MN = -IM*cos(THETA_ARRAY)
-
- *** Compute and Store the repeatedly used Parameters
-
- R_MN_ARRAY_VR = R_MN_ARRAY
- THETA_ARRAY_VR = THETA_ARRAY
- A1MN_VR = A1MN
- A2MN_VR = A2MN
-
- DO m_t=-1,1,1
- DO m_e=-1,1,1
- DO m_b=-1,1,1
- DO n_b=-1,1,1
-
- A1=A1MN(M_T,M_E,M_B,N_B)
- A2=A2MN(M_T,M_E,M_B,N_B)
- A1MN_SQP1_VR(M_T,M_E,M_B,N_B)=DBLE(a1*a1+1.D0)
- A2MN_SQP1_VR(M_T,M_E,M_B,N_B)=DBLE(a2*a2+1.D0)
-
- IF (ABS(a1)/=1.D0) THEN
- AMN_NE1_VR(M_T,M_E,M_B,N_B)=.TRUE.
- ELSE
- AMN_NE1_VR(M_T,M_E,M_B,N_B)=.FALSE.
- END IF
- IF (ABS(a2)/=1.D0) THEN
- AMN_NE2_VR(M_T,M_E,M_B,N_B)=.TRUE.
- ELSE
- AMN_NE2_VR(M_T,M_E,M_B,N_B)=.FALSE.
- END IF
-
- END DO
- END DO
- END DO
- END DO
-
-
- NUMA=4
-
-
- *------------------case VR combination ends--------------------
- *-------------------------------------------------------------
- E_L = -1
- E_U = 1
-
-
-
- cx IF (RHO_MIN .GT. -LOG(1.D-5)/KRHO_MAX/1.5D0) THEN
- IF (RHO_MIN .GT. RHO_MAX/1.5D0) then
- * USE approximation
-
- YVFK =0.D0
- P=1
- COUNT=-1
- ECOUNT=0
-
- DO
- TERM_sum =0.D0
- if (p>P_max) exit
-
-
-
- END DO
-
- *--------------------------------------
- *======================================
- * NEED USE EXACT FORM OF REACTION
- ELSE
-
-
- p_max=20000
- YVFK =0.D0
-
- *------------------------------------------------------------
- *============================================================
-
- * ------------begin pole loop for exponential term summation ---------
- P=1
- COUNT=-1
- ECOUNT=0
-
- CHOICE =1
- DO
-
-
- CALL FIND_KRHO(KRHO,K,K_Z,P,NUM_KRHO,P_STEP,STORE_KRHO
- & ,D, NUM_GREEN_ARRAY)
-
-
-
- TERM_sum =0.D0
- if (p>P_max) exit
- if(1) then
- * ----------- 1. VL combination begins----------------
-
- TERM_sum = TERM_sum+
- & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
- & ,z_n,l_n,b_n,k,d,w,epsilon,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_L
- & ,RHO_MIN, RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
- & R_MN_ARRAY_VL, THETA_ARRAY_VL, A1MN_VL, A2MN_VL,
- & AMN_NE1_VL, AMN_NE2_VL, A1MN_SQP1_VL, A2MN_SQP1_VL )
-
- end if
-
- *-----------------1. VL combination ends------------
-
- if (1) then
- *-----------------2. VR combination begins----------
-
- TERM_sum = TERM_sum+
- & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
- & ,z_n,l_n2,b_n,k,d,w,epsilon,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_U
- & ,RHO_MIN,RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
- & R_MN_ARRAY_VR, THETA_ARRAY_VR, A1MN_VR, A2MN_VR,
- & AMN_NE1_VR, AMN_NE2_VR, A1MN_SQP1_VR, A2MN_SQP1_VR )
-
- *-----------------2. VR combination ends------------
- end if
-
-
- YVFK = YVFK + TERM_sum
-
-
- IF (( ABS(TERM_sum/YVFK) .LT. RELERR_RES_HV) .or.
- & (TERM_sum==0 .and. YVFK==0) ) THEN
- IF (p .EQ. COUNT+p_step) THEN
- * SINCE THE ODD MODES YIELD ZERO RESULTS WHEN EITHER TEST OR EXPANSION
- * FUNCTION IS AT THE MIDDLE OF THE PARALLEL PLATES. WE NEED TO MAKE SURE
- * THAT THE RELATIVE ERROR BOUND IS SATISFIED AT LEAST 3 TIMES
- * LATER COMMENT: NOTE THAT THIS IS POSSIBLY TAKEN CARE OF BY ONLY SUMMING
- * THE EVEN RESIDUES FOR TEST AND EXPANSION FUNCTIONS SITTING AT THE CENTER
- * OF THE PARALLEL PLATE WAVEGUIDE
- IF(ECOUNT==1) EXIT
- ECOUNT=ECOUNT+1
- ELSE
- ECOUNT=0
- ENDIF
- COUNT=p
- END IF
- * GO BACK TO CALCULATE NEW RESIDUE VALUE FOR EXPONENTIALS
-
-
- p=p+p_step
-
- END DO
- c print *, "p=", p, YVFK
- *-------------------------------END of Exponential term summation-------------
- *
- *
- *------------------------------BEGIN pole loop for ILHI term summation-------------------
- P=1
- COUNT=-1
- ECOUNT=0
-
- CHOICE =2
- DO
-
- CALL FIND_KRHO(KRHO,K,K_Z,P,NUM_KRHO,P_STEP,STORE_KRHO
- & ,D, NUM_GREEN_ARRAY)
-
-
- TERM_sum =0.D0
- if (p>P_max) exit
- if(1) then
- * ----------- 1. VL combination begins----------------
-
- TERM_sum = TERM_sum+
- & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
- & ,z_n,l_n,b_n,k,d,w,epsilon,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_L
- & ,RHO_MIN,RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
- & R_MN_ARRAY_VL, THETA_ARRAY_VL, A1MN_VL, A2MN_VL,
- & AMN_NE1_VL, AMN_NE2_VL, A1MN_SQP1_VL, A2MN_SQP1_VL )
-
-
- end if
-
- *-----------------1. VL combination ends------------
-
- if (1) then
- *-----------------2. VR combination begins----------
-
- TERM_sum = TERM_sum+
- & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
- & ,z_n,l_n2,b_n,k,d,w,epsilon,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_U
- & ,RHO_MIN,RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
- & R_MN_ARRAY_VR, THETA_ARRAY_VR, A1MN_VR, A2MN_VR,
- & AMN_NE1_VR, AMN_NE2_VR, A1MN_SQP1_VR, A2MN_SQP1_VR )
-
- *-----------------2. VR combination ends------------
- end if
-
-
- YVFK = YVFK + TERM_sum
-
-
- IF (( ABS(TERM_sum/YVFK) .LT. RELERR_RES_HV) .or.
- & (TERM_sum==0 .and. YVFK==0) ) THEN
- IF (p .EQ. COUNT+p_step) THEN
- * SINCE THE ODD MODES YIELD ZERO RESULTS WHEN EITHER TEST OR EXPANSION
- * FUNCTION IS AT THE MIDDLE OF THE PARALLEL PLATES. WE NEED TO MAKE SURE
- * THAT THE RELATIVE ERROR BOUND IS SATISFIED AT LEAST 3 TIMES
- * LATER COMMENT: NOTE THAT THIS IS POSSIBLY TAKEN CARE OF BY ONLY SUMMING
- * THE EVEN RESIDUES FOR TEST AND EXPANSION FUNCTIONS SITTING AT THE CENTER
- * OF THE PARALLEL PLATE WAVEGUIDE
- IF(ECOUNT==1) EXIT
- ECOUNT=ECOUNT+1
- ELSE
- ECOUNT=0
- ENDIF
- COUNT=p
- END IF
- * GO BACK TO CALCULATE NEW RESIDUE VALUE FOR EXPONENTIALS
-
-
- p=p+p_step
-
- END DO
- c print *, "p=", p, YVFK
- *-------------------------------END pole loop for ILHI term summation-------------
- *------------------------------BEGIN Brach point calculation-------------------
-
- KRHO = 0
- k_z = k
-
- TERM_sum =0.D0
-
- CHOICE =3
-
- * ----------- 1. VL combination----------------
- if (1) then
-
- TERM_sum = TERM_sum+
- & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
- & ,z_n,l_n,b_n,k,d,w,epsilon,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_L
- & ,RHO_MIN,RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
- & R_MN_ARRAY_VL, THETA_ARRAY_VL, A1MN_VL, A2MN_VL,
- & AMN_NE1_VL, AMN_NE2_VL, A1MN_SQP1_VL, A2MN_SQP1_VL )
-
- *-----------------1. VL combination ends------------
- end if
-
- if (1) then
- *-----------------2. VR combination begins----------
-
- TERM_sum = TERM_sum+
- & YVYK(x_m,y_m,z_m,alpha_mn,l_m,b_m,x_n,y_n
- & ,z_n,l_n2,b_n,k,d,w,epsilon,STORE_KRHO,NUM_KRHO,T_L,T_U,0,E_U
- & ,RHO_MIN,RHO_MAX, curr_flag_m, CHOICE, P, KRHO, K_Z,
- & R_MN_ARRAY_VR, THETA_ARRAY_VR, A1MN_VR, A2MN_VR,
- & AMN_NE1_VR, AMN_NE2_VR, A1MN_SQP1_VR, A2MN_SQP1_VR )
-
- *-----------------2. VR combination ends------------
- end if
-
-
- YVFK = YVFK + TERM_sum
- c print *, TERM_sum
- *-------------------------------END of Branch Point summation
-
-
- *
-
-
- END IF
- * above endif is the judge for exact/approximation
-
- c print *, YVFK
-
-
- END FUNCTION YVFK