IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) c10================================================================ c10 RAY TRACING PROGRAM FOR VLF WAVE (X-MODE) IN DE-MODEL c10 Dipole model is selected in REF by Call dipole c10 IGRF model is selected in REF by Call IGRF c10 kizami H is shown in Adam_update c10 Result of ray tracing (OUT1) is written in fort.21 c10 Instantaneous H is written in fort.25 c10======== initial conditions ================================== c fq00: wave frequency in kHz c rr00: starting radius distance from the earth center in km c th00: geomagnetic colatitude in deg c ph00: geomagnetic longitude in deg c es00: angle between the initial wave normal and geomagnetic c meiridian plane, in deg (Epsilon) c dl00: angle between the projected direction of the initial c wave normal onto the geomagnetic meridial plane and c the radial direction, in deg (Delta) rr00 = 6370.0+91.0d0 ph00 = 0.0d0 ameg = 1000.0d0 fq00 = 0.003d0*ameg dl00 = 0.0d0 es00 = 0.0d0 print*,'VLF_DE_dipole_F.f',', rr00=', rr00,', fq00=',fq00, &', ph00=', ph00,', es00=', es00,', dl00=',dl00 DO 100 I=1,5 WRITE(21,*) write(25,*) th00 = dble(i-1)*1.0d0+40.0d0 print*, 'th00=', th00 CALL INIT(fq00,rr00,th00,ph00,es00,dl00) CALL OUT0 CALL ADAMS 100 CONTINUE STOP END ************************************************************************ SUBROUTINE INIT(fq00,rr00,th00,ph00,es00,dl00) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QINIT /R0,T0,P0,E0,D0,F0 save inited data inited / 0 / PAI = dASIN(1.0D0)*2.0D0 RAD = PAI/180.0D0 DEG = 1.0D0/RAD RE = 6370.0D0 RREF = 7370.0D0 RREFP = RREF RBTM = 9370.0D0 CVLCTY = 3.00D5 XME = 9.109534D-31 XMP = 1.6726485D-27 QE = 1.6021892D-19 EP = 8.854185D-12 XK = 1.380662D-23 GM = 3.98603D14 PF0(1) = QE**2/4.0D0/PAI**2/EP/ XME PF0(2) = QE**2/4.0D0/PAI**2/EP/ XMP PF0(3) = QE**2/4.0D0/PAI**2/EP/(XMP*4.0D0) PF0(4) = QE**2/4.0D0/PAI**2/EP/(XMP*16.0D0) GF0(1) =-QE/2.0D0/PAI/ XME *1.0D-3 GF0(2) = QE/2.0D0/PAI/ XMP *1.0D-3 GF0(3) = QE/2.0D0/PAI/(XMP* 4.0D0)*1.0D-3 GF0(4) = QE/2.0D0/PAI/(XMP*16.0D0)*1.0D-3 c XN0 = 32000.0D0 XN0 = 50000.0D0 c XN0 = 4000.0D0 THERM = 1000.0D0 SHINV(2) = XMP* (GM/RREF**2)/XK/THERM*1.0D-3 SHINV(3) = XMP* 4.0D0*(GM/RREF**2)/XK/THERM*1.0D-3 SHINV(4) = XMP*16.0D0*(GM/RREF**2)/XK/THERM*1.0D-3 SHINVP = SHINV(2) ETA(2) = 0.152D0 ETA(3) = 0.82D0 ETA(4) = 0.025D0 c Following RRLI: location of the lower edge of the ionosphere RRLI = 6370.0d0 + 90.0d0 c The following parameters areused in subroutine DENSLI HWLI = 140.0D0 XLPP = 4.0D0 ERPP = 0.2D0 DRR1 = 5.0D0 DTH1 = 0.05D0 DPH1 = 0.2D0 DRR2 = 2.5D0 DTH2 = 0.01D0 DPH2 = 0.1D0 c Following RRMAX & RRMIN: Upper & lower limit of ray tracing(Kmj c Following ERRR and ERRA; Initial values for the limiting error c and Hmin: Limit of integration increment RRMAX = 4.0*RE RRMIN = 6370.0d0 ERRR = 1.0D-5 ERRA = 1.0D-5 HMIN = 1.0D-10 c Folowing DISOUT: Every this count, the output data to be printed DISOUT = 50 c Following NLPMAX: Maximum number of integration loops NLPMAX = 30000 INDOUT = 20 INDMSH = 51 IEOF = -1 C...FREQ.(KHZ),GC-DIST.(KM),GM-COLAT.,GM-LONG.,EPS,DEL(DEG) F0 = fq00 R0 = rr00 T0 = th00 P0 = ph00 E0 = es00 D0 = dl00 RETURN END ************************************************************************ SUBROUTINE ADAMS ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) *** CONSTANT COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY *** INPUT PARAMETER (NSPEC=4) COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QINIT /R0,T0,P0,E0,D0,F0 COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY *** OUTPUT COMMON /QCOND /L,LL,NN,NLOOP,A *** Variables below are only used in the ADAMS method. COMMON /QADAM /Y(8),P(8),F(8), & H,HOLD,WT(8),PHI(8,14),PSIOLD(12), & PSI(12),BETA(12),SIGMA(13),G(13), & ERKM1,ERK,K,KOLD, & KNEW,KS,IPHASE,IFAIL REAL*8 SUML(8) CALL OUT0 *** Initialization of Y Y(1) = R0 Y(2) = T0 Y(3) = P0 EPS = E0 DEL = D0 FKC = F0 Y(4) = COS(EPS*RAD)*COS(DEL*RAD) Y(5) = COS(EPS*RAD)*SIN(DEL*RAD) Y(6) = SIN(EPS*RAD) Y(7) = 0.0D0 Y(8) = 0.0D0 L = 0 LL = 0 NN = 0 NLOOP = 0 ERR0 = MAX(ERRR,ERRA) ERRR = ERRR/ERR0 ERRA = ERRA/ERR0 * Derivative for Y CALL FUNCT(Y,F,1) yp20 = F(2) do L=1,8 WT(L)=ERRR*dABS(Y(L)+ERRA) SUML(L) = (F(L)/WT(L))**2 enddo SUM = 0.0D0 do L=1,8 SUM = SUM+SUML(L) enddo SUM = DSQRT(SUM) H = 0.00001D0*dSQRT(ERR0/SUM) HOLD = 0.0D0 K = 1 KOLD = 0 IPHASE = 1 IFAIL = 0 do L=1,8 PHI(L,1) = F(L) PHI(L,2) = 0.0D0 enddo CALL OUT1(Y) DO NLOOP=0,NLPMAX-1 * Adams-Bashforth Predictor CALL ADAM_B(NLOOP) * Adams-Moulton Corrector CALL ADAM_M(ERR0,NLOOP) *VOCL LOOP,SCALAR IF(LL.EQ.0) THEN IF(IFAIL.NE.0) THEN L=11 ELSE IF(L.GE.21) THEN LL=L ELSE IF(H.LT.HMIN) THEN L =23 LL=23 ELSE IF(Y(1).LT.RRMIN) THEN L =24 LL=24 ELSE IF(Y(1).GT.RRMAX) THEN L =25 LL=25 C ELSE IF((T0-90.0)*F(2).GT.0.0D0) THEN C L =26 C LL=26 C ELSE IF(yp20*yp(2).lt.0.0D0) THEN C L =27 C LL=27 ELSE L =0 NN=NN+1 ENDIF ENDIF ENDIF CALL OUT1(Y) IF (LL.NE.0) THEN write(*,*) "Stop Condition: ", LL GOTO 999 ENDIF do I=1,7 WT(I)=ERRR*ABS(Y(I)+ERRA) enddo ENDDO 999 CONTINUE write(*,*) 'Number of LOOP=',NLOOP,', Number of calculation=',A write(*,*) 'Last interval=',H,', IFAIL=',IFAIL write(21,*) write(25,*) RETURN END ************************************************************************ SUBROUTINE ADAM_B(NLOOP) * Calculation of Adams-Bashforth Predictor ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QADAM /Y(8),P(8),F(8), & H,HOLD,WT(8),PHI(8,14),PSIOLD(12), & PSI(12),BETA(12),SIGMA(13),G(13), & ERKM1,ERK,K,KOLD, & KNEW,KS,IPHASE,IFAIL real*8 ALPHA(12),V(12),W(12) save ALPHA,V,W G(1) = 1.0D0 G(2) = 0.5D0 SIGMA(1) = 1.0D0 * KS is the maximum order * In case step size is changed, KS->0. IF (H .NE.HOLD) KS = 0 IF (KS.LE.KOLD) KS = KS + 1 IF (K.GE.KS) THEN *** PSI do I=KS,K-1 PSIOLD(I) = PSI(I) enddo PSI(KS) = H*KS do I=KS+1,K PSI(I) = PSIOLD(I-1)+H enddo *** ALPHA ALPHA(KS) = 1.0D0/DBLE(KS) do I=KS+1,K ALPHA(I) = H/(PSIOLD(I-1)+H) enddo *** BETA BETA(KS) = 1.0D0 do I=KS+1,K BETA(I) = BETA(I-1)*PSI(I-1)/PSIOLD(I-1) enddo *** SIGMA SIGMA(KS+1) = 1.0D0 do I=KS+1,K SIGMA(I+1) = I*ALPHA(I)*SIGMA(I) enddo *** G (Integral of C) if (KS.EQ.0) write(0,*) 'Error! KS=0' if (KS.EQ.1) then do I=1,K V(I) = 1.0D0/DBLE(I*(I+1)) W(I) = V(I) enddo else if (K.gt.KOLD) then V(K) = 1.0D0/DBLE(K*(K+1)) do J=1,KS-2 V(K-J) = V(K-J)-ALPHA(J+1)*V(K-J+1) enddo endif do IQ=1,K+1-KS V(IQ) = V(IQ)-ALPHA(KS)*V(IQ+1) W(IQ) = V(IQ) enddo G(KS+1) = W(1) endif do I=KS+1,K do IQ=1,K+1-I W(IQ) = W(IQ)-ALPHA(I)*W(IQ+1) enddo G(I+1) = W(1) enddo ENDIF *** PREDICTOR do L=1,8 do I=KS+1,K PHI(L,I) = BETA(I)*PHI(L,I) enddo PHI(L,K+2) = PHI(L,K+1) PHI(L,K+1) = 0.0D0 SUM = 0.0D0 do I=K,1, -1 SUM = SUM + G(I)*PHI(L,I) PHI(L,I) = PHI(L,I) + PHI(L,I+1) enddo P(L) = Y(L) + H * SUM enddo *** Derivative CALL FUNCT(P,F,0) RETURN END ************************************************************************ SUBROUTINE ADAM_M(ERR0,NLOOP) * Calculation of Adams-Moulton Corrector ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QADAM /Y(8),P(8),F(8), & H,HOLD,WT(8),PHI(8,14),PSIOLD(12), & PSI(12),BETA(12),SIGMA(13),G(13), & ERKM1,ERK,K,KOLD, & KNEW,KS,IPHASE,IFAIL DIMENSION GSTR(13) DATA GSTR/0.500D0 ,0.0833D0 ,0.0417D0 ,0.0264D0 ,0.0188D0 , & 0.0143D0 ,0.0114D0 ,0.00936D0,0.00789D0,0.00679D0, & 0.00592D0,0.00524D0,0.00468D0/ KNEW = K *** ERK, ERKM1, ERKM2 ERK = 0.0D0 ERKM1 = 0.0D0 ERKM2 = 0.0D0 do L=1,8 TEMP1 = 1.0D0/WT(L) TEMP2 = F(L)-PHI(L,1) ERK = ERK+(TEMP2 *TEMP1)**2 IF (K.EQ.2) THEN ERKM1 = ERKM1+((PHI(L,K )+TEMP2)*TEMP1)**2 ELSE IF (K.GT.2) THEN ERKM1 = ERKM1+((PHI(L,K )+TEMP2)*TEMP1)**2 ERKM2 = ERKM2+((PHI(L,K-1)+TEMP2)*TEMP1)**2 ENDIF enddo TEMP = DABS(H)*DSQRT(ERK) ERK = TEMP*SIGMA(K+1)*GSTR(K) IF (k.EQ.2) THEN ERKM1 = DABS(H)*SIGMA(K )*GSTR(K-1)*DSQRT(ERKM1) IF(ERKM1.LE.0.5D0*ERK) KNEW = K-1 ELSE IF (K.GT.2) THEN ERKM1 = DABS(H)*SIGMA(K )*GSTR(K-1)*DSQRT(ERKM1) ERKM2 = DABS(H)*SIGMA(K-1)*GSTR(K-2)*DSQRT(ERKM2) IF(DMAX1(ERKM1,ERKM2).LE.ERK) KNEW = K-1 ENDIF *** Validity check ERR = TEMP*(G(K)-G(K+1)) IF(ERR.LE.ERR0) THEN CALL ADAM_UPDATE(ERR0) ELSE CALL ADAM_RETRY(ERR0) ENDIF RETURN END ************************************************************************ SUBROUTINE ADAM_UPDATE(ERR0) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QADAM /Y(8),P(8),F(8), & H,HOLD,WT(8),PHI(8,14),PSIOLD(12), & PSI(12),BETA(12),SIGMA(13),G(13), & ERKM1,ERK,K,KOLD, & KNEW,KS,IPHASE,IFAIL DIMENSION GSTR(13),TWO(13) DATA GSTR/0.500D0 ,0.0833D0 ,0.0417D0 ,0.0264D0 ,0.0188D0 , & 0.0143D0 ,0.0114D0 ,0.00936D0,0.00789D0,0.00679D0, & 0.00592D0,0.00524D0,0.00468D0/ DATA TWO/ 2.0D0, 4.0D0, 8.0D0, 16.0D0, 32.0D0, & 64.0D0, 128.0D0, 256.0D0, 512.0D0,1024.0D0, & 2048.0D0,4096.0D0,8192.0D0/ do L=1,8 Y(L) = P(L) + H * G(K+1) * (F(L)-PHI(L,1)) enddo IFAIL = 0 KOLD = K HOLD = H *** Derivative CALL FUNCT(Y,F,1) do L=1,8 PHI(L,K+1) = F(L)-PHI(L,1) PHI(L,K+2) = PHI(L,K+1)-PHI(L,K+2) enddo do I=1,K do L=1,8 PHI(L,I) = PHI(L,I)+PHI(L,K+1) enddo enddo ERKP1 = 0.0D0 IF(KNEW.EQ.K-1.OR.K.EQ.12) IPHASE=0 IF(IPHASE.EQ.1) THEN K = K+1 ERK = ERKP1 GOTO 500 ENDIF IF(KNEW.EQ.K-1) THEN K = K-1 ERK = ERKM1 GOTO 500 ENDIF IF(K+1.LE.KS) THEN do L=1,8 ERKP1 = ERKP1+(PHI(L,K+2)/WT(L))**2 enddo ERKP1 = DABS(H)*GSTR(K+1)*DSQRT(ERKP1) IF(K.LE.1) THEN IF(ERKP1.LT.0.5D0*ERK) THEN K = K+1 ERK = ERKP1 ENDIF GOTO 500 ENDIF IF(ERKM1.LE.DMIN1(ERK,ERKP1)) THEN K = K-1 ERK = ERKM1 ELSE IF(ERKP1.LT.ERK .AND. K.NE.12) THEN K = K+1 ERK = ERKP1 ENDIF ENDIF 500 CONTINUE *** H ->> IF (ERK*TWO(K+1).LT.0.5D0*ERR0 .OR. IPHASE.EQ.1) THEN HNEW = 1.01D0*H GOTO 1000 ENDIF *** <<- H IF (ERK.GT.0.5D0*ERR0) THEN R = (0.5D0*ERR0/ERK)**(1.0D0/DBLE(K+1)) HNEW = DMIN1(0.9D0,DMAX1(0.5D0,R)) * H GOTO 1000 ENDIF *** -> H <- HNEW = H 1000 CONTINUE H = HNEW WRITE(25,800) Y(1)-6370.0, HNEW, ERK 800 FORMAT(F18.2, F18.6, F18.8) RETURN END ************************************************************************ SUBROUTINE ADAM_RETRY(ERR0) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) COMMON /QADAM /Y(8),P(8),F(8), & H,HOLD,WT(8),PHI(8,14),PSIOLD(12), & PSI(12),BETA(12),SIGMA(13),G(13), & ERKM1,ERK,K,KOLD, & KNEW,KS,IPHASE,IFAIL IPHASE = 0 do I=1,K do L=1,8 PHI(L,I) = (PHI(L,I)-PHI(L,I+1))/BETA(I) enddo enddo do I = 2,K PSI(I-1) = PSI(I)-H enddo IFAIL = IFAIL+1 * Euler method IF (IFAIL.GE.3) KNEW = 1 K = KNEW IF (IFAIL.GT.3 .AND. 0.5D0*ERR0.LT.0.25D0*ERK) THEN H = DSQRT(0.5D0*ERR0/ERK) * H ELSE H = 0.5D0 * H ENDIF *** Derivative CALL FUNCT(Y,F,1) RETURN END ************************************************************************ SUBROUTINE OUT0 ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QINIT /R0,T0,P0,E0,D0,F0 COMMON /QOUT /XX,YY,ZZ,OO,PP COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0 XX = R0*dSIN(T0*RAD)*dCOS(P0*RAD) YY = R0*dSIN(T0*RAD)*dSIN(P0*RAD) ZZ = R0*dCOS(T0*RAD) OO = 0.0D0 PP = 0.0D0 RETURN END ************************************************************************ SUBROUTINE OUT1(Y) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) DIMENSION Y(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY COMMON /QCOND /L,LL,NN,NLOOP,A COMMON /QOUT /XX,YY,ZZ,OO,PP COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0 COMMON /QREF /PMU,PMUS,PMU0, & PSI , & COPSI ,SIPSI , & COPSIS,SIPSIS, & PX(NSPEC),PY(NSPEC),PZ(NSPEC), & PA1,PA2,PA3, & PK1,PK2,PK4,PK5, & PA ,PB ,PC ,PD COMMON /QDENS /XNS (NSPEC), & DNSDRR(NSPEC), & DNSDTH(NSPEC), & DNSDPH(NSPEC) real*8 RR,TH,PH,DL RR = Y(1) TH = Y(2) PH = Y(3) DL = DEL XXX = RR*dSIN(TH*RAD)*dCOS(PH*RAD) YYY = RR*dSIN(TH*RAD)*dSIN(PH*RAD) ZZZ = RR*dCOS(TH*RAD) DDD = dSQRT((XXX-XX)**2+(YYY-YY)**2+(ZZZ-ZZ)**2) OO = OO+DDD XX = XXX YY = YYY ZZ = ZZZ IF(L.NE.0 .OR. LL.NE.0) GOTO 1020 IF(OO.GE.PP) THEN write(21,1015) & REAL(90.-TH), & REAL(RR-6370.), & REAL(PH), & REAL(DL), & PSI, & 8.98*dsqrt(XNS(1)), & REAL(PMU), & REAL(Y(7)) c & Above time(ms)is the propagation delay c & ((dcos((90.0-TH)*RAD)**2+0.5)/1.5)* c & (XNS(1)+XNS(2)+XNS(3)+XNS(4)) 1015 FORMAT(F8.2,F9.2,F9.5,F7.1,F10.1,F12.3,F11.2,F12.4) PP=PP+DISOUT A=A+1 END IF 1020 CONTINUE RETURN END ************************************************************************ SUBROUTINE FUNCT(PATH,DPATH,IMDFY) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8),dpath(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QCOND /L,LL,NN,NLOOP,a COMMON /QREF /PMU,PMUS,PMU0, & PSI , & COPSI ,SIPSI , & COPSIS,SIPSIS, & PX(NSPEC),PY(NSPEC),PZ(NSPEC), & PA1,PA2,PA3, & PK1,PK2,PK4,PK5, & PA ,PB ,PC ,PD COMMON /QIGRF /YR ,YT ,YP , & DYRDRR,DYTDRR,DYPDRR, & DYRDTH,DYTDTH,DYPDTH, & DYRDPH,DYTDPH,DYPDPH, & BB , & DBBDRR,DBBDTH,DBBDPH COMMON /QDENS /XNS (NSPEC), & DNSDRR(NSPEC), & DNSDTH(NSPEC), & DNSDPH(NSPEC) real*8 DMUDRR,DMUDTH,DMUDPH, & DK2DF ,DA3DF ,DENOM CALL REF(PATH) IF(IMDFY.EQ.1) THEN PATH(4) = PATH(4)*PMU/PMU0 PATH(5) = PATH(5)*PMU/PMU0 PATH(6) = PATH(6)*PMU/PMU0 DEL = dATAN2(PATH(5),PATH(4)) EPS = dATAN (PATH(6)/dSQRT(PATH(4)**2+PATH(5)**2)) DEL = DEL*DEG EPS = EPS*DEG END IF DMUDRR = 0.0D0 DMUDTH = 0.0D0 DMUDPH = 0.0D0 DK2DF = 0.0D0 DA3DF = 0.0D0 DENOM = (4.0D0*PA*PMUS+2.0D0*PB)*PMU DO 20 I=1,NSPEC DADX = - COPSIS DBDX = -(PK1/PZ(I)-PK2)*(1.0D0+COPSIS) DCDX = - PK5 PK2A3Y = 2.0D0*(PK2-PA3*PY(I)) DADX = DADX+ SIPSIS/PZ(I) DBDX = DBDX- PK2A3Y*SIPSIS/PZ(I) DCDX = DCDX+PK1*PK2A3Y /PZ(I) DMUDX = -(DADX*PMUS**2+DBDX*PMUS+DCDX)/DENOM DK2DY = -2.0D0 *PX(I)* PY(I) /PZ(I)**2 DK5DY = 2.0D0*PA3*PX(I)*(PZ(I)+2.0D0)/PZ(I)**2 DK5DY = DK5DY+2.0D0*DK2DY*PK2 DADY = DK2DY*SIPSIS DBDY = -PK1*DK2DY*(1.0D0+COPSIS)-DK5DY*SIPSIS DCDY = PK1*DK5DY DMUDY = -(DADY*PMUS**2+DBDY*PMUS+DCDY)/DENOM DXDRR = PX(I)*DNSDRR(I) DXDTH = PX(I)*DNSDTH(I) DXDPH = PX(I)*DNSDPH(I) DYDRR = PY(I)*DBBDRR DYDTH = PY(I)*DBBDTH DYDPH = PY(I)*DBBDPH DMUDRR = DMUDRR+DMUDX*DXDRR+DMUDY*DYDRR DMUDTH = DMUDTH+DMUDX*DXDTH+DMUDY*DYDTH DMUDPH = DMUDPH+DMUDX*DXDPH+DMUDY*DYDPH DK2DFI = 2.0D0 *PX(I) /PZ(I)**2 DA3DFI = (2.0D0-PZ(I))*PX(I)*PY(I)/PZ(I)**2 DK2DF = DK2DF+DK2DFI/FKC DA3DF = DA3DF+DA3DFI/FKC 20 CONTINUE DADPS = PA1+PA2 DBDPS = PK4-PK5 DMUDPS = (DADPS*PMUS+DBDPS)*PMUS DMUDPS = -2.0D0*COPSI*DMUDPS/DENOM DPSDRR = DYRDRR*PATH(4) & +DYTDRR*PATH(5) & +DYPDRR*PATH(6) DPSDTH = DYRDTH*PATH(4) & +DYTDTH*PATH(5) & +DYPDTH*PATH(6) DPSDPH = DYRDPH*PATH(4) & +DYTDPH*PATH(5) & +DYPDPH*PATH(6) DPSDRR = -DPSDRR/PMU DPSDTH = -DPSDTH/PMU DPSDPH = -DPSDPH/PMU if (XNS(1).NE.0.0d0) then DMUDRR = (DMUDRR+DMUDPS*DPSDRR)/PMU DMUDTH = (DMUDTH+DMUDPS*DPSDTH)/PMU DMUDPH = (DMUDPH+DMUDPS*DPSDPH)/PMU DMUDVR = DMUDPS*(PATH(4)*COPSI/PMU-YR) DMUDVT = DMUDPS*(PATH(5)*COPSI/PMU-YT) DMUDVP = DMUDPS*(PATH(6)*COPSI/PMU-YP) else DMUDRR = 0.0d0 DMUDTH = 0.0d0 DMUDPH = 0.0d0 DMUDVR = 0.0d0 DMUDVT = 0.0d0 DMUDVP = 0.0d0 endif DK1DF = 2.0D0* PA1/FKC DK4DF = PK1*DK2DF+PK2*DK1DF DK5DF = 2.0D0*(PK2*DK2DF-PA3*DA3DF) DCDF = PK1*DK5DF +PK5*DK1DF DADF = DK1DF* COPSIS +DK2DF*SIPSIS DBDF = -DK4DF*(1.0D0+COPSIS)-DK5DF *SIPSIS DMUDF = -((DADF*PMUS+DBDF)*PMUS+DCDF)/DENOM if (XNS(1).EQ.0.0d0) dmudf=0.0d0 GMU = PMU+FKC*DMUDF IF (PATH(1).LT.6460.0) THEN DMUDVR = 0.0d0 DMUDVR = 0.0d0 DMUDVP = 0.0d0 DMUDRR = 0.0d0 DMUDTH = 0.0d0 DMUDPH = 0.0d0 END IF DPATH(1) = PATH(4)-DMUDVR DPATH(2) = PATH(5)-DMUDVT DPATH(3) = PATH(6)-DMUDVP DPATH(1) = DPATH(1)/PMUS DPATH(2) = DPATH(2)/PMUS/PATH(1) DPATH(3) = DPATH(3)/PMUS/PATH(1)/dSIN(PATH(2)*RAD) DPATH(4) = DMUDRR & +PATH(5)*DPATH(2) & +PATH(6)*DPATH(3)*dSIN (PATH(2)*RAD) DPATH(5) = DMUDTH/PATH(1) & -PATH(5)*DPATH(1)/PATH(1) & +PATH(6)*DPATH(3)*dCOS (PATH(2)*RAD) DPATH(6) = DMUDPH/PATH(1)/dSIN(PATH(2)*RAD) & -PATH(6)*DPATH(1)/PATH(1) & -PATH(6)*DPATH(2)/dTAN(PATH(2)*RAD) DPATH(2) = DPATH(2)*DEG DPATH(3) = DPATH(3)*DEG DPATH(7) = 1.0D0/PMU/CVLCTY DPATH(8) = 1.0D0/CVLCTY if (XNS(1).eq.0.0d0) alf=0.0d0 ALF = -ATAN(DMUDPS*SIPSI/PMU)*DEG RAY = PSI+ALF RETURN END ************************************************************************ SUBROUTINE REF(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QCOND /L,LL,NN,NLOOP,a COMMON /QREF /PMU,PMUS,PMU0, & PSI , & COPSI ,SIPSI , & COPSIS,SIPSIS, & PX(NSPEC),PY(NSPEC),PZ(NSPEC), & PA1,PA2,PA3, & PK1,PK2,PK4,PK5, & PA ,PB ,PC ,PD COMMON /QIGRF /YR ,YT ,YP , & DYRDRR,DYTDRR,DYPDRR, & DYRDTH,DYTDTH,DYPDTH, & DYRDPH,DYTDPH,DYPDPH, & FH , & DBBDRR,DBBDTH,DBBDPH COMMON /QDENS /XNS (NSPEC), & DNSDRR(NSPEC), & DNSDTH(NSPEC), & DNSDPH(NSPEC) CALL DIPOLE(PATH) CALL DENS(PATH) DO 10 I=1,NSPEC PX(I) = PF0(I)*XNS(I)/FKC**2 PY(I) = GF0(I)*FH /FKC PZ(I) = PY(I)**2-1.0D0 10 CONTINUE PMU0 = dSQRT(PATH(4)**2+PATH(5)**2+PATH(6)**2) COPSI = (YR*PATH(4) & +YT*PATH(5) & +YP*PATH(6))/PMU0 COPSIS = COPSI**2 SIPSIS = dABS(1.0D0-COPSIS) SIPSI = dSQRT(SIPSIS) if (dabs(copsi).gt.1.0) copsi=1.0d0 PSI = dACOS(COPSI)*DEG PA1 = 0.0D0 PA2 = 0.0D0 PA3 = 0.0D0 DO 20 I=1,NSPEC PA1 = PA1+PX(I) PA2 = PA2+PX(I) /PZ(I) PA3 = PA3+PX(I)*PY(I)/PZ(I) 20 CONTINUE PK1 = 1.0D0-PA1 PK2 = 1.0D0+PA2 PK4 = PK1*PK2 PK5 = PK2**2-PA3**2 PA = PK1*COPSIS+PK2*SIPSIS PB =-PK4*(COPSIS+1.0D0)-PK5*SIPSIS PC = PK1*PK5 PD = DABS(PB**2-4.0D0*PA*PC) if (pd.lt.0.0d0) go to 1040 PMUS= (-PB-SQRT(PD))/2.0/PA c Above is Extraordinary mode(whistler mode) if (pmus.lt.0.0d0) go to 1050 PMU = dSQRT(PMUS) 1040 CONTINUE IF (PATH(1).LT.6370.0d0+90.0d0) THEN PMU=1.0d0 gmu=1.0d0 END IF 1050 CONTINUE RETURN END ************************************************************************ SUBROUTINE DENS(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QDENS /XNS (NSPEC), & DNSDRR(NSPEC), & DNSDTH(NSPEC), & DNSDPH(NSPEC) COMMON /QDNSDE/XNDE (NSPEC), & DNDDRR(NSPEC) COMMON /QDNSLI/XNLI ,DNLDRR COMMON /QDNSAO/XNAO ,DNADRR,DNADTH,DNADPH COMMON /QDNSOP/XNOP ,DNODRR,DNODTH,DNODPH CALL DENSDE(PATH) CALL DENSLI(PATH) CALL DENSAO(PATH) CALL DENSOP(PATH) DO 10 IS=1,NSPEC XNS(IS) = XNDE(IS)*XNLI*XNAO*XNOP *... 1/N * dn/dr DNSDRR(IS) = DNDDRR(IS)+DNLDRR+DNADRR+DNODRR DNSDTH(IS) = DNADTH+DNODTH DNSDPH(IS) = DNADPH+DNODPH 10 CONTINUE RETURN END ************************************************************************ SUBROUTINE DENSAO(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QIGRF /YR ,YT ,YP , & DYRDRR,DYTDRR,DYPDRR, & DYRDTH,DYTDTH,DYPDTH, & DYRDPH,DYTDPH,DYPDPH, & BB , & DBBDRR,DBBDTH,DBBDPH COMMON /QDNSAO/XNAO ,DNADRR,DNADTH,DNADPH COMMON /QFLPRM/TH0 , & DT0DRR,DT0DTH,DT0DPH, & PH0 , & DP0DRR,DP0DTH,DP0DPH, & BB0 , & DB0DRR,DB0DTH,DB0DPH, & XLV , & DLVDRR,DLVDTH,DLVDPH XNAO = 1.0D0 DNADRR = 0.0D0 DNADTH = 0.0D0 DNADPH = 0.0D0 RETURN END ************************************************************************ SUBROUTINE DENSDE(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QDNSDE/XNDE (NSPEC), & DNDDRR(NSPEC) GPH = RREF*(PATH(1)-RREF)/PATH(1) DGPHDR = (RREF/PATH(1))**2 ETEXP2 = ETA(2)*EXP(-GPH*SHINV(2)) ETEXP3 = ETA(3)*EXP(-GPH*SHINV(3)) ETEXP4 = ETA(4)*EXP(-GPH*SHINV(4)) SUM1 = ETEXP2 & + ETEXP3 & + ETEXP4 SUM2 = ETEXP2*SHINV(2) & + ETEXP3*SHINV(3) & + ETEXP4*SHINV(4) XNDE (1) = dSQRT(SUM1) XNDE (2) = ETEXP2/XNDE(1) XNDE (3) = ETEXP3/XNDE(1) XNDE (4) = ETEXP4/XNDE(1) DNDDRR(1) =-DGPHDR*SUM2/SUM1/2.0D0 DNDDRR(2) =-DGPHDR*SHINV(2)-DNDDRR(1) DNDDRR(3) =-DGPHDR*SHINV(3)-DNDDRR(1) DNDDRR(4) =-DGPHDR*SHINV(4)-DNDDRR(1) RETURN END ************************************************************************ SUBROUTINE DENSLI(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QDNSLI/XNLI, DNLDRR IF ((PATH(1)-RRLI).GE.0) THEN ALTNRM = (PATH(1)-RRLI)/HWLI XNLI = 1.0D0-dEXP(-ALTNRM**2) DNLDRR = dEXP(-ALTNRM**2)*ALTNRM*2.0D0/HWLI/XNLI ELSE XNLI=0.0 DNLDRR=0.0 END IF RETURN END ************************************************************************ SUBROUTINE DENSOP(PATH) ************************************************************************ IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NSPEC=4) real*8 path(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY COMMON /QPARAM/RREF,RBTM,RREFP, & PF0(NSPEC),GF0(NSPEC), & XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP, & RRLI,HWLI, & XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2, & RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT, & NPATH,NLPMAX,INDOUT,INDMSH,IEOF, & IYEAR,IMNTH,IDATE COMMON /QDNSOP/XNOP ,DNODRR,DNODTH,DNODPH XNOP = XN0 dNODRR = 0.0D0 DNODTH = 0.0d0 dnodph = 0.0d0 RETURN END ********************************************************************** SUBROUTINE DIPOLE(PATH) ********************************************************************** IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) DIMENSION PATH(7) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QIGRF /YR ,YT ,YP , & DYRDRR,DYTDRR,DYPDRR, & DYRDTH,DYTDTH,DYPDTH, & DYRDPH,DYTDPH,DYPDPH, & BB , & DBBDRR,DBBDTH,DBBDPH DATA GMDM/8.07D6/ RR = PATH(1) TH = PATH(2) PH = PATH(3) COSTH = COS(TH*RAD) SINTH = SIN(TH*RAD) COSTHS = 1.0D0+3.0D0*COSTH**2 YR =-2.0D0*COSTH/SQRT(COSTHS) YT =- SINTH/SQRT(COSTHS) YP = 0.0D0 DYRDRR = 0.0D0 DYTDRR = 0.0D0 DYPDRR = 0.0D0 DYRDTH =-2.0D0/COSTHS*YT DYTDTH = 2.0D0/COSTHS*YR DYPDTH = 0.0D0 DYRDPH = 0.0D0 DYTDPH = 0.0D0 DYPDPH = 0.0D0 BB = GMDM*SQRT(COSTHS)/RR**3 DBBDRR =-3.0D0/RR DBBDTH =-3.0D0*SINTH*COSTH/COSTHS DBBDPH = 0.0D0 RETURN END