C********************************************************************* C Two-Dimensional Time-Domain Multiresolution-Analysis Program C********************************************************************* C C Written by M.Fujii, July 6, 1998 C Completed, July 9, 1998 C Revised for field interpolation for plot, July 15, 1998 C Revised for PML ABC, July 23, 1998 C Revised for field singularity correction, Aug. 6, 1998 C Minor correction, Aug. 14, 1998 C Revised for Impedance Boundary Condition, Aug. 30, 1998 C Revised for thin PEC, Sep. 4, 1998 C Revised for 0 degree edge correction, Sep. 8, 1998 C Revised for Impedance Boundary Condition, Oct. 16, 1998 C (Arbitrary real impedance condition C using inverse scattering for wavelet coefficients) C C Features: C -- TE wave (propagats on xz-plane. C Structures are uniform in y-direction.) C -- Rectangular uniform mesh C -- Homogeneous lossy dielectric material C -- PEC and PMC boundary with 3rd-order Lagrange interpolation C -- PML absorbing boundary condition (only in z-direction)(ref.[1]) C -- Field plot interpolation from scaling function coefficients C -- Singularity field correction for 90 degree corners (ref.[2]) C -- Arbitrary real impedance boundary condition for TEM wave C -- Thin PEC with edge node (ZU-side only) C -- Singularity field correction for 0 degree edges (ref.[2]) C C References: C [1] J.P.Berenger, "A Perfectry Matched Layer for the Absorption C of Electomagnetic Waves," J.of Computational Physics, C vol.114, pp.185-200, 1994. C C [2] G.Mur, "The Modeling of Singularities in the Finite- C difference Approximation of the Time-Domain Electromagnetic- C Field Equations," IEEE MTT, vol.29, pp.1073, 1981. C C********************************************************************** PROGRAM TDMRA2D_HOMOGENEOUS_1_06 C********************************************************************** C PARAMETERS C C IN,KN : MAXIMUM CALCULATION REGION IN X AND Z DIRECTIONS C NB : NUMBER OF BASIS FUNCTIONS (4) C NP : NUMBER OF POINTS FOR TIME SIGNAL DETECTION C NN : MAXIMUM NUMBER OF TIME SERIES DATA C NF : NUMBER OF FIELDS C ((TOTAL,COEFFS,INTERPOLATION) X (EY,HX,HZ) = 9) C NS : NUMBER OF FIELD DISTRIBUTION SNAPSHOTS (MAX 99) C NPEC0 : MAXIMUM NUMBER OF PEC BOUNDARY CONDITIONS C NTPEC0 : MAXIMUM NUMBER OF THIN PEC C NECNR0 : MAX NUMBER OF PERFECT ELECTRIC CONDUCTOR CORNER NODES C NPMC0 : MAXIMUM NUMBER OF PMC BOUNDARY CONDITIONS C NQITP0 : MAX NUMBER OF REGIONS FOR QUADRATIC FILD INTERPOLATION C NEX0 : MAX NUMBER OF EXCITATION POSITIONS C NPMLZ0 : NUMBER OF PML ABSORBERS IN Z-DIRECTION C NPMLX0 : NUMBER OF PML ABSORBERS IN X-DIRECTION C MPML : NUMBER OF PML LAYERS X 2 C NABC0 : MAXIMUM NUMBER OF ABSORBING BOUNDARY CONDITIONS C NIMP0 : MAXIMUM NUMBER OF IMPEDANCE BOUNDARY CONDITIONS C PARAMETER (IN=400, KN=1400) PARAMETER (NB=4) PARAMETER (NP=3) PARAMETER (NN=10000) PARAMETER (NF=9) PARAMETER (NS=1) PARAMETER (NPEC0=10) PARAMETER (NTPEC0=10) PARAMETER (NECNR0=10) PARAMETER (NPMC0=10) PARAMETER (NQITP0=10) PARAMETER (NEX0=10) PARAMETER (NPMLZ0=4) PARAMETER (NPMLX0=2) PARAMETER (MPML=2*8) PARAMETER (NABC0=NPMLZ0+NPMLX0) PARAMETER (NIMP0=10) C C CALCULATION REGION PARAMETERS C DEFINITION : C |----+----+----+----+----+- - ---+----+----+----+----| X(I) C IL ILP ILP2 KLP4 ... IUM4 IUM2 IUM IU C <---------> C YEE'S CELL C INTEGER + IL, ILP, ILP2, ILP4, IU, IUM, IUM2, IUM4, ID, + KL, KLP, KLP2, KLP4, KU, KUM, KUM2, KUM4, KD C C MAIN DIMENSION C V : FIELD COMPONENTS C SJ : CURRENT SORCE FOR HAAR BASIS C SD : ACTUAL CURRENT SORCE DISTRIBUTION C A : WAVELET DECOMPOSITION AND RECONSTRUCTION MATRIX C MASK : FOR ELIMINATING FIELD OUTSIDE THE CALCULATION REGION C REAL + V (0:IN,0:KN,NB), SJ(0:IN,0:KN,NB), + SD(0:IN+1,0:KN+1), A(NB,NB), MASK(0:IN,0:KN) C C UPDATE EQUATION COEFFICIENTS C REAL + C0, CE, CMX, CMZ C C CALCULATION CONDITIONS C REAL + DT, UNIT, CVACUUM, EP0, MU0, ER, MR INTEGER + NMAX, NSTORE, NSTORE2, NPRNT, NMOD, NPD C C OUTPUT DATA C REAL + FLDOUT(0:IN+1,0:KN+1,NF,NS), TMOUT(0:NN,0:NP) C C FOURIER TRANSFORM C REAL + FTOUT(0:NN,0:NP), SPARA(0:NN,NP-1) C C FIELD DISTRIBUTION FILE NAME C CHARACTER*8 + FL(NF,NS) C C PORT POSITIONS C INTEGER + IPORTL(NP), IPORTU(NP), KPORTL(NP), KPORTU(NP) C C PEC BOUNDARY POSITIONS C INTEGER + IPECL(NPEC0), IPECU(NPEC0), KPECL(NPEC0), KPECU(NPEC0) CHARACTER*2 + PECDRC(NPEC0) C C THIN PEC CONDUCTORS C INTEGER + ITPECL(NTPEC0), ITPECU(NTPEC0), + KTPECL(NTPEC0), KTPECU(NTPEC0) CHARACTER*2 + TPECDRC(NTPEC0) INTEGER + TPECEDGL(NTPEC0), TPECEDGU(NTPEC0) REAL + EDGTMP(0:2,-2:2) C C FOR 0-DEGREE EDGE NODE STATIC FIELD CORRECTION C REAL + EUVY01M(NB,NTPEC0), EUVY10(NB,NTPEC0), EUVY01(NB,NTPEC0), + EUVY00(NB,NTPEC0), + EUHX0HM(NB,NTPEC0), EUHZH0(NB,NTPEC0), EUHX0H(NB,NTPEC0), + EUIX0HM(NTPEC0), EUIZH0(NTPEC0), EUIX0H(NTPEC0), + ELVY01M(NB,NTPEC0), ELVY1M0(NB,NTPEC0), ELVY01(NB,NTPEC0), + ELVY00(NB,NTPEC0), + ELHX0HM(NB,NTPEC0), ELHZHM0(NB,NTPEC0), ELHX0H(NB,NTPEC0), + ELIX0HM(NTPEC0), ELIZHM0(NTPEC0), ELIX0H(NTPEC0) C C PEC 90-DEGREE CORNER NODES C INTEGER + IECNR(NECNR0), KECNR(NECNR0) CHARACTER*2 + PECNR(NECNR0) REAL + EYTMP(0:2,0:2), + VY01 (NB,NECNR0), VY10 (NB,NECNR0), + VY01M(NB,NECNR0), VY1M0(NB,NECNR0), VY00 (NB,NECNR0), + HX0H (NB,NECNR0), HZH0 (NB,NECNR0), + HX0HM(NB,NECNR0), HZHM0(NB,NECNR0), + AIX0H (NECNR0), AIZH0 (NECNR0), + AIX0HM(NECNR0), AIZHM0(NECNR0) C C LAGRANGE INTERPOLATION COEFFICIENTS FOR PEC BOUNDARIES C REAL + PECL1, PECL2 C C PMC BOUNDARY POSITIONS C INTEGER + IPMCL(NPMC0), IPMCU(NPMC0), KPMCL(NPMC0), KPMCU(NPMC0) CHARACTER*2 + PMCDRC(NPMC0) C C ABSORBING BOUNDARY POSITIONS C INTEGER + IABCL(NABC0), IABCU(NABC0), KABCL(NABC0), KABCU(NABC0) CHARACTER*1 + ABCDRC(NABC0) CHARACTER*3 + SIDEL(NABC0), SIDEU(NABC0) C C PML-Z (PERFECTLY MATCHED LAYER IN Z-DIRECTION) C REAL + PMLZ(0:IN+1,-1:MPML+1,NB,NPMLZ0), + CM0Z, CM1Z, CM0X (0:MPML+1), CM1X (0:MPML+1), + CE0YX, CE1YX, CE0YZ(0:MPML+1), CE1YZ(0:MPML+1), + SGEZ(0:MPML+1), SGMZ(0:MPML+1) REAL + RA !APPARENT REFRECTION COEFFICIENT AT PML INTERFACE C C IMPEDANCE BOUNDARY POSITIONS AND DIMENSIONS FOR PREVIOUS EY-FIELD C INTEGER + IIMPL(NIMP0), IIMPU(NIMP0), KIMPL(NIMP0), KIMPU(NIMP0) CHARACTER*1 + IMPDRC(NIMP0) REAL + CIX0(NIMP0), CIX1(NIMP0), CIZ0(NIMP0), CIZ1(NIMP0), + CIX0I(NIMP0), CIX1I(NIMP0), CIZ0I(NIMP0), CIZ1I(NIMP0), + ZIMP(NIMP0), ZIMPI(NIMP0), + PRVK3(0:IN,NB,NIMP0), PRVK1(0:IN,NB,NIMP0) C C QUADRATIC INTERPOLATION FOR FIELD PLOT C INTEGER + IQITPL(NQITP0), IQITPU(NQITP0), + KQITPL(NQITP0), KQITPU(NQITP0) C C EXCITATION POSITION C INTEGER + IEXCTL(NEX0), IEXCTU(NEX0), + KEXCTL(NEX0), KEXCTU(NEX0) C C PRINTING FIELD WITH NUMERAL PICTURES C REAL + FLDPRNT(0:IN+1,0:KN+1) C C TEMPORARY WORKING ARRAYS FOR LU BASIS C REAL + B(NB) INTEGER + LL, LU, UL, UU C C LAGRANGE COEFFICIENTS FOR FIELD INTERPOLATION AND EXTRAPOLATION C REAL + FLI1, FLI2, FLI3, + FLI11, FLI12, FLI13, FLI22, FLI23, FLI33, + FLE1, FLE2, FLE3, + FLE11, FLE12, FLE21, FLE22, FLE31, FLE32 C C TEMPORALY DIMENSIONS FOR CHECKING PML C c REAL c + TMPEYX(0:IN+1,-1:MPML+1,NB), TMPEYZ(0:IN+1,-1:MPML+1,NB), c + TMPHX (0:IN+1,-1:MPML+1,NB), TMPHZ (0:IN+1,-1:MPML+1,NB) C C INITIAL DATA C DATA IL,KL /0,0/ DATA PI, CVACUUM /3.14159265, 2.998E8/ DATA EP0, MU0 /8.854E-12, 1.257E-6/ C C 2D HAAR BASIS SAMPLING (RECONSTRUCTING) AND DECOMPOSING MATRIX C DATA A / 0.5, 0.5, 0.5, 0.5, + 0.5,-0.5, 0.5,-0.5, + 0.5, 0.5,-0.5,-0.5, + 0.5,-0.5,-0.5, 0.5 / C DATA LL, LU, UL, UU / 1, 2, 3, 4 / C C SINGULARITY-CORRECTION COEFFICIENTS C DATA C90 /0.839947366/ !FOR 90-DEGREE CORNER DATA C00 /0.707106781/ !FOR 0-DEGREE EDGE C DATA C00 /1.0/ !FOR TEST C C 3RD-ORDER LAGRANGE INTERPOLATION COEFFICIENTS C FOR PEC AND PMC BOUNDARIES C PECL1 = 5.0/8.0 PECL2 = -1.0/16.0 C C 2ND-ORDER LAGRANGE INTERPOLATION COEFFICIENTS FOR PEC CORNER NODES C ECL0 = 3.0/8.0 ECL1 = 3.0/4.0 ECL2 = -1.0/8.0 ECL01 = ECL0*ECL1 ECL02 = ECL0*ECL2 ECL11 = ECL1*ECL1 ECL12 = ECL1*ECL2 ECL22 = ECL2*ECL2 C C LAGRANGE INTERPOLATION COEFFICIENTS FOR QUADRATIC FIELD INTERPOLATION C FLI1 = 15.0/16.0 FLI2 = 5.0/32.0 FLI3 = -3.0/32.0 FLI11 = FLI1*FLI1 FLI12 = FLI1*FLI2 FLI13 = FLI1*FLI3 FLI22 = FLI2*FLI2 FLI23 = FLI2*FLI3 FLI33 = FLI3*FLI3 C C LAGRANGE COEFFICIENTS FOR QUADRATIC FILED EXTRAPOLATION AT BOUNDARIES C FLE1 = 45.0/32.0 FLE2 = -9.0/16.0 FLE3 = 5.0/32.0 FLE11 = FLE1*3.0/4.0 FLE12 = FLE1/4.0 FLE21 = FLE2*3.0/4.0 FLE22 = FLE2/4.0 FLE31 = FLE3*3.0/4.0 FLE32 = FLE3/4.0 C WRITE (6,6000) 6000 FORMAT ( + /' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + /' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + /' !!!! !!!!' + /' !!!! TD-MRA 2D PROGRAM !!!!' + /' !!!! !!!!' + /' !!!! 2D Time-Domain Multiresolution Analysis !!!!' + /' !!!! using Haar Wavelets !!!!' + /' !!!! !!!!' + /' !!!! Written by M.Fujii, July 6, 1998 !!!!' + /' !!!! Revised for field plot, July 15, 1998 !!!!' + /' !!!! Revised for PML ABC , July 22, 1998 !!!!' + /' !!!! Revised for singularity correction, !!!!' + /' !!!! Aug. 6, 1998 !!!!' + /' !!!! !!!!' + /' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + /' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') C C ******************************** C INPUT THE CALCULATION CONDITIONS C ******************************** C C ------------------------------------------------- C (1) SET UNIT IN LENGTH AND OTHER BASIC PARAMETERS C ------------------------------------------------- C UNIT = 1.0E-3 !mm GHZ = 1.0E9*UNIT/CVACUUM Z0 = 120*PI !FOR VARIABLE CONVERSION C C -------------------------- C (2) SET CALCULATION REGION C -------------------------- C ZMAX0 = 100.0 !LENGTH OF THE WAVEGUIDE C C ------------------------------------------------------------------ C (3) SET NUMBERS OF GRIDS (EVEN NUMBERS, 2 X NUMBER OF YEE'S CELLS) C --------------------------============--==========---------------- C ID = 26 KD = 104 C C GRID SPACINGS (!!! 1/2 OF YEE'S CELL !!!) C ^^^^^^^^^^^^^^^^^^^^^^^^^ DX = 22.86/2.0/6.0/2.0 !WR-90 WAVEGUIDE (WIDTH: 22.86mm) DZ = ZMAX0/FLOAT(KD+1) c dz = dz*1.25 !non-square-grid stability test c dz = dx !square grid C C CALCULATION REGION PARAMETERS C ILP = IL+1 ILP2 = IL+2 ILP4 = IL+4 KLP = KL+1 KLP2 = KL+2 KLP4 = KL+4 IU = IL+ID KU = KL+KD IUM = IU-1 IUM2 = IU-2 IUM4 = IU-4 KUM = KU-1 KUM2 = KU-2 KUM4 = KU-4 C C ---------------------------- C (4) SET DIELECTRIC CONSTANTS C ---------------------------- C ER = 1.0 C MR = 1.0 C C ---------------------------------------------------- C (5) SET CONDUCTIVITY C SG SHOULD BE NORMALIZED BY MULTIPLYING (Z0*UNIT) C ---------------------------------------------------- C C SG = 1.0E-2*Z0*UNIT SG = 0.0 C C -------------------------------------------------------- C (6) SET PEC BOUNDARY POSITIONS AND DIRECTION C THE NUMBERS SHOULD SPECIFY ELECTRIC FIELD (EY) NODES C -------------------------------------------------------- C C Example of notation: C C X C | C | KPECL (=KPECU) C ---------------------------------------- C | | C | | C | upper | | IPECU C | side | | C | PEC | (U) | C | | | C | | | IPECL C | | C | | C | lower side PEC | C | ---------------- | IPECL (=IPECU) C | (L) | C | | C | | C | | C ---------------------------------------- - - - ->Z C KPECL KPECU C C NPEC = 2 !NUMBER OF PEC WALL C IPECL(1) = 0 !Xmin FOR MAIN STRUCTURE IPECU(1) = 0 KPECL(1) = 0 KPECU(1) = KD PECDRC(1) = 'XL' C IPECL(2) = 14 !Xmin FOR REFERENCE STRUCTURE IPECU(2) = 14 KPECL(2) = 0 KPECU(2) = 30 PECDRC(2) = 'XL' C IPECL(3) = 0 !PEC_3 IPECU(3) = 4 KPECL(3) = 50 KPECU(3) = 50 PECDRC(3) = 'ZU' C IPECL(4) = 6 !PEC_4 IPECU(4) = 6 KPECL(4) = 52 KPECU(4) = 52 PECDRC(4) = 'XL' C IPECL(5) = 0 !PEC_5 IPECU(5) = 4 KPECL(5) = 54 KPECU(5) = 54 PECDRC(5) = 'ZL' C C ------------------------------------------------------------------- C (7) SET THIN PEC POSITIONS AND DIRECTION C THE NUMBERS SHOULD SPECIFY ELECTRIC FIELD (EY) NODES C NUMBER OF THIN PEC, NTPEC C EDGE SINGULARITY CORRECTION, NSNG0 = 1: ON, 0: OFF. C THIN PEC DIRECTION, TPECDRC = 'ZU' (ONLY 'ZU' IS AVAILABLE NOW) C EDGE CONDITIONS, TPECEDGL, TPECEDGU = C 2: 0 DEGREE EDGE WITH 2ND-ORDER INTERPOLATION C 1: 0 DEGREE EDGE WITH LINEAR AVERAGING C 0: NOT EDGE, CONECTED TO PEC OR PMC SIDE WALLS (DEFAULT) C ------------------------------------------------------------------- C NTPEC = 1 !NUMBER OF THIN PEC C NSNG0 = 1 !EDGE SINGULARITY CORRECTION C ITPECL(1) = 0 ITPECU(1) = 6 c ITPECL(1) = 6 !FOR TEST c ITPECU(1) = 12 !FOR TEST KTPECL(1) = 52 KTPECU(1) = 52 TPECDRC(1)= 'ZU' !ONLY THIS SIDE IS AVAILABLE NOW TPECEDGL(1) = 0 !NOT EDGE, CONNECTED TO A PEC WALL TPECEDGU(1) = 1 !0 DEGREE EDGE WITH 2ND-ORDER INTERPOLATION C C ---------------------------------------------------- C (8) SET PERFECT ELECTRIC CONDUCTOR (PEC) CORNER NODE C THE NUMBERS SHOULD SPECIFY ELECTRIC FIELD NODES C ---------------------------------------------------- C C NODE SPECIFICATION: C C X(I) C | C | Ey NODE C |/ | | C -X-------------X------------X-- C | | | C | | | C | UL | UU | C | o | o | C | | | C ICNR(*)-X-------------X------------X-- C | | | C | o | o | C | LL | LU | C | | | C | | | C -X-------------X------------X------> Z(K) C | | | C KCNR(*) C C NUMBER OF PEC CORNER NODES C NECNR = 0 C C SINGULARITY CORRECTION (NSNG90 = 1: ON, 0: OFF) C NSNG90 = 1 C C POSISIONS C IECNR(1) = 6 KECNR(1) = 50 PECNR(1) = 'LU' C IECNR(2) = 6 KECNR(2) = 54 PECNR(2) = 'LL' C C ---------------------------------------------------------------- C (9) SET PMC BOUNDARY POSITIONS AND DIRECTION C THE NUMBERS SHOULD SPECIFY MAGNETIC FIELD (HX OR HZ) NODES C PMC CAN OVERRIDE PEC BOUNDARY WHEN SET AT THE SAME POSITION C (EY-FIELDS OUTSIDE PMC WALLS ARE ONLY DUE TO PEC ARGORITHM.) C ---------------------------------------------------------------- C NPMC = 2 C IPMCL(1) = 11 !AT X=XMAX IPMCU(1) = 11 KPMCL(1) = 0 KPMCU(1) = KD PMCDRC(1) = 'XU' C IPMCL(2) = 25 !AT X=XMAX FOR REFERENCE STRUCTURE IPMCU(2) = 25 KPMCL(2) = 0 KPMCU(2) = 30 PMCDRC(2) = 'XU' C C -------------------------------------------------------------- C (10) SET IMPEDANCE BOUNDARY POSITIONS AND DIRECTION C THE NUMBERS SHOULD SPECIFY ELECTRIC FIELD (EY) NODES C (VALID ONLY FOR NORMAL INCIDENT PLANE WAVE) C ZIMP > 0 C (ZIMP SHOULD NOT BE ZERO, OTHERWIZE, ZERO DIVIDE OCCURS.) C -------------------------------------------------------------- C NIMP = 0 !NUMBER OF IMPEDANCE WALL C IIMPL(1) = 0 !AT Z=ZMAX IIMPU(1) = ID KIMPL(1) = KD KIMPU(1) = KD ZIMP(1) = SQRT(MR/ER) !ZIMP = 1.0 IMPDRC(1) = 'U' C IIMPL(2) = 0 !AT Z=ZMIN IIMPU(2) = ID KIMPL(2) = 0 KIMPU(2) = 0 ZIMP(2) = SQRT(MR/ER) !ZIMP = 1.0 IMPDRC(2) = 'L' C IIMPL(3) = 0 !AT X=XMIN IIMPU(3) = 0 KIMPL(3) = 0 KIMPU(3) = KD ZIMP(3) = SQRT(MR/ER) !Z = 1.0 IMPDRC(3) = 'L' C IIMPL(4) = ID !AT X=XMAX IIMPU(4) = ID KIMPL(4) = 0 KIMPU(4) = KD ZIMP(4) = SQRT(MR/ER) !Z = 1.0 IMPDRC(4) = 'U' C C ------------------------------------------------------------------- C (11) SET PML ABSORBING BOUNDARY POSITIONS AND DIRECTION C THE NUMBERS SHOULD SPECIFY ELECTRIC FIELD (EY) NODES C (EVEN NUMBERS) C THIS WILL OVERRIDE PEC BOUNDARY WHEN SET AT THE SAME POSITION. C SIDE WALLS CAN BE SPECIFIED BY SETTING "SIDEL" AND "SIDEU" AS C SIDEL(*) = 'PEC' OR 'PMC'. C ------------------------------------------------------------------- C C APPARENT REFRECTION COEFFICIENT AT THE INTERFACE BETWEEN C PML AND NORMAL REGIONS (TYPICALLY [1.0E-5 -- 1.OE-6]) C RA = 1.0E-5 C C NUMBER OF ABSORBING BOUNDARIES C NABC = 4 C IF (NABC.GT.NABC0) STOP 'DIMENSION EXCEEDED! NABC > NABC0' C IABCL(1) = 0 !Z=Zmax IABCU(1) = 12 KABCL(1) = KD KABCU(1) = KD ABCDRC(1) = 'U' SIDEL(1) = 'PEC' SIDEU(1) = 'PMC' C IABCL(2) = 0 !Z=Zmin IABCU(2) = 12 KABCL(2) = 0 KABCU(2) = 0 ABCDRC(2) = 'L' SIDEL(2) = 'PEC' SIDEU(2) = 'PMC' C IABCL(3) = 14 !Z=Zmax FOR REFERENCE STRUCTURE IABCU(3) = 26 KABCL(3) = 30 KABCU(3) = 30 ABCDRC(3) = 'U' SIDEL(3) = 'PEC' SIDEU(3) = 'PMC' C IABCL(4) = 14 !Z=Zmin FOR REFERENCE STRUCTURE IABCU(4) = 26 KABCL(4) = 0 KABCU(4) = 0 ABCDRC(4) = 'L' SIDEL(4) = 'PEC' SIDEU(4) = 'PMC' C C --------------------------------------------------------------- C (12) SET MASK ZERO TO ELIMINATE FIELDS OUTSIDE THE CALCULATION C REGION (OUTSIDE OF PEC,PMC AND ABC) DEFAULT OF MASK = 0.0 C --------------------------------------------------------------- C DO 300 K=KL,KU !ALL REGION DO 300 I=IL,IU MASK(I,K) = 1.0 300 CONTINUE C K = 52 C DO 310 I=IL,IL+4 C MASK(I,K) = 0.0 C310 CONTINUE DO 320 K=KL,KL+KD !BETWEEN MAIN AND REFERENCE STRUCTURES DO 320 I=IL+12,IL+13 MASK(I,K) = 0.0 320 CONTINUE DO 330 K=31,KD !REST OF THE REFERENCE STRUCTURE DO 330 I=14,26 MASK(I,K) = 0.0 330 CONTINUE C C -------------------------------------------------------------------- C (13) SET REGIONS FOR FIELD PLOT USING QUADRATIC INTERPOLATION C NUMBER SHOULD SPECIFY EY NODES RANGING FROM 2 TO MAX-2. C OUTSIDE OF THIS REGION IS AUTOMATICALY BILINEARLY INTERPOLATED. C FIELD VALUES AT BOUNDARIES ARE QUADRATICALY EXTRAPOLATED. C FIELDS ARE STORED IN "??ITRP##" WITH FIELD ?? AND NUMBER ##. C -------------------------------------------------------------------- C NQITP = 0 C IQITPL(1) = 2 IQITPU(1) = ID-2 KQITPL(1) = 2 KQITPU(1) = KD-2 C C ---------------------------------------------------------- C (14) SET OUTPUT PORT POSITIONS (MAXIMUM: NP POINTS (NP=3)) C THE NUMBERS SHOULD SPECIFY ELECTRIC FIELD (EY) NODES C THE THIRD PORT IS TO BE USED AS A REFERENCE PORT WHEN C EXTRACTING S-PARAMETERS. C ---------------------------------------------------------- C IPORTL(1) = 8 IPORTU(1) = 10 KPORTL(1) = 12 KPORTU(1) = 12 C IPORTL(2) = 8 IPORTU(2) = 10 KPORTL(2) = KD-12 KPORTU(2) = KD-12 C IPORTL(3) = 22 !REFERENCE PORT IPORTU(3) = 24 KPORTL(3) = 12 KPORTU(3) = 12 C C ---------------------------------------------------------------- C (15) SET EXCITATION POSITION (ON ACTUAL GRID) C NUMBERS COULD BE ANY INTEGER BETWEEN (0,0) AND (ID+1,KD+1), C BUT NOT AT ABSORBING BOUNDARY CONDITION (ABC) WALL C ---------------------------------------------------------------- C C NUMBER OF EXCITATION POSITIONS C NEX = 2 C IEXCTL(1) = 0 C IEXCTU(1) = 13 !For PEC symmetry IEXCTU(1) = 12 !For PMC symmetry KEXCTL(1) = 2 KEXCTU(1) = 8 C IEXCTL(2) = 14 !FOR REFERENCE C IEXCTU(2) = 27 !For PEC symmetry IEXCTU(2) = 26 !For PMC symmetry KEXCTL(2) = 2 KEXCTU(2) = 8 C IEXCTL(3) = 0 C IEXCTU(3) = 13 !For PEC symmetry IEXCTU(3) = 12 !For PMC symmetry KEXCTL(3) = KD-8 KEXCTU(3) = KD-2 C C SMOOTH EXCITATION DISTRIBUTION (RAISED-COSINE DISTRIBUTION) C DO 500 L=1,NEX DO 500 K=KEXCTL(L),KEXCTU(L) DO 500 I=IEXCTL(L),IEXCTU(L) SD(I,K) = + SIN(PI*FLOAT(K-KEXCTL(L))/FLOAT(KEXCTU(L)-KEXCTL(L)))**2 C + *SIN(PI*FLOAT(I-IEXCTL(L))/FLOAT(IEXCTU(L)-IEXCTL(L)))**2 C For half-sine distribution in x-dirction (PMC symmetry) + *SIN(0.5*PI*FLOAT(I-IEXCTL(L))/FLOAT(IEXCTU(L)-IEXCTL(L)))**2 500 CONTINUE C C EXCITATION FOR EACH BASIS FUNCTION AT EACH EY NODE C DO 520 K=KL,KU,2 DO 520 I=IL,IU,2 B(LL) = SD(I, K ) B(LU) = SD(I, K+1) B(UL) = SD(I+1,K ) B(UU) = SD(I+1,K+1) DO 540 L=1,NB DO 560 M=1,NB SJ(I,K,L) = SJ(I,K,L)+A(M,L)*B(M) 560 CONTINUE 540 CONTINUE 520 CONTINUE C C --------------------------------------------------- C (16) SET TIME INTERVAL C DTMRGN : MARGIN FOR DT TO MAXIMUM LIMIT OF DT C DT : NORMALIZED TIME DISCRETIZATION C DTREAL : UNNORMALIZED VALUE OF DT IN SECOND C --------------------------------------------------- C DTMRGN = 0.9 C DD = 1.0/SQRT(1.0/(2.0*DX)**2+1.0/(2.0*DZ)**2) DT = DTMRGN*DD/SQRT(ER) DTREAL = DT*UNIT/CVACUUM C C --------------------------------------------------------------------- C (17) SET EXCITING PULSE DURATION (RAISED-COSINE-MODULATED SINE WAVE) C FTGT : TARGET FREQUENCY C NMOD : NUMBER OF HALF SINE WAVE WITHIN EXCITATION DURATION C TT : PULSE DURATION TIME C NPD : NUMBER OF TIME STEPS FOR PULSE DURATION C --------------------------------------------------------------------- C FTGT = 10.0*GHZ NMOD = 2 C TT = 1.0/FTGT NPD = INT(TT/2.0/DT)*NMOD C C ---------------------------------------------------------- C (18) SET TIME STEP C NMAX : NUMBER OF TIME STEPS ( NMAX/NSTORE2 < NN ) C NSTORE2 : TIME SERIES DATA STORE INTERVAL C NSTORE : FIELD DISTRIBUTION DATA STORE INTERVAL C NPRNT : FIELD DATA PRINT INTERVAL C FNQ : NYQUIST FREQUENCY (FOR REFERENCE) C ---------------------------------------------------------- C NMAX = INT(4.42E-10/DTREAL)*100 C NSTORE2 = (NMAX-1)/NN+1 NSTORE = (NMAX-1)/NS+1 NPRNT = (NMAX-1)/20+1 FNQ = 1.0/(2.0*DT*FLOAT(NSTORE2)) IF (NMAX/NSTORE2.GT.NN) + STOP 'DIMENSION EXCEEDED! NMAX/NSTORE2 > NN' C C ----------------------------------- C (19) FOURIER TRANSFORM CONDITIONS C NDFT = 1 (FFT ON), or 0 (OFF) C FCTR: CENTER FREQUENCY (Hz) C FSPN: FREQ. SPAN (Hz) C FRES: FREQ. RESOLUTION (Hz) C ----------------------------------- C NDFT = 1 C FCTR = 11E9 FSPN = 18E9 FRES = 5E7 C C ---------------------------- C PRINT CALCULATION CONDITIONS C ---------------------------- C WRITE (6,6010) ID/2,KD/2,2.0*DX,2.0*DZ,UNIT WRITE (6,6020) DT,DTREAL,DTMRGN,FTGT/GHZ,FTGT,FNQ, + TT,NMOD,NPD, + NMAX,NSTORE,NPRNT,NSTORE2 C 6010 FORMAT (/'[2D TD-MRA CALCULATION CONDITONS]' + /,' A.MODEL GEOMETRIES' + /,' YEE''S GRID COUNTS FOR' + /,' X,Z-DIRECTION = ',2I6 + /,' YEE''S GRID DIMENSIONS FOR' + /,' X-DIRECTION = ',G15.5 + /,' Z-DIRECTION = ',G15.5 + /,' LENGTH UNIT IN M = ',G15.5) 6020 FORMAT (' B.TIME SPACE CONDITIONS' + /,' TIME INTERVAL, DT = ',G15.5 + /,' TIME INTERVAL /SEC, DTREAL = ',G15.5 + /,' TIME STEP MARGIN, DTMRGN = ',G15.5 + /,' TERGET FREQUENCY, FTGT /GHZ = ',G15.5 + /,' TERGET FREQUENCY, FTGT = ',G15.5 + /,' NYQIST FREQUENCY, FNQ = ',G15.5 + /,' PULSE DURATION TIME ,TT = ',G15.5 + /,' NUMBER OF HALF-SINE WAVE IN PULSE = ',I8 + /,' PULSE DURATION /IT, NPD = ',I8 + /,' MAXIMUM NUMBER OF TIME STEPS, NMAX = ',I8 + /,' FIELD DATA STORE INTERVAL, NSTORE = ',I8 + /,' FIELD DATA PRINT INTERVAL, NPRNT = ',I8 + /,' TIME DATA STORE INTERVAL, NSTORE2 = ',I8) C C ************* C PREPROCESSING C ************* C ---------------------------- C UPDATE-EQUATION COEFFICIENTS C ---------------------------- C C0 = (2.0*ER-SG*DT)/(2.0*ER+SG*DT) CE = 2.0*DT/(2.0*ER+SG*DT)/(4.0*DX*DZ) CMX = -DT*DX/(MR*DZ) CMZ = -DT*DZ/(MR*DX) C C ---------------- C PML COEFFICIENTS C ---------------- C C RERATIVE PERMEABILITY (TEMPORARY CONSTANT) C AMU = MR C C MAXIMUM ELECTRIC LOSS FACTOR C SGMAX = -3.0*ER*EP0*CVACUUM*LOG(RA)/(2.0*DZ*UNIT*FLOAT(MPML/2)) C C VARIABLE CONVERSION C SGMAX = SGMAX*Z0*UNIT C C LOSS FACTOR AT NODAL POINTS IN UNIT OF YEE'S GRID INTERVAL C (TAKE AVARAGE VALUE OVER YEE'S GRID INTERVAL, EQ(40) OF REF[1]) C RMAX = 0.5*FLOAT(MPML) DO 700 K=0,MPML+1 R = 0.5*FLOAT(K) SGEZ(K) = SGMAX*(3.0*R*R+0.25)/(3.0*RMAX*RMAX) SGMZ(K) = AMU*SGEZ(K)/ER 700 CONTINUE SGEZ(0) = 0.5*SGEZ(0) !SGEZ(K) IS A QUADRATIC FUNCTION. SGMZ(0) = 0.5*SGMZ(0) C C PML UPDATE EQUATION COEFFICIENTS C CM0Z = 1.0 CM1Z = -DT/AMU*DZ/DX DO 800 K=1,MPML-1,2 CM0X(K) = EXP(-SGMZ(K)/AMU*DT) CM1X(K) = (1.0-EXP(-SGMZ(K)/AMU*DT))/SGMZ(K)*DX/DZ 800 CONTINUE DO 810 K=0,MPML,2 CE0YZ(K) = EXP(-SGEZ(K)/ER*DT) CE1YZ(K) = (1.0-EXP(-SGEZ(K)/ER*DT))/SGEZ(K)/(DZ*DX*4.0) 810 CONTINUE CE0YX = 1.0 CE1YX = -DT/ER/(DX*DZ*4.0) C C ----------------------------------------------- C IMPEDANCE BOUNDARY UPDATE EQUATION COEFFICIENTS C ----------------------------------------------- C ZSPC = SQRT(MR/ER) !SPACE IMPEDANCE DO 850 M=1,NIMP IF (IIMPL(M).EQ.IIMPU(M)) THEN !PERPENDICULAR TO X-AXIS ANM = 4.0*MR*2.0*DX-3.0*ZIMP(M)*DT DNM = 4.0*MR*2.0*DX+3.0*ZIMP(M)*DT CIZ0(M) = ANM/DNM CIZ1(M) = DT/DNM ZIMPI(M) = ZSPC*ZSPC/ZIMP(M) !INVERSE SCATTERING IMPEDANCE ANM = 4.0*MR*2.0*DX-3.0*ZIMPI(M)*DT DNM = 4.0*MR*2.0*DX+3.0*ZIMPI(M)*DT CIZ0I(M) = ANM/DNM !INVERSE SCATTERING COEFF. CIZ1I(M) = DT/DNM !INVERSE SCATTERING COEFF. ELSEIF (KIMPL(M).EQ.KIMPU(M)) THEN !PERPENDICULAR TO Z-AXIS ANM = 4.0*MR*2.0*DZ-3.0*ZIMP(M)*DT DNM = 4.0*MR*2.0*DZ+3.0*ZIMP(M)*DT CIX0(M) = ANM/DNM CIX1(M) = DT/DNM ZIMPI(M) = ZSPC*ZSPC/ZIMP(M) !INVERS SCATERING IMPEDANCE ANM = 4.0*MR*2.0*DZ-3.0*ZIMPI(M)*DT DNM = 4.0*MR*2.0*DZ+3.0*ZIMPI(M)*DT CIX0I(M) = ANM/DNM !INVERSE SCATTERING COEFF. CIX1I(M) = DT/DNM !INVERSE SCATTERING COEFF. ENDIF 850 CONTINUE C C ***************************** C ***************************** C START TIME MARCHING ALGORITHM C ***************************** C ***************************** C PRINT *,'SOLVING UPDATE EQUATIONS ...' C DO 10000 IT=1,NMAX C IF (IT.LT.NPD) THEN ST = SIN(PI*FLOAT(IT-1)/FLOAT(NPD))**2 + *SIN(FLOAT(NMOD)*PI*FLOAT(IT-1)/FLOAT(NPD)) ELSE ST = 0.0 END IF C C ************************************************* C SOLVE UPDATE EQUATIONS FOR MAXWELL CURL EQUATIONS C ************************************************* C C -------------------------------- C EY NODE UPDATE (INTERNAL POINTS) C -------------------------------- C DO 1000 L=1,NB DO 1000 K=KLP2,KUM2,2 DO 1000 I=ILP2,IUM2,2 V(I,K,L) = C0*V(I,K,L) + +CE*(V(I,K+1,L)-V(I,K-1,L)-V(I+1,K,L)+V(I-1,K,L) + -ST*SJ(I,K,L)) 1000 CONTINUE C C ----------------------------------------------------------- C SET MASK OUTSIDE THE ACTUAL CALCULATION REGION FOR EY FIELD C ----------------------------------------------------------- C C ACTUAL CALCULATION REGION : MASK = 1.0 C OUTSIDE OF THE REGION : MASK = 0.0 C DO 1100 L=1,NB DO 1100 K=KL,KU,2 DO 1100 I=IL,IU,2 V(I,K,L) = MASK(I,K)*V(I,K,L) 1100 CONTINUE C C ------------ C PEC BOUNDARY C ------------ C C APPLIED ONLY TO EY NODE C DO 1300 L=1,NPEC C C PEC PERPENDICULAR TO X-AXIS C C LOWER SIDE IF (PECDRC(L).EQ.'XL') THEN I = IPECL(L) IP2 = I+2 IP4 = I+4 DO 1320 K=KPECL(L),KPECU(L),2 V(I,K,1) = 0.5*(PECL1*(V(IP2,K,1)+V(IP2,K,3)) + +PECL2*(V(IP4,K,1)+V(IP4,K,3))) V(I,K,3) = -V(I,K,1) V(I,K,2) = 0.5*(PECL1*(V(IP2,K,2)+V(IP2,K,4)) + +PECL2*(V(IP4,K,2)+V(IP4,K,4))) V(I,K,4) = -V(I,K,2) 1320 CONTINUE C UPPER SIDE ELSEIF (PECDRC(L).EQ.'XU') THEN I = IPECL(L) IM2 = I-2 IM4 = I-4 DO 1340 K=KPECL(L),KPECU(L),2 V(I,K,1) = 0.5*(PECL1*(V(IM2,K,1)-V(IM2,K,3)) + +PECL2*(V(IM4,K,1)-V(IM4,K,3))) V(I,K,3) = V(I,K,1) V(I,K,2) = 0.5*(PECL1*(V(IM2,K,2)-V(IM2,K,4)) + +PECL2*(V(IM4,K,2)-V(IM4,K,4))) V(I,K,4) = V(I,K,2) 1340 CONTINUE ENDIF C C PEC PERPENDICULAR TO Z-AXIS C C LOWER SIDE IF (PECDRC(L).EQ.'ZL') THEN K = KPECL(L) KP2 = K+2 KP4 = K+4 DO 1350 I=IPECL(L),IPECU(L),2 V(I,K,1) = 0.5*(PECL1*(V(I,KP2,1)+V(I,KP2,2)) + +PECL2*(V(I,KP4,1)+V(I,KP4,2))) V(I,K,2) = -V(I,K,1) V(I,K,3) = 0.5*(PECL1*(V(I,KP2,3)+V(I,KP2,4)) + +PECL2*(V(I,KP4,3)+V(I,KP4,4))) V(I,K,4) = -V(I,K,3) 1350 CONTINUE C UPPER SIDE ELSEIF (PECDRC(L).EQ.'ZU') THEN K = KPECL(L) KM2 = K-2 KM4 = K-4 DO 1360 I=IPECL(L),IPECU(L),2 V(I,K,1) = 0.5*(PECL1*(V(I,KM2,1)-V(I,KM2,2)) + +PECL2*(V(I,KM4,1)-V(I,KM4,2))) V(I,K,2) = V(I,K,1) V(I,K,3) = 0.5*(PECL1*(V(I,KM2,3)-V(I,KM2,4)) + +PECL2*(V(I,KM4,3)-V(I,KM4,4))) V(I,K,4) = V(I,K,3) 1360 CONTINUE ENDIF C 1300 CONTINUE C C ------------------ C THIN PEC CONDUCTER C ------------------ C DO 11000 M = 1,NTPEC C IF (TPECDRC(M).EQ.'ZU') THEN C C AT K = KTPECL(M), CONVENTIONAL PEC CONDITION C K = KTPECL(M) KM2 = K-2 KM4 = K-4 DO 11010 I=ITPECL(M),ITPECU(M),2 V(I,K,1) = 0.5*(PECL1*(V(I,KM2,1)-V(I,KM2,2)) + +PECL2*(V(I,KM4,1)-V(I,KM4,2))) V(I,K,2) = V(I,K,1) V(I,K,3) = 0.5*(PECL1*(V(I,KM2,3)-V(I,KM2,4)) + +PECL2*(V(I,KM4,3)-V(I,KM4,4))) V(I,K,4) = V(I,K,3) 11010 CONTINUE C C AT K = KTPECL(M)+2, SELF-CONSISTENT PEC CONDITION C K = KTPECL(M)+2 KP2 = K+2 DO 11020 I=ITPECL(M),ITPECU(M),2 EKU = V(I,K,1) -V(I,K,2) EKP2U = V(I,KP2,1)-V(I,KP2,2) V(I,K,1) = 0.5*((PECL1+1)*EKU+PECL2*EKP2U) V(I,K,2) = V(I,K,1)-EKU EKU = V(I,K,3) -V(I,K,4) EKP2U = V(I,KP2,3)-V(I,KP2,4) V(I,K,3) = 0.5*((PECL1+1)*EKU+PECL2*EKP2U) V(I,K,4) = V(I,K,3)-EKU 11020 CONTINUE C C 0 DEGREE EDGE NODE C DO 11030 K=-2,2 DO 11030 I=0,2 EDGTMP(I,K) = 0.0 11030 CONTINUE C C LOWER EDGE C I = ITPECL(M) C IF (TPECEDGL(M).EQ.1.OR.TPECEDGL(M).EQ.2) THEN C C COMPUTE SURROUNDING UU-SUBNODES C JDIM = TPECEDGL(M) DO 11040 K=-JDIM,JDIM DO 11040 I=0,JDIM DO 11040 IB=1,NB ITMP = ITPECL(M)-2*I KTMP = KTPECL(M)+2*K EDGTMP(I,K) = EDGTMP(I,K)+A(IB,UU)*V(ITMP,KTMP,IB) 11040 CONTINUE C C AT THE CONVENTIONAL CORNER NODE C C INTERPOLATION C IF (TPECEDGL(M).EQ.1) THEN C SIMPLE AVERAGE B(UU) = 0.0 B(UL) = 0.5*EDGTMP(0,-1) B(LU) = 0.5*EDGTMP(1,0) B(LL) = 0.25*(EDGTMP(0,-1)+EDGTMP(1,-1)+EDGTMP(1,0)) ELSEIF (TPECEDGL(M).EQ.2) THEN C 2ND-ORDER INTERPOLATION B(UU) = 0.0 B(UL) = ECL1*EDGTMP(0,-1)+ECL2*EDGTMP(0,-2) B(LU) = ECL1*EDGTMP(1,0) +ECL2*EDGTMP(2,0) B(LL) = ECL0*(B(UL)+B(LU)) + +ECL11*EDGTMP(1,-1) + +ECL12*(EDGTMP(1,-2)+EDGTMP(2,-1)) + +ECL22*EDGTMP(2,-2) ENDIF C C RE-DECOMPOSITION INTO HAAR BASIS COEFFICIENTS C I = ITPECL(M) K = KTPECL(M) DO 11050 IB=1,NB VTMP = 0.0 DO 11060 L=1,NB VTMP = VTMP+A(L,IB)*B(L) 11060 CONTINUE V(I,K,IB) = VTMP 11050 CONTINUE C C AT THE SELF-CONSISTENT CORNER NODE C C RECONSTRUCTION OF FIELD AT LU-SUBNODES C I = ITPECL(M) K = KTPECL(M)+2 EYLU = 0.0 DO 11070 L=1,NB EYLU = EYLU+A(L,LU)*V(I,K,L) 11070 CONTINUE C C INTERPOLATION C IF (TPECEDGL(M).EQ.1) THEN C SIMPLE AVERAGE B(UU) = EDGTMP(0,1) B(LU) = EYLU B(UL) = 0.5*B(UU) B(LL) = 0.25*(EDGTMP(1,0)+EDGTMP(1,1)+EDGTMP(0,1)) ELSEIF (TPECEDGL(M).EQ.2) THEN C 2ND-ORDER INTERPOLATION B(UU) = EDGTMP(0,1) B(LU) = EYLU B(UL) = ECL1*EDGTMP(0,1)+ECL2*EDGTMP(0,2) B(LL) = ECL01*(EDGTMP(1,0)+EDGTMP(0,1)) + +ECL02*(EDGTMP(2,0)+EDGTMP(0,2)) + +ECL11* EDGTMP(1,1) + +ECL12*(EDGTMP(1,2)+EDGTMP(2,1)) + +ECL22*EDGTMP(2,2) ENDIF C C RE-DECOMPOSITION INTO HAAR BASIS COEFFICIENTS C I = ITPECL(M) K = KTPECL(M)+2 DO 11080 IB=1,NB VTMP = 0.0 DO 11090 L=1,NB VTMP = VTMP+A(L,IB)*B(L) 11090 CONTINUE V(I,K,IB) = VTMP 11080 CONTINUE C ENDIF C C UPPER EDGE C IF (TPECEDGU(M).EQ.1.OR.TPECEDGU(M).EQ.2) THEN C C COMPUTE SURROUNDING LU-SUBNODES C JDIM = TPECEDGU(M) DO 11140 K=-JDIM,JDIM DO 11140 I=0,JDIM DO 11140 IB=1,NB ITMP = ITPECU(M)+2*I KTMP = KTPECL(M)+2*K EDGTMP(I,K) = EDGTMP(I,K)+A(IB,LU)*V(ITMP,KTMP,IB) 11140 CONTINUE C C AT THE CONVENTIONAL CORNER NODE C C INTERPOLATION C IF (TPECEDGU(M).EQ.1) THEN C SIMPLE AVERAGE B(LU) = 0.0 B(LL) = 0.5*EDGTMP(0,-1) B(UU) = 0.5*EDGTMP(1,0) B(UL) = 0.25*(EDGTMP(0,-1)+EDGTMP(1,-1)+EDGTMP(1,0)) ELSEIF (TPECEDGU(M).EQ.2) THEN C 2ND-ORDER INTERPOLATION B(LU) = 0.0 B(LL) = ECL1*EDGTMP(0,-1)+ECL2*EDGTMP(0,-2) B(UU) = ECL1*EDGTMP(1,0) +ECL2*EDGTMP(2,0) B(UL) = ECL0*(B(LL)+B(UU)) + +ECL11*EDGTMP(1,-1) + +ECL12*(EDGTMP(1,-2)+EDGTMP(2,-1)) + +ECL22*EDGTMP(2,-2) ENDIF C C RE-DECOMPOSITION INTO HAAR BASIS COEFFICIENTS C I = ITPECU(M) K = KTPECL(M) DO 11150 IB=1,NB VTMP = 0.0 DO 11160 L=1,NB VTMP = VTMP+A(L,IB)*B(L) 11160 CONTINUE V(I,K,IB) = VTMP 11150 CONTINUE C C AT THE SELF-CONSISTENT CORNER NODE C C RECONSTRUCTION OF FIELD AT UU-SUBNODES C I = ITPECU(M) K = KTPECL(M)+2 EYUU = 0.0 DO 11170 L=1,NB EYUU = EYUU+A(L,UU)*V(I,K,L) 11170 CONTINUE C C INTERPOLATION C IF (TPECEDGU(M).EQ.1) THEN C SIMPLE AVERAGE B(LU) = EDGTMP(0,1) B(UU) = EYUU B(LL) = 0.5*B(LU) B(UL) = 0.25*(EDGTMP(1,0)+EDGTMP(1,1)+EDGTMP(0,1)) ELSEIF (TPECEDGU(M).EQ.2) THEN C 2ND-ORDER INTERPOLATION B(LU) = EDGTMP(0,1) B(UU) = EYUU B(LL) = ECL1*EDGTMP(0,1)+ECL2*EDGTMP(0,2) B(UL) = ECL01*(EDGTMP(1,0)+EDGTMP(0,1)) + +ECL02*(EDGTMP(2,0)+EDGTMP(0,2)) + +ECL11* EDGTMP(1,1) + +ECL12*(EDGTMP(1,2)+EDGTMP(2,1)) + +ECL22*EDGTMP(2,2) ENDIF C C RE-DECOMPOSITION INTO HAAR BASIS COEFFICIENTS C I = ITPECU(M) K = KTPECL(M)+2 DO 11180 IB=1,NB VTMP = 0.0 DO 11190 L=1,NB VTMP = VTMP+A(L,IB)*B(L) 11190 CONTINUE V(I,K,IB) = VTMP 11180 CONTINUE C ENDIF C ELSEIF (TPECDRC(M).EQ.'ZL') THEN C --- TO BE IMPLEMENTED --- ELSEIF (TPECDRC(M).EQ.'XU') THEN C --- TO BE IMPLEMENTED --- ELSEIF (TPECDRC(M).EQ.'XL') THEN C --- TO BE IMPLEMENTED --- ENDIF C 11000 CONTINUE C C ---------------------------------- C PEC CORNER NODE (90 DEGREE CORNER) C ---------------------------------- C DO 11300 L=1,NECNR C C COMPUTE GIVEN NODAL VALUES C (RESEARCH NOTE VOL.R4 PP.13 EQ(64)-(66)) C ISGN = 1 KSGN = 1 IF (PECNR(L)(1:1).EQ.'U') ISGN = -1 IF (PECNR(L)(2:2).EQ.'U') KSGN = -1 C IF (PECNR(L).EQ.'LL') THEN NLU = 1 ELSEIF (PECNR(L).EQ.'LU') THEN NLU = 2 ELSEIF (PECNR(L).EQ.'UL') THEN NLU = 3 ELSEIF (PECNR(L).EQ.'UU') THEN NLU = 4 ENDIF C DO 11490 K=0,2 DO 11490 I=0,2 EYTMP(I,K) = 0.0 11490 CONTINUE C DO 11500 I=0,2 DO 11500 K=0,2 DO 11500 IB=1,NB ITMP = IECNR(L)+ISGN*2*I KTMP = KECNR(L)+KSGN*2*K EYTMP(I,K) = EYTMP(I,K)+A(IB,NLU)*V(ITMP,KTMP,IB) 11500 CONTINUE C C COMPUTE INTERPOLATION (ASSUME EYTMP(0,0)=0) C EX = ECL1 *EYTMP(1,0)+ECL2 *EYTMP(2,0) EZ = ECL1 *EYTMP(0,1)+ECL2 *EYTMP(0,2) EXZ = ECL0 *(EX+EZ) + +ECL11*EYTMP(1,1)+ECL12*EYTMP(1,2) + +ECL12*EYTMP(2,1)+ECL22*EYTMP(2,2) C C CONVERT FIELDS INTO HAAR BASIS COEFFICIENTS C IF (PECNR(L).EQ.'LL') THEN B(LL) = 0.0 B(LU) = EZ B(UL) = EX B(UU) = EXZ ELSEIF (PECNR(L).EQ.'LU') THEN B(LL) = EZ B(LU) = 0.0 B(UL) = EXZ B(UU) = EX ELSEIF (PECNR(L).EQ.'UL') THEN B(LL) = EX B(LU) = EXZ B(UL) = 0.0 B(UU) = EZ ELSEIF (PECNR(L).EQ.'UU') THEN B(LL) = EXZ B(LU) = EX B(UL) = EZ B(UU) = 0.0 ENDIF C DO 11540 IB=1,NB VTMP = 0.0 DO 11560 M=1,NB VTMP = VTMP+A(M,IB)*B(M) 11560 CONTINUE V(IL+IECNR(L),KL+KECNR(L),IB) = VTMP 11540 CONTINUE C 11300 CONTINUE C C ------- C PML ABC C ------- C DO 1400 L=1,NABC C IF (KABCL(L).EQ.KABCU(L)) THEN C C Z-DIRECTION (ABC INTERFACE PERPENDICULAR TO Z-AXIS) C C TRANSFER MAGNETIC FIELD VALUES AT NODES NEXT TO THE INTERFACE C INTO PML DIMENSION "PMLZ" C C LOWER SIDE CASE IF (ABCDRC(L).EQ.'L') THEN DO 1410 IB=1,NB DO 1410 I=IL+IABCL(L),IL+IABCU(L),2 PMLZ(I,-1,IB,L) = -V(I,KABCL(L)+1,IB) 1410 CONTINUE C UPPER SIDE CASE ELSEIF (ABCDRC(L).EQ.'U') THEN DO 1420 IB=1,NB DO 1420 I=IL+IABCL(L),IL+IABCU(L),2 PMLZ(I,-1,IB,L) = V(I,KABCL(L)-1,IB) 1420 CONTINUE ENDIF C C EXPONENTIAL TIME STEPPING EQUATIONS FOR INNER EYX AND EYZ NODES C DO 1440 IB=1,NB DO 1440 K=0,MPML-2,2 DO 1440 I=IL+IABCL(L)+2,IL+IABCU(L)-2,2 C FOR Eyx PMLZ(I,K,IB,L) = PMLZ(I,K,IB,L) + +CE1YX*(PMLZ(I+1,K,IB,L)-PMLZ(I-1,K,IB,L)) C FOR Eyz PMLZ(I+1,K+1,IB,L) = CE0YZ(K)*PMLZ(I+1,K+1,IB,L) + +CE1YZ(K)*(PMLZ(I,K+1,IB,L)-PMLZ(I,K-1,IB,L)) 1440 CONTINUE C C PEC WALLS INSIDE PML REGIONS C C UPPER X SIDE C IF (SIDEU(L).EQ.'PEC') THEN I = IABCU(L) IM2 = I-2 IM4 = I-4 DO 1450 K=0,MPML,2 C FOR Eyx PMLZ(I,K,1,L) = 0.5* + (PECL1*(PMLZ(IM2,K,1,L)-PMLZ(IM2,K,3,L)) + +PECL2*(PMLZ(IM4,K,1,L)-PMLZ(IM4,K,3,L))) PMLZ(I,K,3,L) = PMLZ(I,K,1,L) PMLZ(I,K,2,L) = 0.5* + (PECL1*(PMLZ(IM2,K,2,L)-PMLZ(IM2,K,4,L)) + +PECL2*(PMLZ(IM4,K,2,L)-PMLZ(IM4,K,4,L))) PMLZ(I,K,4,L) = PMLZ(I,K,2,L) C FOR Eyz PMLZ(I+1,K+1,1,L) = 0.5* + (PECL1*(PMLZ(IM2+1,K+1,1,L)-PMLZ(IM2+1,K+1,3,L)) + +PECL2*(PMLZ(IM4+1,K+1,1,L)-PMLZ(IM4+1,K+1,3,L))) PMLZ(I+1,K+1,3,L) = PMLZ(I+1,K+1,1,L) PMLZ(I+1,K+1,2,L) = 0.5* + (PECL1*(PMLZ(IM2+1,K+1,2,L)-PMLZ(IM2+1,K+1,4,L)) + +PECL2*(PMLZ(IM4+1,K+1,2,L)-PMLZ(IM4+1,K+1,4,L))) PMLZ(I+1,K+1,4,L) = PMLZ(I+1,K+1,2,L) 1450 CONTINUE ENDIF C C LOWER X SIDE C IF (SIDEL(L).EQ.'PEC') THEN I = IABCL(L) IP2 = I+2 IP4 = I+4 DO 1460 K=0,MPML,2 C FOR Eyx PMLZ(I,K,1,L) = 0.5* + (PECL1*(PMLZ(IP2,K,1,L)+PMLZ(IP2,K,3,L)) + +PECL2*(PMLZ(IP4,K,1,L)+PMLZ(IP4,K,3,L))) PMLZ(I,K,3,L) = -PMLZ(I,K,1,L) PMLZ(I,K,2,L) = 0.5* + (PECL1*(PMLZ(IP2,K,2,L)+PMLZ(IP2,K,4,L)) + +PECL2*(PMLZ(IP4,K,2,L)+PMLZ(IP4,K,4,L))) PMLZ(I,K,4,L) = -PMLZ(I,K,2,L) C FOR Eyz PMLZ(I+1,K+1,1,L) = 0.5* + (PECL1*(PMLZ(IP2+1,K+1,1,L)+PMLZ(IP2+1,K+1,3,L)) + +PECL2*(PMLZ(IP4+1,K+1,1,L)+PMLZ(IP4+1,K+1,3,L))) PMLZ(I+1,K+1,3,L) = -PMLZ(I+1,K+1,1,L) PMLZ(I+1,K+1,2,L) = 0.5* + (PECL1*(PMLZ(IP2+1,K+1,2,L)+PMLZ(IP2+1,K+1,4,L)) + +PECL2*(PMLZ(IP4+1,K+1,2,L)+PMLZ(IP4+1,K+1,4,L))) PMLZ(I+1,K+1,4,L) = -PMLZ(I+1,K+1,2,L) 1460 CONTINUE ENDIF C C Z BACK SIDE (FIXED PEC BOUNDARY) C K = MPML KM2 = K-2 KM4 = K-4 DO 1470 I=IABCL(L),IABCU(L),2 C FOR Eyx PMLZ(I,K,1,L) = 0.5* + (PECL1*(PMLZ(I,KM2,1,L)-PMLZ(I,KM2,2,L)) + +PECL2*(PMLZ(I,KM4,1,L)-PMLZ(I,KM4,2,L))) PMLZ(I,K,2,L) = PMLZ(I,K,1,L) PMLZ(I,K,3,L) = 0.5* + (PECL1*(PMLZ(I,KM2,3,L)-PMLZ(I,KM2,4,L)) + +PECL2*(PMLZ(I,KM4,3,L)-PMLZ(I,KM4,4,L))) PMLZ(I,K,4,L) = PMLZ(I,K,3,L) C FOR Eyz PMLZ(I+1,K+1,1,L) = 0.5* + (PECL1*(PMLZ(I+1,KM2+1,1,L)-PMLZ(I+1,KM2+1,2,L)) + +PECL2*(PMLZ(I+1,KM4+1,1,L)-PMLZ(I+1,KM4+1,2,L))) PMLZ(I+1,K+1,2,L) = PMLZ(I+1,K+1,1,L) PMLZ(I+1,K+1,3,L) = 0.5* + (PECL1*(PMLZ(I+1,KM2+1,3,L)-PMLZ(I+1,KM2+1,4,L)) + +PECL2*(PMLZ(I+1,KM4+1,3,L)-PMLZ(I+1,KM4+1,4,L))) PMLZ(I+1,K+1,4,L) = PMLZ(I+1,K+1,3,L) 1470 CONTINUE C C EXPONENTIAL TIME STEPPING EQUATION FOR HZ NODES IN PML REGION C DO 1480 IB=1,NB DO 1480 K=0,MPML-2,2 DO 1480 I=IL+IABCL(L)+1,IL+IABCU(L)-1,2 PMLZ(I,K,IB,L) = PMLZ(I,K,IB,L) + +CM1Z*(PMLZ(I+2,K+1,IB,L)+PMLZ(I+1,K,IB,L) + -PMLZ(I ,K+1,IB,L)-PMLZ(I-1,K,IB,L)) 1480 CONTINUE C C EXPONENTIAL TIME STEPPING EQUATION FOR HX NODES IN PML REGION C DO 1490 IB=1,NB DO 1490 K=1,MPML-1,2 DO 1490 I=IL+IABCL(L)+2,IL+IABCU(L)-2,2 PMLZ(I,K,IB,L) = + CM0X(K)* PMLZ(I,K,IB,L) + +CM1X(K)*(PMLZ(I+1,K+2,IB,L)+PMLZ(I,K+1,IB,L) + -PMLZ(I+1,K ,IB,L)-PMLZ(I,K-1,IB,L)) 1490 CONTINUE C C PMC WALLS INSIDE PML REGIONS C C UPPER X SIDE C IF (SIDEU(L).EQ.'PMC') THEN I = IABCU(L)-1 IM2 = I-2 IM4 = I-4 DO 1455 K=0,MPML,2 C FOR Hz PMLZ(I,K,1,L) = 0.5* + (PECL1*(PMLZ(IM2,K,1,L)-PMLZ(IM2,K,3,L)) + +PECL2*(PMLZ(IM4,K,1,L)-PMLZ(IM4,K,3,L))) PMLZ(I,K,3,L) = PMLZ(I,K,1,L) PMLZ(I,K,2,L) = 0.5* + (PECL1*(PMLZ(IM2,K,2,L)-PMLZ(IM2,K,4,L)) + +PECL2*(PMLZ(IM4,K,2,L)-PMLZ(IM4,K,4,L))) PMLZ(I,K,4,L) = PMLZ(I,K,2,L) 1455 CONTINUE ENDIF C C LOWER X SIDE C IF (SIDEL(L).EQ.'PMC') THEN I = IABCL(L)+1 IP2 = I+2 IP4 = I+4 DO 1465 K=0,MPML,2 C FOR Hz PMLZ(I,K,1,L) = 0.5* + (PECL1*(PMLZ(IP2,K,1,L)+PMLZ(IP2,K,3,L)) + +PECL2*(PMLZ(IP4,K,1,L)+PMLZ(IP4,K,3,L))) PMLZ(I,K,3,L) = -PMLZ(I,K,1,L) PMLZ(I,K,2,L) = 0.5* + (PECL1*(PMLZ(IP2,K,2,L)+PMLZ(IP2,K,4,L)) + +PECL2*(PMLZ(IP4,K,2,L)+PMLZ(IP4,K,4,L))) PMLZ(I,K,4,L) = -PMLZ(I,K,2,L) 1465 CONTINUE ENDIF C C TRANSFER ELECTRIC FIELDS AT PML INTERFACE TO THE NORMAL DIMENSION C DO 1500 IB=1,NB DO 1500 I=IL+IABCL(L),IL+IABCU(L),2 V(I,KABCL(L),IB) = PMLZ(I ,0,IB,L)+PMLZ(I+1,1,IB,L) 1500 CONTINUE C C CHECK EYX & EYZ FIELD VALUES C c DO 1510 IB=1,NB c DO 1520 K=0,MPML,2 c DO 1520 I=IL+IABCL(L),IL+IABCU(L),2 c TMPEYX(I,K,IB) = PMLZ(I,K,IB,L) !EXTRACT ONLY EYX c1520 CONTINUE c DO 1530 K=1,MPML+1,2 c DO 1530 I=IL+IABCL(L)+1,IL+IABCU(L)+1,2 c TMPEYZ(I,K,IB) = PMLZ(I,K,IB,L) !EXTRACT ONLY EYZ c1530 CONTINUE c IF (IT/NPRNT*NPRNT.EQ.IT) THEN c if (ib.eq.1) then c PRINT *,'EYX FIELD;',IB,'-TH COEFFICIENT; IT=',IT c CALL PRNT(IN+2,MPML+3,TMPEYX(0,-1,IB)) c PRINT *,'EYZ FIELD;',IB,'-TH COEFFICIENT; IT=',IT c CALL PRNT(IN+2,MPML+3,TMPEYZ(0,-1,IB)) c endif c ENDIF c1510 CONTINUE C ELSEIF (IABCL(L).EQ.IABCU(L)) THEN C C X-DIRECTION (ABC INTERFACE PERPENDICULAR TO X-AXIS) C C ------- TO BE IMPLEMENTED ------- C C C ENDIF C 1400 CONTINUE C C ------------------------------------------------ C HX AND HZ NODE STANDARD UPDATE (INTERNAL POINTS) C ------------------------------------------------ C DO 1600 L=1,NB C HX NODE DO 1620 K=KLP,KUM,2 DO 1620 I=IL,IU,2 V(I,K,L) = V(I,K,L)+CMX*(-V(I,K+1,L)+V(I,K-1,L)) 1620 CONTINUE C HZ NODE DO 1640 K=KL,KU,2 DO 1640 I=ILP,IUM,2 V(I,K,L) = V(I,K,L)+CMZ*( V(I+1,K,L)-V(I-1,K,L)) 1640 CONTINUE 1600 CONTINUE C C ---------------------------------------------- C SET MASK OUTSIDE THE ACTUAL CALCULATION REGION C ---------------------------------------------- C C ACTUAL CALCULATION REGION : MASK = 1.0 C OUTSIDE OF THE REGION : MASK = 0.0 C DO 1650 L=1,NB C HX DO 1660 K=KLP,KUM,2 DO 1660 I=IL,IU,2 V(I,K,L) = MASK(I,K)*V(I,K,L) 1660 CONTINUE C HZ DO 1670 K=KL,KU,2 DO 1670 I=ILP,IUM,2 V(I,K,L) = MASK(I,K)*V(I,K,L) 1670 CONTINUE 1650 CONTINUE C C ------------ C PMC BOUNDARY C ------------ C DO 1700 L=1,NPMC C C HZ NODE (PMC PERPENDICULAR TO X-AXIS) C C LOWER SIDE IF (PMCDRC(L).EQ.'XL') THEN I = IPMCL(L) IP2 = I+2 IP4 = I+4 DO 1720 K=KPMCL(L),KPMCU(L),2 V(I,K,1) = 0.5*(PECL1*(V(IP2,K,1)+V(IP2,K,3)) + +PECL2*(V(IP4,K,1)+V(IP4,K,3))) V(I,K,3) = -V(I,K,1) V(I,K,2) = 0.5*(PECL1*(V(IP2,K,2)+V(IP2,K,4)) + +PECL2*(V(IP4,K,2)+V(IP4,K,4))) V(I,K,4) = -V(I,K,2) 1720 CONTINUE C UPPER SIDE ELSEIF (PMCDRC(L).EQ.'XU') THEN I = IPMCL(L) IM2 = I-2 IM4 = I-4 DO 1740 K=KPMCL(L),KPMCU(L),2 V(I,K,1) = 0.5*(PECL1*(V(IM2,K,1)-V(IM2,K,3)) + +PECL2*(V(IM4,K,1)-V(IM4,K,3))) V(I,K,3) = V(I,K,1) V(I,K,2) = 0.5*(PECL1*(V(IM2,K,2)-V(IM2,K,4)) + +PECL2*(V(IM4,K,2)-V(IM4,K,4))) V(I,K,4) = V(I,K,2) 1740 CONTINUE ENDIF C C HX NODE (PMC PERPENDICULAR TO Z-AXIS) C C LOWER SIDE IF (PMCDRC(L).EQ.'ZL') THEN K = KPMCL(L) KP2 = K+2 KP4 = K+4 DO 1750 I=IPMCL(L),IPMCU(L),2 V(I,K,1) = 0.5*(PECL1*(V(I,KP2,1)+V(I,KP2,2)) + +PECL2*(V(I,KP4,1)+V(I,KP4,2))) V(I,K,2) = -V(I,K,1) V(I,K,3) = 0.5*(PECL1*(V(I,KP2,3)+V(I,KP2,4)) + +PECL2*(V(I,KP4,3)+V(I,KP4,4))) V(I,K,4) = -V(I,K,3) 1750 CONTINUE C UPPER SIDE ELSEIF (PMCDRC(L).EQ.'ZU') THEN K = KPMCL(L) KM2 = K-2 KM4 = K-4 DO 1760 I=IPMCL(L),IPMCU(L),2 V(I,K,1) = 0.5*(PECL1*(V(I,KM2,1)-V(I,KM2,2)) + +PECL2*(V(I,KM4,1)-V(I,KM4,2))) V(I,K,2) = V(I,K,1) V(I,K,3) = 0.5*(PECL1*(V(I,KM2,3)-V(I,KM2,4)) + +PECL2*(V(I,KM4,3)-V(I,KM4,4))) V(I,K,4) = V(I,K,3) 1760 CONTINUE ENDIF C 1700 CONTINUE C C ------------------------------------------- C IMPEDANCE BOUNDARY C (VALID ONLY FOR NORMAL INCIDENT PLANE WAVE) C ------------------------------------------- C C EY NODE C DO 2000 L=1,NIMP C C IMPEDANCE BOUNDARY PERPENDICULAR TO X-AXIS C IF (IIMPL(L).EQ.IIMPU(L)) THEN I = IIMPL(L) C LOWER SIDE IF (IMPDRC(L).EQ.'L') THEN DO 2010 IB = 1,NB DO 2010 K=KIMPL(L),KIMPU(L),2 C -- TO BE IMPLEMENTED -- 2010 CONTINUE C UPPER SIDE ELSEIF (IMPDRC(L).EQ.'U') THEN DO 2020 IB = 1,NB DO 2020 K=KIMPL(L),KIMPU(L),2 C -- TO BE IMPLEMENTED -- 2020 CONTINUE ENDIF C C IMPEDANCE BOUNDARY PERPENDICULAR TO Z-AXIS C ELSEIF (KIMPL(L).EQ.KIMPU(L)) THEN K = KIMPL(L) C LOWER SIDE IF (IMPDRC(L).EQ.'L') THEN KP = K+1 KP2 = K+2 KP3 = K+3 DO 2030 IB = 1,3,2 !NORMAL SCATTERING DO 2030 I=IIMPL(L),IIMPU(L),2 V(I,KP,IB) = CIX0(L)*PRVK1(I,IB,L) + +CIX1(L)*(4.0*2.0*DX*V(I,KP2,IB) + +ZIMP(L)*(V(I,KP3,IB)+PRVK3(I,IB,L))) PRVK1(I,IB,L) = V(I,KP,IB) !UPDATE PREVIOUS H-FILED PRVK3(I,IB,L) = V(I,KP3,IB) !UPDATE PREVIOUS H-FILED 2030 CONTINUE DO 2035 IB = 2,4,2 !INVERSE SCATTERING DO 2035 I=IIMPL(L),IIMPU(L),2 V(I,KP,IB) = CIX0I(L)*PRVK1(I,IB,L) + +CIX1I(L)*(4.0*2.0*DX*V(I,KP2,IB) + +ZIMPI(L)*(V(I,KP3,IB)+PRVK3(I,IB,L))) PRVK1(I,IB,L) = V(I,KP,IB) !UPDATE PREVIOUS H-FILED PRVK3(I,IB,L) = V(I,KP3,IB) !UPDATE PREVIOUS H-FILED 2035 CONTINUE C UPPER SIDE ELSEIF (IMPDRC(L).EQ.'U') THEN KM = K-1 KM2 = K-2 KM3 = K-3 DO 2040 IB = 1,3,2 !NORMAL SCATTERING DO 2040 I=IIMPL(L),IIMPU(L),2 V(I,KM,IB) = CIX0(L)*PRVK1(I,IB,L) + -CIX1(L)*(4.0*2.0*DX*V(I,KM2,IB) + -ZIMP(L)*(V(I,KM3,IB)+PRVK3(I,IB,L))) PRVK1(I,IB,L) = V(I,KM,IB) !UPDATE PREVIOUS H-FILED PRVK3(I,IB,L) = V(I,KM3,IB) !UPDATE PREVIOUS H-FILED 2040 CONTINUE DO 2045 IB = 2,4,2 !INVERSE SCATTERING DO 2045 I=IIMPL(L),IIMPU(L),2 V(I,KM,IB) = CIX0I(L)*PRVK1(I,IB,L) + -CIX1I(L)*(4.0*2.0*DX*V(I,KM2,IB) + -ZIMPI(L)*(V(I,KM3,IB)+PRVK3(I,IB,L))) PRVK1(I,IB,L) = V(I,KM,IB) !UPDATE PREVIOUS H-FILED PRVK3(I,IB,L) = V(I,KM3,IB) !UPDATE PREVIOUS H-FILED 2045 CONTINUE ENDIF ENDIF C 2000 CONTINUE C C ------------------------------------------------------------- C SINGULARITY CORRECTION AT PEC 90-DEGREE CORNER NODE (REF.[2]) C ------------------------------------------------------------- C IF (NSNG90.EQ.1) THEN C DO 1800 L=1,NECNR I = IECNR(L) K = KECNR(L) C IF (PECNR(L).EQ.'LL') THEN C C CALCULATE EY-FIELD IN LU-BASIS AT (I,K+2),(I+2,K) AND (I,K) C DO 1770 IB=1,NB VY01(IB,L) = 0.0 VY10(IB,L) = 0.0 VY00(IB,L) = 0.0 DO 1770 M=1,NB VY01(IB,L) = VY01(IB,L)+A(M,IB)*V(I,K+2,M) VY10(IB,L) = VY10(IB,L)+A(M,IB)*V(I+2,K,M) VY00(IB,L) = VY00(IB,L)+A(M,IB)*V(I,K,M) 1770 CONTINUE C C UPDATE H-FILEDS FOR LU-BASIS C DO 1780 IB=1,NB HX0H(IB,L) = HX0H(IB,L)+CMX*(-VY01(IB,L)+VY00(IB,L)) HZH0(IB,L) = HZH0(IB,L)+CMZ*( VY10(IB,L)-VY00(IB,L)) 1780 CONTINUE C C UPDATE H-FIELDS FOR LL-SUBNODE WITH SINGULARITY-CORRECTION C AIX0H(L) = AIX0H(L)-CMX*C90*VY01(LL,L) AIZH0(L) = AIZH0(L)+CMZ*C90*VY10(LL,L) C C REPLACE H-FIELD AT LU SUBNODE WITH CORRECTED VALUES C HX0H(LL,L) = AIX0H(L) HZH0(LL,L) = AIZH0(L) C C RE-DECOMPOSE H-FIELD VALUES INTO HAAR BASIS COEFFICIENTS C DO 1790 IB=1,NB VXTMP = 0.0 VZTMP = 0.0 DO 1795 M=1,NB VXTMP = VXTMP+A(M,IB)*HX0H(M,L) VZTMP = VZTMP+A(M,IB)*HZH0(M,L) 1795 CONTINUE V(I,K+1,IB) = VXTMP V(I+1,K,IB) = VZTMP 1790 CONTINUE C ELSEIF (PECNR(L).EQ.'LU') THEN C C CALCULATE EY-FIELD IN LU-BASIS AT (I,K-2),(I+2,K) AND (I,K) C DO 1807 IB=1,NB VY01M(IB,L) = 0.0 VY10 (IB,L) = 0.0 VY00 (IB,L) = 0.0 DO 1810 M=1,NB VY01M(IB,L) = VY01M(IB,L)+A(M,IB)*V(I,K-2,M) VY10 (IB,L) = VY10 (IB,L)+A(M,IB)*V(I+2,K,M) VY00 (IB,L) = VY00 (IB,L)+A(M,IB)*V(I,K,M) 1810 CONTINUE 1807 CONTINUE C C UPDATE STANDARD H-FILEDS IN LU-BASIS C DO 1817 IB=1,NB HX0HM(IB,L) = HX0HM(IB,L)+CMX*(-VY00(IB,L)+VY01M(IB,L)) HZH0 (IB,L) = HZH0 (IB,L)+CMZ*( VY10(IB,L)-VY00(IB,L)) 1817 CONTINUE C C UPDATE H-FIELDS AT LU-SUBNODE WITH SINGULARITY-CORRECTION C AIX0HM(L) = AIX0HM(L)+CMX*C90*VY01M(LU,L) AIZH0 (L) = AIZH0 (L)+CMZ*C90*VY10 (LU,L) C C REPLACE H-FIELD AT LU-SUBNODE WITH CORRECTED VALUES C HX0HM(LU,L) = AIX0HM(L) HZH0 (LU,L) = AIZH0 (L) C C RE-DECOMPOSE H-FIELD VALUES INTO HAAR BASIS COEFFICIENTS C DO 1830 IB=1,NB VXTMP = 0.0 VZTMP = 0.0 DO 1835 M=1,NB VXTMP = VXTMP+A(M,IB)*HX0HM(M,L) VZTMP = VZTMP+A(M,IB)*HZH0 (M,L) 1835 CONTINUE V(I,K-1,IB) = VXTMP V(I+1,K,IB) = VZTMP 1830 CONTINUE C ELSEIF (PECNR(L).EQ.'UL') THEN C C CALCULATE EY-FIELD IN LU-BASIS AT (I,K+2),(I-2,K) AND (I,K) C DO 1850 IB=1,NB VY01 (IB,L) = 0.0 VY1M0(IB,L) = 0.0 VY00 (IB,L) = 0.0 DO 1850 M=1,NB VY01 (IB,L) = VY01 (IB,L)+A(M,IB)*V(I,K+2,M) VY1M0(IB,L) = VY1M0(IB,L)+A(M,IB)*V(I-2,K,M) VY00 (IB,L) = VY00 (IB,L)+A(M,IB)*V(I,K,M) 1850 CONTINUE C C UPDATE H-FILEDS FOR UL-BASIS C DO 1860 IB=1,NB HX0H (IB,L) = HX0H (IB,L)+CMX*(-VY01(IB,L)+VY00 (IB,L)) HZHM0(IB,L) = HZHM0(IB,L)+CMZ*( VY00(IB,L)-VY1M0(IB,L)) 1860 CONTINUE C C UPDATE H-FIELDS FOR UL-SUBNODE WITH SINGULARITY-CORRECTION C AIX0H (L) = AIX0H (L)-CMX*C90*VY01 (UL,L) AIZHM0(L) = AIZHM0(L)-CMZ*C90*VY1M0(UL,L) C C REPLACE H-FIELD AT LU SUBNODE WITH CORRECTED VALUES C HX0H (UL,L) = AIX0H (L) HZHM0(UL,L) = AIZHM0(L) C C RE-DECOMPOSE H-FIELD VALUES INTO HAAR BASIS COEFFICIENTS C DO 1870 IB=1,NB VXTMP = 0.0 VZTMP = 0.0 DO 1875 M=1,NB VXTMP = VXTMP+A(M,IB)*HX0H (M,L) VZTMP = VZTMP+A(M,IB)*HZHM0(M,L) 1875 CONTINUE V(I,K+1,IB) = VXTMP V(I-1,K,IB) = VZTMP 1870 CONTINUE C ELSEIF (PECNR(L).EQ.'UU') THEN C C CALCULATE EY-FIELD IN LU-BASIS AT (I,K-2),(I-2,K) AND (I,K) C DO 1880 IB=1,NB VY01M(IB,L) = 0.0 VY1M0(IB,L) = 0.0 VY00 (IB,L) = 0.0 DO 1880 M=1,NB VY01M(IB,L) = VY01M(IB,L)+A(M,IB)*V(I,K-2,M) VY1M0(IB,L) = VY1M0(IB,L)+A(M,IB)*V(I-2,K,M) VY00 (IB,L) = VY00 (IB,L)+A(M,IB)*V(I,K,M) 1880 CONTINUE C C UPDATE H-FILEDS FOR LU-BASIS C DO 1890 IB=1,NB HX0HM(IB,L) = HX0HM(IB,L)+CMX*(-VY00(IB,L)+VY01M(IB,L)) HZHM0(IB,L) = HZHM0(IB,L)+CMZ*( VY00(IB,L)-VY1M0(IB,L)) 1890 CONTINUE C C UPDATE H-FIELDS FOR UL-SUBNODE WITH SINGULARITY-CORRECTION C AIX0HM(L) = AIX0HM(L)+CMX*C90*VY01M(UU,L) AIZHM0(L) = AIZHM0(L)-CMZ*C90*VY1M0(UU,L) C C REPLACE H-FIELD AT LU SUBNODE WITH CORRECTED VALUES C HX0HM(UU,L) = AIX0HM(L) HZHM0(UU,L) = AIZHM0(L) C C RE-DECOMPOSE H-FIELD VALUES INTO HAAR BASIS COEFFICIENTS C DO 1900 IB=1,NB VXTMP = 0.0 VZTMP = 0.0 DO 1905 M=1,NB VXTMP = VXTMP+A(M,IB)*HX0HM(M,L) VZTMP = VZTMP+A(M,IB)*HZHM0(M,L) 1905 CONTINUE V(I,K-1,IB) = VXTMP V(I-1,K,IB) = VZTMP 1900 CONTINUE C ENDIF C 1800 CONTINUE C ENDIF C C ------------------------------------------------------- C SINGULARITY CORRECTION FOR 0-DEGREE EDGE NODE (REF.[2]) C (ONLY FOR ZU-SIDE THIN PEC) C ------------------------------------------------------- C IF (NSNG0.EQ.1) THEN C DO 2100 L=1,NTPEC C C X-UPPER SIDE (SINGULARITY AT LU-SUBNODE) C IF (TPECEDGU(L).EQ.1.OR.TPECEDGU(L).EQ.2) THEN C I = ITPECU(L) K = KTPECL(L) KP = K+1 KP2 = K+2 KM = K-1 KM2 = K-2 IP = I+1 IP2 = I+2 C DO 2110 IB=1,NB EUVY01M(IB,L) = 0.0 EUVY10 (IB,L) = 0.0 EUVY00 (IB,L) = 0.0 EUVY01 (IB,L) = 0.0 DO 2120 M=1,NB EUVY01M(IB,L) = EUVY01M(IB,L)+A(M,IB)*V(I ,KM2,M) EUVY10 (IB,L) = EUVY10 (IB,L)+A(M,IB)*V(IP2,K ,M) EUVY00 (IB,L) = EUVY00 (IB,L)+A(M,IB)*V(I ,K ,M) EUVY01 (IB,L) = EUVY01 (IB,L)+A(M,IB)*V(I ,KP2,M) 2120 CONTINUE 2110 CONTINUE C C UPDATE STANDARD H-FILEDS IN LU-BASIS C DO 2130 IB=1,NB EUHX0HM(IB,L) = EUHX0HM(IB,L) + +CMX*(-EUVY00(IB,L)+EUVY01M(IB,L)) EUHZH0 (IB,L) = EUHZH0 (IB,L) + +CMZ*( EUVY10(IB,L)-EUVY00(IB,L)) EUHX0H (IB,L) = EUHX0H (IB,L) + +CMX*(-EUVY01(IB,L)+EUVY00(IB,L)) 2130 CONTINUE C C UPDATE H-FIELDS AT LU-SUBNODE WITH SINGULARITY-CORRECTION C EUIX0HM(L) = EUIX0HM(L)+CMX*C00*EUVY01M(LU,L) EUIZH0 (L) = EUIZH0 (L)+CMZ*C00*EUVY10 (LU,L) EUIX0H (L) = EUIX0H (L)-CMX*C00*EUVY01 (LU,L) C C REPLACE H-FIELD AT LU-SUBNODE WITH CORRECTED VALUES C EUHX0HM(LU,L) = EUIX0HM(L) EUHZH0 (LU,L) = EUIZH0 (L) EUHX0H (LU,L) = EUIX0H (L) C C RE-DECOMPOSE H-FIELD VALUES INTO HAAR BASIS COEFFICIENTS C DO 2140 IB=1,NB VXTMPA = 0.0 VZTMP = 0.0 VXTMPB = 0.0 DO 2150 M=1,NB VXTMPA = VXTMPA+A(M,IB)*EUHX0HM(M,L) VZTMP = VZTMP +A(M,IB)*EUHZH0 (M,L) VXTMPB = VXTMPB+A(M,IB)*EUHX0H (M,L) 2150 CONTINUE V(I ,KM,IB) = VXTMPA V(IP,K ,IB) = VZTMP V(I ,KP,IB) = VXTMPB 2140 CONTINUE C ENDIF C C X-LOWER SIDE (SINGULARITY AT UU-SUBNODE) C IF (TPECEDGL(L).EQ.1.OR.TPECEDGL(L).EQ.2) THEN C I = ITPECL(L) K = KTPECL(L) C KP = K+1 KP2 = K+2 KM = K-1 KM2 = K-2 IM = I-1 IM2 = I-2 C DO 2210 IB=1,NB ELVY01M(IB,L) = 0.0 ELVY1M0(IB,L) = 0.0 ELVY00 (IB,L) = 0.0 ELVY01 (IB,L) = 0.0 DO 2220 M=1,NB ELVY01M(IB,L) = ELVY01M(IB,L)+A(M,IB)*V(I ,KM2,M) ELVY1M0(IB,L) = ELVY1M0(IB,L)+A(M,IB)*V(IM2,K ,M) ELVY00 (IB,L) = ELVY00 (IB,L)+A(M,IB)*V(I ,K ,M) ELVY01 (IB,L) = ELVY01 (IB,L)+A(M,IB)*V(I ,KP2,M) 2220 CONTINUE 2210 CONTINUE C C UPDATE STANDARD H-FILEDS IN LU-BASIS C DO 2230 IB=1,NB ELHX0HM(IB,L) = ELHX0HM(IB,L) + +CMX*(-ELVY00(IB,L)+ELVY01M(IB,L)) ELHZHM0(IB,L) = ELHZHM0(IB,L) + +CMZ*( ELVY00(IB,L)-ELVY1M0(IB,L)) ELHX0H (IB,L) = ELHX0H (IB,L) + +CMX*(-ELVY01(IB,L)+ELVY00(IB,L)) 2230 CONTINUE C C UPDATE H-FIELDS AT UU-SUBNODE WITH SINGULARITY-CORRECTION C ELIX0HM(L) = ELIX0HM(L)+CMX*C00*ELVY01M(UU,L) ELIZHM0(L) = ELIZHM0(L)-CMZ*C00*ELVY1M0(UU,L) ELIX0H (L) = ELIX0H (L)-CMX*C00*ELVY01 (UU,L) C C REPLACE H-FIELD AT UU-SUBNODE WITH CORRECTED VALUES C ELHX0HM(UU,L) = ELIX0HM(L) ELHZHM0(UU,L) = ELIZHM0(L) ELHX0H (UU,L) = ELIX0H (L) C C RE-DECOMPOSE H-FIELD VALUES INTO HAAR BASIS COEFFICIENTS C DO 2240 IB=1,NB VXTMPA = 0.0 VZTMP = 0.0 VXTMPB = 0.0 DO 2250 M=1,NB VXTMPA = VXTMPA+A(M,IB)*ELHX0HM(M,L) VZTMP = VZTMP +A(M,IB)*ELHZHM0(M,L) VXTMPB = VXTMPB+A(M,IB)*ELHX0H (M,L) 2250 CONTINUE V(I ,KM,IB) = VXTMPA V(IM,K ,IB) = VZTMP V(I ,KP,IB) = VXTMPB 2240 CONTINUE C ENDIF C 2100 CONTINUE C ENDIF C C **************** C END OF EQUATIONS C **************** C C ************************************ C FIELD DISTRIBUTION SNAPSHOTS STORAGE C ************************************ C IF (IT/NSTORE*NSTORE.EQ.IT) THEN C IS = IT/NSTORE C C -------- C EY FIELD C -------- C C EY - TOTAL FIELD C DO 3000 K=KL,KU,2 DO 3000 I=IL,IU,2 DO 3000 L=1,NB FLDOUT(I ,K ,1,IS) = FLDOUT(I ,K ,1,IS)+A(L,LL)*V(I,K,L) FLDOUT(I ,K+1,1,IS) = FLDOUT(I ,K+1,1,IS)+A(L,LU)*V(I,K,L) FLDOUT(I+1,K ,1,IS) = FLDOUT(I+1,K ,1,IS)+A(L,UL)*V(I,K,L) FLDOUT(I+1,K+1,1,IS) = FLDOUT(I+1,K+1,1,IS)+A(L,UU)*V(I,K,L) 3000 CONTINUE C C EY - 2D HAAR BASIS COEFFICIENTS C DO 3100 K=KL,KU,2 DO 3100 I=IL,IU,2 FLDOUT(I/2 ,K/2 ,2,IS) = V(I,K,1) !PP FLDOUT(I/2 ,K/2+KU/2+1,2,IS) = V(I,K,2) !PS FLDOUT(I/2+IU/2+1,K/2 ,2,IS) = V(I,K,3) !SP FLDOUT(I/2+IU/2+1,K/2+KU/2+1,2,IS) = V(I,K,4) !SS 3100 CONTINUE C C EY - FIELDS INTERPOLATED FROM SCALING COEFFICIENTS C C (EY - BILINEAR INTERPOLATION FOR INNER POINTS) C DO 3110 K=KL,KU-2,2 DO 3110 I=IL,IU-2,2 FLDOUT(I+1,K+1,3,IS) = + 0.5*(9.0*V(I ,K ,1)+3.0*V(I+2,K ,1) + +3.0*V(I ,K+2,1)+ V(I+2,K+2,1))/16.0 FLDOUT(I+2,K+1,3,IS) = + 0.5*(3.0*V(I ,K ,1)+9.0*V(I+2,K ,1) + + V(I ,K+2,1)+3.0*V(I+2,K+2,1))/16.0 FLDOUT(I+1,K+2,3,IS) = + 0.5*(3.0*V(I ,K ,1)+ V(I+2,K ,1) + +9.0*V(I ,K+2,1)+3.0*V(I+2,K+2,1))/16.0 FLDOUT(I+2,K+2,3,IS) = + 0.5*( V(I ,K ,1)+3.0*V(I+2,K ,1) + +3.0*V(I ,K+2,1)+9.0*V(I+2,K+2,1))/16.0 3110 CONTINUE C C (EY - QUADRATIC EXTRAPOLATION FOR BOUNDARY FIELD VALUES) C SIDES K = KL DO 3120 I=IL,IU-2,2 FLDOUT(I+1,K ,3,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I+2,K ,1) + +FLE21*V(I ,K+2,1)+FLE22*V(I+2,K+2,1) + +FLE31*V(I ,K+4,1)+FLE32*V(I+2,K+4,1)) FLDOUT(I+2,K ,3,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I+2,K ,1) + +FLE22*V(I ,K+2,1)+FLE21*V(I+2,K+2,1) + +FLE32*V(I ,K+4,1)+FLE31*V(I+2,K+4,1)) 3120 CONTINUE K = KU DO 3130 I=IL,IU-2,2 FLDOUT(I+1,K+1,3,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I+2,K ,1) + +FLE21*V(I ,K-2,1)+FLE22*V(I+2,K-2,1) + +FLE31*V(I ,K-4,1)+FLE32*V(I+2,K-4,1)) FLDOUT(I+2,K+1,3,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I+2,K ,1) + +FLE22*V(I ,K-2,1)+FLE21*V(I+2,K-2,1) + +FLE32*V(I ,K-4,1)+FLE31*V(I+2,K-4,1)) 3130 CONTINUE I = IL DO 3140 K=KL,KU-2,2 FLDOUT(I ,K+1,3,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I ,K+2,1) + +FLE21*V(I+2,K ,1)+FLE22*V(I+2,K+2,1) + +FLE31*V(I+4,K ,1)+FLE32*V(I+4,K+2,1)) FLDOUT(I ,K+2,3,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I ,K+2,1) + +FLE22*V(I+2,K ,1)+FLE21*V(I+2,K+2,1) + +FLE32*V(I+4,K ,1)+FLE31*V(I+4,K+2,1)) 3140 CONTINUE I = IU DO 3150 K=KL,KU-2,2 FLDOUT(I+1,K+1,3,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I ,K+2,1) + +FLE21*V(I-2,K ,1)+FLE22*V(I-2,K+2,1) + +FLE31*V(I-4,K ,1)+FLE32*V(I-4,K+2,1)) FLDOUT(I+1,K+2,3,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I ,K+2,1) + +FLE22*V(I-2,K ,1)+FLE21*V(I-2,K+2,1) + +FLE32*V(I-4,K ,1)+FLE31*V(I-4,K+2,1)) 3150 CONTINUE C CORNER POINTS FLDOUT(IL,KL,3,IS) = 0.5*( + 2.0*FLDOUT(ILP ,KL ,3,IS)-FLDOUT(ILP2,KL ,3,IS) + +2.0*FLDOUT(IL ,KLP,3,IS)-FLDOUT(IL ,KLP2,3,IS)) FLDOUT(IU+1,KL,3,IS) = 0.5*( + 2.0*FLDOUT(IU ,KL ,3,IS)-FLDOUT(IU-1,KL ,3,IS) + +2.0*FLDOUT(IU+1,KLP,3,IS)-FLDOUT(IU+1,KLP2,3,IS)) FLDOUT(IL,KU+1,3,IS) = 0.5*( + 2.0*FLDOUT(ILP,KU+1,3,IS)-FLDOUT(ILP2,KU+1,3,IS) + +2.0*FLDOUT(IL ,KU ,3,IS)-FLDOUT(IL ,KU-1,3,IS)) FLDOUT(IU+1,KU+1,3,IS) = 0.5*( + 2.0*FLDOUT(IU ,KU+1,3,IS)-FLDOUT(IU-1,KU+1,3,IS) + +2.0*FLDOUT(IU+1,KU ,3,IS)-FLDOUT(IU+1,KU-1,3,IS)) C C (EY - QUADRATIC POLINOMIAL INTERPOLATION ONLY FOR SPECIFIED AREA) C DO 3180 N=1,NQITP DO 3180 K=KL+KQITPL(N),KL+KQITPU(N),2 DO 3180 I=IL+IQITPL(N),IL+IQITPU(N),2 FLDOUT(I ,K ,3,IS) = 0.5*( + FLI22*V(I-2,K-2,1)+FLI12*V(I-2,K ,1)+FLI23*V(I-2,K+2,1) + +FLI12*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI13*V(I ,K+2,1) + +FLI23*V(I+2,K-2,1)+FLI13*V(I+2,K ,1)+FLI33*V(I+2,K+2,1)) FLDOUT(I+1,K ,3,IS) = 0.5*( + FLI23*V(I-2,K-2,1)+FLI13*V(I-2,K ,1)+FLI33*V(I-2,K+2,1) + +FLI12*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI13*V(I ,K+2,1) + +FLI22*V(I+2,K-2,1)+FLI12*V(I+2,K ,1)+FLI23*V(I+2,K+2,1)) FLDOUT(I ,K+1,3,IS) = 0.5*( + FLI23*V(I-2,K-2,1)+FLI12*V(I-2,K ,1)+FLI22*V(I-2,K+2,1) + +FLI13*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI12*V(I ,K+2,1) + +FLI33*V(I+2,K-2,1)+FLI13*V(I+2,K ,1)+FLI23*V(I+2,K+2,1)) FLDOUT(I+1,K+1,3,IS) = 0.5*( + FLI33*V(I-2,K-2,1)+FLI13*V(I-2,K ,1)+FLI23*V(I-2,K+2,1) + +FLI13*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI12*V(I ,K+2,1) + +FLI23*V(I+2,K-2,1)+FLI12*V(I+2,K ,1)+FLI22*V(I+2,K+2,1)) 3180 CONTINUE C C -------- C HX FIELD C -------- C C HX - TOTAL FIELD C DO 3200 K=KLP,KUM,2 DO 3200 I=IL,IU,2 DO 3200 L=1,NB FLDOUT(I ,K ,4,IS) = FLDOUT(I ,K ,4,IS)+A(L,LL)*V(I,K,L) FLDOUT(I ,K+1,4,IS) = FLDOUT(I ,K+1,4,IS)+A(L,LU)*V(I,K,L) FLDOUT(I+1,K ,4,IS) = FLDOUT(I+1,K ,4,IS)+A(L,UL)*V(I,K,L) FLDOUT(I+1,K+1,4,IS) = FLDOUT(I+1,K+1,4,IS)+A(L,UU)*V(I,K,L) 3200 CONTINUE C C HX - 2D HAAR BASIS COEFFICIENTS C DO 3300 K=KLP,KUM,2 DO 3300 I=IL,IU,2 FLDOUT(I/2 ,K/2 ,5,IS) = V(I,K,1) !PP FLDOUT(I/2 ,K/2+KU/2+1,5,IS) = V(I,K,2) !PS FLDOUT(I/2+IU/2+1,K/2 ,5,IS) = V(I,K,3) !SP FLDOUT(I/2+IU/2+1,K/2+KU/2+1,5,IS) = V(I,K,4) !SS 3300 CONTINUE C C HX - FIELDS INTERPOLATED FROM SCALING COEFFICIENTS C C (HX - BILINEAR INTERPOLATION) C DO 3310 K=KLP,KUM-2,2 DO 3310 I=IL,IU-2,2 FLDOUT(I+1,K+1,6,IS) = + 0.5*(9.0*V(I ,K ,1)+3.0*V(I+2,K ,1) + +3.0*V(I ,K+2,1)+ V(I+2,K+2,1))/16.0 FLDOUT(I+2,K+1,6,IS) = + 0.5*(3.0*V(I ,K ,1)+9.0*V(I+2,K ,1) + + V(I ,K+2,1)+3.0*V(I+2,K+2,1))/16.0 FLDOUT(I+1,K+2,6,IS) = + 0.5*(3.0*V(I ,K ,1)+ V(I+2,K ,1) + +9.0*V(I ,K+2,1)+3.0*V(I+2,K+2,1))/16.0 FLDOUT(I+2,K+2,6,IS) = + 0.5*( V(I ,K ,1)+3.0*V(I+2,K ,1) + +3.0*V(I ,K+2,1)+9.0*V(I+2,K+2,1))/16.0 3310 CONTINUE C C (HX - QUADRATIC EXTRAPOLATION FOR BOUNDARY FIELD VALUES) C SIDES K = KLP DO 3320 I=IL,IU-2,2 FLDOUT(I+1,K ,6,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I+2,K ,1) + +FLE21*V(I ,K+2,1)+FLE22*V(I+2,K+2,1) + +FLE31*V(I ,K+4,1)+FLE32*V(I+2,K+4,1)) FLDOUT(I+2,K ,6,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I+2,K ,1) + +FLE22*V(I ,K+2,1)+FLE21*V(I+2,K+2,1) + +FLE32*V(I ,K+4,1)+FLE31*V(I+2,K+4,1)) 3320 CONTINUE K = KUM DO 3330 I=IL,IU-2,2 FLDOUT(I+1,K+1,6,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I+2,K ,1) + +FLE21*V(I ,K-2,1)+FLE22*V(I+2,K-2,1) + +FLE31*V(I ,K-4,1)+FLE32*V(I+2,K-4,1)) FLDOUT(I+2,K+1,6,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I+2,K ,1) + +FLE22*V(I ,K-2,1)+FLE21*V(I+2,K-2,1) + +FLE32*V(I ,K-4,1)+FLE31*V(I+2,K-4,1)) 3330 CONTINUE I = IL DO 3340 K=KLP,KUM-2,2 FLDOUT(I ,K+1,6,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I ,K+2,1) + +FLE21*V(I+2,K ,1)+FLE22*V(I+2,K+2,1) + +FLE31*V(I+4,K ,1)+FLE32*V(I+4,K+2,1)) FLDOUT(I ,K+2,6,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I ,K+2,1) + +FLE22*V(I+2,K ,1)+FLE21*V(I+2,K+2,1) + +FLE32*V(I+4,K ,1)+FLE31*V(I+4,K+2,1)) 3340 CONTINUE I = IU DO 3350 K=KLP,KUM-2,2 FLDOUT(I+1,K+1,6,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I ,K+2,1) + +FLE21*V(I-2,K ,1)+FLE22*V(I-2,K+2,1) + +FLE31*V(I-4,K ,1)+FLE32*V(I-4,K+2,1)) FLDOUT(I+1,K+2,6,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I ,K+2,1) + +FLE22*V(I-2,K ,1)+FLE21*V(I-2,K+2,1) + +FLE32*V(I-4,K ,1)+FLE31*V(I-4,K+2,1)) 3350 CONTINUE C CORNER POINTS FLDOUT(IL,KLP,6,IS) = 0.5*( + 2.0*FLDOUT(ILP ,KLP ,6,IS)-FLDOUT(ILP2,KLP ,6,IS) + +2.0*FLDOUT(IL ,KLP2,6,IS)-FLDOUT(IL ,KLP2+1,6,IS)) FLDOUT(IU+1,KLP,6,IS) = 0.5*( + 2.0*FLDOUT(IU ,KLP ,6,IS)-FLDOUT(IU-1,KLP ,6,IS) + +2.0*FLDOUT(IU+1,KLP2,6,IS)-FLDOUT(IU+1,KLP2+1,6,IS)) FLDOUT(IL,KU,6,IS) = 0.5*( + 2.0*FLDOUT(ILP ,KU ,6,IS)-FLDOUT(ILP2,KU ,6,IS) + +2.0*FLDOUT(IL ,KUM ,6,IS)-FLDOUT(IL ,KUM2 ,6,IS)) FLDOUT(IU+1,KU,6,IS) = 0.5*( + 2.0*FLDOUT(IU ,KU ,6,IS)-FLDOUT(IU-1,KU ,6,IS) + +2.0*FLDOUT(IU+1,KUM ,6,IS)-FLDOUT(IU+1,KUM2 ,6,IS)) C C (HX - QUADRATIC POLINOMIAL INTERPOLATION ONLY AT SPECIFIED AREA) C DO 3380 N=1,NQITP DO 3380 K=KL+KQITPL(N)+1,KL+KQITPU(N)-1,2 DO 3380 I=IL+IQITPL(N),IL+IQITPU(N),2 FLDOUT(I ,K ,6,IS) = 0.5*( + FLI22*V(I-2,K-2,1)+FLI12*V(I-2,K ,1)+FLI23*V(I-2,K+2,1) + +FLI12*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI13*V(I ,K+2,1) + +FLI23*V(I+2,K-2,1)+FLI13*V(I+2,K ,1)+FLI33*V(I+2,K+2,1)) FLDOUT(I+1,K ,6,IS) = 0.5*( + FLI23*V(I-2,K-2,1)+FLI13*V(I-2,K ,1)+FLI33*V(I-2,K+2,1) + +FLI12*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI13*V(I ,K+2,1) + +FLI22*V(I+2,K-2,1)+FLI12*V(I+2,K ,1)+FLI23*V(I+2,K+2,1)) FLDOUT(I ,K+1,6,IS) = 0.5*( + FLI23*V(I-2,K-2,1)+FLI12*V(I-2,K ,1)+FLI22*V(I-2,K+2,1) + +FLI13*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI12*V(I ,K+2,1) + +FLI33*V(I+2,K-2,1)+FLI13*V(I+2,K ,1)+FLI23*V(I+2,K+2,1)) FLDOUT(I+1,K+1,6,IS) = 0.5*( + FLI33*V(I-2,K-2,1)+FLI13*V(I-2,K ,1)+FLI23*V(I-2,K+2,1) + +FLI13*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI12*V(I ,K+2,1) + +FLI23*V(I+2,K-2,1)+FLI12*V(I+2,K ,1)+FLI22*V(I+2,K+2,1)) 3380 CONTINUE C C -------- C HZ FIELD C -------- C C HZ - TOTAL FIELD C DO 3400 K=KL,KU,2 DO 3400 I=ILP,IUM,2 DO 3400 L=1,NB FLDOUT(I ,K ,7,IS) = FLDOUT(I ,K ,7,IS)+A(L,LL)*V(I,K,L) FLDOUT(I ,K+1,7,IS) = FLDOUT(I ,K+1,7,IS)+A(L,LU)*V(I,K,L) FLDOUT(I+1,K ,7,IS) = FLDOUT(I+1,K ,7,IS)+A(L,UL)*V(I,K,L) FLDOUT(I+1,K+1,7,IS) = FLDOUT(I+1,K+1,7,IS)+A(L,UU)*V(I,K,L) 3400 CONTINUE C C HZ - 2D HAAR BASIS COEFFICIENTS C DO 3500 K=KL,KU,2 DO 3500 I=ILP,IUM,2 FLDOUT(I/2 ,K/2 ,8,IS) = V(I,K,1) !PP FLDOUT(I/2 ,K/2+KU/2+1,8,IS) = V(I,K,2) !PS FLDOUT(I/2+IU/2+1,K/2 ,8,IS) = V(I,K,3) !SP FLDOUT(I/2+IU/2+1,K/2+KU/2+1,8,IS) = V(I,K,4) !SS 3500 CONTINUE C C HZ - FIELDS INTERPOLATED FROM SCALING COEFFICIENTS C C (HZ - BILINEAR INTERPOLATION) C DO 3510 K=KL,KU-2,2 DO 3510 I=ILP,IUM-2,2 FLDOUT(I+1,K+1,9,IS) = + 0.5*(9.0*V(I ,K ,1)+3.0*V(I+2,K ,1) + +3.0*V(I ,K+2,1)+ V(I+2,K+2,1))/16.0 FLDOUT(I+2,K+1,9,IS) = + 0.5*(3.0*V(I ,K ,1)+9.0*V(I+2,K ,1) + + V(I ,K+2,1)+3.0*V(I+2,K+2,1))/16.0 FLDOUT(I+1,K+2,9,IS) = + 0.5*(3.0*V(I ,K ,1)+ V(I+2,K ,1) + +9.0*V(I ,K+2,1)+3.0*V(I+2,K+2,1))/16.0 FLDOUT(I+2,K+2,9,IS) = + 0.5*( V(I ,K ,1)+3.0*V(I+2,K ,1) + +3.0*V(I ,K+2,1)+9.0*V(I+2,K+2,1))/16.0 3510 CONTINUE C C (HZ - QUADRATIC EXTRAPOLATION FOR BOUNDARY FIELD VALUES) C SIDES K = KL DO 3520 I=ILP,IUM-2,2 FLDOUT(I+1,K ,9,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I+2,K ,1) + +FLE21*V(I ,K+2,1)+FLE22*V(I+2,K+2,1) + +FLE31*V(I ,K+4,1)+FLE32*V(I+2,K+4,1)) FLDOUT(I+2,K ,9,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I+2,K ,1) + +FLE22*V(I ,K+2,1)+FLE21*V(I+2,K+2,1) + +FLE32*V(I ,K+4,1)+FLE31*V(I+2,K+4,1)) 3520 CONTINUE K = KU DO 3530 I=ILP,IUM-2,2 FLDOUT(I+1,K+1,9,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I+2,K ,1) + +FLE21*V(I ,K-2,1)+FLE22*V(I+2,K-2,1) + +FLE31*V(I ,K-4,1)+FLE32*V(I+2,K-4,1)) FLDOUT(I+2,K+1,9,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I+2,K ,1) + +FLE22*V(I ,K-2,1)+FLE21*V(I+2,K-2,1) + +FLE32*V(I ,K-4,1)+FLE31*V(I+2,K-4,1)) 3530 CONTINUE I = ILP DO 3540 K=KL,KU-2,2 FLDOUT(I ,K+1,9,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I ,K+2,1) + +FLE21*V(I+2,K ,1)+FLE22*V(I+2,K+2,1) + +FLE31*V(I+4,K ,1)+FLE32*V(I+4,K+2,1)) FLDOUT(I ,K+2,9,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I ,K+2,1) + +FLE22*V(I+2,K ,1)+FLE21*V(I+2,K+2,1) + +FLE32*V(I+4,K ,1)+FLE31*V(I+4,K+2,1)) 3540 CONTINUE I = IUM DO 3550 K=KL,KU-2,2 FLDOUT(I+1,K+1,9,IS) = 0.5*( + FLE11*V(I ,K ,1)+FLE12*V(I ,K+2,1) + +FLE21*V(I-2,K ,1)+FLE22*V(I-2,K+2,1) + +FLE31*V(I-4,K ,1)+FLE32*V(I-4,K+2,1)) FLDOUT(I+1,K+2,9,IS) = 0.5*( + FLE12*V(I ,K ,1)+FLE11*V(I ,K+2,1) + +FLE22*V(I-2,K ,1)+FLE21*V(I-2,K+2,1) + +FLE32*V(I-4,K ,1)+FLE31*V(I-4,K+2,1)) 3550 CONTINUE C CORNER POINTS FLDOUT(ILP,KL,9,IS) = 0.5*( + 2.0*FLDOUT(ILP2,KL ,9,IS)-FLDOUT(ILP2+1,KL ,9,IS) + +2.0*FLDOUT(ILP ,KLP ,9,IS)-FLDOUT(ILP ,KLP2,9,IS)) FLDOUT(IU,KL,9,IS) = 0.5*( + 2.0*FLDOUT(IUM ,KL ,9,IS)-FLDOUT(IUM2 ,KL ,9,IS) + +2.0*FLDOUT(IU ,KLP ,9,IS)-FLDOUT(IU ,KLP2,9,IS)) FLDOUT(ILP,KU+1,9,IS) = 0.5*( + 2.0*FLDOUT(ILP2,KU+1,9,IS)-FLDOUT(ILP2+1,KU+1,9,IS) + +2.0*FLDOUT(ILP ,KU ,9,IS)-FLDOUT(ILP ,KU-1,9,IS)) FLDOUT(IU,KU+1,9,IS) = 0.5*( + 2.0*FLDOUT(IUM ,KU+1,9,IS)-FLDOUT(IUM2 ,KU+1,9,IS) + +2.0*FLDOUT(IU ,KU ,9,IS)-FLDOUT(IU ,KU-1,9,IS)) C C (HZ - QUADRATIC POLINOMIAL INTERPOLATION ONLY AT SPECIFIED AREA) C DO 3580 N=1,NQITP DO 3580 K=KL+KQITPL(N),KL+KQITPU(N),2 DO 3580 I=IL+IQITPL(N)+1,IL+IQITPU(N)-1,2 FLDOUT(I ,K ,9,IS) = 0.5*( + FLI22*V(I-2,K-2,1)+FLI12*V(I-2,K ,1)+FLI23*V(I-2,K+2,1) + +FLI12*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI13*V(I ,K+2,1) + +FLI23*V(I+2,K-2,1)+FLI13*V(I+2,K ,1)+FLI33*V(I+2,K+2,1)) FLDOUT(I+1,K ,9,IS) = 0.5*( + FLI23*V(I-2,K-2,1)+FLI13*V(I-2,K ,1)+FLI33*V(I-2,K+2,1) + +FLI12*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI13*V(I ,K+2,1) + +FLI22*V(I+2,K-2,1)+FLI12*V(I+2,K ,1)+FLI23*V(I+2,K+2,1)) FLDOUT(I ,K+1,9,IS) = 0.5*( + FLI23*V(I-2,K-2,1)+FLI12*V(I-2,K ,1)+FLI22*V(I-2,K+2,1) + +FLI13*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI12*V(I ,K+2,1) + +FLI33*V(I+2,K-2,1)+FLI13*V(I+2,K ,1)+FLI23*V(I+2,K+2,1)) FLDOUT(I+1,K+1,9,IS) = 0.5*( + FLI33*V(I-2,K-2,1)+FLI13*V(I-2,K ,1)+FLI23*V(I-2,K+2,1) + +FLI13*V(I ,K-2,1)+FLI11*V(I ,K ,1)+FLI12*V(I ,K+2,1) + +FLI23*V(I+2,K-2,1)+FLI12*V(I+2,K ,1)+FLI22*V(I+2,K+2,1)) 3580 CONTINUE C ENDIF C C ************************************** C EY FIELD DISTRIBUTION PRINT ON DISPLAY C ************************************** C IF (IT/NPRNT*NPRNT.EQ.IT) THEN C DO 4000 K=KL,KU+1 DO 4000 I=IL,IU+1 FLDPRNT(I,K) = 0.0 4000 CONTINUE C DO 4050 K=KL,KU,2 DO 4050 I=IL,IU,2 DO 4060 L=1,NB FLDPRNT(I ,K ) = FLDPRNT(I ,K )+A(L,LL)*V(I,K,L) FLDPRNT(I ,K+1) = FLDPRNT(I ,K+1)+A(L,LU)*V(I,K,L) FLDPRNT(I+1,K ) = FLDPRNT(I+1,K )+A(L,UL)*V(I,K,L) FLDPRNT(I+1,K+1) = FLDPRNT(I+1,K+1)+A(L,UU)*V(I,K,L) 4060 CONTINUE 4050 CONTINUE C VMAX = 0.0 DO 4100 K=KL,KU+1 DO 4100 I=IL,IU+1 VMAX = AMAX1(VMAX,ABS(FLDPRNT(I,K))) 4100 CONTINUE IF (VMAX.EQ.0.0) THEN VMAX = 1.0 ENDIF WRITE (6,6050) IT,VMAX 6050 FORMAT (//' IT = ',I8,/' MAX EY = ',G15.7/) C DO 4200 I=IU+1,IL,-1 WRITE (6,6100) (INT(10.0*ABS(FLDPRNT(I,K)/VMAX)),K=KL,KU+1) 4200 CONTINUE 6100 FORMAT (1H ,500I1) C ENDIF C C ************************************************ C TIME SERIES DATA STORAGE (BY DETECTING EY VALUE) C ************************************************ C IF (IT/NSTORE2*NSTORE2.EQ.IT) THEN C ITS = INT(IT/NSTORE2) C DO 4990 IP=0,NP TMOUT(ITS,IP) = 0.0 4990 CONTINUE C DO 5000 IP=1,NP VTMP = 0.0 DO 5100 I=IL+IPORTL(IP),IL+IPORTU(IP),2 DO 5100 K=KL+KPORTL(IP),KL+KPORTU(IP),2 VTMP = VTMP+V(I,K,1)+V(I,K,2)+V(I,K,3)+V(I,K,4) 5100 CONTINUE TMOUT(ITS,IP) = -VTMP + /FLOAT(ABS(IPORTU(IP)-IPORTL(IP))/2+1) + /FLOAT(ABS(KPORTU(IP)-KPORTL(IP))/2+1) 5000 CONTINUE C C -------------------------- C STORE RAW EXCITATION PULSE C -------------------------- C TMOUT(ITS,0) = ST C ENDIF C C ****************************** C ****************************** C END OF TIME MARCHING ALGORITHM C ****************************** C ****************************** C 10000 CONTINUE C C ***************************** C POSTPROCESSING OF OUTPUT DATA C ***************************** C ------------------------------------------------------------ C FIELD DISTRIBUTION NORMALIZATION AND DATA STORAGE INTO FILES C ------------------------------------------------------------ C PRINT * PRINT *,'POST-PROCESSING ...' C C ------------------------ C FIELD DATA NORMALIZATION C ------------------------ C C FINDING THE MAXIMUM VALUE C EMAX = -1.0E10 HMAX = -1.0E10 C DO 7000 JS=1,NS DO 7000 K=KL,KU+1 DO 7000 I=IL,IU+1 DO 7100 M=1,3 !EY TOTAL(1), COEFFS(2), INTERP(3) EMAX = AMAX1(EMAX,ABS(FLDOUT(I,K,M,JS))) 7100 CONTINUE DO 7200 M=4,6 !HX TOTAL(4), COEFFS(5), INTERP(6) HMAX = AMAX1(HMAX,ABS(FLDOUT(I,K,M,JS)/(2.0*DX))) 7200 CONTINUE DO 7220 M=7,9 !HZ TOTAL(7), COEFFS(8), INTERP(9) HMAX = AMAX1(HMAX,ABS(FLDOUT(I,K,M,JS)/(2.0*DZ))) 7220 CONTINUE 7000 CONTINUE C WRITE (6,6200) EMAX, HMAX 6200 FORMAT (1H ,5X,'EMAX = ',G15.5/1H ,5X,'HMAX = ',G15.5/) C C disable normalization C emax = 1.0 hmax = 1.0 C C NORMALIZATION C DO 7300 JS=1,NS DO 7300 K=KL,KU+1 DO 7300 I=IL,IU+1 DO 7400 M=1,3 !EY TOTAL(1), COEFFS(2), INTERP(3) FLDOUT(I,K,M,JS) = FLDOUT(I,K,M,JS)/EMAX 7400 CONTINUE DO 7500 M=4,6 !HX TOTAL(4), COEFFS(5), INTERP(6) FLDOUT(I,K,M,JS) = FLDOUT(I,K,M,JS)/(2.0*DX*HMAX) 7500 CONTINUE DO 7520 M=7,9 !HZ TOTAL(7), COEFFS(8), INTERP(9) FLDOUT(I,K,M,JS) = FLDOUT(I,K,M,JS)/(2.0*DZ*HMAX) 7520 CONTINUE 7300 CONTINUE C C CREATE FILE NAME C DO 7600 JS=1,NS NDIGIT1 = 48+(JS-JS/10*10) NDIGIT2 = 48+JS/10 FL(1,JS) = 'EYTOTL'//CHAR(NDIGIT2)//CHAR(NDIGIT1) FL(2,JS) = 'EYCOEF'//CHAR(NDIGIT2)//CHAR(NDIGIT1) FL(3,JS) = 'EYITRP'//CHAR(NDIGIT2)//CHAR(NDIGIT1) FL(4,JS) = 'HXTOTL'//CHAR(NDIGIT2)//CHAR(NDIGIT1) FL(5,JS) = 'HXCOEF'//CHAR(NDIGIT2)//CHAR(NDIGIT1) FL(6,JS) = 'HXITRP'//CHAR(NDIGIT2)//CHAR(NDIGIT1) FL(7,JS) = 'HZTOTL'//CHAR(NDIGIT2)//CHAR(NDIGIT1) FL(8,JS) = 'HZCOEF'//CHAR(NDIGIT2)//CHAR(NDIGIT1) FL(9,JS) = 'HZITRP'//CHAR(NDIGIT2)//CHAR(NDIGIT1) 7600 CONTINUE C C ------------------ C FIELD DATA STORAGE C ------------------ C PRINT *,'STORING DATA ...' C DO 7700 JF=1,NF !NUMBER OF FIELDS DO 7700 JS=1,INT(NMAX/NSTORE) !NUMBER OF SNAPSHOTS C OPEN (11,FILE=FL(JF,JS),FORM='FORMATTED') WRITE (11,6300) FL(JF,JS) C IF (JF.LE.3.OR.JF.GE.7) KOFFSET = 0 !EY OR HZ FIELD IF (JF.GE.4.AND.JF.LE.6) KOFFSET = 1 !HX FIELD DO 7800 K=KL+KOFFSET,KU+1-KOFFSET Z = DZ*FLOAT(K) IF (JF.LE.6) IOFFSET = 0 !EY OR HX FIELD IF (JF.GE.7) IOFFSET = 1 !HZ FIELD DO 7900 I=IL+IOFFSET,IU+1-IOFFSET X = DX*FLOAT(I) WRITE (11,6400) Z,X,FLDOUT(I,K,JF,JS) 7900 CONTINUE WRITE (11,*) 7800 CONTINUE C CLOSE(11) C 7700 CONTINUE C 6300 FORMAT ('# ',A8,' FIELD ON ZX PLANE' + /'# Z',14X,'X',14X,'FIELD VALUE') 6400 FORMAT (' ',3(G16.8,3X)) C C ---------------------- C STORE TIME SERIES DATA C ---------------------- C NDATA = INT(NMAX/NSTORE2) TIMESTEP = DT*FLOAT(NSTORE2)*UNIT/CVACUUM C OPEN (10,FILE='TIMEDATA',FORM='FORMATTED') WRITE(10,6500) 6500 FORMAT('# TD-MRA TIME SERIES DATA (IN MKSA UNIT)' + /'# FOR S-PARAMETER EXTRACTION' + /'#',18X,'TIME(SEC)',9X,'EXCITATION',17X,'V1',17X,'V2',17X,'V3') DO 8000 N=0,NDATA TIME = TIMESTEP*FLOAT(N) WRITE(10,6520) N,TIME,(TMOUT(N,I),I=0,NP) 8000 CONTINUE 6510 FORMAT(1H ,5(3X,G16.8)) 6520 FORMAT(1H ,I8,5(3X,G16.8)) CLOSE (10) C C ----------------------------------------------------------------- C SPECTRUM ANALYSIS BY DFT C C TIME AND FREQUENCY SIGNAL ASSIGNMENTS: C TMOUT(N,0) , FTOUT(M,0) : PURE EXCITATION PULSE C TMOUT(N,1 -- NP-1), FTOUT(M,1 -- NP-1) : SIGNAL AT PORTS C TMOUT(N,NP) , FTOUT(M,NP) : REFERENCE SIGNAL C ---------------------------------------------------------------- C IF (NDFT.EQ.1) THEN C PRINT *,'PROCESSING SPECTRUM ANARYSIS BY DFT ...' C DO 9000 N=0,NMAX/NSTORE2 TMOUT(N,1) = TMOUT(N,1)-TMOUT(N,NP) !SUBTRACT REFERENCE 9000 CONTINUE C OMG1 = 2.0*PI*(FCTR-0.5*FSPN) OMGR = 2.0*PI*FRES DTS = DTREAL*FLOAT(NSTORE2) OMG1DTS = 2.0*PI*(FCTR-0.5*FSPN)*DTS OMGRDTS = 2.0*PI*FRES*DTS C C DFT C DO 9010 IP=0,NP DO 9020 M=0,INT(FSPN/FRES) OMGDTS = OMG1DTS+OMGRDTS*FLOAT(M) R = 0.0 X = 0.0 DO 9030 N=0,NMAX/NSTORE2 R = R+TMOUT(N,IP)*COS(OMGDTS*FLOAT(N)) !DFT_REAL X = X-TMOUT(N,IP)*SIN(OMGDTS*FLOAT(N)) !DFT_IMAG 9030 CONTINUE R = R*DTS X = X*DTS FTOUT(M,IP) = SQRT(R*R+X*X) 9020 CONTINUE 9010 CONTINUE C C CALCULATE S-PARAMETERS FROM POWER SPECTRUM DENSITIES C DO 9040 IP=1,NP-1 DO 9040 M=0,INT(FSPN/FRES) SPARA(M,IP) = FTOUT(M,IP)/FTOUT(M,NP) 9040 CONTINUE C C STORE POWER SPECTRUM DENSITIES AND S-PARAMETERS C OPEN (12,FILE='FREQDATA') WRITE (12,6600) OPEN (13,FILE='FRQSPARA') WRITE (13,6620) F1 = FCTR-0.5*FSPN DO 9050 M=0,INT(FSPN/FRES) FREQ = F1+FRES*FLOAT(M) WRITE (12,6610) FREQ,(20.0*LOG10(FTOUT(M,IP)),IP=0,NP) WRITE (13,6630) FREQ,(20.0*LOG10(SPARA(M,IP)),IP=1,NP-1) 9050 CONTINUE CLOSE (12) CLOSE (13) C ENDIF C 6600 FORMAT ('# FREQENCY POWER SPECTRUM DENSITIES (dB)', + /'# ', + /'# ',10X,'FREQ.',5X,'EXCITATION', + 10X,'INPUT',9X,'OUTPUT',6X,'REFERENCE') 6610 FORMAT (1H ,5(3X,G15.7)) 6620 FORMAT ('# S-PARAMETERS (dB)', + /'# ', + /'# ',10X,'FREQ.',12X,'S11',12X,'S21') 6630 FORMAT (1H ,3(3X,G15.7)) C STOP END C********************************************************************* SUBROUTINE PRNT (N,M,A) C********************************************************************* C C Supplemental subroutine for checking values in an array A with C dimensions N x M. C DIMENSION A(N,M) C VMAX = 0.0 DO 1100 K=1,M DO 1100 I=1,N VMAX = AMAX1(VMAX,ABS(A(I,K))) 1100 CONTINUE IF (VMAX.EQ.0.0) THEN VMAX = 1.0 ENDIF WRITE (6,6050) N,M,VMAX 6050 FORMAT (//' ',I5,' X ',I5,' MATRIX, WITH MAX VALUE = ',G15.7/) C DO 1200 I=1,N WRITE (6,6100) (INT(10.0*ABS(A(I,K))/VMAX),K=1,M) 1200 CONTINUE 6100 FORMAT (1H ,1X,200I1) C RETURN END