IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) PARAMETER (NX=121,NY=41,NZ=431) DIMENSION ENE(NX,NY,NZ) COMMON/ ENECOM/ENE c10================================================================ c10 RAY TRACING PROGRAM FOR HF WAVE O- and X-mode IN IRI-MODEL c10 Dipole model is selected in REF by Call dipole c10 IGRF model is selected in REF by Call IGRF c10 O-mode is for MODEPR= 0 and X-mode is for MODEPR = 1 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 Maximum number of hops is assigned by NOHOPS here c10======== initial conditions ================================== C IRI data are processed from IRI using e.g. IRI_UT03_5, C in which ionospheric data at UT03hour for a certain altitude, C latitude and longitude range assigned in the iri_main.f. C Processed data are stored in 'iri_UT03_5.dat' which is taken C to the present ray tracing program, and rearranged as ENE(I,J,K). OPEN(20,FILE='iri_UT03_5.dat') DO 95 J=1,NY DO 96 I=1,NX READ(20,*) (ENE(I,J,K), K=1,NZ) 96 CONTINUE 95 CONTINUE CLOSE(20) XW=1.0 IXW=1 YW=1.0 IYW=1 ZW=1.0 IZW=1 J=21 DO 98 I=1,NX,10 DO 99 K=1,NZ WRITE(27,105) I,J,K,30+(I-1)*IXW,(70.+(K-1)*ZW), &ENE(I,J,K) 105 FORMAT(I5, I5, I5, I5, F10.0, F12.0, F12.7) 99 CONTINUE WRITE(27,*) 98 CONTINUE WRITE(27,*) PAI = dASIN(1.0D0)*2.0D0 RAD = PAI/180.0D0 DEG = 1.0D0/RAD c Mode of propagation (MODEPR): 0 for O-mode, 1 for X-mode c fq00:(kHz) wave frequency in kHz c rr00:(km) starting position in radius distance(km) c th00: (deg) initial colatitude c ph00: (deg) initial longitude c bt00: (deg) initial eastward deflection angle of c K-vector from the geomagnetic meridian plane c al00: (deg) initial zenith angle of K-vector (alfa) c es00: (deg) angle between initial K-vector and the c projection of initial K-vector onto the c geomagnetic meridian plan (epsilon) c dl00: (deg) zenith angle of the initial K-vector (delta) c When bt00=0, then dl00=al00 c MULHOPM: Number of Hops c NOHOPS: Muximum number of hops to be pre-assigned here MODEPR = 0 NOHOPS = 3 rr00= 6370.0d0 ameg = 1000.0d0 fq00 = 8.0d0*ameg th00 = 55.0 ph00 = 135.0 bt00 = 0.0 if (MODEPR.EQ.0) then print*,'HF_IRI_dipole_O-MODE' print*,' fq00=', fq00,', rr00=',rr00,', th00=', th00, &', ph00=', ph00,', bt00=',bt00 else print*,'HF_IRI_dipole_X-MODE' print*,' fq00=', fq00,', rr00=',rr00,', th00=', th00, &', ph00=', ph00,', bt00=',bt00 end if DO 100 I=1,3 al00 = (i-1)*10.0+50.0 es00 = DASIN(sin(al00*RAD)*sin(bt00*RAD))*DEG dl00 = DASIN(sin(al00*RAD)*cos(bt00*RAD)/cos(es00*RAD))*DEG print*,' al00=', al00,', es00=', es00,', dl00=',dl00 WRITE(21,*) WRITE(25,*) CALL INIT(fq00,rr00,th00,ph00,es00,dl00,bt00,MODEPR,NOHOPS) CALL OUT0 CALL ADAMS 100 CONTINUE STOP END ************************************************************************ SUBROUTINE INIT(fq00,rr00,th00,ph00,es00,dl00,bt00,MODEPR,NOHOPS) ************************************************************************ 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,Bt0,Ph0,Th0,NHOPS,MODPROP save inited data inited / 0 / MODPROP = MODEPR NHOPS = NOHOPS PAI = dASIN(1.0D0)*2.0D0 RAD = PAI/180.0D0 DEG = 1.0D0/RAD RE = 6370.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 GF0(1) =-QE/2.0D0/PAI/ XME *1.0D-3 c RRMAX: upper limit altitude for ray tracing in km RRMIN = 6370.0d0 c RRMAX: upper limit altitude for ray tracing in km RRMAX = RE + 500.0 c Initial values for the limiting error; ERRR and ERRA, c and limit of integration increment(Hmin) ERRR = 1.0D-5 ERRA = 1.0D-5 HMIN = 1.0D-10 DISOUT = 3 c DISOUT: Every this count, the output data to be printed NLPMAX = 60000 INDOUT = 20 INDMSH = 51 IEOF = -1 C...FREQ.(KHZ),GC-DIST.(KM),GM-COLAT.,GM-LONG.,EPS,DEL(DEG) F0 = fq00 Bt0= bt00 Ph0= ph00 Th0= th00 R0 = rr00 T0 = th00 P0 = ph00 E0 = es00 D0 = dl00 if(inited.eq.0) then inited = 1 endif 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,Bt0,Ph0,Th0,NHOPS,MODPROP 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 MULHOP = 0 *** 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) *** Main Loop 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.6370.0d0) THEN Y(1)=6371.0 Y(2)=Y(2) Y(3)=Y(3) Y(4)=-Y(4) Y(5)=Y(5) Y(6)=Y(6) MULHOP=MULHOP+1 IF ((MULHOP.GT.0).AND.(MULHOP.LT.NHOPS)) THEN c print*,'MULHOP=',MULHOP GOTO 40 ELSE IF (MULHOP.EQ.NHOPS) THEN GOTO 999 ENDIF 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 40 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,', NHOPS=',NHOPS 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,Bt0,Ph0,Th0,NHOPS,MODPROP 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 /QINIT/R0,T0,P0, & E0,D0,F0,Bt0,Ph0,Th0,NHOPS,MODPROP 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 PHB=DATAN(DSIN((TH-Th0)*RAD)*DTAN(Bt0*RAD))*DEG write(30,500) Y(1),Y(2),Y(3),Y(4),Y(5),Y(6) 500 format(6F10.3) write(21,1015) & REAL(TH), & REAL(RR-6370.0), & REAL(PH), & Ph0+PHB, & PH-(Ph0 + PHB), & REAL(DL), c & PSI, & 8.98*dsqrt(XNS(1)), & REAL(PMU) 1015 FORMAT(F8.2,F9.2,F9.3,F9.3,F9.3,F8.1,F10.1,F12.6) 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) go to 702 DMUDRR = 0.0d0 DMUDTH = 0.0d0 DMUDPH = 0.0d0 DMUDVR = 0.0d0 DMUDVT = 0.0d0 DMUDVP = 0.0d0 go to 703 702 continue 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) 703 continue 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 DMUDVT = 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 c fq00: kHz c rr00: km c th00: deg of gm-colat. c ph00: deg of gm-long. c es00: deg c dl00: deg ************************************************************************ 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) COMMON /QINIT /R0,T0,P0, & E0,D0,F0,Bt0,Ph0,Th0,NHOPS,MODPROP CALL DIPOLE(PATH) CALL DENS(PATH) DO 10 I=1,NSPEC PX(I) = PF0(I)*XNS(I)/FKC**2 PY(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 if (modprop.EQ.0) then PMUS= (-PB+SQRT(PD))/2.0/PA c Above is Ordinary mode else PMUS= (-PB-SQRT(PD))/2.0/PA c Above is Extraordinary mode end if 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) PARAMETER (NX=121,NY=41,NZ=431) DIMENSION PATH(8) DIMENSION ENE(NX,NY,NZ) 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 /ENECOM /ENE ZW=1.0 XW=1.0 YW=1.0 PHI00=115.0 c 高度（ALT)、余緯度（TH)、経度（PH)情報を読み込み C Altitude(ALT), CoLatitude(TH) and Longitude(PH) data C are read in. ALT=PATH(1)-6370. C TH=90.0-PATH(2) TH=PATH(2) PH=PATH(3) c1 print*,' ALT=', ALT, ' TH=',TH, ' PH=',PH c Start of Main program c 高度70kmより下と1200kmより上の時はすべて０を返す C Electron density is set zero above an altitude of C 1200km and below 70km. ALTMIN=70. ALTMAX=500. IF((ALT.LT.ALTMIN).OR.(ALTMAX.LT.ALT))THEN DO 250 IS=1,NSPEC XNS (IS) = 0. DNSDRR(IS) = 0. DNSDTH(IS) = 0. DNSDPH(IS) = 0. 250 CONTINUE GO TO 700 ELSE c 高度60km〜1200kmの間の時IRIデータより電子密度、勾配を求める C In the altitude range from 60 to 1200km, electron density C is calculated for given grid points from IRI model. c 高度、余緯度、経度の位置判定 C Judgement of Altitude, Colatitude and Longitude of the ray path c 高度の位置の判定 C Judgement of Altitude of the ray path DO IZ=1,(NZ-1) IF (((70.+(IZ-1)*ZW).LT.ALT).AND.(ALT.LT.(70.+IZ*ZW)))THEN NZ1=IZ NZ2=IZ+1 RR2=(ALT-(70.+(IZ-1)*ZW))/ZW ZZ=2 ELSE IF((70.+(IZ-1)*ZW).EQ.ALT)THEN NZ1=IZ NZ2=IZ+1 ZZ=1 ENDIF ENDDO c 緯度の位置の判定 C Judgement of Colaitude of the ray path DO IX=1,(NX-1) IF(((30.+(IX-1)*XW).LT.TH).AND.(TH.LT.(30.+IX*XW)))THEN NX1=IX NX2=IX+1 TH2=(TH-(30.+(IX-1)*XW))/XW XX=2 ELSE IF(TH.EQ.(30.+(IX-1)*XW))THEN NX1=IX NX2=IX+1 XX=1 ENDIF ENDDO c 経度位置の判定 C Judgement of Longitude of the ray path DO IY=1,(NY-1) IF(((PHI00+(IY-1)*YW).LT.PH).AND.(PH.LT.(PHI00+IY*YW))) THEN NY1=IY NY2=IY+1 PH2=(PH-(PHI00+(IY-1)*YW))/YW YY=2 ELSE IF(PH.EQ.(PHI00+(IY-1)*YW)) THEN NY1=IY NY2=IY+1 YY=1 END IF ENDDO c1 print*, ' NX1=',NX1,' NX2=',NX2,' NY1=',NY1,' NY2=',NY2,' NZ1=',NZ1,' NZ c1 12=',NZ2 c1 print*,' XX=',XX,' YY=',YY,' ZZ=',ZZ,' TH2=',TH2,' RR2=',RR2,' PH2=',PH2 c1 Stop c 高度、余緯度、経度の位置判定終わり C End of judgement of Altitude, Colatitude and Longitude c 電子密度、勾配を計算するために各格子の電子密度を計算する C c わかりやすいように基準点を変数に置き換える A1=ENE(NX1,NY1,NZ1) A2=ENE(NX2,NY1,NZ1) A3=ENE(NX1,NY2,NZ1) A4=ENE(NX2,NY2,NZ1) B1=ENE(NX1,NY1,NZ2) B2=ENE(NX2,NY1,NZ2) B3=ENE(NX1,NY2,NZ2) B4=ENE(NX2,NY2,NZ2) C Return c 内挿の準備点を求める c 高度、緯度面 C Altitude-Latitude plane C1=A1+(B1-A1)*RR2 C2=A1+(A2-A1)*TH2 C3=B1+(B2-B1)*TH2 C4=A2+(B2-A2)*RR2 D1=A3+(B3-A3)*RR2 D2=A3+(A4-A3)*TH2 D3=B3+(B4-B3)*TH2 D4=A4+(B4-A4)*RR2 E1=A1+(A3-A1)*PH2 E2=A2+(A4-A2)*PH2 E3=B1+(B3-B1)*PH2 E4=B2+(B4-B2)*PH2 c1 print*,' A1=',A1,' A2=',A2,' B1=',B1,' B2=',B2 c1 print*,' C1=',C1,' C2=',C2,' C3=',C3,' C4=',C4 c 内挿の準備点から面の内挿点を求める c 高度、緯度面 C Altitude-Latitude plane IF((ZZ.EQ.1).AND.(XX.EQ.1))THEN F1=A1 ELSE IF((ZZ.EQ.1).AND.(XX.EQ.2))THEN F1=C2 ELSE IF((ZZ.EQ.2).AND.(XX.EQ.1))THEN F1=C1 ELSE IF((ZZ.EQ.2).AND.(XX.EQ.2))THEN F1=C2+(C3-C2)*RR2 ENDIF C 高度、経度面 C Altitude-Longitude plane IF((ZZ.EQ.1).AND.(YY.EQ.1))THEN G1=A1 ELSE IF((ZZ.EQ.1).AND.(YY.EQ.2))THEN G1=E1 ELSE IF((ZZ.EQ.2).AND.(YY.EQ.1))THEN G1=D2 ELSE IF((ZZ.EQ.2).AND.(YY.EQ.2))THEN G1=E1+(E3-E1)*RR2 ENDIF C 緯度、経度面 C Latitude-Londitude plane IF((YY.EQ.1).AND.(XX.EQ.1))THEN H1=A1 ELSE IF((YY.EQ.1).AND.(XX.EQ.2))THEN H1=C2 ELSE IF((YY.EQ.2).AND.(XX.EQ.1))THEN H1=E1 ELSE IF((YY.EQ.2).AND.(XX.EQ.2))THEN H1=E1+(E2-E1)*TH2 ENDIF c 各格子の電子密度の計算終わり C So far, the electron density on the grid points were C calculated from IRI model. c 求めた各格子の電子密度より、電子密度、勾配を計算する C In the following, electron density and its space gradients C are calculated from the given electron density on the C grid point. c 高度・緯度・経度が格子上にある場合 C The case where the Altitude, Latitude and Longitude are C all on the grid point. IF((ZZ.EQ.1).AND.(XX.EQ.1).AND.(YY.EQ.1))THEN XNS (1) = A1 DNSDRR(1) = 2.*(B1-A1)/ZW/(B1+A1) DNSDTH(1) = 2.*(A2-A1)/XW/RAD/(A2+A1) DNSDPH(1) = 2.*(A3-A1)/YW/RAD/(A3+A3) c 高度・経度が格子上にあり緯度が格子上にない場合 C The case where the Altitude and Longitude are on the C grid point but the Latitude is not on the grid point. ELSE IF((ZZ.EQ.1).AND.(XX.EQ.2).AND.(YY.EQ.1))THEN XNS (1) = C2 DNSDRR(1) = 2.*(C3-C2)/ZW/(C3+C2) DNSDTH(1) = 2.*(A2-A1)/XW/RAD/(A2+A1) DNSDPH(1) = 2.*(D2-C2)/YW/RAD/(D2+C2) C 高度・緯度が格子上にあり、経度が格子上にない場合 C The case where the Altitude and Latitude are on the C grid point but the Longitude is not. ELSE IF((ZZ.EQ.1).AND.(XX.EQ.1).AND.(YY.EQ.2))THEN XNS (1) = E1 DNSDRR(1) = 2.*(E3-E1)/ZW/(E3+E1) DNSDTH(1) = 2.*(E2-E1)/XW/RAD/(E2+E1) DNSDPH(1) = 2.*(A3-A1)/YW/RAD/(A3+A1) C 緯度と経度が格子上にあり、高さが格子上にない場合 C The case where the Latitude and Longitude are on the C grid point but Altitude is not. ELSE IF((ZZ.EQ.2).AND.(XX.EQ.1).AND.(YY.EQ.1))THEN XNS (1) = E1 DNSDRR(1) = 2.*(B1-A1)/ZW/(B1+A1) DNSDTH(1) = 2.*(C4-C1)/XW/RAD/(C4+C1) DNSDPH(1) = 2.*(D1-C1)/YW/RAD/(D1+C1) C 緯度が格子上にあり、高さ・経度とも格子上にない場合 C The case where the Latitude is on the grid point, C but Altitude and Longitude are not. ELSE IF((ZZ.EQ.2).AND.(XX.EQ.1).AND.(YY.EQ.2))THEN XNS (1) = E1+(E3-E1)*RR2 DNSDRR(1) = 2.*(E3-E1)/ZW/(E3+E1) DNSDTH(1) = (((E2-E1)/(E2+E1))+((E4-E3)/(E4+E3))) 1/XW/RAD DNSDPH(1) = 2.*(D1-C1)/YW/RAD/(D1+C1) C 経度が格子上にあり、高さ・緯度ともに格子上にない場合 C The case where the Longitude is on the grid point but C both Altitude and Latitude are not on the grid point. ELSE IF((ZZ.EQ.2).AND.(XX.EQ.2).AND.(YY.EQ.1))THEN XNS (1) = C1+(C4-C1)*TH2 DNSDRR(1) = 2.*(C3-C2)/ZW/(C3+C2) DNSDTH(1) = 2.*(C4-C1)/XW/RAD/(C4+C1) DNSDPH(1) = 2.*(D2-C2)/YW/RAD/(D2+C2) C 高さが格子上にあり、経度・緯度ともに格子上にない場合 C The case where the Altitude is on the grid point, but both C Longitude and Latitude are not on the grid point. ELSE IF((ZZ.EQ.1).AND.(XX.EQ.2).AND.(YY.EQ.2))THEN XNS (1) = C2+(D2-C2)*PH2 DNSDRR(1) = 2.*(D2-C2)/YW/(D2+C2) DNSDTH(1) = 2.*(E2-E1)/XW/RAD/(E2+E1) DNSDPH(1) = 2.*(D2-C2)/YW/RAD/(D2+C2) C 高度・経度・緯度ともに格子上にない場合 C The case where the Altitude, Longitude and Latitude are C not on the grid point. ELSE IF((ZZ.EQ.2).AND.(XX.EQ.2).AND.(YY.EQ.2))THEN XNS (1) = (E1+(E3-E1)*RR2)+((E2+(E4-E2)*RR2)- 1(E1+(E3-E1)*RR2))*TH2 DNSDRR(1) = 2.*((C3+(D3-C3)*PH2)-(C2+(D2-C2)*PH2))/ 1ZW/((C3+(D3-C3)*PH2)+(C2+(D2-C2)*PH2)) DNSDTH(1) = 2.*((E2+(E4-E2)*RR2)-(E1+(E3-E1)*RR2))/ 1XW/RAD/((E2+(E4-E2)*RR2)+(E1+(E3-E1)*RR2)) DNSDPH(1) = 2.*((D2+(D3-D2)*RR2)-(C2+(C3-C2)*RR2))/ 1YW/RAD/((D2+(D3-D2)*RR2)+(C2+(C3-C2)*RR2)) ENDIF c Set the followings to zero except electron DO 610 IS=2,NSPEC XNS (IS) = 0. DNSDRR(IS) = 0. DNSDTH(IS) = 0. DNSDPH(IS) = 0. 610 CONTINUE c End in calculation of electron density and its gradients in three dimension 700 ENDIF RETURN END ********************************************************************** SUBROUTINE DIPOLE(PATH) ********************************************************************** IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) real*8 path(8) COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY COMMON /QIGRF /YR ,YT ,YP , & DYRDRR,DYTDRR,DYPDRR, & DYRDTH,DYTDTH,DYPDTH, & DYRDPH,DYTDPH,DYPDPH, & FH , & DBBDRR,DBBDTH,DBBDPH c1 RR = PATH(1) TH = PATH(2) c1 PH = PATH(3) COSTH = dCOS(TH*RAD) c1 SINTH = SIN(TH*RAD) TANTH = dSIN(TH*RAD)/dCOS(TH*RAD) TANTHS = TANTH**2 COSTHS = 1.0D0+3.0D0*COSTH**2 FH = -870*(6370/PATH(1))**3*dSQRT(COSTHS) YR = TANTH/dABS(TANTH)/dSQRT(1.0+0.25*TANTHS) YT = dSQRT(1-YR**2) YP = 0.0 RETURN END