C C modified by Gary Burgess C date: Aug. 18,1998 C update: Feb. 4,1999 C update: Nov. 24,1999 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C INVERSION OF SATELLITE H'(F) DATA TO N(H) PROFILES. C C C DEPOSITED IN NATIONAL SPACE SCIENCES DATA CENTER 1973. C C C METHODS USED DISCUSSED IN FOLLOWING REPORTS: C C C J.E. JACKSON, THE ANALYSIS OF TOPSIDE IONOGRAMS, C C C SEPT 1967, GSFC REPT X-615-67-452. C C C J.E. JACKSON, THE REDUCTION OF TOPSIDE IONOGRAMS TO C C C ELECTRON DENSITY PROFILES, JUNE 1969, PROC. IEEE. C C C J.E. JACKSON, THE P'(F) TO N(H) INVERSION PROBLEM IN C C C IONOSPHERIC SOUNDINGS, MAY 1971, GSFC REPT C C C X-625-71-186. C C C J.E. JACKSON AND C.J. MCQUILLAN, ANALYSIS OF TOPSIDE C C C IONOGRAMS WITH ALLOWANCE FOR SATELLITE MOTION AND C C C UNKNOWN FXS, 1973, GSFC REPT X-625-73- . C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IMPLICIT REAL*8(A-H,O-Y) REAL*4 DLAT,DLONG,DATE,HT COMPLEX*16 TITX,TITY,TYPE LOGICAL*1 TITL(80) real*4 xmid,ymid COMMON /SWJTCH/ ISAT,LOCSAT,IFXS,NCASES,ICODE(4),N,METHOD,IORD,IP COMMON /CNT/ZA(2500,3),ZXCO(50),LA,LXCO,ZA2(80,3),LA2 COMMON /MID/XMID,YMID COMMON /JKK/J common/com_icode/iicode C C REWIND THE TEMPORARY SAVE DATA FILE (UNIT 9)=LIB.DATA(SAVEIT) C C REWIND 9 C GET NAME OF INPUT X-TRACE FILE AND OPEN VARIOUS INPUT AND OUTPUT FILES. CALL GET_FILE C iicode: 1 = FIXED SATELLITE POSITION AND KNOWN FXS; C 2 = FIXED SATELLITE POSITION AND UNKNOWN FXS; C 3 = VARIABLE SATELLITE POSITION AND KNOWN FXS; C 4 = VARIABLE SATELLITE POSITION AND UNKNOWN FXS. icode(1)=iicode icode(2)=0 icode(3)=0 icode(4)=0 LA=0 LA2=0 LXCO=0 C READ PLOT AND CONTOUR CHARACTERISTICS FROM UNIT 4. CALCULATE CENTER C AND SCALING FACTOR(ZY AND ZX3) FOR POLAR CONTOURS. CALL READ4(ISP,IX,IY,ISVE,ISH,ITYP,ZHMAX,ZX3,ZY,TITX,TITY,TYPE, 1 TITL,ZSANG,ZEANG,LBL,IHDR,IHEM) C READ TITLE FOR ENTIRE JOB AND IDENTIFICATION, PROCESSING SWITCHES C (ICODE, SEE BELOW), AND H'(F) DATA FOR ONE IONOGRAM AT A TIME. C DETERMINE INITIAL SATELLITE POSITION BY INTERPOLATING WORLD MAP C DATA IF REQUESTED. WRITE OUT INPUT DATA. DETERMINE NCASES, THE NO. C OF DIFFERENT WAYS IN WHICH THE IONOGRAM IS TO BE PROCESSED. RETURN C TO STATEMENT 300 IF AN END-OF-FILE IS READ. 100 CALL NEWREAD(DLAT,DLONG,DATE,HT,LQ,&300) DO 200 ICT=1,NCASES J=ICODE(ICT) C ICODE(ICT) DETERMINES PROCEDURE TO BE USED FOR THIS ANALYSIS OF C THIS IONOGRAM, AS FOLLOWS: C 1 = FIXED SATELLITE POSITION AND KNOWN FXS; C 2 = FIXED SATELLITE POSITION AND UNKNOWN FXS; C 3 = VARIABLE SATELLITE POSITION AND KNOWN FXS; C 4 = VARIABLE SATELLITE POSITION AND UNKNOWN FXS. IF (ICT .NE. 1) THEN TYPE*,' CALLING PRINT ',J CALL PRINT !(J) ENDIF LOCSAT=0 IFXS=0 IF (J .EQ. 2 .OR. J .EQ. 4) IFXS=1 IF (J .EQ. 3 .OR. J .EQ. 4) LOCSAT=1 C ANALYZE IONOGRAM ACCORDING TO REQUIRED CONDITIONS: DETERMINE FXS C IF UNKNOWN (IFXS=1); ALLOW SATELLITE POSITION TO VARY IF REQUESTED C (LOCSAT=1). LAMINATION PROCEDURE FOR INVERSION OF H'(F) DATA TO C N(H) DATA IS CONTAINED IN SUBROUTINE XLAM. RETURN TO STATEMENT 100 C IF PROBLEMS ARE ENCOUNTERED IN COMPUTING FXS. CALL NEWINVERT(DLAT,DLONG,DATE,HT,LQ,&100) C COMPUTE ORDINARY TRACE VIRTUAL HEIGHTS FROM N(H) PROFILE DERIVED C FROM EXTRAORDINARY TRACE H'(F) DATA, IF REQUESTED (IORD=1). IF C PROBLEMS ARISE IN THIS COMPUTATION, DELETE IT FROM THE OUTPUT. IF(IORD.EQ.1) CALL NEWOTRACE(DLAT,DLONG,DATE,LQ) C INTERPOLATE RESULTANT N(H) DATA TO OBTAIN FINAL OUTPUT: HEIGHTS C AT FIXED DENSITIES AND DENSITIES AT FIXED HEIGHTS. RETURN TO C STATEMENT 100 IF PROBLEMS ARISE IN INTERPOLATION (I.E.TERMINATE C PROCESSING OF THIS IONOGRAM). CALL INTERP(&100) C IF REQUESTED, PLOT IONOGRAM INFO AND ELECTRON DENSITY VS. ALTITUDE. C CALCULATE AND SAVE VALUES FOR CONTOURING. IF(IP.NE.0) CALL PLOT1(IHDR,TITL,ISP,ZHMAX,ISVE,ITYP,IX,IY,DLAT,DL 1ONG) C CALCULATE FIELD LINES FOR DRAWING ON CONTOURS. CALL FLINE(DLAT,DLONG,DATE,HT,LQ,ZXCO,LXCO,IX,IHEM) 200 CONTINUE GO TO 100 C DRAW CONTOURS(FLAT OR CURVED EARTH OR BOTH),WITH X-AXIS LATITUDE,LONG, C MINUTES OR SECONDS, Y-AXIS ALTITUDE OR DISTANCE FROM SATELLITE. VALUES C MAY BE ELECTRON DENSITY,FN,FN/FH,FX,FZ,OR FT. FIELD LINES, SATELLITE C AND F-PEAK POINTS ARE OPTIONAL,AS ARE HEADERS. NOT ALL OPTIONS ARE C COMPATIBLE. 300 CONTINUE C 300 IF(IP.EQ.1) CALL CONT(ZA,LA,ZA2,LA2,ISP,ZX3,ZY,IHDR,TITX,TITY,TYPE C 1 ,TITL,ISH,ZXCO,LXCO,ZHMAX,ZSANG,ZEANG,LBL,IX,IY) STOP END SUBROUTINE FLINE(DLAT,DLONG,DATE,HT,LQ,XCO,LXCO,IX,IHEM) C FLINE CALCULATES FIELD LINES TO BE PLOTTED ON CONTOURS. ACTUAL C PLOTTING IS DONE IN CONT. COMMON /FLIN/ FLINES(2,3,180),NLINES DIMENSION XCOS(3,30),XCO(50),XCO1(30) DATA INIT /0/,N/1/ C IF N IS NEGATIVE, WE HAVE PREVIOUSLY HAD A REQUEST FOR FIELD LINES C TO BE CALCULATED AT ALL IONOGRAMS. IF(N.LT.0) GO TO 30 C IF N=0, WE HAVE PREVIOUSLY REQUESTED NO FIELD LINES. IF(N.EQ.0) RETURN IF(INIT.GT.0) GO TO 10 C NOW WE INITIALIZE THIS ROUTINE. NLINES=0 INIT=1 IX1=IX+1 FACT=1. C WE NOW READ IN THE FIELD LINE REQUEST, N<0 ALL IONOGRAMS, C N=0 NO IONOGRAMS, OTHERWISE N SPECIFIED IONOGRAMS. READ(4,1100) N 1100 FORMAT(I2) IF(N.EQ.0) RETURN IF(N.LT.0) GO TO 30 GO TO(100,110,120,100),IX1 1000 FORMAT(13(3F2.0,2X)) 120 FACT=1/60. C NOW WE READ IN THE X-COEFFICIENTS OF SELECTED IONOGRAMS, THEN C WE HAVE TESTS BASED ON THE SPECIFIED (BY IX, SEE READ4) KIND OF C X-COEFFICIENT WE ARE DEALING WITH. 110 READ(4,1000) ((XCOS(J,I),J=1,3),I=1,N) DO 5 I=1,N XCO1(I)=XCOS(1,I)*3600+XCOS(2,I)*60+XCOS(3,I) 5 XCO1(I)=XCO1(I)*FACT GO TO 10 100 READ(4,2000) (XCO1(I),I=1,N) 2000 FORMAT(13F6.2) WRITE(16,2010) 2010 FORMAT(1X,'FIELD LINE LATITUDES:') WRITE(16,2015) (XCO1(I),I=1,N) 2015 FORMAT(1X,13(F6.2,2X)) C NOW WE CHECK TO SEE IF THIS IONOGRAM IS A SPECIFIED ONE. 10 DO 20 I=1,N IF(ABS(XCO1(I)-XCO(LXCO)).LT..015) GO TO 30 20 CONTINUE RETURN C IF THIS IS A SPECIFIED IONOGRAM, OR IF ALL WERE REQUESTED, CALCULATE C THE FIELD LINES. 30 NLINES=NLINES+1 FLINES(1,1,NLINES)=DLAT FLINES(1,2,NLINES)=DLONG FLINES(1,3,NLINES)=HT HT6=HT/6 C NOW WE CALCULATE THE FIELD LINE IN SIX PIECES. DO 40 I=1,6 N1=NLINES+I-1 CALL TRACE(FLINES(1,1,N1),FLINES(1,2,N1),FLINES(1,3,N1),FLINES(2,1 X ,N1),FLINES(2,2,N1),FLINES(2,3,N1),HT6,DATE,IHEM) C NOW STORE THE CALCULATED VALUES, AND INCREMENT THE COUNTER. DO 40 J=1,3 40 FLINES(1,J,N1+1)=FLINES(2,J,N1) NLINES=NLINES+5 RETURN END SUBROUTINE INTERP(*) C SUBROUTINE INTERP. INTERPOLATE THE N(H) DATA OBTAINED FROM THE C SCALED X TRACE TO DETERMINE HEIGHTS AT FIXED DENSITIES AND C DENSITIES AT FIXED HEIGHTS. WRITE OUT INTERPOLATED RESULTS. IMPLICIT REAL*8(A-H,O-Z) C REAL*4 DLAT,DLONG,DATE COMMON /CALC/ AA(50),BB(50),F(50),HP(50),FH(50),TH(70),DEN(70) COMMON /SWJTCH/ ISAT,LOCSAT,IFXS,NCASES,ICODE(4),N,METHOD,IORD,IP COMMON /PL/ TH1(200),ED1(200),TH2(200),ED2(200),M,M1 COMMON /OCALC/ FO(50),HPO(50) COMMON /IO/ IYR,IDY,ITI,LMT,TITLE(40),DIP,FHS DIMENSION DENREF(56),D(2),Z(2) c DATA THINT/50.0/ data thint/20.0/ DATA DENREF / 1.0E2,2.0E2,2.5E2,3.0E2,4.0E2,5.0E2,6.0E2,7.0E2, D 8.0E2,9.0E2,1.0E3,1.25E3,1.5E3,2.0E3,2.5E3,3.0E3, D 4.0E3,5.0E3,6.0E3,7.0E3,8.0E3,9.0E3,1.0E4,1.2E4, D 1.4E4,1.6E4,1.8E4,2.0E4,2.3E4,2.7E4,3.0E4,3.5E4, D 4.0E4,5.0E4,6.0E4,7.0E4,8.0E4,9.0E4,1.0E5,1.25E5, D 1.5E5,2.0E5,2.5E5,3.0E5,3.5E5,4.0E5,4.5E5,5.0E5, D 6.0E5,7.0E5,8.0E5,9.0E5,1.0E6,1.25E6,1.5E6,2.0E6 / C WRITE OUT RESULTS FROM SCALED POINTS. WRITE(16,1) WRITE(16,2) (TH(I),DEN(I),I=1,N) C DENSITIES AT FIXED HEIGHTS OBTAINED FROM LAMINATION CONTAINING C DESIRED HEIGHT AND FROM EQUATION OF LAMINATION PT=3500.0 1400 IF (PT .LT. TH(1)) GO TO 1500 PT=PT-THINT c type*,' pt th(1) = ',pt,th(1) IF (PT) 1600,1600,1400 1500 CONTINUE type*,' ' ED1(1)=DEN(1) TH1(1)=TH(1) D2LN=DLOG(DEN(1)) PT1=TH(1) K=2 DO 1700 I=2,N D1LN=D2LN D2LN=DLOG(DEN(I)) PT2=TH(I) 1800 IF(PT2.GE.PT1) GO TO 1600 IF(PT2.GT.PT) GO TO 1900 IF (PT2 .LT. PT) GO TO 2000 PT=PT-THINT IF (PT .LE. 0.0) GO TO 1600 2000 TH1(K)=PT IF (METHOD .EQ. 2 .AND. (BB(I)**2) .GT. 0.01) GO TO 2100 ED1(K)=D2LN+(PT-PT2)*(D1LN-D2LN)/(PT1-PT2) c type*,'1, K,ED1 =',k,ed1(k) GO TO 2200 2100 CONTINUE ED1(K)=D1LN-(AA(I)+DSQRT(AA(I)**2-4.*BB(I)*(PT1-PT)))/(2.*BB(I)) c type*,'2, K,ED1 =',k,ed1(k) 2200 K=K+1 PT=PT-THINT IF (PT) 1600,1600,1800 1900 PT1=PT2 1700 CONTINUE TH1(K)=TH(N) ED1(K)=DEN(N) M=K-1 DO 800 I=2,M 800 ED1(I)=DEXP(ED1(I)) C WRITE OUT INTERPOLATED RESULTS. WRITE(16,3) WRITE(16,2) (TH1(I),ED1(I),I=1,K) WRITE(20,3) WRITE(20,12) (th1(i),ed1(i),i=1,k) C HEIGHTS AT FIXED DENSITIES OBTAINED FROM LAMINATION CONTAINING C DESIRED DENSITY AND FROM EQUATION OF LAMINATION. N=N-1 M1=1 K=1 DO 200 L=1,N K=K-1 DO 300 I=1,2 D(I)=DEN(L+I-1) 300 Z(I)=TH(L+I-1) D1LN=DLOG(D(1)) 400 K=K+1 D1=DENREF(K) DO 500 I=1,2 ED2(M1)=D1 IF (D1-D(I)) 600,700,500 700 TH2(M1)=Z(I) M1=M1+1 IF(I.EQ.1) GO TO 400 500 CONTINUE GO TO 200 600 IF(I.EQ.1) GO TO 400 V=DLOG(D1)-D1LN TH2(M1)=Z(1)+V*(AA(L+1)+BB(L+1)*V) M1=M1+1 IF(K.GE.56) GO TO 900 GO TO 400 200 CONTINUE 900 M1=M1-1 N = N+ 1 C WRITE OUT INTERPOLATED RESULTS. WRITE (16,4) WRITE(16,2) (TH2(I),ED2(I),I=1,M1) WRITE (21,4) WRITE(21,12) (TH2(I),ED2(I),I=1,M1) IF(IORD.EQ.0) RETURN WRITE(16,8) WRITE(16,9) (FO(I),HPO(I),I=1,N) RETURN 1600 WRITE (16,7) RETURN 1 ENTRY INREF WRITE (16,5) WRITE (16,6) (DENREF(K),K=1,56) RETURN 1 FORMAT(' ALTITUDES(KM) AND DENSITIES(EL/CC) CORRESPONDING TO SCALE FD H''(F) DATA') 2 FORMAT(1H ,3X,4('ALT.',3X,'DENSITY',5X)/4(2X,0PF6.1,1X,1PE10.3)) 12 FORMAT(1H ,3X,('ALT.',3X,'DENSITY',5X)/(2X,0PF6.1,1X,1PE10.3)) 3 FORMAT(' ALTITUDES AND DENSITIES AT SELECTED HEIGHTS') 4 FORMAT(' ALTITUDES AND DENSITIES AT SELECTED DENSITIES') 5 FORMAT(1x,'REFERENCE DENSITIES FOR INTERPOLATION') 6 FORMAT(2X,1P8E10.2) 7 FORMAT(1x,'ERROR IN DATA. FIXED-HEIGHT INTERPOLATION NOT POSSIBLE F.'/1H ,'FURTHER PROCESSING OF THIS IONOGRAM TERMINATED.') 8 FORMAT(' CALCULATED ORDINARY TRACE'/' ',5(3X,'F',6X,'HP',2X)) 9 FORMAT(5(1X,F6.3,1X,F6.1)) END SUBROUTINE NEWOTRACE(DLAT,DLONG,DATE,LQ) C SUBROUTINE OTRACE: LAMINATION PROCEDURE FOR COMPUTING O-TRACE C H'(F) DATA FROM N(H) PROFILE DERIVED FROM X-TRACE. INTEGRATE C USING GAUSSIAN QUADRATURE WITH 3 COEFFICIENTS (EXCEPTION: 16 C ARE USED ON LAST TWO LAMINATIONS). THE CHANGE OF VARIABLE C T**2 = 1-X C IS USED TO KEEP THE INTEGRAND FINITE AT ALL POINTS. IMPLICIT REAL*8(A-H,O-Z) REAL*4 DLAT,DLONG,HGT,DATE,BN,BE,BV,BT COMMON /CALC/ AA(50),BB(50),F(50),HP(50),FH(50),TH(70),DEN(70) COMMON /MISC/ DP,FH1,ASUM,BSUM,FA,CPH,X1,X2,GMT1 COMMON /SWJTCH/ ISAT,LOCSAT,IFXS,NCASES,ICODE(4),N,METHOD,IORD,IP COMMON /OCALC/ FO(50),HPO(50) COMMON /IO/ IYR,IDY,ITI,LMT,TITLE(40),DIP,FHS DIMENSION SS(19),W(19) DATA SS/.77459667D0,0.D0,-.77459667D0,.98940094D0,.94457502D0, D .86563120D0,.75540441D0,.61787624D0,.45801678D0, D .28160355D0,.09501251D0,-.09501251D0,-.28160355D0, D -.45801678D0,-.61787624D0,-.75540441D0,-.86563120D0, D -.94457502D0,-.98940094D0 / DATA W/.55555556D0,.88888889D0,.55555556D0,.02715246D0, D .06225352D0,.09515851D0,.12462897D0,.14959599D0, D .16915652D0,.18260342D0,.18945061D0,.18945061D0, D .18260342D0,.16915652D0,.14959599D0,.12462897D0, D .09515851D0,.06225352D0,.02715246D0 / DTF=8.98D-3 DO 100 I=1,N 100 FO(I)=DTF*DSQRT(DEN(I)) HPO(1)=0. FS=F(1) C SAVE INITIAL VALUES OF PARAMETERS. DLT=DLAT DLG=DLONG TH1=TH(1) DEN1=DEN(1) AA2=1./AA(2) FACT=DEN1*DEXP(-TH1*AA2) M=2 DO 1300 I=2,N DP=0. FA=FO(I) IF(M.EQ.3 .OR. LOCSAT.EQ.0) GO TO 1000 CALL FRTIME(FS,TIME) C CONVERT TIME FROM DECIMAL HR TO HHMMSS. IH=TIME RM=(TIME-IH)*60. IM=RM IS=(RM-IM)*60 C REPLACE ORBINT WITH CALL TO ROUTINE TO CALCULATE ORBITAL ELEMENTS. DLAT1=DLAT DLONG1=DLONG HGT1=HGT CALL JWORLDMAP(ISAT,IYR,IDY,IH,IM,IS,DLAT1,DLONG1,HGT1,0) DLAT=DLAT1 DLONG=DLONG1 HGT=HGT1 C CALL ORBINT(1,TIME,DLAT,DLONG,HGT) TH(1)=HGT DEN(1)=FACT*DEXP(TH(1)*AA2) IF(TH(1).GT.TH(2)) GO TO 2800 WRITE(16,1) TH(1),FA,DLAT,DLONG,TH(2) TH(1)=TH(2) M=3 GO TO 1000 2800 DO 3100 L=1,I CALL FIELDG(DLAT,DLONG,HGT,DATE,10,LQ,BN,BE,BV,BT) FH(L)=BT*2.8D-5 3100 HGT=TH(L+1) 1000 FA2=1.0/DEN(I) FA1=1.0/FA DO 300 J=M,I IL=1 IU=3 IF(J.LT.(I-1)) GO TO 400 IL=4 IU=19 400 X1=DEN(J-1)*FA2 X2=DEN(J)*FA2 X1LN=DLOG(X1) X2LN=DLOG(X2) Y1=FH(J-1)*FA1 Y2=FH(J)*FA1 C YR=Y AT REFLECTION POINT. CHANGE FROM VARIABLE X TO T.USE B9 AND C C9 TO CHANGE Y WITHIN EACH LAMINATION. T1=DSQRT(1.D0-X1) T2=0. IF(J.NE.I) T2=DSQRT(1.D0-X2) S1 = 0.0 S2=0.0 B9 = DLOG(Y2/Y1)/3.0 V=X2LN-X1LN RB=BB(J)/AA(J) RA = 1.0+RB*V C9 =(DEXP(B9)-1.0)/(V*RA) T1=0.5*(T1+T2) T2=T1-T2 C GAUSSIAN INTEGRATION DO 900 K=IL,IU T=T2*SS(K)+T1 X=(1.0-T**2) RX=1.0/(1.0-X) IF(CPH.GE.0.0009) GO TO 1200 C UX IS VALUE OF INTEGRAND AT THE REFLECTION POINT 2400 UX=T*DSQRT(RX)/X GO TO 1100 1200 A=X2LN-DLOG(X) Y = Y2/((1.0+C9*A*(1.0+RB*(2.0*V-A)))**3) YL = Y * CPH YT2=Y**2-YL**2 R=0.5*YT2*RX/YL SPSI=DSIN(DATAN(R)) A=(1.0+SPSI)*DSQRT(1.0+R**2) SO=1.+YL/A A=X/SO IF (A.GE.1.0D0) GO TO 1400 IF(A.GE.0.9999999D0) GO TO 2400 UX=T*(1.+.5*A*(1.-SO)*(1.-(1.+X)*SPSI*RX)/SO)/(X*DSQRT(1.-A)) 1100 S1=S1+UX*W(K) 900 S2=S2+UX*W(K)*(DLOG(X)-X1LN) BSUM=2.0*T2 ASUM=BSUM*S1 BSUM=2.0*BSUM*S2 300 DP=DP-AA(J)*ASUM-BB(J)*BSUM 1300 HPO(I)=DP 1500 DLAT=DLT DLONG=DLG TH(1)=TH1 DEN(1)=DEN1 FH(1)=FH1 RETURN C RESET IORD TO OMIT ORDINARY TRACE FROM OUTPUT. 1400 WRITE(16,2) J,I,X,SO,A IORD=0 GO TO 1500 1 FORMAT('0SATELLITE (',F6.1,' KM) BELOW FIRST LAMINATION AT',F7.3, F ' MHZ ON O-TRACE.'/' FIXED POSITION ASSUMED: LAT =',F8.2, F ' LONG =',F8.2,' HT =',F7.1) 2 FORMAT('0FOR OTRACE LAM',I3,' POINT',I3,' X =',1PD11.3,' 1-Y** F2/(1-X) =',D11.3/' RATIO =',D11.3,'>1. CALCULATIONS TERMINATED.') END SUBROUTINE NEWREAD(DLAT,DLONG,DATE,HT,LQ,*) C SUBROUTINE READ. READ TITLE FOR JOB. FOR EACH IONOGRAM, READ ID C (NO. OF POINTS, YR,DAY #,GMT OF DATA START); PROCESSING SWITCHES: C MAP=1 IF SATELLITE POSITION TO BE INTERPOLATED FROM INPUT WORLD C MAP DATA; MAP=0 IF THIS DATA IS INPUT DIRECTLY. C ISAT = SATELLITE ID: 1=AL-1; 2=AL-2; 3=ISIS-1; 4=ISIS-2. C IF ZERO OR BLANK, VARIABLE POSITION ANALYSIS NOT PERFORMED C ICODE IDENTIFIES CONDITIONS APPLIED IN EACH ANALYSIS (SEE MAIN) C IORD=1 IF O-TRACE VIRTUAL HEIGHTS ARE TO BE COMPUTED FROM C X-TRACE N(H) PROFILE. IF ZERO, OMIT THIS CALCULATION. C IP=1 IF PLOTTING DESIRED. IP=0 OR BLANK OTHERWISE. C METH IDENTIFIES TYPE OF LAMINATION USED IN ANALYSIS (1=LINEAR C IN LOG(N); 2(OR BLANK)=PARABOLIC IN LOG(N)). C AND SCALED H'(F) DATA. IMPLICIT REAL*8(A-H,O-Z) REAL*8 APPF(8)/2.,2.165,2.39,2.75,3.,5.,5.52,6./, REALF(8)/2.,2.1, R 2.3,2.7,3.,5.,5.5,6./, SATID(5)/'AL-1: ','AL-2: ','ISIS-1: ' R ,'ISIS-2: ',' '/, RAD/57.29577951D0/ REAL*4 DLAT,DLONG,DATE,HT,BN,BE,BV,BT LOGICAL*4 TEST LOGICAL*1 STRING(80) COMMON /CALC/ AA(50),BB(50),F(50),HP(50),FH(50),TH(70),DEN(70) COMMON /MISC/ DP,FH1,ASUM,BSUM,FA,CPH,X1,X2,GMT1 COMMON /TIME/IH,IM,IS COMMON /IO/ IYR,IDY,ITI,LMT,TITLE(40),DIP,FHS,MON,IDAY COMMON /SWJTCH/ ISAT,LOCSAT,IFXS,NCASES,ICODE(4),N,METHOD,IORD,IP COMMON /JKK/JK DIMENSION ITLO(3),DDEC(27),DET(27),DNO(27),DD(27,4),EQM(27,4) DATA ISKP/0/, K/0/ C OPEN(UNIT=3,NAME='NHINPUT.DAT',TYPE='OLD',FORM='FORMATTED') IF (ISKP .NE. 0) GO TO 100 C LIST ALL USER INPUT DATA CARDS. I1=0 I=0 J=1 WRITE(16,22)J 1400 READ(3,23,END=1300) STRING I1=I1+1 IF(I1.LT.29) GO TO 1400 I=I+1 IF(MOD(I,57).NE.0) GO TO 1500 J=J+1 WRITE(16,22)J 1500 WRITE(16,24) STRING GO TO 1400 1300 REWIND 3 C SET UP AND TEST FIELDG. INITIALIZE PARAMETERS FOR THIS JOB. LQ = 0 AA(1)=0.0 BB(1)=0.0 BB(2)=0.0 BB(3)=0. CALL FIELDG(-60.,-180.,100.,1960.,10,LQ,BN,BE,BV,BT) WRITE (16,11) BN,BE,BV,BT C INPUT: DECLINATION FROM NAUTICAL ALMANAC, EQUATION OF TIME FROM EP C DATA AT NOON GREENWICH, IN DECIMAL UNITS, FOR 4 SUCCESSIVE YEARS. DO 25 J=1,4 25 READ(3,18) (DNO(I),DD(I,J),EQM(I,J),I=1,27) DO 20 I=1,27 20 DNO(I)=DNO(I)+0.5 RADI=1.0D0/RAD STM=1.0D0/60.0D0 HTD=1.0D0/24.0D0 CALL INREF READ (3,1) TITLE ISKP=1 100 READ(3,2,END=200) MAP,ISAT,IORD,IP,METH C WRITE(*,2) MAP,ISAT,IORD,IP,METH CALL READ_SCALED(N,IYR,IDY,IH,IM,rIS,F,HP) is=ris DATE=FLOAT(IYR)+FLOAT(IDY)/365.0+1900.0 ITI=IH*10000 + IM*100 + IS c GMT1=DFLOAT(IH)+DFLOAT(IM)/60.0+DFLOAT(IS)/3600.0 GMT1=DFLOAT(IH)+DFLOAT(IM)/60.0+rIS/3600.0 IF (MAP .LT. 1) GO TO 300 TIME=GMT1 C REPLACE CALL TO ORBINT WITH CALL TO WORLDMAP PROGRAM TO CALCULATE OBBITAL C ELEMENTS DLAT1=DLAT DLONG1=DLONG HT1=HT c type*,' worldmap from newread not used' c type*,' enter lat, long, hgt' c accept*, dlat1,dlong1,ht1 CALL JWORLDMAP(ISAT,IYR,IDY,IH,IM,IS,DLAT1,DLONG1,HT1,1) DLAT=DLAT1 DLONG=DLONG1 HT=HT1 C CALL ORBINT(0,TIME,DLAT,DLONG,HT) IF(MAP.EQ.2) READ(3,3) FHS GO TO 400 300 READ(3,3) HT,DIP,FHS,DLAT,DLONG C 400 READ(3,3) (HP(I),F(I),I=1,N) 400 CONTINUE IF (ISAT .NE. 3) GO TO 800 C CORRECT FREQUENCIES SCALED FOR ISIS-1 IONOGRAMS. DO 900 I=1,N FREQ=F(I) IF((FREQ.GT.2..AND.FREQ.LT.3.) .OR. (FREQ.GT.5..AND.FREQ.LT.6.)) I F(I)=TAB(APPF,REALF,8,1,FREQ) 900 CONTINUE C CALCULATE LOCAL MEAN TIME 800 TLO=GMT1+DLONG/15.D0+24.D0 DO 700 I=1,3 ITLO(I)=TLO 700 TLO=(TLO-ITLO(I))*60.0 LMT=ITLO(1)*1.E4+ITLO(2)*1.0E2+ITLO(3) IF (LMT .GT. 240000) LMT=LMT-240000 C COMPUTE MONTH AND DAY FROM DAY OF YEAR. COMPUTE SOLAR ZENITH ANGLE C USING DECLINATION AND EQN OF TIME FOR CORRESPONDING YEAR IN PERIOD C 1968 TO 1971. CALL DAYMON(IYR,IDY,MON,IDAY) I=K K=MOD(IYR,4)+1 IF(K.EQ.I) GO TO 1200 DO 1100 I=1,27 DDEC(I)=DD(I,K) 1100 DET(I)=EQM(I,K) 1200 CPH=DLAT*RADI DDATE=GMT1*HTD+DFLOAT(IDY) ET=TAB(DNO,DET,27,4,DDATE)*STM HA=(DABS(180.0D0-DLONG-15.0D0*(GMT1-ET)))*RADI DEC=TAB(DNO,DDEC,27,4,DDATE)*RADI COSZ=DSIN(DEC)*DSIN(CPH)+DCOS(DEC)*DCOS(CPH)*DCOS(HA) ZENI=DACOS(COSZ)*RAD HA=HA*RAD DEC=DEC*RAD ET=ET*60.0D0 C COMPUTE DIP AND FH AT ORIGINAL SATELLITE POSITION CALL FIELDG(DLAT,DLONG,HT,DATE,10,LQ,BN,BE,BV,BT) DIP=ATAN2(BV,SQRT(BT**2-BV**2))*RAD C CPH=COS(PROPAGATION ANGLE). VERTICAL PROPAGATION IS ASSUMED. CPH=DCOS((90.0-DABS(DIP))/RAD) FH(1)=BT*2.8000E-05 C EXT. RAY EXIT FREQ. MUST BE GREATER THAN FHS FXS=FH(1)+0.001 IF (F(1).GE.FXS) GO TO 710 F(1)=FXS WRITE (16,17) 710 CONTINUE FH1=FH(1) IF(MAP.EQ.2) F(1)=F(1)-FHS+FH1 FHS = FH(1) TH(1) = HT C DETERMINE NCASES, THE NUMBER OF WAYS IN WHICH THIS IONOGRAM IS TO C BE PROCESSED. IF (ICODE(1) .EQ. 0) GO TO 1000 C EXCLUDE VARIABLE SATELLITE ANALYSIS IF ORBIT DATA OR SATELLITE ID C NOT SPECIFIED, AND CALL FOR FIXED POSITION ANALYSIS IF NOT C PREVIOUSLY PERFORMED. TEST=MAP.LT.1 .OR. ISAT.EQ.0 DO 600 I=1,4 J=I IF (TEST .AND. ICODE(I).GT.2) ICODE(I)=ICODE(I)-2 IF (ICODE(I) .EQ. 0 .OR. (TEST .AND. J.GT.2)) GO TO 500 600 CONTINUE J=5 GO TO 500 1000 ICODE(1)=1 J=2 500 NCASES=J-1 J=ICODE(1) C INITIALIZE METHOD FOR TYPE OF LAMINATION DESIRED. WRITE OUT INPUT. C TYPE*,' ENTRY START' ENTRY PRINT !(JK) C TYPE*,' JK IN PRINT ',JK C TYPE*,' ENTRY END ' J=JK METHOD = 2 IF(METH.EQ.1)METHOD=1 WRITE (16,4) TITLE IF(METHOD.EQ.1)WRITE(16,19) IF(METHOD.EQ.2)WRITE(16,21) WRITE(16,5) IDY,MON,IDAY,IYR,ITI,LMT,DEC,ET,ZENI WRITE (16,6) ITI,HT,DLAT,DLONG,DIP,FHS I2=ISAT IF(ISAT.EQ.0) I2=5 WRITE(16,16) SATID(I2) IF (J .EQ. 1 ) WRITE (16,12) IF (J .EQ. 2 ) WRITE (16,13) IF (J .EQ. 3 ) WRITE (16,14) IF (J .EQ. 4 ) WRITE (16,15) IF (ISAT .EQ. 3) WRITE (16,9) WRITE (16,7) WRITE (16,8) (F(I),HP(I),I=1,N) RETURN 200 RETURN 1 1 FORMAT (20A4) 2 FORMAT(2I10,10X,I5,I2,I3) 3 FORMAT (8F10.4) 4 FORMAT(1x,20A4/1H ,20A4/' GEOMAGNETIC FIELD MODEL GSFC 9/65, NMAX F=10 (HENDRICKS AND CAIN, 1966).') 5 FORMAT(' INPUT DATA'/' DAY',I4,' (',I2,'/',I2,') OF 19',I2,' AT', F I7,' GMT LOCAL MEAN TIME',I7/' DECLINATION',F7.2,' DEG EQN OF FTIME',F7.2,' MIN ZENITH ANGLE',F8.2,' DEG') 6 FORMAT(' SATELLITE POSITION AT ',I6,' GMT'/' HT=',F7.1,' LAT=', F F7.2,' LONG=',F7.2,' DIP=',F7.2,' FHS=',F5.3) 7 FORMAT(1x,'X TRACE'/1H ,5(3X,'F',6X,'HP',2X)) 8 FORMAT(5(1X,F6.3,1X,F6.1)) 9 FORMAT(1H ,8X,'CORRECTED SCALED FREQUENCIES FOLLOW.') 11 FORMAT(1H ,'SEPTEMBER 65 COEFFICIENTS'/1H ,2X,'LAT',5X,'LONG',5X, F 'HT',4X,'DATE',6X,'BN',11X,'BE',11X,'BV',11X,'BT'/1H , F '-60.00',2X,'-180.00',2X,'100.0',2X,'1960.0',4F13.4) 12 FORMAT(1H ,8X,'FIXED SATELLITE POSITION AND KNOWN FXS ASSUMED.') 13 FORMAT(1H ,8X,'FIXED SATELLITE POSITION AND UNKNOWN FXS ASSUMED.') 14 FORMAT(1H ,8X,'VARIABLE SATELLITE POSITION AND KNOWN FXS ASSUMED.' F ) 15 FORMAT(1H ,8X,'VARIABLE SATELLITE POSITION AND UNKNOWN FXS ASSUMED F.') 16 FORMAT(' ',A8) 17 FORMAT(1H ,'FREQUENCY READ TOO LOW ON FIRST DATA POINT.') 18 FORMAT(4(F4.0,2F8.5)) 19 FORMAT(' LAMINATIONS LINEAR IN LOG(N)') 21 FORMAT(' LAMINATIONS PARABOLIC IN LOG(N)') 22 FORMAT(//' ***** LISTING OF ALL INPUT DATA CARDS: PAGE',I3, F ' *****'//) 23 FORMAT(80A1) 24 FORMAT(' ',80A1) END SUBROUTINE READ4(ISP,IX,IY,ISVE,ISH,ITYP,HMAX,X3,Y,TITX, 1 TITY,TYPE,TITL,SANG,EANG,LBL,IHDR,IHEM) C READ4 IS A SUBROUTINE TO INITIALIZE ALL OF THE NEW PLOT AND CONTOUR C CONTROLLING PARAMETERS. THESE ARE EITHER READ FROM FT04F001 OR C CALCULATED. LOGICAL*1 TITL(80) COMPLEX*16 TITX,TITY,TITLES(13),TYPE C COMMON /OPTION/ IS A COMMON BLOCK USED BY THE WOLFPLOT CONTOUR PACKAGE COMMON /OPTION/IMIN,IMAX,JMIN,JMAX,FINTER,NUMV,CONTUR(50) C TITLES HOLDS THE VARIOUS TITLES ASSOCIATED WITH DIFFERENT CHOICES OF C INPUT PARAMETERS. DATA TITLES/' LATITUDE (DEG) ',' TIME (MINUTES) ', 1' TIME (SECONDS) ',' LONGITUDE (DEG)',' ALTITUDE (KM) ', 2' DISTANCE (KM) ','ELECTRON DENSITY','PLASMA FREQUENCY', 3' FN/FH ',' X-MODE CUTOFF ',' Z-MODE CUTOFF ', 4' UPPER HYBRID ',' FH '/ C HMAX IS THE MAXIMUM ALTITUDE SHOWN ON THE ELECTRON DENSITY PROFILES C AND THE CONTOURS. IT DEFAULTS TO 4000 KM, AND SHOULD BE A MULTIPLE OF C 500 KM, NO LARGER THAN 4000 KM. C FOR VAX OPEN FILE C OPEN(UNIT=4,NAME='NH2.DAT',TYPE='OLD',FORM='FORMATTED',READONLY) HMAX=4000. FINTER=0. NUMV=8 IMIN=1 IMAX=80 JMIN=1 JMAX=80 C THE FOLLOWING DEFAULT CONTOUR LEVELS ARE USED ONLY FOR ELECTRON C DENSITY CONTOURS. CONTUR(1)=1.E2 CONTUR(2)=3.E2 CONTUR(3)=1.E3 CONTUR(4)=3.E3 CONTUR(5)=1.E4 CONTUR(6)=3.E4 CONTUR(7)=1.E5 CONTUR(8)=3.E5 C WE NOW READ IN THE MAXIMUM ALTITUDE, AND NUMBER OF CONTOUR LEVELS. READ(4,1100) HMAX1,NUMV1 1100 FORMAT(F6.0,I2) C WE NOW CHECK FOR NEW CONTOUR LEVELS, IF THE DEFAULT IS NOT BEING USED. IF(NUMV1.LE.0) GO TO 1 NUMV=NUMV1 READ(4,1200) (CONTUR(I),I=1,NUMV) 1200 FORMAT(8E10.3) C WE NOW SET THE MAXIMUM ALTITUDE, IF THE DEFAULT IS NOT BEING USED. 1 IF(HMAX1.LE.0) GO TO 2 HMAX=HMAX1 2 CONTINUE C WE NOW READ IN A SERIES OF FLAGS: C ISP CONTROLS THE USE OF SATELLITE AND F-PEAK POINTS ON THE C CONTOURS(0=NONE,1=SATELLITE,2=F-PEAK,3=BOTH). C IX CONTROLS THE X AXIS ON CONTOURS (0=LATITUDE, 1=TIME(SECS), C 2=TIME(MINS), 3=LONGITUDE). C IY CONTROLS Y AXIS(0=ALTITUDE, 1=DISTANCE FROM SATELLITE). C ISVE IS NOT YET IMPLEMENTED. C IHDR CONTROLS THE PLOTTING OF TITLES(0=YES, 1=NO). C ISH CONTROLS THE CONTOUR SHAPE(0=NONE, 1=FLAT EARTH, C 2=CURVED EARTH, 3=PLOT BOTH CONTOURS). C ITYP CONTROLS THE VALUE BEING CONTOURED (0=ELECTRON C DENSITY, 1=FN, 2=FN/FH, 3=FX, 4=FZ, 5=FT, 6=FH). C LBL CONTROLS THE LABELING OF CONTOUR LINES (SEE WOLFPLOT C CONTOURING MANUAL FOR DETAILS). READ(4,1000) ISP,IX,IY,ISVE,IHDR,ISH,ITYP,LBL 1000 FORMAT(7I1,I2) C FOR CONTOURS OTHER THAN ELECTRON DENSITY FOR WHICH NON-DEFAULT C CONTOUR LEVELS HAVE NOT BEEN SET, SET NEW DEFAULTS. IF(ITYP.EQ.0.OR.NUMV1.NE.0) GO TO 3 CONTUR(1)=.1 CONTUR(2)=.5 CONTUR(3)=1. CONTUR(4)=2. CONTUR(5)=3. CONTUR(6)=4. CONTUR(7)=5. CONTUR(8)=6. C IF HEADERS ARE TO BE USED, READ IN A USER SPECIFIED TITLE. 3 IF(IHDR.EQ.0) READ(4,2000) (TITL(I),I=1,80) 2000 FORMAT(80A1) TITX=TITLES(IX+1) TITY=TITLES(IY+5) TYPE=TITLES(ITYP+7) C NOW THAT THE TITLES HAVE BEEN SELCTED, CHECK FOR CURVED-EARTH CONTOURS C IF THEY ARE PRESENT, READ IN THE START AND END ANGLES DESIRED FOR THE C CONTOURS (IN DEGREES). IF(ISH.GT.1) READ(4,3000) SANG,EANG 3000 FORMAT(2F6.2) C C SET NORTH-SOUTH HEMISHERE FLAG,IHEM C IHEM=1 - SOUTHERN HEMISPHERE C IHEM=-1 - NORTHERN HEMISPHERE C IF(ISH .GT. 1 .AND. SANG .LT. 0.) IHEM=1 IF(ISH .GT. 1 .AND. SANG .GT. 0.) IHEM=-1 C C C PRINT OUT CONTROLLING VALUES C WRITE(16,4000) 4000 FORMAT(1X,'INITIAL CONDITIONS:') WRITE(16,4002)HMAX 4002 FORMAT(1X,'MAX ALTITUDE(KM) ',F6.0) WRITE(16,4003) 4003 FORMAT(1X,'CONTOUR LEVELS ') WRITE(16,4004) (CONTUR(I),I=1,NUMV) 4004 FORMAT(1X,E10.3) JUMP=ISP+1 GOTO (4010,4020,4030,4040),JUMP 4010 WRITE(16,4015) 4015 FORMAT(1X,'NO CONTOUR CONTROLS') GOTO 4100 4020 WRITE(16,4025) 4025 FORMAT(1X,'SATELLITE CONTOUR CONTROLS') GOTO 4100 4030 WRITE(16,4035) 4035 FORMAT(1X,'F-PEAK CONTOUR CONTROLS') GOTO 4100 4040 WRITE(16,4045) 4045 FORMAT(1X,'BOTH SATELLITE AND F-PEAK CONTOUR CONTROLS') 4100 JUMP=IX+1 GOTO(4110,4120,4130,4140),JUMP 4110 WRITE(16,4115) 4115 FORMAT(1X,'X-AXIS - LATITUDE') GOTO 4200 4120 WRITE(16,4125) 4125 FORMAT(1X,'X-AXIS - TIME(SECS)') GOTO 4200 4130 WRITE(16,4135) 4135 FORMAT(1X,'X-AXIS - TIME(MINS)') GOTO 4200 4140 WRITE(16,4145) 4145 FORMAT(1X,'X-AXIS - LONGITUDE') 4200 JUMP=IY+1 GOTO(4210,4220),JUMP 4210 WRITE(16,4215) 4215 FORMAT(1X,'Y-AXIS - ALTITUDE') GOTO 4300 4220 WRITE(16,4225) 4225 FORMAT(1X,'Y-AXIS - DISTANCE FROM SATELLITE') 4300 JUMP=ISH+1 GOTO(4310,4320,4330,4340),JUMP 4310 WRITE(16,4315) 4315 FORMAT(1X,'NO CONTOURS') GOTO 4400 4320 WRITE(16,4325) 4325 FORMAT(1X,'FLAT EARTH CONTOURS') GOTO 4400 4330 WRITE(16,4335) 4335 FORMAT(1X,'CURVED EARTH CONTOURS') GOTO 4400 4340 WRITE(16,4345) 4345 FORMAT(1X,'BOTH CONTOURS PLOTTED') 4400 JUMP=ITYP+1 GOTO(4410,4420,4430,4440,4450,4460,4470),JUMP 4410 WRITE(16,4415) 4415 FORMAT(1X,'ELECTRON DENSITY CONTOURED') GOTO 4500 4420 WRITE(16,4425) 4425 FORMAT(1X,'FN CONTOURED') GOTO 4500 4430 WRITE(16,4435) 4435 FORMAT(1X,'FN/FH CONTOURED') GOTO 4500 4440 WRITE(16,4445) 4445 FORMAT(1X,'FX CONTOURED') GOTO 4500 4450 WRITE(16,4455) 4455 FORMAT(1X,'FZ CONTOURED') GOTO 4500 4460 WRITE(16,4465) 4465 FORMAT(1X,'FT CONTOURED') GO TO 4500 4470 WRITE(16,4475) 4475 FORMAT(1X,'FH CONTOURED') 4500 WRITE(16,4510)SANG,EANG 4510 FORMAT(1X,'START ANGLE ',F6.2,4X,'END ANGLE ',F6.2) C NOW CALCULATE THE MAXIMUM SCALING FACTOR AND CORRESPONDING CENTER C OF EARTH POINT. THIS WILL GIVE THE LARGEST POSSIBLE CURVED EARTH C CONTOUR WITH X AND Y SCALES EQUAL. X1=1793.5/((6371.2+HMAX)*SIN(ABS(EANG-SANG)/(2*57.29578))) X2=2525.5/((6371.2+HMAX)-6371.2*COS((EANG-SANG)/(2*57.29578))) X3=AMIN1(X1,X2)*0.9 Y=355.-X3*6371.2*COS((EANG-SANG)/(2*57.29578)) RETURN END SUBROUTINE FIELDG (DLAT,DLONG,ALT,TM,NMX,L,X,Y,Z,F) EQUIVALENCE (SHMIT(1,1),TG(1,1)) COMMON/FLDCOM/ST,CT,SPH,CPH,R,BT,BP,BR,B,NMAX COMMON /COEFFS/TG(18,18) DIMENSION G(18,18),GT(18,18),SHMIT(18,18),AID(15) IF (A.EQ.6378.165) IF(L)19,1,2 DATA MAXN,TZERO,J,K/10,1960.,0,0/,G/ X0.0,-30424.9,-1536.1,1300.9,957.6,-227.7,49.8,70.9,4.8,9.9,8*0.,57 X74.8,-2161.6,3000.2,-1987.0,802.8,359.5,60.7,-57.2,6.7,2.9,8*0.,-1 X949.8,204.3,1585.3,1290.4,502.6,231.3,4.5,5.6,-8.8,7.4,8*0.,-431.0 X,230.8,-130.0,871.2,-394.0,-31.2,-241.7,7.5,-13.8,-15.6,8*0.,152.0 X,-268.4,2.9,-250.5,271.4,-157.3,-1.2,-24.4,-3.3,11.4,8*0.,8.6,121. X2,-116.0,-110.4,79.9,-65.2,0.5,-1.5,7.1,11.1,8*0.,-11.9,102.8,60.9 X,-27.2,-12.4,-11.6,-109.1,14.1,-5.6,1.,8*0.,-54.,-24.4,-9.1,2.2,27 X.6,-21.1,-20.1,5.8,11.7, 9*0.,6.9,-12.2,5.8,-17.0,2.6,23.6,-2.5,-1 X6.0,6.4,1.6,8*0.,-22.0,15.6,5.1,-3.5,-1.8,9.6,12.1,0.2,-2.5,1.5,15 X2*0./ DATA GT/ X0.0,20.59,-29.07,2.66,-0.86,2.55,-0.70,11*0.,-3.94,6.02,1.21,-10.0 X3,1.94,-0.08,0.99,11*0.,-13.69,-15.78,-0.70,1.63,-1.17,1.53,0.85,1 X1*0.,6.49,2.93,-9.24,-1.30,-0.54,-0.43,2.11,11*0.,1.77,-1.54,3.18, X-5.48,-4.17,-0.72,1.57,11*0.,3.04,2.88,-1.86,1.25,0.80,1.64,-0.09, X11*0.,-1.39,0.12,1.53,-0.73,-0.06,0.45,0.06,209*0. / A=6378.165 FLAT=1.-1./298.3 A2=A**2 A4=A**4 B2=(A*FLAT)**2 A2B2=A2*(1.-FLAT**2) A4B4=A4*(1.-FLAT**4) IF (L)14,14,2 1 IF (TM-TLAST) 17,19,17 2 READ (3,3) J,K,TZERO,AID ! WAS UNIT 5 3 FORMAT (2I1,1X,F6.1,15A4,A3) L=0 WRITE (16,4) J,K,TZERO,AID 4 FORMAT (2I3,5X'EPOCH=',F7.1,5X15A4,A3) MAXN=0 TEMP=0. 5 READ (3,6) N,M,GNM,HNM,GTNM,HTNM ! WAS UNIT 5 6 FORMAT(2I3,6F11.4) IF (N.LE.0) GOTO7 MAXN=(MAX0(N,MAXN)) G(N,M)=GNM GT(N,M)=GTNM TEMP= MAX1(TEMP, ABS(GTNM)) IF (M.EQ.1) GOTO5 G(M-1,N)=HNM GT(M-1,N)=HTNM GO TO 5 7 WRITE (16,8) 8 FORMAT ('0 N M',6X'G',10X'H',9X'GT',9X'HT'//) DO 12 N=2,MAXN DO 12 M=1,N MI=M-1 IF (M.EQ.1) GOTO10 WRITE (16,9) N,M,G(N,M),G(MI,N),GT(N,M),GT(MI,N) 9 FORMAT (2I3,4F11.4) GO TO 12 10 WRITE (16,11) N,M,G(N,M),GT(N,M) 11 FORMAT (2I3,F11.4,11X,F11.4) 12 CONTINUE WRITE (16,13) 13 FORMAT (1H1) IF (TEMP.EQ.0.) L=-1 14 IF (K.NE.0) GOTO17 SHMIT(1,1)=-1. DO 15 N=2,MAXN SHMIT(N,1)=SHMIT(N-1,1)*FLOAT(2*N-3)/FLOAT(N-1) SHMIT(1,N)=0. JJ=2 DO 15 M=2,N SHMIT(N,M)=SHMIT(N,M-1)*SQRT(FLOAT((N-M+1)*JJ)/FLOAT(N+M-2)) SHMIT(M-1,N)=SHMIT(N,M) 15 JJ=1 DO 16 N=2,MAXN DO 16 M=1,N G(N,M)=G(N,M)*SHMIT(N,M) GT(N,M)=GT(N,M)*SHMIT(N,M) IF (M.EQ.1) GOTO16 G(M-1,N)=G(M-1,N)*SHMIT(M-1,N) GT(M-1,N)=GT(M-1,N)*SHMIT(M-1,N) 16 CONTINUE 17 T=TM-TZERO DO 18 N=1,MAXN DO 18 M=1,N TG(N,M)=G(N,M)+T*GT(N,M) IF (M.EQ.1) GOTO18 TG(M-1,N)=G(M-1,N)+T*GT(M-1,N) 18 CONTINUE TLAST=TM 19 SINLA= SIN(DLAT/57.29578) RLONG=DLONG/57.29578 CPH= COS(RLONG) SPH= SIN(RLONG) IF (J.EQ.0) GOTO20 R=ALT+6371.2 CT=SINLA GO TO 21 20 SINLA2=SINLA**2 COSLA2=1.-SINLA2 DEN2=A2-A2B2*SINLA2 DEN= SQRT(DEN2) FAC=(((ALT*DEN)+A2)/((ALT*DEN)+B2))**2 CT=SINLA/ SQRT(FAC*COSLA2+SINLA2) R= SQRT(ALT*(ALT+2.*DEN)+(A4-A4B4*SINLA2)/DEN2) 21 ST= SQRT(1.-CT**2) NMAX=MIN0(NMX,MAXN) CALL FIELD Y=BP F=B IF (J) 22,23,22 22 X=-BT Z=-BR RETURN C TRANSFORMS FIELD TO GEODETIC DIRECTIONS 23 SIND=SINLA*ST- SQRT(COSLA2)*CT COSD= SQRT(1.0-SIND**2) X=-BT*COSD-BR*SIND Z=BT*SIND-BR*COSD RETURN END SUBROUTINE FRTIME(F1,TIME) C SUBROUTINE FRTIME. COMPUTE GMT (DECIMAL HRS) AT WHICH FA OCCURRED C FROM THE SWEEP RATE OF THE SATELLITE. INITIAL FREQUENCY F1 C OCCURRED AT TIME GMT1 (IN DECIMAL HOURS). C ISAT IDENTIFIES THE SATELLITE: 1=AL-1; 2=AL-2; 3=ISIS-1; 4=ISIS-2. C REFERENCE TIMES AT FREQUENCIES FBR ARE STORED, FOR EACH SATELLITE, C IN THE ARRAY REF. THESE ARE BASED ON THE FOLLOWING SWEEP RATES: C ALOUETTE-1: SWEEP RATE = 1.0 MHZ/SEC (0.5-12 MHZ). C ALOUETTE-2: SWEEP RATE = 0.125 MHZ/SEC (0.2-2 MHZ); =1 (2-15 MHZ). C ISIS-1: SWEEP RATE = 0.31 MHZ/SEC (.25-2 MHZ); = 0.88 (2-5 MHZ); C = 1.03 (5-20 MHZ). C ISIS-2: SWEEP RATE = 0.375 MHZ/SEC(.25-2 MHZ); = 1.125 (2-5 MHZ); C = 1.50 (5-20 MHZ). IMPLICIT REAL*8(A-H,O-Z) COMMON /MISC/ DP,FH1,ASUM,BSUM,FA,CPH,X1,X2,GMT1 COMMON /SWJTCH/ ISAT,LOCSAT,IFXS,NCASES,ICODE(4),N,METHOD,IORD,IP REAL*8 FBR(4)/0.,2.,5.,20./, REF(16)/0.,2.,5.,20.,0.,16.,19.,34., R 0.,6.45,9.86,24.41,0.,5.34,8.,18./, REFT(4) K=4*(ISAT-1) DO 100 I=1,4 100 REFT(I)=REF(K+I) DELT=TAB(FBR,REFT,4,1,FA)-TAB(FBR,REFT,4,1,F1) TIME=GMT1+DELT/3.6D3 RETURN END SUBROUTINE CONT(A,LA,A2,LA2,ISP,X3,Y,IHDR,TITX,TITY,TYPE,TITL, 1 ISH,XCO,LXCO,HMAX,SANG,EANG,LBL,IX,IY) C SUBROUTINE CONT DOES THE ACTUAL CONTOURING OF DATA. IT IS PASSED C MOST OF ITS DATA BY PLOT1 AND CONTROL DATA FROM READ4. COMMON /FLIN/FLINES(2,3,180),NLINES DIMENSION A(2500,3),ANGLE(1),Z(80,80),A2(80,3),XCO(LXCO) C DIMENSION B2(80,3) DATA UNDE /Z77777777/ LOGICAL*1 TITL(80) COMPLEX*16 TITX,TITY,TYPE COMMON /OPTION/IMIN,IMAX,JMIN,JMAX,FINTER,NUMV,CONTUR(50) COMMON /MID/XMID,YMID C SET THE DEFAULT STORAGE VALUE AS THE UNDEFINED VALUE FOR CONTOURING. c CALL UNDEF(UNDE) C C ATTEMPT TO ENTER FUDGE FACTOR TO MOVE ASTERISK C C C NOW WE CALCULATE THE NUMBERS OF ROWS AND COLUMNS IN THE ARRAY. C=LXCO L=HMAX/500 D=L*10 C CALCULATE THE LOCATION OF FIELD LINE INFORMATION FOR THIS X-AXIS C TYPE(LATITUDE OR LONGITUDE ONLY). IFL=IX/3+1 C IF ISH=0(NO CONTOURS REQUESTED) GO TO END. IF(ISH.EQ.0) GO TO 200 C CONVERT INPUT ARRAY TO ROWS AND COLUMNS. do jq=1,2500 write(12,*) (a(jq,jw),jw=1,3) enddo CALL CONVRT(A,LA,2500,XCO(1),XCO(LXCO),50.,HMAX,1.,C,1.,D) X0=XCO(1)-.1*(XCO(LXCO)-XCO(1)) X1=XCO(LXCO)+.1*(XCO(LXCO)-XCO(1)) L1=HMAX/500. C C CHANGE MADE HERE FOR TICK MARKS 500 USED TO BE 250 C C IF ONLY CURVED EARTH CONTOURS, GO TO APPROPRIATE LOCATION. IF(ISH.EQ.2) GO TO 100 C IF HEADERS REQUESTED, PLOT THEM. c IF(IHDR.NE.0) GO TO 30 c CALL HORLIN(TITL,80,XMID,2850,0,4) c CALL HORLIN(TYPE,16,XMID,2850,0,2) c CALL VERLIN(TITY,-16,300.5,YMID,-12,0) c CALL HORLIN(TITX,16,XMID,320.5,0,-4) c 30 CALL OGRID(X0,X1,L,'F5.0)',1,0.,HMAX,L1,'F5.0)',2,0) C IF SATELLITE OR F-PEAK POINTS REQUESTED, PLOT THEM. c IF(ISP.GT.0) CALL DATAPT(A2,LA2,80,1,'*') IMAX=LXCO JMAX=HMAX/50 C IF FIELD LINE PLOTTING REQUESTED AND LEGAL, PLOT THEM. IF(NLINES.LE.0.OR.IX.EQ.1.OR.IX.EQ.2.OR.IY.EQ.1) GO TO 6 c DO 5 I=1,NLINES c 5 CALL PLOT(FLINES(1,IFL,I),FLINES(1,3,I),2,' ') C CLEAR OUT Z ARRAY, AND PUT VALUES TO BE CONTOURED AT THE GRID POINTS C NEAREST INPUT VALUES. 6 JMAX=HMAX/50 type*,' imin imax jmin jmax ',imin,imax,jmin,jmax DO 10 I=IMIN,IMAX DO 10 J=JMIN,JMAX 10 Z(I,J)=UNDE type*,' unde ',unde type*,' la ',la DO 20 I=1,LA J=A(I,1)+.5 C type*,' j ',j K=A(I,2)+.5 C type*,' k ',k Z(J,K)=A(I,3) C print*,' z j k ',z(j,k),j,k 20 write(13,*) j,k,z(j,k) C DUMP THE Z ARRAY. c CALL ZDUMP(Z,80,80) X2=IMIN-(IMAX-IMIN)*.1 X4=IMAX+(IMAX-IMIN)*.1 Y1=JMIN Y2=JMAX C RESCALE FOR CONTOURS, AND PLOT THEM. c CALL SCALE(X2,X4,Y1,Y2,0) c CALL LABELI(LBL) c CALL TRACER(Z,80,80) c CALL FRMADV C IF ONLY FLAT EARTH REQUESTED, GO TO END. IF(ISH.EQ.1) GO TO 200 100 CONTINUE C NOW HANDLE CURVED-EARTH CONTOURS C COORDINATES MUST ALL BE REVERSED. DO 110 I=1,2500 ST1=A(I,1) A(I,1)=A(I,2) 110 A(I,2)=ST1 DO 120 I=1,80 ST1=A2(I,1) A2(I,1)=A2(I,2) 120 A2(I,2)=ST1 C C BUMP ALL ALTITUDES UP BY 50 KM FOR COSMETIC REASONS C 7/31/81 P.K. PARKER C c DO 122 I=1,80 c 122 A2(I,1)=A2(I,1)+50. C C ARRAY B2 IS A2 'FUDGED' TO MOVE THE ASTERISK SLIGHTLY LEFT C I=IMIN IMIN=JMIN JMIN=I I=IMAX IMAX=JMAX JMAX=I C NOW THAT EVERYTHING HAS BEEN REVERSED FOR THE POLAR CONTOURS, SET UP C THE POLAR GRID. C IF HEADERS HAVE BEEN REQUESTED, PLOT THEM. c IF(IHDR.NE.0) GO TO 40 c CALL HORLIN(TITL,80,XMID,2850.5,0,4) c CALL HORLIN(TYPE,16,XMID,2850.5,0,2) c CALL VERLIN(TITY,-16,300.5,YMID,-12,0) c CALL HORLIN(TITX,16,XMID,320.5,0,-4) c 40 CONTINUE C SET UP FOR POLAR (CURVED EARTH) CONTOURS. c CALL POLAR(2) AMAX=ABS(EANG-SANG)/2 AMIN=-AMAX RMINO=X3*6371.2 RMAXO=X3*(6371.2+HMAX) c CALL PSTGRD(RMINO,AMIN,RMAXO,AMAX,2148.5,Y) ANGLE(1)=AMIN-3 L=(EANG-SANG)/5+.5 c CALL PLABEL(ANGLE,1) c CALL OGRID(0.,HMAX,L1,'F5.0)',2,SANG,EANG,L,'F5.0)',1,0) C IF REQUESTED AND LEGAL, PLOT THE FIELD LINES. c IF(NLINES.LE.0.OR.IX.EQ.1.OR.IX.EQ.2.OR.IY.EQ.1) GO TO 7 c DO 130 I=1,NLINES c 130 CALL PLOT(FLINES(1,3,I),FLINES(1,IFL,I),2,' ') C IF REQUESTED, PLOT SATELLITE OR F-PEAK POINTS. c 7 IF(ISP.GT.0) CALL DATAPT(A2,LA2,80,1,'*') C NOW CLEAR Z ARRAY, THEN FILL RECTANGULAR GRID. DO 140 I=IMIN,IMAX DO 140 J=JMIN,JMAX 140 Z(I,J)=UNDE DO 150 I=1,LA C type*,' la second ',la J=A(I,1)+.5 C type*,' j second ',j K=A(I,2)+.5 C type*,' k second ',k 150 Z(J,K)=A(I,3) C DUMP Z-ARRAY. c CALL ZDUMP(Z,80,80) C CALCULATE END VALUES OF X COORDINATES AT START AND END ANGLES. do ij=1,80 WRITE(11,*) (z(ij,ik),ik=1,80) enddo SLOPE=(JMAX-JMIN)/(XCO(LXCO)-XCO(1)) X2=JMIN-SLOPE*(XCO(1)-SANG) X4=JMAX+SLOPE*(EANG-XCO(LXCO)) Y1=IMIN Y2=IMAX C NOW DO ACTUAL CURVED EARTH CONTOURS, AND END THE PLOT PACKAGE. c CALL SCALE(Y1,Y2,X2,X4,0) c CALL LABELI(LBL) c CALL TRACER(Z,80,80) c CALL FRMADV 200 continue c 200 CALL ENDPLT RETURN END SUBROUTINE NEWINVERT(DLAT,DLONG,DATE,HT,LQ,*) C SUBROUTINE INVERT. COMPUTE FX AT SATELLITE IF REQUIRED (IFXS=1). C INVERT H'(F) DATA TO OBTAIN N(H) PROFILE USING THE LAMINATION C METHOD CONTAINED IN SUBROUTINE XLAM. ALLOW THE SATELLITE POSITION C TO VARY IF REQUIRED (LOCSAT=1). IMPLICIT REAL*8(A-H,O-Z) REAL*4 DLAT,DLONG,HT,DATE,HGT,BN,BE,BV,BT REAL*8 HT2I(5),RAD/57.29577951D0/ COMMON /CALC/ AA(50),BB(50),F(50),HP(50),FH(50),TH(70),DEN(70) COMMON /MISC/ DP,FH1,ASUM,BSUM,FA,CPH,X1,X2,GMT1 COMMON /SWJTCH/ ISAT,LOCSAT,IFXS,NCASES,ICODE(4),N,METHOD,IORD,IP COMMON /IO/ IYR,IDY,ITI,LMT,TITLE(40),DIP,FHS common/times/times(50),iyrs(50),imos(50),idays(50),idoys(50) DLT=DLAT DLG=DLONG DEN(1)=(F(1)*(F(1)-FH(1)))*12400.00 F1=F(1) KK=2 IF(IFXS.EQ.0) GO TO 1600 C IF FXSAT IS CALLED, KK=4 AND THE CALCULATION OF THE FIRST TWO C LAMINATIONS IS NOT REPEATED. C DEN1 = INITIAL DENSITY AT ORIGINAL SATELLITE POSITION. HT = HEIGHT KK=4 CALL FXSAT(DLAT,DLONG,DATE,LQ,&100) AA2=1.0/AA(2) FACT=DEN(1)*DEXP(-HT*AA2) 1600 DEN1=DEN(1) FS=F(1) NL = N M=2 IFLAG=0 DO 300 I=KK,N DP=0.0 IF (I .GT. NL) GO TO 300 IF (F(I) .GE. (F(I-1)+0.001)) GO TO 500 C TERMINATE IF F(I) IS NOT MONOTONICALLY INCREASING. WRITE(16,2) F(I) GO TO 100 500 FA= F(I) IF (I .LT. 3 .OR. M .EQ. 3 .OR. LOCSAT .EQ. 0) GO TO 800 C VARY SATELLITE POSITION AS FOLLOWS: COMPUTE GMT IN SEC AT WHICH C F(I) OCCURRED FROM NOMINAL SWEEP RATE OF THE SATELLITE (SUBROUTINE C FRTIME). COMPUTE CORRESPONDING SATELLITE POSITION BY INTERPOLATING C PREVIOUSLY INPUT WORLD MAP DATA (SUBROUTINE ORBINT). CALCULATE FH C AT EACH INPUT H' USING NEW LAT AND LONG IN THE CALL TO FIELD G. C M=3 IF FIRST LAMINATION COLLAPSES. c 900 CALL FRTIME(FS,TIME) 900 time=times(i) iyr=iyrs(i) idy=idoys(i) C CONVERT TIME FROM DEC HOUR TO HHMMSS. C TYPE*,' TIME ',TIME IH=TIME RM=(TIME-IH)*60. IM=RM rs=(rm-im)*60. IS=(RM-IM)*60 C REPLACE CALL TO ORBINT WITH CALL TO PROGRAM TO CALCULATE ORBITAL ELEMENTS DLAT1=DLAT DLONG1=DLONG HGT1=HGT CALL JWORLDMAP(ISAT,IYR,IDY,IH,IM,IS,DLAT1,DLONG1,HGT1,0) c TYPE*,' WORLDMAP from newinvert routnie NOT CALLED. ' c TYPE*,' ENTER LAT AND LONG AND ALT' c ACCEPT*,DLAT1,DLONG1,HGT1 DLAT=DLAT1 DLONG=DLONG1 HGT=HGT1 C CALL ORBINT(1,TIME,DLAT,DLONG,HGT) TH(1)=HGT IF(TH(1).GE.TH(2)) GO TO 1000 WRITE(16,6) TH(1),FA WRITE(16,7) DLAT,DLONG,TH(2) TH(1)=TH(2) M=3 C FIXED SATELLITE ASSUMED FOR REMAINDER OF ANALYSIS. GO TO 800 1000 IFH=I-1 DO 1100 L=1,IFH CALL FIELDG(DLAT,DLONG,HGT,DATE,10,LQ,BN,BE,BV,BT) FH(L)=BT*2.8D-5 IF (I.EQ.NL.AND. L.EQ.1) DIPS=ATAN2(BV,SQRT(BT**2-BV**2))*RAD 1100 HGT=TH(L+1) 1700 DEN(1)=FACT*DEXP(TH(1)*AA2) C INITIAL FH(I) BASED UPON 0.1 PER CENT INCREASE IN DENSITY 800 FH(I)=FA-DEN(I-1)*8.07258870D-5/FA IF(TH(I-1).GT.1.D3) GO TO 600 AAA=1.15*FH(I-1) IF ((AAA-FH(I)) .LT. 0.0) FH(I)=AAA 600 CALL XLAM(M,I,DLAT,DLONG,DATE,LQ,1,&100) IF (I .NE. 2 .OR. LOCSAT .EQ. 0) GO TO 300 C FOR VARIABLE SATELLITE POSITION AND KNOWN FXS, ITERATE ON FIRST C LAMINATION TH(2) UNTIL CHANGE IS LESS THAN OR EQUAL TO 0.01 KM. C IF CALCULATION DOES NOT CONVERGE AFTER FIVE ITERATIONS, PRINT C MESSAGE AND ASSUME FIXED SATELLITE POSITION FOR ANALYSIS. AA2=1.0/AA(2) FACT=DEN1*DEXP(-HT*AA2) IFLAG=IFLAG+1 DP=0.0 HT2I(IFLAG)=TH(2) IF(IFLAG.GT.1) GO TO 1300 HT2=TH(2) GO TO 900 1300 HT1=HT2 HT2=TH(2) IF ((HT2-HT1)**2 .LE. 1.D-4) GO TO 300 IF(IFLAG.NE.5) GO TO 1700 DLAT=DLT DLONG=DLG WRITE(16,3) WRITE (16,7) DLAT,DLONG,HT LOCSAT=0 300 CONTINUE IF (LOCSAT.EQ.0) GO TO 1500 I=IFLAG-1 IF(IFXS.EQ.0) WRITE(16,4)I,(HT2I(L),L=1,IFLAG) WRITE(16,5) F(NL),TH(1),DLAT,DLONG,DIPS,FH(1) C RESET VALUES AT ORIGINAL SATELLITE POSITION. 1500 FH(1)=FH1 F(1)=F1 TH(1)=HT DEN(1)=DEN1 N = NL DLAT=DLT DLONG=DLG RETURN C USE RETURN 1 IF PROBLEMS ENCOUNTERED IN XLAM OR FXSAT, OR IF C INPUT DATA IS NOT MONOTONICALLY INCREASING. 100 RETURN 1 2 FORMAT(' DATA FOR F=',F6.3,' NOT MONOTONICALLY INCREASING. PROCES FSING TERMINATED.') 3 FORMAT(' CALCULATION OF FIRST LAMINATION WITH VARIABLE POSITION DI FD NOT CONVERGE.') 4 FORMAT(1x,I2,' ITERATIONS REQUIRED TO COMPUTE INITIAL SLOPE.'/ F 1H ,'SUCCESSIVE HEIGHTS OF SECOND POINT:',4F8.2) 5 FORMAT(1x,'SATELLITE POSITION AT FX = ',F6.3,' MHZ.'/' HT=',F6.1 F ,' LAT=',F7.2,' LONG=',F7.2,' DIP=',F7.2,' FHS=',F5.3) 6 FORMAT(' SATELLITE (',F6.1,' KM) BELOW FIRST LAMINATION AT',F7.3, F ' MHZ ON X-TRACE.') 7 FORMAT(' FIXED POSITION ASSUMED: LAT =',F8.2,' LONG =',F8.2, F ' HT =',F7.1) END FUNCTION TAB(XT,YT,N,K,X) IMPLICIT REAL*8(A-H,O-Z) DIMENSION XT(100),YT(100),X1(6),Y1(6) C THIS ROUTINE USES AITKEN S METHOD OF INTERPOLATION. FOR A GIVEN C INDEPENDENT VARIABLE X, AN INTERPOLATION OF ORDER K IS PERFORMED. C THE CALLING SEQUENCE IS Y = TAB(XT,YT,N,K,X) WHERE: C XT = ARRAY OF INDEPENDENT VALUES IN A MONOTONIC SEQUENCE. C YT = ARRAY OF DEPENDENT VALUES. C N = NUMBER OF INDEPENDENT VALUES. C K = ORDER OF INTERPOLATION (MUST BE LESS THAN 6). C X = INDEPENDENT VALUE FOR WHICH INTERPOLATION IS DESIRED. NPTS=K+1 C CHECK FOR ASCENDING OR DESCENDING SEQUENCE IF (XT(1)-XT(N))5,5,1 C SEQUENCE IS DESCENDING, CHECK IF EXTRAPOLATION IS NEEDED 1 IF(X-XT(1))2,2,10 2 IF(X-XT(N))30,3,3 3 DO 4 I=1,N IF(X-XT(I))4,4,50 4 CONTINUE C SEQUENCE IS ASCENDING, CHECK IF EXTRAPOLATION IS NEEDED 5 IF(X-XT(1))10,10,20 C SET UP FOR EXTRAPOLATION AT BEGINNING OF TABLE 10 DO 15 I=1,NPTS X1(I)=XT(I) 15 Y1(I)=YT(I) GO TO 100 20 IF(X-XT(N))40,30,30 C SET UP FOR EXTRAPOLATION AT END OF TABLE 30 DO 35 I=1,NPTS NN=N-NPTS+I X1(I)=XT(NN) 35 Y1(I)=YT(NN) GO TO 100 C INTERPOLATION NEEDED, FIND THE TWO TABLE VALUES THAT THE C INDEPENDENT VARIABLE,X,LIES BETWEEN 40 DO 45 I=1,N IF(X-XT(I))50,45,45 45 CONTINUE I=N C CHOOSE K+1 INDEPENDENT TABLE VALUES CLOSEST TO X 50 IF(MOD(K,2))51,51,60 51 IF(DABS(X-XT(I-1))-DABS(XT(I)-X)) 52,56,56 52 IF(I-2-K/2)10,53,53 53 IF(I-1+K/2-N)54,54,30 54 DO 55 J=1,NPTS NN=I-2-K/2+J X1(J)=XT(NN) 55 Y1(J)=YT(NN) GO TO 100 56 IF(I-1-K/2)10,57,57 57 IF(I+K/2-N)58,58,30 58 DO 59 J=1,NPTS NN=I-K/2+J-1 X1(J)=XT(NN) 59 Y1(J)=YT(NN) GO TO 100 60 IF(I-1+(K+1)/2-N)61,61,30 61 IF(I-2-(K-1)/2)10,62,62 62 DO 63 J=1,NPTS NN=I-1-(K+1-2*J)/2 X1(J)=XT(NN) 63 Y1(J)=YT(NN) 100 DO 110 II=1,K SAVEY1=Y1(1) L=II+1 DO 110 JJ=L,NPTS NN=JJ-II LL=NN+1 110 Y1(NN)=(SAVEY1*(X1(JJ)-X)-Y1(LL)*(X1(II)-X))/(X1(JJ)-X1(II)) TAB=Y1(1) RETURN END SUBROUTINE TRACE(SLAT,SLONG,SHGT,RLAT,RLONG,RHGT,DIST,DATE,IHEM) COMMON /STEPP/LCT,Y,YOLD,BN,BE,BV,B,ST,SGN,D DATA RAD/57.2958/ DIMENSION Y(3),YOLD(3) Y(1)=SHGT Y(2)=SLAT Y(3)=SLONG L=0 SGN=IHEM D=DIST/10. ST=COS(SLAT/RAD) DO 10 LCT=1,14 CALL FIELDG(Y(2),Y(3),Y(1),DATE,10,L,BN,BE,BV,B) 10 CALL STEP RLAT=Y(2) RLONG=Y(3) RHGT=Y(1) RETURN END SUBROUTINE DAYMON (IYEAR,IDAYNO,MON,IDAY) C ROUTINE TO COMPUTE MONTH DAY, GIVEN YEAR AND DAY OF YEAR (IDAYNO) DIMENSION NMON(4),K(4) DATA K/28,29,30,31/ YEAR=IYEAR DO 10 I=1,4 10 NMON(I)=0 IDM=0 IDT=0 C NMON(I) COUNTS THE NUMBERS OF DIFFERENT-LENGTH MONTHS UP TO THAT I C WHICH IDAYNO OCCURS. I=1 FOR 28 DAY,2 FOR 29,3 FOR 30,4 FOR 31. C IDM IS SET TO THE NUMBER OF DAYS IN THE ITH MONTH. IDT IS THE TOTA C NUMBER OF DAYS, UP TO THE ITH MONTH, CONSIDERED SO FAR. DO 20 I=1,12 IF ( I .NE. 2) GO TO 30 ID=2 IF (AMOD(YEAR,4.) .NE. 0.0) ID=1 C IF ID=1, THIS IS NOT A LEAP YEAR. NMON(ID)=1 IDM=K(ID) GO TO 60 30 IF (I .EQ. 4 .OR. I .EQ. 6 .OR. I .EQ. 9 .OR. I .EQ. 11) GO TO 50 C THIS MONTH HAS 31 DAYS. IDM=31 NMON(4)=NMON(4)+1 GO TO 60 C THIS MONTH HAS 30 DAYS. 50 IDM=30 NMON(3)=NMON(3)+1 60 IDT=IDT+IDM IF (IDAYNO-IDT) 70,70,20 70 IDAY=IDAYNO MON=0 DO 80 J=1,4 IDAY=IDAY-NMON(J)*K(J) 80 MON=MON+NMON(J) IDAY=IDAY+IDM RETURN 20 CONTINUE RETURN END SUBROUTINE FIELD COMMON/FLDCOM/ST,CT,SPH,CPH,R,BT,BP,BR,B,NMAX COMMON /COEFFS/G(18,18) DIMENSION P(18,18),DP(18,18),CONST(18,18),SP(18),CP(18),FN(18),FM( 118) IF (P(1,1).EQ.1.0) GO TO 3 1 P(1,1)=1. DP(1,1)=0. SP(1)=0. CP(1)=1. DO 2 N=2,18 FN(N)=N DO 2 M=1,N FM(M)=M-1 2 CONST(N,M)=FLOAT((N-2)**2-(M-1)**2)/FLOAT((2*N-3)*(2*N-5)) 3 SP(2)=SPH CP(2)=CPH DO 4 M=3,NMAX SP(M)=SP(2)*CP(M-1)+CP(2)*SP(M-1) 4 CP(M)=CP(2)*CP(M-1)-SP(2)*SP(M-1) AOR=6371.2/R AR=AOR**2 BT=0. BP=0. BR=0. DO 8 N=2,NMAX AR=AOR*AR DO 8 M=1,N IF(N.EQ.M) GO TO 5 6 P(N,M)=CT*P(N-1,M)-CONST(N,M)*P(N-2,M) DP(N,M)=CT*DP(N-1,M)-ST*P(N-1,M)-CONST(N,M)*DP(N-2,M) GO TO 7 5 P(N,N)=ST*P(N-1,N-1) DP(N,N)=ST*DP(N-1,N-1)+CT*P(N-1,N-1) 7 PAR=P(N,M)*AR IF (M.EQ.1) GO TO 9 TEMP=G(N,M)*CP(M)+G(M-1,N)*SP(M) BP=BP-(G(N,M)*SP(M)-G(M-1,N)*CP(M))*FM(M)*PAR GO TO 10 9 TEMP=G(N,M)*CP(M) BP=BP-(G(N,M)*SP(M))*FM(M)*PAR 10 BT=BT+TEMP*DP(N,M)*AR 8 BR=BR-TEMP*FN(N)*PAR BP=BP/ST B= SQRT(BT*BT+BP*BP+BR*BR) RETURN END SUBROUTINE XLAM(M,I,DLAT,DLONG,DATE,LQ,IFLAG,*) C SUBROUTINE XLAM: LAMINATION PROCEDURE FOR INVERTING X-TRACE H'(F) C DATA TO OBTAIN N(H) PROFILE. FOR PARABOLIC IN LOG(N) LAMINATIONS, C THE EQUATION USED FOR THE JTH INTERVAL IS: C H = H(J-1) + AA(J)*LN(N/N(J-1)) + BB(J)*LN(N/N(J-1))**2 C WHERE AA(I),BB(I) ARE THE CORRESPONDING SLOPES. FOR LINEAR IN C LOG(N) LAMINATIONS,BB(I)=0. INTEGRATE BY USING GAUSSIAN QUADRATURE C WITH 3 COEFFICIENTS (EXCEPTION: 7 COEFFICIENTS USED WHEN Y2>0.93). C THE CHANGE OF VARIABLE: C T**2 = 1 -X/(1-YR) C IS USED TO KEEP THE INTEGRAND FINITE AT ALL POINTS. C IFLAG=1 FOR USUAL COMPUTATIONS; =2 FOR COMPUTATIONS OF FXSMAX; C =3 FOR COMPUTATIONS OF FXS WHICH YIELDS AA(3)=AA(2). C PROCESSES ONE H'(F) POINT PER CALL. RESULTS RETURNED VIA COMMON. IMPLICIT REAL*8(A-H,O-Z) REAL*4 DLAT,DLONG,HGT,DATE,BN,BE,BV,BT COMMON /CALC/ AA(50),BB(50),F(50),HP(50),FH(50),TH(70),DEN(70) COMMON /MISC/ DP,FH1,ASUM,BSUM,FA,CPH,X1,X2,GMT1 COMMON /SWJTCH/ ISAT,LOCSAT,IFXS,NCASES,ICODE(4),N,METHOD,IORD,IP DIMENSION SS(10),W(10),DIFF(20) DATA SS/.77459667D0,0.D0,-.77459667D0,.94910791D0,.74153119D0, D .40584515D0,0.D0,-.40584515D0,-.74153119D0,-.94910791D0/ DATA W/.55555556D0,.88888889D0,.55555556D0,.12948497D0,.27970539D0 D ,.38183005D0,.41795918D0,.38183005D0,.27970539D0,.12948497D0/ DATA I1,I2/1,2/ IF(IFXS.EQ.1 .AND. I.EQ.4) AAI=AA(3) MM=M IQ=0 ILIM=2 IF (IFLAG.NE.1) ILIM=3 FA2=1.0/(1.24D4*FA**2) FA1=1.0/FA GIRO=FH(I) TH(I) =TH(I-1) HGT=TH(I-1) 200 HGTA = TH(I) YR=FH(I)*FA1 YRM1=1.0-YR RYR=1.0/YRM1 DEN(I) = (F(I)*(F(I)-FH(I)))*12400.00 IF(DEN(I).GT.(DEN(I-1)*1.0005)) GO TO 2600 J=I-1 IF(IFLAG.NE.3) WRITE(16,1)I,DEN(I),J,DEN(J) GO TO 2500 2600 DO 300 J=M,I Y1=FH(J-1)*FA1 Y2=FH(J)*FA1 IL=1 IU=3 IF (Y2.LT. 0.93) GO TO 100 IL=4 IU=10 100 X1=DEN(J-1)*FA2 X2=DEN(J)*FA2 X1LN=DLOG(X1) X2LN=DLOG(X2) C YR=Y AT REFLECTION POINT. CHANGE FROM VARIABLE X TO T.USE B9 AND C C9 TO CHANGE Y WITHIN EACH LAMINATION. A=X1*RYR IF(A.LT.1.D0) GO TO 2800 WRITE(16,3) I1,A,J,I GO TO 2500 2800 T1=DSQRT(1.D0-A) S1 = 0.0 S2=0.0 B9 = DLOG(Y2/Y1)/3.0 V=X2LN-X1LN C CALCULATION FOR LAST LAMINATION. IF (J .LT. I) GO TO 400 T2 = 0.0 IF (METHOD .GT. 1 .AND. I .GT. ILIM .AND. IQ .GE. 2) GO TO 500 C LAMINATIONS LINEAR IN LOG(N). 800 ALIM = 1.0 RA = 1.0 RB = 0.0 GO TO 600 C LAMINATIONS PARABOLIC IN LOG(N). 500 AA(I) = AAI BB(I) = (TH(I)-TH(I-1)-AA(I)*V)/(V*V) AA(I+1)=AA(I)+(2.0)*BB(I)*V ALIM = AA(I+1)/AA(I) GO TO 700 400 A=X2*RYR IF(A.LT.1.D0) GO TO 2900 WRITE(16,3) I2,A,J,I GO TO 2500 C CALCULATION FOR ALL LAMINATIONS EXCEPT THE LAST. 2900 T2=DSQRT(1.D0-A) IF (METHOD .EQ. 1) GO TO 800 ALIM = 1.0 700 RB = BB(J)/AA(J) RA = 1.0+RB*V 600 C9 =(DEXP(B9)-1.0)/(V*RA) T1=0.5*(T1+T2) T2=T1-T2 C GAUSSIAN INTEGRATION DO 900 K=IL,IU T=T2*SS(K)+T1 X=(1.0-T**2)*YRM1 RX=1.0/(1.0-X) A=X2LN-DLOG(X) Y = Y2/((1.0+C9*A*(1.0+RB*(2.0*V-A)))**3) IF (T .GT. 0.0002) GO TO 1000 C UX IS VALUE OF INTEGRAND AT THE REFLECTION POINT 2400 UX=.707107*(2.0-YR)*RYR**2/(DSQRT(1.0+CPH**2)*DSQRT(1.0+3.0*C9*YR U *ALIM*RYR)) GO TO 1100 1000 YL = Y * CPH YT2=Y**2-YL**2 IF (CPH .LT. 0.0009) GO TO 1200 R=0.5*YT2*RX/YL SPSI=DSIN(DATAN(R)) A=(1.0+SPSI)*DSQRT(1.0+R**2) SX = 1.0-A*YL GO TO 1300 C PROPAGATION ANGLE NEAR 90 DEGREES (TRANSVERSE PROPAGATION). 1200 SPSI=1.0D0 SX=1.0-RX*Y**2 1300 A=X/SX IF (A.GE.1.0D0) GO TO 1400 IF(A.GE.0.9999999D0) GO TO 2400 C REJECTS BAD DATA CARDS UX=T*(1.+.5*A*(1.-SX)*(1.+(1.+X)*SPSI*RX)/SX)/(X*DSQRT(1.-A)) 1100 S1=S1+UX*W(K) 900 S2=S2+UX*W(K)*(DLOG(X)-X1LN) BSUM=2.0*T2*YRM1 ASUM=BSUM*S1 BSUM=2.0*BSUM*S2 IF (J .GE. I) GO TO 1500 300 DP=DP-AA(J)*ASUM-BB(J)*BSUM 1500 IF (IFLAG .EQ. 2) RETURN C ITERATE ON LAST LAMINATION USING IMPROVED TECHNIQUE. IF (I .GT. 2) GO TO 1600 AA(2)=-HP(2)/ASUM TH(2)=TH(1)+AA(2)*V AAI = AA(2) GO TO 1800 1600 IF (METHOD .EQ. 2 .AND. IQ .GE. 2) GO TO 1900 AA(I)=(DP-HP(I))/ASUM TH(I)=TH(I-1)+AA(I)*V BB(I)=0.0 HGT=TH(I) IF (METHOD-1) 1800,1800,2000 1900 BB(I)=(DP-AA(I)*ASUM-HP(I))/BSUM TH(I)=TH(I-1)+V*(AA(I)+BB(I)*V) 1800 HGT =TH(I) C CONVERGENCE TEST TO INSURE THAT SUCCESSIVE HT. ARE WITHIN 0.01 KM IF ((HGT-HGTA)**2 .LE. 0.0001) GO TO 2100 2000 IQ=IQ+1 C ITERATE ONLY FOR LAST LAMINATION (M=I),10 ITERATIONS MAX. NO. IF(IQ.LE.9) GO TO 3000 WRITE(16,2) I,HGT,HGTA GO TO 2100 3000 M=I CALL FIELDG(DLAT,DLONG,HGT,DATE,10,LQ,BN,BE,BV,BT) DIFF(IQ)=FH(I)-BT*2.8000E-05 IF (IQ .GT. 1) GO TO 2200 FH(I)=BT*2.8000E-05 TEST=GIRO IF(FH(I).GT.GIRO) TEST=FH(I) GO TO 200 2200 FH(I)=GIRO-DIFF(1)*(FH(I)-GIRO)/(DIFF(IQ)-DIFF(1)) IF (FH(I) .GT. FH(I-1)) GO TO 2300 FH(I)=FH(I-1) + 0.001 GO TO 200 2300 IF (FH(I) .LT. TEST) GO TO 200 FH(I)=TEST-0.0001 GO TO 200 2100 AAI=AAI+2.0*BB(I)*V DEN(I)=(F(I)*(F(I)-FH(I)))*12400.0 IF(DEN(I).GT.(DEN(I-1)*1.0005)) GO TO 2700 J=I-1 IF(IFLAG.NE.3) WRITE(16,1)I,DEN(I),J,DEN(J) GO TO 2500 2700 M=MM RETURN C USE RETURN 1 IF BAD DATA FOUND. 1400 WRITE(16,4) J,I,X,SX,A 2500 M=MM IF(IFXS.EQ.0) WRITE(16,5) RETURN 1 1 FORMAT(' DEN(',I2,') = ',1PD10.3,' LESS THAN DEN(',I2,') = ', F 1PD10.3/' DENSITIES NOT MONOTONICALLY INCREASING. PROCESSING OF FIONOGRAM TERMINATED.') 2 FORMAT(' ITERATION DID NOT CONVERGE FOR I =',I3,'. LAST TWO HEIGHT FS:',2F8.2) 3 FORMAT(' X',I1,'/(1-YR) =',1PD11.3,' > 1. (LAMINATION',I3,' POINT' F ,I3,'). PROCESSING TERMINATED.') 4 FORMAT(' FOR LAMINATION',I3,' POINT',I3,' X =',1PD11.3,' 1-Y** F2/(1-X) =',D11.3/' RATIO =',D11.3,'>1. CALCULATIONS TERMINATED.') 5 FORMAT(' PROCESSING MAY BE POSSIBLE USING FXSAT SUBROUTINE (ICODE( F1)=2).') END SUBROUTINE STEP COMMON/STEPP/LCT,Y,YOLD,BN,BE,BV,B,ST,SGN,D DATA RAD/57.29578/ DIMENSION Y(3),YOLD(3),YP(3,4) C FIELD LINE INTEGRATION PROGRAM, SIMILAR TO FAP VERSION OF GL AIDE C SHARE PGM NO. 413). USES BOOTSTRAP STARTING PROCEDURE, THEN ADAMS C FORMULA, SET LCT = 1 TO START, THEN ADVANCE LCT BY 1 EACH TIME THR C Y(3) ARE INITIAL COORDINATES - GEOCENTRIC DISTANCE, LATITUDE AND L C IN RADIANS. BR,ETC, ARE FIELD COMPONENTS AT Y(3). ST = SIN THETA O C LATITUDE AT Y(3). STEP GENERATES NEW Y(3) AND RETAINS PREVIOUS POI C YOLD(3). IF SGN = +1., INTEGRATION PROCEEDS OPPOSITE TO FLD DIRECT C I.E., NORTH TO SOUTH IN EARTH,S FIELD. SGN = -1. REVERSES DIRECTIO C D = STEP SIZE ALONG FLD LINE. Y(1) AND D IN SAME UNITS. 3 STEPS AR C IN 7 ITERATIONS, THEN ONE STEP EACH ADDITIONAL ITERATION, YP(3,4) C DIFFERENTIALS AT Y(3) AND 3 PREVIOUS POINTS. YP(1,4)=SGN*BV/B BH=B*(Y(1)+6371.2) YP(2,4)=-SGN*RAD*BN/BH YP(3,4)=-SGN*RAD*BE/(BH*ST) IF(LCT.GT.7) GO TO 9 C BOOTSTRAPPING PROCEDURE FOR FIRST 7 ITERATIONS, 3 STEP DO 8 I=1,3 GO TO (1,2,3,4,5,6,7),LCT 1 D2=D/2.0 D6=D/6.0 D12=D/12.0 D24=D/24.0 YP(I,1)=YP(I,4) YOLD(I)=Y(I) Y(I)=YOLD(I)+D*YP(I,1) GO TO 8 2 YP(I,2)=YP(I,4) Y(I)=YOLD(I)+D2*(YP(I,2)+YP(I,1)) GO TO 8 3 Y(I)=YOLD(I)+D6*(2.*YP(I,4)+YP(I,2)+3.*YP(I,1)) GO TO 8 4 YP(I,2)=YP(I,4) YOLD(I)=Y(I) Y(I)=YOLD(I)+D2*(3.*YP(I,2)-YP(I,1)) GO TO 8 5 Y(I)=YOLD(I)+D12*(5.*YP(I,4)+8.*YP(I,2)-YP(I,1)) GO TO 8 6 YP(I,3)=YP(I,4) YOLD(I)=Y(I) Y(I)=YOLD(I)+D12*(23.*YP(I,3)-16.*YP(I,2)+5.*YP(I,1)) GO TO 8 7 Y(I)=YOLD(I)+D24*(9.*YP(I,4)+19.*YP(I,3)-5.*YP(I,2)+YP(I,1)) 8 CONTINUE RETURN C 8TH ITERATION AND BEYOND -ADAMS 4-POINT FORMULA 9 DO 10 I=1,3 YOLD(I)=Y(I) Y(I)=YOLD(I)+D24*(55.*YP(I,4)-59.*YP(I,3)+37.*YP(I,2)-9.*YP(I,1)) DO 10 J=1,3 10 YP(I,J)=YP(I,J+1) RETURN END SUBROUTINE FXSAT(DLAT,DLONG,DATE,LQ,*) C SUBROUTINE FXSAT. DETERMINE FX AT SATELLITE TO GIVE IDENTICAL C SLOPES FOR THE FIRST TWO LAMINATIONS, ASSUMING LINEAR IN LOG(N) C LAMINATIONS. IMPLICIT REAL*8(A-H,O-Z) REAL*4 DLAT,DLONG,DATE,HGT,BN,BE,BV,BT LOGICAL LOWN COMMON /CALC/ AA(50),BB(50),F(50),HP(50),FH(50),TH(70),DEN(70) COMMON /MISC/ DP,FH1,ASUM,BSUM,FA,CPH,X1,X2,GMT1 COMMON /SWJTCH/ ISAT,LOCSAT,IFXS,NCASES,ICODE(4),N,METHOD,IORD,IP DIMENSION FXS(10),DIFF(10),PCNT(10) DATA PCNT / .90,.92,.94,.96,.97,.98,.99,.993,.996,.998 / C BEGIN BY FINDING MAXIMUM PERMISSIBLE VALUE OF FXS BASED ON A C NEARLY VERTICAL PROFILE AT THE SATELLITE. FOR A REGULAR IONOGRAM, C FIRST TRY TEN SELECTED PERCENTAGES OF F(2), AND FIND THE TWO C WHICH YIELD CLOSEST AGREEMENT BETWEEN RESULTING FH(2) AND ACTUAL C FH(2). SEARCH STOPS BEFORE TEN VALUES HAVE BEEN TRIED IF THESE C DIFFERENCES CHANGE SIGN OR BECOME NON-MONOTONIC. NO ITERATIONS ARE C DONE IN LAMINATION CALCULATIONS. FOR LOW DENSITY IONOGRAMS C (FHS > .89*F(2)), USE TEN VALUES OF FXS EQUALLY SPACED BETWEEN FHS C AND F(2). METH=METHOD METHOD=1 IF(F(3).LT.(F(2)+0.001)) GO TO 2000 FHS=FH1 L=2 I=2 FA=F(2) F1MAX=FA-0.0011 TEST=(FA-FHS)/11.0 IF(TEST.LT.0.0011) GO TO 1500 LOWN=.FALSE. IF (FHS.GE.(0.89*FA)) LOWN=.TRUE. DO 100 M=1,10 FXS(M)=PCNT(M)*FA IF (LOWN) FXS(M)=FHS+M*TEST IF(FXS(M).GT.F1MAX) GO TO 800 FH(2)=FA-1.001*FXS(M)*(FXS(M)-FHS)/FA DEN(1)=FXS(M)*(FXS(M)-FHS)*12400.0 DEN(2)=FA*(FA-FH(2))*12400.0 CALL XLAM(L,I,DLAT,DLONG,DATE,LQ,2,&2000) HGT=TH(1)-(HP(2)*DLOG(X2/X1))/ASUM CALL FIELDG(DLAT,DLONG,HGT,DATE,10,LQ,BN,BE,BV,BT) DIFF(M)=(BT*2.8000E-05) - FH(2) TH(2)=HGT K=M DA=DABS(DIFF(M)) IF (M-2) 300,400,500 300 MA=1 FC=DA GO TO 100 500 IF(DIFF(M).LT.DIFF(M-1)) GO TO 800 C FIND THE TWO SMALLEST VALUES OF DA (FC,FB), AND CORRESPONDING C VALUES OF M (MA,MB). IF (DA.GE.FB .OR. DA.EQ.FC) GO TO 600 400 IF (DA.LT.FC) GO TO 700 FB=DA MB=M GO TO 600 700 FB=FC FC=DA MB=MA MA=M 600 IF((DIFF(M)*DIFF(M-1)) .LT. 0.) GO TO 800 100 CONTINUE 800 FC=FXS(MA) FB=FXS(MB) DA=DIFF(MA) DB=DIFF(MB) NL=K C FIND MAXIMUM PERMISSIBLE VALUE OF FXS (=FB) USING FALSE POSITION C RULE AND STARTING WITH FXS(MA) AND FXS(MB) FROM THE DO 100 LOOP. C SEARCH STOPS BEFORE TEN VALUES HAVE BEEN TRIED IF THE ABSOLUTE C VALUE OF THE DIFFERENCE BETWEEN COMPUTED AND ACTUAL FH(2)<.00001. DO 900 M=1,10 FXS(M)=FC-DA*(FB-FC)/(DB-DA) IF(FXS(M).GT.F1MAX) FXS(M)=F1MAX FH(2)=FA-1.001*FXS(M)*(FXS(M)-FHS)/FA DEN(1)=FXS(M)*(FXS(M)-FHS)*12400.0 DEN(2)=FA*(FA-FH(2))*12400.0 CALL XLAM(L,I,DLAT,DLONG,DATE,LQ,2,&2000) HGT=TH(1)-(HP(2)*DLOG(X2/X1))/ASUM CALL FIELDG(DLAT,DLONG,HGT,DATE,10,LQ,BN,BE,BV,BT) DIFF(M)=(BT*2.8000E-05) - FH(2) TH(2)=HGT FB=FXS(M) DB=DIFF(M) K=M IF (DABS(DB).LT. 0.00001) GO TO 1000 900 CONTINUE C COMPUTE FXS=F(1) WHICH GIVES THE SAME SLOPES FOR THE FIRST TWO C LAMINATIONS (AA(2)=AA(3)). INITIAL VALUE OF F(1) IS THE AVERAGE OF C MAXIMUM (FB) AND MINIMUM (FH(1)) PERMISSIBLE VALUES. FB IS THE C REFERENCE POINT FOR FALSE POSITION RULE USED TO REFINE SUCCESSIVE C VALUES OF FXS. PROCEDURE IS APPLIED TO DEL2=(1-AA(3)/AA(2)) FOR C LOW DENSITY IONOGRAMS (DENSITY CORRESPONDING TO FB < 500 EL/CC); C OTHERWISE, DEL2=(1-AA(3)/AA(2))**2. CONVERGENCE IS ACHIEVED WHEN C (1.0-AA(3)/AA(2))**2 < .0001. A MAXIMUM OF 15 VALUES OF FXS TRIED. 1000 F(1)=FB DELF=FB-FHS LOWN=.FALSE. IF((1.24D4*FB*(FB-FH(1))) .LE. 500.) LOWN=.TRUE. K=K+NL MA=1 DELA = 1.0 NL=N DO 1100 M=1,15 FX1=F(1) DEN(1)=12400.*FX1*(FX1-FH(1)) DO 1300 I=2,3 IBAD=0 DP=0.0 FA=F(I) FH(I)=FA-DEN(I-1)*8.07258870D-5/FA DEN(I)=12400.*FA*(FA-FH(I)) CALL XLAM(L,I,DLAT,DLONG,DATE,LQ,3,&1700) GO TO 1300 1700 IBAD=1 1300 CONTINUE IF (IBAD .EQ. 0) GO TO 1800 IF(M.EQ.15) GO TO 200 IF (M .GT. MA) GO TO 1400 C RESET REFERENCE. MA=MA+1 F(1)=F(1)-0.001 GO TO 1100 C RESTART ON FALSE POSITION METHOD. 1400 F(1)=FHS+DFLOAT(M)*DELF/16. GO TO 1100 1800 RATIO=AA(3)/AA(2) DELA=1.-RATIO R=RATIO IF(.NOT.LOWN) R=R**2 IF (M .GT. MA) GO TO 1900 DEL1=1.-R F1=F(1) F(1)=0.5*(FH(1)+FB) GO TO 1100 1900 DEL2=1.-R F2=F(1) IF(DELA**2.LT.0.0001) GO TO 1200 F(1)=F1-(DEL1*(F2-F1))/(DEL2-DEL1) 1100 CONTINUE 1200 WRITE(16,1)FB,K,FX1,RATIO,M METHOD=METH N=NL RETURN 200 WRITE(16,5) GO TO 2100 1500 WRITE(16,4) GO TO 2100 2000 WRITE(16,3) 2100 METHOD=METH WRITE(16,2) RETURN 1 1 FORMAT(' FXSMAX=',F6.3,' MHZ.',I4,' ITERATIONS (MAX=20).'/ F ' COMPUTED FXS=',F6.3,' MHZ SLOPE2/SLOPE1=',F7.4,I4, F ' ITERATIONS (MAX=15).') 2 FORMAT(' PROBLEM OCCURRED DURING CALL TO FXSAT SUBROUTINE. PROCES FSING OF THIS IONOGRAM TERMINATED.') 3 FORMAT(' DATA FOR F=',F6.3,' NOT MONOTONICALLY INCREASING.') 4 FORMAT(' FH(2) TOO CLOSE TO FHS.') 5 FORMAT(' CALCULATION OF FXS DID NOT CONVERGE WITH 15 ITERATIONS.') END C THIS ROUTINE IS A FAST DOUBLE-PRECISION POINTER BUBBLE SORT. C IT DOES NOT DO THE ACTUAL EXCHANGING OF ELEMENTS. TO DO SO, C THE MAIN ROUTINE SHOULD HAVE A(I)=AR(IPTR(I)). SUBROUTINE PTRSRT(AR,N,IPTR) IMPLICIT REAL*8(A) DIMENSION AR(N),IPTR(N) DO 10 K=1,N 10 IPTR(K)=K N1=N-1 DO 20 K=1,N1 K1=N-K IFLG=0 DO 30 L=1,K1 IF(AR(IPTR(L)).LE.AR(IPTR(L+1))) GO TO 30 IFLG=1 I=IPTR(L) IPTR(L)=IPTR(L+1) IPTR(L+1)=I 30 CONTINUE C NOTE THAT THE ROUTINE STOPS IF NO EXCHANGES ARE MADE ON C ANY GIVEN PASS. MAXIMUM SORT STEPS IS N(N-1)/2. THE C POINTERS WIL SORT THE ARRAY TO BE INCREASING. IF(IFLG.EQ.0) RETURN 20 CONTINUE RETURN END SUBROUTINE ZLIMIT(Z,ISIZE1,JSIZE,ZMIN,ZMAX) C C NAME: SUBROUTINE ZLIMIT C C DESCRIPTION: ZLIMIT FINDS THE MAXIMUM (ZMAX) AND MINIMUM (ZMIN) OF C THE Z ARRAY. C DIMENSION Z(ISIZE1,JSIZE) COMMON/OPTION/IMIN,IMAX,JMIN,JMAX COMMON/CON$/DUM$(33),UNDEF C INCLUDE 'OVRLAY.FTN' DATA XMIN,XMAX/2*0.0/ FIRST=1.0 SECOND=0.0 DO 1 I=IMIN,IMAX DO 1 J=JMIN,JMAX ZZ=Z(I,J) IF (ZZ.EQ.UNDEF) GO TO 1 XMIN=AMIN1(XMIN,ZZ)*SECOND+ZZ*FIRST XMAX=AMAX1(XMAX,ZZ)*SECOND+ZZ*FIRST FIRST=0.0 SECOND=1.0 1 CONTINUE ZMIN=XMIN ZMAX=XMAX RETURN END C THIS SUBROUTINE PLOTS ELECTRON DENSITY VS ALTITUDE FOR C EACH INPUT IONOGRAM. IT ALSO PLOTS THE XTRACE AND THE C CALCULATED O-TRACE (RANGE VS FREQUENCY). FINALLY, IT C COLLECTS THE DATA FOR USE IN THE CONTOURING ROUTINE, C CONT. SUBROUTINE PLOT1(IHDR,TITL,ISP,ZHMAX,ISVE,ITYP,IX,IY,ZDLAT,ZDLONG) IMPLICIT REAL*8 (A-H,O-Y) dimension plas_freq(70) LOGICAL*1 TIT2(80),TITL(80) LOGICAL*4 TIT4(20) COMPLEX*16 TIT1(5) real*4 xmid,ymid EQUIVALENCE (TIT1(1),TIT2(1)),(TIT4(1),TIT2(1)) COMMON/MID/XMID,YMID C TIT1 IS USED TO BUILD THE SKELETON TITLE FOR THE ELECTRON DENSITY C PLOTS. WE THEN USE TIT2 (WHICH IS EQUIVALENCED WITH IT) TO FILL C IN THE HOLES. DATA TIT1/'DAY ( / ) ','OF 19 AT ',' GMT LMT' X,' ',' '/ C ZTH4,ZED4,ZF,ZFO,ZHP,ZHPO ARE SINGLE PRECISION ARRAYS DESIGNED C TO HOLD THE DATA FOR THE TWO PLOTS MADE FROM EACH IONOGRAM. DIMENSION IPTR(210),TH3(210),ED3(210),ZTH4(210),ZED4(210),ZF(50), X ZFO(50),FH1(70),FN(70),ZHPO(50),ZHP(50) COMMON /CALC/ AA(50),BB(50),F(50),HP(50),FH(50),TH(70),DEN(70) C COMMON /CNT/ HOLDS THE DATA TO BE CONTOURED, AND PASSES IT C BACK TO MAIN. COMMON /CNT/ZA(2500,3),ZXCO(50),LA,LXCO,ZA2(80,3),LA2 COMMON /IO/IYR,IDY,ITI,LMT,TITLE(40),DIP,FHS,MON,IDAY COMMON /TIME/IH,IM,IS COMMON /OCALC/ FO(50),HPO(50) COMMON /SWJTCH/ ISAT,LOCSAT,IFXS,NCASES,ICODE(4),N,METHOD,IORD,IP COMMON /PL/ TH1(200),ED1(200),TH2(200),ED2(200),M,M1 C IF THIS IS THE FIRST CALL TO THIS ROUTINE,WE MUST INITIALIZE. DATA INIT /0/ c open(unit=10,file='plas_freq.dat',status='new',form='formatted') IF(INIT.NE.0) GO TO 10 c CALL PLOTST(00001,1) c XL=500. c YB=375. c XR=3942. c YT=2850. c CALL SETGRD(500.,375.,3942.,2850.,1) c XMID=((XR-XL)/2.)+XL c YMID=((YT-YB)/2.)+YB 10 CONTINUE C WE NOW FIT DAY, MONTH, TIME, ETC INTO THE SKELETON TITLE. c CALL EDIT(REAL(IDY),'I3)',TIT2(5),NN,IBL,1,IERR) c CALL EDIT(REAL(MON),'I2)',TIT2(10),NN,IBL,1,IERR) c CALL EDIT(REAL(IDAY),'I2)',TIT2(13),NN,IBL,1,IERR) c CALL EDIT(REAL(IYR),'I2)',TIT2(22),NN,IBL,1,IERR) c CALL EDIT(REAL(ITI),'I6)',TIT2(28),NN,IBL,1,IERR) c CALL EDIT(REAL(LMT),'I6)',TIT2(39),NN,IBL,1,IERR) C NOW SIGNAL THAT THIS ROUTINE HAS BEEN CALLED. INIT=1 DO 15 I=1,N !20 ZF(I)=F(I) ZHP(I)=HP(I) ZFO(I)=FO(I) 15 ZHPO(I)=HPO(I) WRITE(7,*) ' X TRACE ' WRITE(7,*) ' FREQ X_APPRGE(HP) ' WRITE(8,*) ' O TRACE ' WRITE(8,*) ' FREQ O_APPRGE(HP) ' do ka=1,N ! 20 if(zf(ka).ne.0.OR.zhp(ka).ne.0.) write(7,*) zf(ka),zhp(ka) if(zfo(ka).ne.0.OR.zhpo(ka).ne.0.) write(8,*) zfo(ka),zhpo(ka) if(zfo(ka).ne.0.OR.zhpo(ka).ne.0.) write(60,*) zfo(ka),zhpo(ka) enddo C WE NOW MUST PLOT THE X AND O TRACES. c CALL SETSIZ(1.5) c CALL GRID(0.,10.,40,'F3.0)',4,4400.,0.,22,'F5.0)',1,0) C IF HEADERS HAVE BEEN REQUESTED (SEE READ4), PLOT THEM. c IF(IHDR.NE.0) GO TO 17 c CALL VERLIN('VIRTUAL HEIGHT(KM)',-18,300.,YMID,-8,0) c CALL HORLIN('FREQUENCY(MHZ)',14,XMID,190.5,0,-2) c 17 CALL PLOT(ZF,ZHP,N,'X') C IF THE O-TRACE HAS BEEN CALCULATED, PLOT IT c IF(IORD.EQ.1) CALL PLOT(ZFO,ZHPO,N,'O') c CALL FRMADV C NEXT THE ELECTRON DENSITY PLOT. C CALCULATE THE NUMBER OF INTERVALS ON THE Y-AXIS FROM THE MAX ALT. L=ZHMAX/500+.5 c CALL OGRID(1.,7.,6,'F2.0)',1,0.,ZHMAX,L,'F5.0)',1,0) c CALL SCALE(10.,1.00E+07,0.,ZHMAX,1) C IF HEADERS HAVE BEEN REQUESTED, PLOT THEM. c IF(IHDR.NE.0) GO TO 18 c CALL HORLIN('LOG ELECTRON DENSITY(CM-3)',26,XMID,190.5,0,-2) c CALL VERLIN('ALTITUDE(KM)',-12,300.,YMID,-8,0) c CALL HORLIN(TITL,70,XMID,2850.,0,4) c CALL HORLIN(TIT2,48,XMID,2850.,0,2) C NOW WE GATHER UP ALL THE ELECTRON DENSITY AND TRUE HEIGHT ARRAYS. 18 DO 20 I=1,N TH3(I)=TH(I) 20 ED3(I)=DEN(I) DO 30 I=1,M TH3(N+I)=TH1(I) 30 ED3(N+I)=ED1(I) DO 40 I=1,M1 TH3(M+N+I)=TH2(I) 40 ED3(M+N+I)=ED2(I) M2=M+N+M1 C NOW THAT THE THREE ELECTRON DENSITY ARRAYS HAVE BEEN COALESCED, C WE MUST SORT THEM. CALL PTRSRT(ED3,M2,IPTR) DO 50 I=1,M2 ZED4(I)=ED3(IPTR(I)) 50 ZTH4(I)=TH3(IPTR(I)) write(9,*) ' ELECTRON DENSITY ' write(9,*) ' ELEC_DENS NE_HGT ' write(10,*) ' PLASMA FREQUENCY ' write(10,*) ' PLAS_FREQ PLAS_FREQ_HGT ' do kb=1,m2 plas_freq(kb)=(.00898*sqrt(zed4(kb))) write(9,*) zed4(kb),zth4(kb) write(63,*) zed4(kb),zth4(kb) write(10,932) plas_freq(kb),zth4(kb) write(62,932) plas_freq(kb),zth4(kb) enddo 932 format(2f12.6) 1010 FORMAT(5I12,F16.7) C NOW FOR THE ACTUAL PLOTTING. c CALL PLOT(ZED4,ZTH4,M2,' ') c CALL FRMADV C NOW WE INTERPOLATE FH VALUES AT EACH 50 KM MARKER BETWEEN THE F-PEAK C AND THE SATELLITE. J=TH(1)/50 K=TH(N)/50+1 L=N DO 110 I=K,J 120 IF(I.LE.(TH(L-1)/50)) GO TO 110 L=L-1 GO TO 120 110 FH1(J-I+2)=FH(L)*(6371.2+TH(L))**3/(6371.2+50*I)**3 FH1(1)=FH(1) FH1(M)=FH(N) C NOW WE PLACE THE ELECTRON DENSITIES AT THE SAME LOCATIONS IN FN. DO 135 I=1,M 135 FN(I)=ED1(I) C IF WE ARE GOING TO CONTOUR ELECTRON DENSITY, SKIP CALCULATIONS. IF(ITYP.EQ.0) GO TO 130 C NOW CALCULATE THE PLASMA FREQUENCY AND PUT IT IN FN. c DO 140 I=1,M FN(I)=DSQRT(ED1(I)*80.6)/1000. c type*,' plasma freq from newplot1 ',i,fn(i) 140 continue c 140 write(10,1981) fn(i),th1(i) 1981 format(2f12.6) GO TO (130,150,160,170,180,190),ITYP C IF FN/FH HAS BEEN REQUESTED, CALCULATE IT. 150 DO 155 I=1,M 155 FN(I)=FN(I)/FH1(I) GO TO 130 C IF FX HAS BEEN REQUESTED, CALCULATE IT. 160 DO 165 I=1,M 165 FN(I)=.5*(FH1(I)+DSQRT(FH1(I)**2+4*FN(I)**2)) GO TO 130 C IF FZ HAS BEEN REQUESTED, CALCULATE IT. 170 DO 175 I=1,M 175 FN(I)=.5*(-FH1(I)+DSQRT(FH1(I)**2+4*FN(I)**2)) GO TO 130 C IF FT HAS BEEN REQUESTED, CALCULATE IT. 180 DO 185 I=1,M 185 FN(I)=DSQRT(FH1(I)**2+FN(I)**2) GO TO 130 C IF FH HAS BEEN REQUESTED PUT IT IN FN ARRAY. 190 DO 195 I=1,M 195 FN(I)=FH1(I) 130 CONTINUE C NOW TO GATHER THEM WITH THE VALUES FOR OTHER IONOGRAMS TO C BE CONTOURED. LXCO=LXCO+1 IX1=IX+1 C WE DECIDE WHAT TO PUT IN ZXCO BASED ON WHAT X-COORDINATE HAS BEEN C SPECIFIED. GO TO(200,210,220,230),IX1 200 ZXCO(LXCO)=ZDLAT GO TO 240 210 ZXCO(LXCO)=IS+IM*60 GO TO 240 220 ZXCO(LXCO)=IM+IS/60. GO TO 240 230 ZXCO(LXCO)=ZDLONG 240 CONTINUE C C THIS CODE ADDED IN ATTEMPT TO IMPLEMENT THE 'SAVE DATA' OPTION C UNIT 9 IS LIB.DATA(SAVEIT) C SAVING FN,FH,ZDLAT,ZDLONG,M,TH1,IYR,IDY,IH,IM,IS C P.K. PARKER 9/11/81 C IF(ISVE .EQ. 0) GOTO 250 WRITE(9,300) M WRITE(9,310) ZDLAT,ZDLONG 300 FORMAT(I6) 310 FORMAT(F14.8,4X,F14.8) DO 320 I=1,M 320 WRITE(9,330) TH1(I),ED1(I),FH1(I) 330 FORMAT(3(F20.8,4X)) WRITE(9,340) IYR,IDY WRITE(9,350) IH,IM,IS 340 FORMAT(2I10) 350 FORMAT(3I10) C C END SAVE OPTION INSERT C C NOW WE STORE THE INFO IN THE ARRAY TO BE CONTOURED. 250 DO 60 I=1,M ZA(LA+I,1)=ZXCO(LXCO) ZA(LA+I,2)=TH1(I) IF(IY.EQ.1) ZA(LA+I,2)=TH1(1)-ZA(LA+I,2) ZA(LA+I,3)=FN(I) c print*,' za ',za(la+i,1),za(la+i,2),za(la+i,3) 60 CONTINUE LA=LA+M C NOW WE STORE SATELLITE OR F-PEAK POINTS, AS REQUESTED. IF(ISP.EQ.0) RETURN IF(ISP.EQ.2) GO TO 70 LA2=LA2+1 ZA2(LA2,1)=ZXCO(LXCO) ZA2(LA2,2)=TH1(1) IF(IY.EQ.1) ZA2(LA2,2)=TH1(1)-ZA2(LA2,2) ZA2(LA2,3)=FN(1) 70 IF(ISP.EQ.1) RETURN LA2=LA2+1 ZA2(LA2,1)=ZXCO(LXCO) ZA2(LA2,2)=TH1(M) IF(IY.EQ.1) ZA2(LA2,2)=TH1(1)-ZA2(LA2,2) ZA2(LA2,3)=FN(N) RETURN END SUBROUTINE CONVRT(DATA,NDIM,N,DXMIN,DXMAX,DYMIN,DYMAX, * XMIN,XMAX,YMIN,YMAX) C THIS SUBROUTINE CONVERTS DATA FROM INPUT UNITS TO UNITS SCALED TO C FIT WITHIN THE INDICES OF THE Z-ARRAY C C DIMENSION DATA(N,3) ! WAS DATA(NDIM,3) C do i=1,N C WRITE(13,*) data(i,1),DATA(I,2),DATA(I,3) C enddo C print*,' ndim, n ',ndim,n C print*,' dxmin,dxmax,dymin,dymax',dxmin,dxmax,dymin,dymax C print*,' xmin,xmax,ymin,ymax',xmin,xmax,ymin,ymax CEI=(XMAX-XMIN)/(DXMAX-DXMIN) C print*,' cei ',cei C cej1=ymax-ymin C print*,' cej1 ',cej1 C cej2=dymax-dymin C print*,' cej2 ',cej2 C cej=(ymax-ymin)/(dymax-dymin) CEJ=(YMAX-YMIN)/(DYMAX-DYMIN) C print*,' cej ',cej DO 1 I=1,ndim !N C type*,' data(i,1) before ',data(i,1) DATA(I,1)=CEI*(DATA(I,1)-DXMIN)+XMIN C type*,' data(i,1) after ',data(i,1) C type*,' data(i,2) before ',data(i,2) DATA(I,2)=CEJ*(DATA(I,2)-DYMIN)+YMIN C type*,' data(i,2) after ',data(i,2) 1 CONTINUE RETURN END SUBROUTINE JWORLDMAP(JSAT,lyear,lday,lhour,min,lsec,DLAT,DLONG,HGT, 1 jrou) C PROGRAM EFCALC13 C C SIMPLE EPHEMERIS CALCULATION C USES BROWER ELEMENTS UP TO DAY 283, 1979 AND C NORAD ELEMENTS FROM DAY 284 ON. C PRODUCES LAT,LONG,HGT,GEOMAGNETIC TIME (TO HALF HOUR) C GROMAGNETIC LAT,GEOMAGNETIC LONG,FH,DIP,CHI, C SL(IN SUNLIGHT), NSL(NOT IN SUNLIGHT) C FORMAT($ EFCALC UPDATED OCTOBER 26/72$) C PRINTS OUT AT INTERVALS THROUGHOUT TIME RANGE C FOR ISIS 1, ISIS2, ALOUETTE 1 OR ALOUETTE 2 C UPDATE: CONVERTED from CP6 Fortran to Sun Fortran (1992) on June 10 C 1993 by Michael G. Fischer (UCalgary) C ALL subroutines and functions are now part of this source code, C except for the standard intrinsic Fortran functions, no other sources C are required. C Four data files are required: ISIS1,ISIS1N & ISIS2,ISIS2N C TO COMPILE: ALL the required source code is listed within this file C (as of June 14, 1993) to compile and link the program (with Sun FortranC 1992) just go: C THERE are two versions of this source code: one requires compilation C with the automatic double precision option, the other doesn't, both C version produce identical results (which are identical to the CP6 C version) C efcalc13.f <- requires automatic double precision option C efcalc13d.f <- doesn't (intended for machines without ADP option C f77 -r8 efcalc13.f <- the -r8 engages automatic double precision C f77 efcalc13d.f <- non ADP version C All the files (as of June 14, 1993) can be found in efcalc13.tar.Z C IMPLICIT REAL*8 (A-H,O-Z) C OPEN(UNIT=77,NAME='TEMPOUT.DAT',TYPE='OLD', C 1 FORM='FORMATTED') C OPEN(UNIT=75,NAME='ISIS_OUT.DAT',TYPE='NEW',FORM='FORMATTED') ISAT=JSAT ERR=0.01 L=0 RE=6371.2 PI=3.14159 DR=PI/180. c 665 WRITE (*,666) ' 1 FOR ALOUETTE 1' c WRITE(*,666) ' 2 FOR ALOUETTE 2' c WRITE(*,666) ' 3 FOR ISIS 1' c WRITE(*,666) ' 4 FOR ISIS 2' 666 FORMAT (1X,'Please enter the ',A) c READ (*,10) ISAT C IF(ISAT.EQ.1) WRITE(75,*) ' ISIS 1' C IF(ISAT.EQ.2) WRITE(75,*) ' ISIS 2' c IF (ISAT.GE.5 .OR. ISAT.EQ.0) GOTO 665 IF (ISAT.EQ.1) THEN ! ALOUETTE 1 OPEN(UNIT=32,FILE='AL1.DAT',STATUS='OLD') C OPEN(UNIT=32,FILE='AL1B.DAT',STATUS='OLD') ELSEIF (ISAT.EQ.2) THEN ! ALOUETTE 2 OPEN(UNIT=32,FILE='AL2.DAT',STATUS='OLD') ELSEIF (ISAT.EQ.3) THEN ! ISIS 1 ISAT = 5 OPEN (UNIT=32,FILE='ISIS1.DAT',STATUS='OLD') OPEN (UNIT=33,FILE='ISIS1N.DAT',STATUS='OLD') ELSEIF (ISAT.EQ.4) THEN ! ISIS 2 ISAT = 6 OPEN (UNIT=32,FILE='ISIS2.DAT',STATUS='OLD') OPEN (UNIT=33,FILE='ISIS2N.DAT',STATUS='OLD') ENDIF if(jrou.eq.1) then write(20,100) write(21,100) write(22,100) write(64,100) endif C WRITE(75,100) 100 FORMAT(' YR DAY TIME LMT LAT LONG HGT GMLTM GMLAT GML &ONG INVLAT FH DIP CHI SUN L') 1 CONTINUE C READ(77,1010,END=999) LYEAR,LDAY,LHOUR,MIN,LSEC 1010 FORMAT(13X,I2,I3,1X,3I2) c WRITE (*,666) 'year LAST 2 DIGITS (ZERO TO END)' c READ(*,10,END=200) LYEAR c IF(LYEAR.EQ.00) GOTO 999 c WRITE (*,666) 'day' c READ(*,10,END=200) LDAY c IF(LDAY.GT.366) THEN c TYPE*,' DAY NUMBER BEYOND 366 ' c STOP c ENDIF 10 FORMAT(I3) c WRITE (*,666) 'hour' c READ(*,10) LHOUR c IF(LHOUR.GT.23) THEN c TYPE*,' HOUR BEYOND 23 ' c STOP c ENDIF c WRITE (*,666) 'min' c READ(*,10) MIN c IF(MIN.GT.59) THEN c TYPE*,' MINUTE BEYOND 59 ' c STOP c ENDIF c WRITE (*,666) 'sec' c READ(*,10) LSEC c IF(LSEC.GT.59) THEN c TYPE*,' SECOND BEYOND 59 ' c STOP c ENDIF c WRITE(*,100) IF(LYEAR.LT.79)GO TO 50 IF(LYEAR.EQ.79.AND.LDAY.LT.284)GO TO 50 CALL ORBITB(LYEAR,LDAY,LHOUR,MIN,LSEC,DLAT,DLONG,HGT,FH, CDIP,LDIP,CHI) GO TO 51 50 IF(L.NE.0)GO TO 58 CALL ORBITX(79,284,0,0,0,DLAT,DLONG,HGT,FH,DIP,LDIP,CHI,L,ISAT) 58 CALL ORBITX(LYEAR,LDAY,LHOUR,MIN,LSEC,DLAT,DLONG, CHGT,FH,DIP,LDIP,CHI,L,ISAT) 51 IF(CHI.LT.90.) THEN ILL=3H SL GO TO 19 END IF DA=(HGT+RE)*SIN(PI-CHI*DR) IF(DA.GT.(RE+30.)) THEN ILL=3H SL GO TO 19 END IF ILL=3HNSL 19 LHI=CHI+0.5 CALL LMT(LHOUR,MIN,LSEC,DLONG,JTIME) JHOUR=LHOUR JMIN=MIN CALL INVAR(DLAT,DLONG,HGT,ERR,BB,FL) IF(FL.LT.1.0)FL=1.0 XINLAT=DATAN(DSQRT(DBLE(FL)-1.0))*57.2957795 XINVLAT=SIGN(XINLAT,DIP) CALL GMLTIME(JHOUR,JMIN,DLAT,DLONG,GLAT,GLONG,ITIME) C WRITE(*,20)LYEAR,LDAY,LHOUR,MIN,LSEC,JTIME,DLAT,DLONG,HGT, C CITIME,GLAT,GLONG,XINVLAT,FH,LDIP,LHI,ILL,FL C WRITE(75,20)LYEAR,LDAY,LHOUR,MIN,LSEC,JTIME,DLAT,DLONG,HGT, C CITIME,GLAT,GLONG,XINVLAT,FH,LDIP,LHI,ILL,FL if(jrou.eq.1) then WRITE(20,20)LYEAR,LDAY,LHOUR,MIN,LSEC,JTIME,DLAT,DLONG,HGT, CITIME,GLAT,GLONG,XINVLAT,FH,LDIP,LHI,ILL,FL WRITE(21,20)LYEAR,LDAY,LHOUR,MIN,LSEC,JTIME,DLAT,DLONG,HGT, CITIME,GLAT,GLONG,XINVLAT,FH,LDIP,LHI,ILL,FL WRITE(22,20)LYEAR,LDAY,LHOUR,MIN,LSEC,JTIME,DLAT,DLONG,HGT, CITIME,GLAT,GLONG,XINVLAT,FH,LDIP,LHI,ILL,FL WRITE(64,20)LYEAR,LDAY,LHOUR,MIN,LSEC,JTIME,DLAT,DLONG,HGT, CITIME,GLAT,GLONG,XINVLAT,FH,LDIP,LHI,ILL,FL endif 20 FORMAT(I3,I4,I3,I2,1H:,I2,I5,F7.2,F8.2,F8.2,I5,F7.2,F8.2,F7.2,F6.3 C,2I4,1X,A3,F8.2) c f5.2 c GO TO 1 c999 CONTINUE CLOSE(UNIT=32) CLOSE(UNIT=33) RETURN END C CHECK STATEMENT 367; ISAT=? C UPDATE: Converted from CP6 Fortran to Sun Fortran (1992) on June 9, 1993 C by Michael G. Fischer (UCalgary) C All Hollerith constants were changed to their ASCII equivalents SUBROUTINE ORBITX(LYEAR,LDAY,LHOUR,MIN,LSEC,DLAT,DLONG,HGT,FH, C DIP,LDIP,CHI,LSTART,ISAT) C C THIS IS THE VERSION OF ORBITX THAT USES THE SATT LIBRARY. C IT MUST BE COMPILED IN AUTOMATIC DOUBLE PRECISION. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION OSCELM(6), G(5,2) SAVE COMMON /ORBDATA/ ELMEAN(6,400), EPOCH(400), NSET COMMON /CCCBRW/ GP(5) C C COMMON /ORBDATA/ IS USED ONLY IN SUBROUTINE ORBITX . C IT IS USED TO PUT THE DATA IN A CORE AREA WHICH WILL NOT BE OVERLAID. C COMMON /CCCBWR/ CONTAINS THE EARTH CONSTANTS USED BY SUBROUTINE BRWR. C DATA G /1.146730234E+4, 5.41095E-4, 2.285E-6, 7.96125E-7, 2.32E-7, C 1.146781276E+4, 5.41240E-4, 2.560E-6, 6.90000E-7, 0.60E-7/ C C OSCELM: OSCULATING ORBITAL ELEMENTS C G: TWO DIFFERENT SETS OF EARTH CONSTANTS - OLD AND NEW C THESE CONSTANTS ARE EXPRESSED IN UNITS OF EARTH RADII AND DAYS. C IF(LSTART.NE.0) GO TO 13 C C READ (BROUWER MEAN) ORBITAL ELEMENTS C 1 LSTART=1 REWIND 32 READ(32,700) NSET IF(NSET.LE.400) GO TO 6 WRITE(*,4) NSET 4 FORMAT(1H1,I6,' SET OF ORBIT ELEMENTS IS GT 400, ORBITX ERROR') STOP 6 DO 10 I=1,NSET READ (32,701) IMO, IDA, IYR, IHR, IMIN, SEC READ (32,702) ELMEAN(1,I), ELMEAN(2,I), ELMEAN(3,I), C ELMEAN(6,I), ELMEAN(5,I), ELMEAN(4,I) 10 EPOCH(I) = JDS(IDA,IMO,IYR) + FDAY(IHR,IMIN,SEC) RETURN C 13 CONTINUE C C SEARCH TABLE FOR THE EPOCH NEAREST TO THE GIVEN TIME. C PERMISSIBLE TIMES RANGE FROM 7.5 DAYS BEFORE THE FIRST EPOCH C TO 7.5 DAYS AFTER THE LAST EPOCH. C SEC = LSEC TIME = JDS(LDAY,1,LYEAR) + FDAY(LHOUR,MIN,SEC) IF (TIME .LT. EPOCH(1)-7.5) THEN WRITE (*,28) LYEAR,LDAY STOP END IF 28 FORMAT (' ORBITX: 19',I2,' DAY',I4, C ' IS BEFORE LAUNCH OF SATELLITE') DO 32 I=2,NSET 32 IF (EPOCH(I) .GT. TIME) GO TO 35 IF (EPOCH(NSET)+7.5 .GT. TIME) THEN I = NSET GO TO 35 END IF WRITE (*,33) LYEAR,LDAY 33 FORMAT (' ORBITX: 19',I2,' DAY',I4,' : ORBITAL ELEMENTS ARE NOT CYET AVAILABLE FOR THIS DATE') STOP 35 IF (EPOCH(I)-TIME .GT. TIME-EPOCH(I-1)) I = I-1 DT = TIME - EPOCH(I) C C CHOOSE A SET OF EARTH CONSTANTS. C THE CONSTANTS WERE CHANGED JULY 1, 1970 (JULIAN SPACE DATE = 4669). C FOR ISIS1 (ISAT=5), THE NEW CONSTANTS WERE USED FROM THE BEGINNING. J = 1 IF (EPOCH(I) .GT. 4669.) J = 2 IF (ISAT .EQ. 5) J = 2 DO 15 II=1,5 15 GP(II) = G(II,J) C C IF THE ORBITAL ELEMENTS ARE THE SAME AS IN THE PREVIOUS CALL, C AND IF ORBITX HAS NOT MEANWHILE BEEN OVERLAID, C ENTRY BRWR2 IS USED TO SAVE TIME. C IF (I .EQ. IPREV .AND. ITEST .EQ. 86428642) GO TO 40 CALL BRWR(ELMEAN(1,I), DT, OSCELM) GO TO 42 40 CALL BRWR2(ELMEAN(1,I),DT,OSCELM) 42 IPREV = I ITEST = 86428642 C C THE SUBROUTINE BRWR HAS GENERATED OSCULATING ELEMENTS C FROM BROUWER MEAN ELEMENTS. C THE FOLLOWING FOUR SUBROUTINES FIND THE POSITION OF THE SATELLITE. RAGG = RAG(TIME) CALL ELMEFX(OSCELM, RAGG, X,Y,Z) CALL CTNPOLR(X,Y,Z,R,DLONG,CLAT) CALL GCTGD(R,CLAT,HGT,DLAT) C C SOLAR ZENITH ANGLE IS CALCULATED IN GEOCENTRIC COORDINATES. C MAXIMUM ERROR DUE TO THIS IS .19 DEGREES. CALL XYZSUN(TIME,RAGG,XSUN,YSUN,ZSUN) CHI = ACOS ((X*XSUN + Y*YSUN + Z*ZSUN) / R) C C CHANGE UNITS TO DEGREES AND KM C RE = 6378.166 RD = 180. / 3.141592653589793 DLAT = DLAT*RD DLONG = DLONG*RD IF(DLONG.GT.180.)DLONG=DLONG-360. CHI = CHI*RD HGT = HGT*RE C C MAGNETIC FIELD TM = 1900. + LYEAR + LDAY/365. L = 0 CALL FIELDGW(DLAT,DLONG,HGT,TM,10,L,X,Y,Z,B,FH,DIP,LDIP) RETURN 700 FORMAT(14I5) 701 FORMAT(5I2,F6.2) 702 FORMAT(6E13.5) END SUBROUTINE FIELDGW (DLAT,DLONG,ALT,TM,NMX,L,X,Y,Z,F,FH,DIP,LDIP) C THIS G ARRAY CONTAINS BOTH G AND H VALUES C THE COEFFICIENTS IN THE DATA STATEMENT ARE GSFC POGO(3/68) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION G(10,10),GT(10,10),SHMIT(10,10),AID(8) COMMON/FLD/TG(10,10) EQUIVALENCE (SHMIT,TG) SAVE DATA G/0.0000,-30469.6182,-1545.0059,1342.4801, 956.6111, & -241.3722, 47.4416, 76.3438, 7.5544, 7.3062, & 5800.1381, -2135.0002, 2987.7382, -1996.3616, 810.2744, & 363.2780, 57.0676, -53.3420, 0.8783, 6.3956, & -1966.5059, 167.1517, 1536.5008, 1302.8855, 507.5039, & 240.8420, 0.1275, 3.0896, -6.8812, 4.2724, & -459.1049, 225.7102, -163.7941, 911.0671, -350.0185, & -41.9546, -260.9984, 18.1722, -1.9995, -18.1036, & 130.9535, -271.0978, 0.5221, -158.7739, 271.1802, & -149.8464, -0.5149, -34.9527, -1.5733, 17.5593, & 24.0910, 121.2381, -103.7740, -101.5531, 102.1834, & -59.7087, -21.9609, -7.7097, 17.6242, 4.5229, & -5.9805, 108.8910, 65.4443, -54.4120, -1.4088, & -25.4429, -124.0177, 10.5723, -2.9160, 7.1386, & -60.3545, -30.8683, -10.7653, 1.3955, 17.8596, & -16.0346, -28.4959, -16.5459, 13.4169, 10.8215, & 7.9880, -12.2720, 10.2653, -4.7203, -1.4643, & 25.7864, -9.4226, -15.8364, 21.8852, 0.2932, & -16.4910, 21.7182, 4.8647, -0.7574, -0.4670, & 2.2628, 14.8247, -4.9400, -11.1155, 6.1647 / DATA GT/0.0000,26.2326,-23.2563,-8.1870,-0.2983, C 3.7970, -0.3454, -0.9405, 0.4405, 0.4571, C -6.0804, 4.6814, 2.2698, -9.2639, -1.4749, C -0.5420, 0.9442, 0.0383, 0.3801, 0.1981, C -8.9823, -8.8134, 11.7852, -0.7963, -5.9657, C 1.1791, 2.0947, 0.0781, 0.4881, -0.4975, C 10.7141, 2.6285, 0.7655, -10.9927, -6.3001, C 1.8080, 5.6534, -0.7792, -1.5293, 0.9545, C 4.2687, -0.8457, 2.8781, -14.2171, -3.3459, C -3.5372, 0.6064, 2.5089, -0.6244, -1.5044, C -0.9694, 1.0785, -4.4750, 0.7108, -5.2754, C 1.1564, 2.0906, 0.4394, -0.7129, -0.5128, C -0.9555, -0.5072, 0.5407, 1.9777, 0.8297, C 0.1223, -3.9239, 0.1895, 1.5310, -0.9824, C -0.5259, 0.8302, 0.4609, 1.1859, 1.1768, C -0.5557, 1.9179, 1.3661, 0.3665, -1.0291, C 0.3579, -0.2168, -0.5430, -1.5834, 0.0360, C 0.0874, -0.1225, 1.0628, -2.5528, 0.4840, C -0.8203, -0.9487, 0.2185, -0.2745, -0.7359, C 0.7973, -0.7113, 0.6129, 2.0906, -1.5406/ IF(A.NE.6378.165) GO TO 31 IF(L) 19,1,2 31 A=6378.165 44 FLAT=1.-1./298.3 45 A2=A**2 46 A4=A**4 47 B2=(A*FLAT)**2 48 A2B2=A2*(1.-FLAT**2) 49 A4B4=A4*(1.-FLAT**4) 50 MAXN=10 J=0 K=0 TZERO=1960.0 IF (L) 14,14,2 51 1 IF (TM-TLAST) 17,19,17 52 2 READ(105,3)J,K, TZERO,(AID(I),I=1,8) 3 FORMAT(1H ,2I1,1X,F6.1,8A8) L=0 55 WRITE(*,4)J,K,TZERO,(AID(I),I=1,8) 4 FORMAT(1H ,2I3,5X'EPOCH=',F7.1,5X8A8) MAXN=0 58 TEMP=0. 59 5 READ(105,6)N,M,GNM,HNM,GTNM,HTNM 6 FORMAT(1H ,2I3,6F11.4) IF (N.LE.0) GO TO 7 36 MAXN=(MAX0(N,MAXN)) 63 G(N,M)=GNM 64 GT(N,M)=GTNM 65 TEMP=DMAX1(TEMP,DABS(GTNM)) 66 IF (M.EQ.1) GO TO 5 G(M-1,N)=HNM GT(M-1,N)=HTNM GO TO 5 70 7 WRITE(*,8) 8 FORMAT ('0 N M',6X'G',10X'H',9X'GT',9X'HT'//) DO 12 N=2,MAXN 73 DO 12 M=1,N 74 MI=M-1 75 IF (M.EQ.1) GO TO 10 WRITE(*,9)N,M,G(N,M),G(MI,N),GT(N,M),GT(MI,N) 9 FORMAT(1H ,2I3,4F11.4) GO TO 12 79 10 WRITE(*,11)N,M,G(N,M),GT(N,M) 11 FORMAT(1H ,2I3,F11.4,11X,F11.4) 12 CONTINUE 82 WRITE(*,13) 13 FORMAT (1H1) 84 IF(TEMP.EQ.0.) L=-1 14 IF(K.NE.0) GO TO 17 SHMIT(1,1)=-1. DO 15 N=2,MAXN 88 SHMIT(N,1)=SHMIT(N-1,1)*DBLE(2*N-3)/DBLE(N-1) 89 SHMIT(1,N)=0. 90 JJ=2 91 DO 15 M=2,N 92 SHMIT(N,M)=SHMIT(N,M-1)*DSQRT(DBLE((N-M+1)*JJ)/DBLE(N+M-2)) 93 SHMIT(M-1,N)=SHMIT(N,M) 94 15 JJ=1 95 DO 16 N=2,MAXN 96 DO 16 M=1,N 97 G(N,M)=G(N,M)*SHMIT(N,M) 98 GT(N,M)=GT(N,M)*SHMIT(N,M) 99 IF (M.EQ.1) GO TO 16 G(M-1,N)=G(M-1,N)*SHMIT(M-1,N) GT(M-1,N)=GT(M-1,N)*SHMIT(M-1,N) 16 CONTINUE 17 T=TM-TZERO 103 DO 18 N=1,MAXN 104 DO 18 M=1,N 105 TG(N,M)=G(N,M)+T*GT(N,M) 106 IF (M.EQ.1) GO TO 18 TG(M-1,N)=G(M-1,N)+T*GT(M-1,N) 108 18 CONTINUE 109 TLAST=TM 110 19 SINLA=DSIN(DLAT/57.2957795) 111 RLONG=DLONG/57.2957795 112 CPH=DCOS(RLONG) 113 SPH=DSIN(RLONG) 114 IF (J.EQ.0) GO TO 20 R=ALT+6371.2 116 CT=SINLA 117 GO TO 21 118 20 SINLA2=SINLA**2 119 COSLA2=1.-SINLA2 120 DEN2=A2-A2B2*SINLA2 121 DEN=DSQRT(DEN2) 122 FAC=(((ALT*DEN)+A2)/((ALT*DEN)+B2))**2 123 CT=SINLA/SQRT(FAC*COSLA2+SINLA2) 124 R=DSQRT(ALT*(ALT+2.*DEN)+(A4-A4B4*SINLA2)/DEN2) 125 21 ST=DSQRT(1.-CT**2) 126 NMAX=MIN0(NMX,MAXN) 127 CALL FIELDX(ST,CT,SPH,CPH,R,NMAX,BT,BP,BR,B) Y=BP F=B IF (J) 22,23,22 129 22 X=-BT 130 Z=-BR 131 RETURN 132 C TRANSFORMS FIELD TO GEODETIC DIRECTIONS 133 23 SIND=SINLA*ST-DSQRT(COSLA2)*CT 134 COSD=DSQRT(1.0-SIND**2) 135 X=-BT*COSD-BR*SIND 136 Z=BT*SIND-BR*COSD 137 FH=F*2.8E-5 HOR=DSQRT(X*X+Y*Y) DIP=DATAN(Z/HOR) satan1=0.5 LDIP=DIP*57.2957795 + DSIGN(satan1,DIP) RETURN 138 END 139 SUBROUTINE FIELDX(ST,CT,SPH,CPH,R,NMAX,BT,BP,BR,B) IMPLICIT REAL*8 (A-H,O-Z) COMMON/FLD/G(10,10) SAVE DIMENSION P(10,10),DP(10,10),CONST(10,10),SP(10),CP(10),FN(10) 1,FM(10) IF (P(1,1).EQ.1.0) GO TO 3 1 P(1,1)=1. 6 DP(1,1)=0. 7 SP(1)=0. 8 CP(1)=1. 9 DO 2 N=2,10 FN(N)=N 11 DO 2 M=1,N 12 FM(M)=M-1 13 2 CONST(N,M)=DBLE((N-2)**2-(M-1)**2)/DBLE((2*N-3)*(2*N-5)) 14 3 SP(2)=SPH 16 CP(2)=CPH 17 DO 4 M=3,NMAX 18 SP(M)=SP(2)*CP(M-1)+CP(2)*SP(M-1) 19 4 CP(M)=CP(2)*CP(M-1)-SP(2)*SP(M-1) 20 AOR=6371.2/R 21 AR=AOR**2 22 BT=0. 23 BP=0. 24 BR=0. 25 DO 8 N=2,NMAX 26 AR=AOR*AR 27 DO 8 M=1,N 28 IF (N-M) 6,5,6 29 5 P(N,N)=ST*P(N-1,N-1) 30 DP(N,N)=ST*DP(N-1,N-1)+CT*P(N-1,N-1) 31 GO TO 7 32 6 P(N,M)=CT*P(N-1,M)-CONST(N,M)*P(N-2,M) 33 DP(N,M)=CT*DP(N-1,M)-ST*P(N-1,M)-CONST(N,M)*DP(N-2,M) 34 7 PAR=P(N,M)*AR 35 IF(M.EQ.1) GO TO 9 TEMP=G(N,M)*CP(M)+G(M-1,N)*SP(M) 361 BP=BP-(G(N,M)*SP(M)-G(M-1,N)*CP(M))*FM(M)*PAR 37 GO TO 10 371 9 TEMP=G(N,M)*CP(M) 38 BP=BP-(G(N,M)*SP(M))*FM(M)*PAR 381 10 BT=BT+TEMP*DP(N,M)*AR 39 8 BR=BR-TEMP*FN(N)*PAR 391 BP=BP/ST B=SQRT(BT*BT+BP*BP+BR*BR) 41 RETURN 42 END 43 SUBROUTINE ORBITB(LYEAR,LDAY,LHOUR,MIN,LSEC,DLAT,DLONG,HGT,FH, C DIP,LDIP,CHI) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IE(54), IDENT(4),E(27),A(6) SAVE EQUIVALENCE (E,IE) DATA IDENT /0, 6509801, 6900901, 7102401/ DATA PI,RE /3.141592653589793, 6378.166/ C RD = 180. / PI C SEC=LSEC TIME=JDS(LDAY,1,LYEAR)+FDAY(LHOUR,MIN,SEC) CALL SDC1(TIME,E,IFLAG) IF(IFLAG.NE.0)THEN WRITE(*,50)LYEAR,LDAY,LHOUR,MIN,LSEC STOP END IF 50 FORMAT(51H CHECK ISIS1N OR ISIS2N FOR ERROR. EPOCH CLOSEST TO 1 ,I2,I3,3I2) C CALL ELMENT (TIME-E(4),E(10),E(16),E(22),A) RAGG = RAG(TIME) CALL ELMEFX (A,RAGG,X,Y,Z) CALL CTNPOLR (X,Y,Z,R,DLONG,ALAT) CALL GCGD (R,ALAT,HGT,DLAT) HGT = HGT * RE DLAT = DLAT * RD DLONG = DLONG * RD IF (DLONG .GT. 180.) DLONG = DLONG - 360. C C SOLAR ZENITH ANGLE IS CALCULATED IN GEOCENTRIC COORDINATES. C MAXIMUM ERROR DUE TO THIS IS .19 DEGREES. CALL XYZSUN(TIME,RAGG,XSUN,YSUN,ZSUN) CHI = ACOS ((X*XSUN + Y*YSUN + Z*ZSUN) / R) CHI = CHI * RD C C MAGNETIC FIELD TM = 1900. + LYEAR + LDAY/365. L = 0 CALL FIELDGW(DLAT,DLONG,HGT,TM,10,L,X,Y,Z,B,FH,DIP,LDIP) RETURN END SUBROUTINE LMT(JHOUR,JMIN,JSEC,DLONG,JTIME) IMPLICIT REAL*8 (A-H,O-Z) SAVE satan1=1440.0 DTIME=DBLE(JHOUR)*60.0+DBLE(JMIN)+DBLE(JSEC)/60.0+ CDLONG*4.0 IF(DTIME) 3,4,4 3 DTIME=DTIME+satan1 4 DTIME=DMOD(DTIME,satan1) JTIME=(DTIME+0.5)/60.0 JJMIN=DTIME-60.0*DBLE(JTIME)+0.5 JTIME=100*JTIME+JJMIN RETURN END SUBROUTINE GMLTIME(JHOUR,JMIN,DLAT,DLONG,GLAT,GLONG,ITIME) C C CALCULATES GEOMAGNETIC LATITUDE AND LONGITUDE AND GEOMAG TIME C NEEDS U.T. , LAT AND LONGITUDE . C GEOMAGNETIC TIME IS PACKED HRMN IN ITIME . C IMPLICIT REAL*8 (A-H,O-Z) SAVE satan1 = 90.0-DLAT satan2 = DLONG-291.0 satan3 = 11.70 satan4 = 1440.0 CALL SPHERE (satan1,satan2,satan3,GLAT,GLONG) C CALL SPHERE(90.-DLAT,DLONG-291.,11.7,GLAT,GLONG) GLAT=90.-GLAT GLONG=GLONG-180. DTIME=DBLE(JHOUR)*60.+DBLE(JMIN)+GLONG*4.-276.0 IF(DTIME) 3,4,4 3 DTIME=DTIME+satan4 4 DTIME=DMOD(DTIME,satan4) ITIME=(DTIME+0.5)/60. JJMIN=DTIME-60.*DBLE(ITIME)+0.5 ITIME=100*ITIME+JJMIN RETURN END C SPHERE SUBROUTINE SPHERE(P1,A1,B,P2,A2) IMPLICIT REAL*8 (A-H,O-Z) PI=180. R = 57.29577950 SIN P1 = DSIN (P1/R) COS P1 = DCOS (P1/R) SIN A1 = DSIN (A1/R) COS A1 = DCOS (A1/R) SIN B = DSIN (B/R) COS B = DCOS (B/R) COS P2 = SIN P1 * COS A1 * SIN B + COS P1 * COS B SIN P2 = DSQRT (DABS (1. - COS P2 * COS P2)) IF(COS P2)2,1,2 1 P2 = 90. GO TO 3 2 RATIO=SIN P2/COS P2 P2=R*DATAN(RATIO) IF(SIN P2.GT.0..AND.COS P2.LT.0.)P2=PI+P2 IF(SIN P2.LT.0..AND.COS P2.LT.0.)P2=P2-PI IF (P2 .LT. 0.) P2 = P2 + 360. IF(P2)3,4,3 3 SIN A2 = -SIN P1 * SIN A1 / SIN P2 IF(COS B)6,5,6 4 A2 = 0. RETURN 5 COS A2 = COS P1 / SIN P2 GO TO 7 6 COS A2 = (SIN B * COS P2 - SIN P1 * COS A1)/(COS B * SIN P2) 7 RATIO=SIN A2/COS A2 A2=R*DATAN(RATIO) IF(SIN A2.GT.0..AND.COS A2.LT.0.)A2=PI+A2 IF(SIN A2.LT.0..AND.COS A2.LT.0.)A2=A2-PI IF (A2 .LT. 0.) A2 = A2 + 360. RETURN END SUBROUTINE NEWMAG(R,ST ,PHI,BR,BT,BP,B,THET) C NEWMAG VERSION OF DECEMBER 1965 C 1965.0 GSFC (9/65) C UPDATE: Converted from CP6 Fortran to Sun Fortran (1992) on June 9 1993 C by Michael G. Fischer (UCalgary) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION G(10,10),BM(10) SAVE DATA (G(N, 1),N=1,10)/ -0 C, 3.032193E 04, 2.522093E 03, -3.285459E 03, -4.170639E 03 C, 1.692819E 03, -6.684202E 02, -1.900312E 03, -2.405232E 02 C, C-9.358495E 02/ DATA (G(N, 2),N=1,10)/ -5.755070E 03 C, 2.131549E 03, -5.206994E 03, 6.237642E 03, -4.496227E 03 C, -3.650850E 03, -1.241578E 03, 2.029996E 03, -4.463745E 02 C, C-3.659410E 02/ DATA (G(N, 3),N=1,10)/ 3.495705E 03 C, -1.085898E 02, -1.369823E 03, -2.514676E 03, -1.943789E 03 C, -1.836598E 03, -1.313045E 02, -1.626874E 02, 4.917246E 02 C, C-8.068787E 02/ DATA (G(N, 4),N=1,10)/ 1.220352E 03 C, -4.753192E 02, 1.392784E 02, -6.836385E 02, 8.297622E 02 C, 1.568303E 02, 2.302497E 03, -1.540896E 02, 5.700617E 02 C, C1.292881E 03/ DATA (G(N, 5),N=1,10)/ -7.922399E 02 C, 1.080333E 03, -3.941087E 01, 2.055201E 02, -1.853181E 02 C, 3.569555E 02, -3.656370E 01, 3.012583E 02, 8.903696E 01 C, C-6.436587E 02/ DATA (G(N, 6),N=1,10)/ -2.424140E 02 C, -1.041800E 03, 5.898179E 02, 2.310479E 02, -5.887414E 01 C, 4.001436E 01, -1.209943E-02, 9.459898E 00, -1.050984E 02 C, C-3.745591E 02/ DATA (G(N, 7),N=1,10)/ 3.563806E 02 C, -1.545264E 03, -6.828717E 02, 1.681499E 02, 2.971388E 01 C, 6.276772E 00, 7.309118E 01, -3.402882E 01, 3.871370E 01 C, C-1.670375E 01/ DATA (G(N, 8),N=1,10)/ 1.915876E 03 C, 7.079673E 02, 1.857451E 02, -2.732077E 01, -1.705171E 02 C, 5.115862E 01, 1.302727E 01, -3.776955E 00, -2.940332E 01 C, C3.510623E-01/ DATA (G(N, 9),N=1,10)/ -4.633602E 02 C, 6.821298E 02, -2.394838E 02, 4.549622E 02, -3.794850E 01 C, -1.617146E 02, 6.268821E 00, 1.004341E 01, -4.002399E 00 C, C-4.152194E 00/ DATA (G(N,10),N=1,10)/ 2.803131E 03 C, -1.698787E 03, -4.244406E 02, 1.998351E 02, 6.192396E 01 C, -1.668931E 02, -9.080082E 01, -5.963821E-01, 1.524572E 00 C, C-9.238670E-01/ DATA J/0/NMAX/10/SUM/5.431251E04/ DATA (BM(I),I=1,10)/ 9.933492E 04 C, 9.933492E 04, 3.746322E 04, 2.457753E 04, 1.329481E 04 C, 6.468820E 03, 3.385349E 03, 1.616258E 03, 7.409154E 02 C, C3.641040E 02/ IF(SUM) 1000,1005,1000 1000 P22 = 0. C**** RELATIVE ERROR IN B BERR=0.0001 AR = 0. DO 1001 L=1,NMAX DO 1001 M=1,NMAX AR = AR + 1. 1001 P22 = P22 + AR*G(M,L) SUM=(SUM-P22)/SUM IF( ABS (SUM) - 1.E-4) 1004,1004,1002 C**** NOTE FOLLOWING PRINT AND STOP STATEMENTS 1002 PRINT 1003,SUM 1003 FORMAT(21H DATA WRONG IN NEWMAG E15.6) STOP 1004 SUM=0. 1005 P22=ST P21=SQRT (1.-P22*P22) AR=R IF(THET-1.570796327) 154,154,156 156 P21=-P21 154 IF(J) 152,152,151 151 SSQ=P22*P22 AR=AR+(14.288-SSQ*(21.3677+.108*SSQ))/6371.2 152 AR=1./AR C N= 2 159 DP22=P21 SP2=SIN (PHI) CP2=COS (PHI) DP21=-P22 AOR=AR*AR*AR C2=G(2,2)*CP2+G(1,2)*SP2 BR=-(AOR+AOR)*(G(2,1)*P21+C2*P22) BT=AOR*(G(2,1)*DP21+C2*DP22) BP=AOR*(G(1,2)*CP2-G(2,2)*SP2)*P22 IF(NMAX-3) 1,3,3 3 AOR=AOR*AR ERR=BERR*SQRT ((BP/P22)**2+BR**2+BT**2) IF(BM(3)*AOR-ERR) 1,1,103 103 SP3=(SP2+SP2)*CP2 CP3=(CP2+SP2)*(CP2-SP2) P31=P21*P21-0.333333333 P32=P21*P22 P33=P22*P22 DP31=-P32-P32 DP32=P21*P21-P33 DP33=-DP31 C2=G(3,2)*CP2+G(1,3)*SP2 C3=G(3,3)*CP3+G(2,3)*SP3 BR=BR-3.0*AOR*(G(3,1)*P31+C2*P32+C3*P33) BT=BT+AOR*(G(3,1)*DP31+C2*DP32+C3*DP33) BP=BP-AOR*((G(3,2)*SP2-G(1,3)*CP2)*P32+2.0*(G(3,3)*SP3-G(2,3)*CP3) 1*P33) C N= 4 IF(NMAX-4) 1,4,4 4 AOR=AOR*AR IF(BM(4)*AOR-ERR) 1,1,104 104 SP4=SP2*CP3+CP2*SP3 CP4=CP2*CP3-SP2*SP3 P41=P21*P31-0.26666666*P21 DP41=P21*DP31+DP21*P31-0.26666666*DP21 P42=P21*P32-0.20000000*P22 DP42=P21*DP32+DP21*P32-0.20000000*DP22 P43=P21*P33 DP43=P21*DP33+DP21*P33 P44=P22*P33 DP44=3.0*P43 C2=G(4,2)*CP2+G(1,4)*SP2 C3=G(4,3)*CP3+G(2,4)*SP3 C4=G(4,4)*CP4+G(3,4)*SP4 BR=BR-4.0*AOR*(G(4,1)*P41+C2*P42+C3*P43+C4*P44) BT=BT+AOR*(G(4,1)*DP41+C2*DP42+C3*DP43+C4*DP44) BP=BP-AOR*((G(4,2)*SP2-G(1,4)*CP2)*P42+2.0*(G(4,3)*SP3-G(2,4)*CP3) 1*P43+3.0*(G(4,4)*SP4-G(3,4)*CP4)*P44) IF(NMAX-5) 1,5,5 5 AOR=AOR*AR IF(BM(5)*AOR-ERR) 1,1,105 105 SP5=(SP3+SP3)*CP3 CP5=(CP3+SP3)*(CP3-SP3) P51=P21*P41-0.25714285*P31 DP51=P21*DP41+DP21*P41-0.25714285*DP31 P52=P21*P42-0.22857142*P32 DP52=P21*DP42+DP21*P42-0.22857142*DP32 P53=P21*P43-0.14285714*P33 DP53=P21*DP43+DP21*P43-0.14285714*DP33 P54=P21*P44 DP54=P21*DP44+DP21*P44 P55=P22*P44 DP55=4.0*P54 C2=G(5,2)*CP2+G(1,5)*SP2 C3=G(5,3)*CP3+G(2,5)*SP3 C4=G(5,4)*CP4+G(3,5)*SP4 C5=G(5,5)*CP5+G(4,5)*SP5 BR=BR-5.0*AOR*(G(5,1)*P51+C2*P52+C3*P53+C4*P54+C5*P55) BT=BT+AOR*(G(5,1)*DP51+C2*DP52+C3*DP53+C4*DP54+C5*DP55) BP=BP-AOR*((G(5,2)*SP2-G(1,5)*CP2)*P52+2.0*(G(5,3)*SP3-G(2,5)*CP3) 1*P53+3.0*(G(5,4)*SP4-G(3,5)*CP4)*P54+4.0*(G(5,5)*SP5-G(4,5)*CP5)*P 255) C N= 6 IF(NMAX-6) 1,6,6 6 AOR=AOR*AR IF(BM(6)*AOR-ERR) 1,1,106 106 SP6=SP2*CP5+CP2*SP5 CP6=CP2*CP5-SP2*SP5 P61=P21*P51-0.25396825*P41 DP61=P21*DP51+DP21*P51-0.25396825*DP41 P62=P21*P52-0.23809523*P42 DP62=P21*DP52+DP21*P52-0.23809523*DP42 P63=P21*P53-0.19047619*P43 DP63=P21*DP53+DP21*P53-0.19047619*DP43 P64=P21*P54-0.11111111*P44 DP64=P21*DP54+DP21*P54-0.11111111*DP44 P65=P21*P55 DP65=P21*DP55+DP21*P55 P66=P22*P55 DP66=5.0*P65 C2=G(6,2)*CP2+G(1,6)*SP2 C3=G(6,3)*CP3+G(2,6)*SP3 C4=G(6,4)*CP4+G(3,6)*SP4 C5=G(6,5)*CP5+G(4,6)*SP5 C6=G(6,6)*CP6+G(5,6)*SP6 BR=BR-6.0*AOR*(G(6,1)*P61+C2*P62+C3*P63+C4*P64+C5*P65+C6*P66) BT=BT+AOR*(G(6,1)*DP61+C2*DP62+C3*DP63+C4*DP64+C5*DP65+C6*DP66) BP=BP-AOR*((G(6,2)*SP2-G(1,6)*CP2)*P62+2.0*(G(6,3)*SP3-G(2,6)*CP3) 1*P63+3.0*(G(6,4)*SP4-G(3,6)*CP4)*P64+4.0*(G(6,5)*SP5-G(4,6)*CP5)*P 265+5.0*(G(6,6)*SP6-G(5,6)*CP6)*P66) IF(NMAX-7) 1,7,7 7 AOR=AOR*AR IF(BM(7)*AOR-ERR) 1,1,107 107 SP7=(SP4+SP4)*CP4 CP7=(CP4+SP4)*(CP4-SP4) P71=P21*P61-0.25252525*P51 DP71=P21*DP61+DP21*P61-0.25252525*DP51 P72=P21*P62-0.24242424*P52 DP72=P21*DP62+DP21*P62-0.24242424*DP52 P73=P21*P63-0.21212121*P53 DP73=P21*DP63+DP21*P63-0.21212121*DP53 P74=P21*P64-0.16161616*P54 DP74=P21*DP64+DP21*P64-0.16161616*DP54 P75=P21*P65-0.09090909*P55 DP75=P21*DP65+DP21*P65-0.09090909*DP55 P76=P21*P66 DP76=P21*DP66+DP21*P66 P77=P22*P66 DP77=6.0*P76 C2=G(7,2)*CP2+G(1,7)*SP2 C3=G(7,3)*CP3+G(2,7)*SP3 C4=G(7,4)*CP4+G(3,7)*SP4 C5=G(7,5)*CP5+G(4,7)*SP5 C6=G(7,6)*CP6+G(5,7)*SP6 C7=G(7,7)*CP7+G(6,7)*SP7 BR=BR-7.0*AOR*(G(7,1)*P71+C2*P72+C3*P73+C4*P74+C5*P75+C6*P76+C7*P7 17) BT=BT+AOR*(G(7,1)*DP71+C2*DP72+C3*DP73+C4*DP74+C5*DP75+C6*DP76+C7* 1DP77) BP=BP-AOR*((G(7,2)*SP2-G(1,7)*CP2)*P72+2.0*(G(7,3)*SP3-G(2,7)*CP3) 1*P73+3.0*(G(7,4)*SP4-G(3,7)*CP4)*P74+4.0*(G(7,5)*SP5-G(4,7)*CP5)*P 275+5.0*(G(7,6)*SP6-G(5,7)*CP6)*P76+6.0*(G(7,7)*SP7-G(6,7)*CP7)*P77 3) IF(NMAX-8) 1,8,8 8 AOR=AOR*AR IF(BM(8)*AOR-ERR) 1,1,108 108 SP8=SP2*CP7+CP2*SP7 CP8=CP2*CP7-SP2*SP7 P81=P21*P71-0.25174825*P61 DP81=P21*DP71+DP21*P71-0.25174825*DP61 P82=P21*P72-0.24475524*P62 DP82=P21*DP72+DP21*P72-0.24475524*DP62 P83=P21*P73-0.22377622*P63 DP83=P21*DP73+DP21*P73-0.22377622*DP63 P84=P21*P74-0.18881118*P64 DP84=P21*DP74+DP21*P74-0.18881118*DP64 P85=P21*P75-0.13986013*P65 DP85=P21*DP75+DP21*P75-0.13986013*DP65 P86=P21*P76-0.07692307*P66 DP86=P21*DP76+DP21*P76-0.07692307*DP66 P87=P21*P77 DP87=P21*DP77+DP21*P77 P88=P22*P77 DP88=7.0*P87 C2=G(8,2)*CP2+G(1,8)*SP2 C3=G(8,3)*CP3+G(2,8)*SP3 C4=G(8,4)*CP4+G(3,8)*SP4 C5=G(8,5)*CP5+G(4,8)*SP5 C6=G(8,6)*CP6+G(5,8)*SP6 C7=G(8,7)*CP7+G(6,8)*SP7 C8=G(8,8)*CP8+G(7,8)*SP8 BR=BR-8.0*AOR*(G(8,1)*P81+C2*P82+C3*P83+C4*P84+C5*P85+C6*P86+C7*P8 17+C8*P88) BT=BT+AOR*(G(8,1)*DP81+C2*DP82+C3*DP83+C4*DP84+C5*DP85+C6*DP86+C7* 1DP87+C8*DP88) BP=BP-AOR*((G(8,2)*SP2-G(1,8)*CP2)*P82+2.0*(G(8,3)*SP3-G(2,8)*CP3) 1*P83+3.0*(G(8,4)*SP4-G(3,8)*CP4)*P84+4.0*(G(8,5)*SP5-G(4,8)*CP5)*P 285+5.0*(G(8,6)*SP6-G(5,8)*CP6)*P86+6.0*(G(8,7)*SP7-G(6,8)*CP7)*P87 3+7.0*(G(8,8)*SP8-G(7,8)*CP8)*P88) IF(NMAX-9) 1,9,9 9 AOR=AOR*AR IF(BM(9)*AOR-ERR) 1,1,109 109 SP9=(SP5+SP5)*CP5 CP9=(CP5+SP5)*(CP5-SP5) P91=P21*P81-0.25128205*P71 DP91=P21*DP81+DP21*P81-0.25128205*DP71 P92=P21*P82-0.24615384*P72 DP92=P21*DP82+DP21*P82-0.24615384*DP72 P93=P21*P83-0.23076923*P73 DP93=P21*DP83+DP21*P83-0.23076923*DP73 P94=P21*P84-0.20512820*P74 DP94=P21*DP84+DP21*P84-0.20512820*DP74 P95=P21*P85-0.16923076*P75 DP95=P21*DP85+DP21*P85-0.16923076*DP75 P96=P21*P86-0.12307692*P76 DP96=P21*DP86+DP21*P86-0.12307692*DP76 P97=P21*P87-0.06666666*P77 DP97=P21*DP87+DP21*P87-0.06666666*DP77 P98=P21*P88 DP98=P21*DP88+DP21*P88 P99=P22*P88 DP99=8.0*P98 C2=G(9,2)*CP2+G(1,9)*SP2 C3=G(9,3)*CP3+G(2,9)*SP3 C4=G(9,4)*CP4+G(3,9)*SP4 C5=G(9,5)*CP5+G(4,9)*SP5 C6=G(9,6)*CP6+G(5,9)*SP6 C7=G(9,7)*CP7+G(6,9)*SP7 C8=G(9,8)*CP8+G(7,9)*SP8 C9=G(9,9)*CP9+G(8,9)*SP9 BR=BR-9.0*AOR*(G(9,1)*P91+C2*P92+C3*P93+C4*P94+C5*P95+C6*P96+C7*P9 17+C8*P98+C9*P99) BT=BT+AOR*(G(9,1)*DP91+C2*DP92+C3*DP93+C4*DP94+C5*DP95+C6*DP96+C7* 1DP97+C8*DP98+C9*DP99) BP=BP-AOR*((G(9,2)*SP2-G(1,9)*CP2)*P92+2.0*(G(9,3)*SP3-G(2,9)*CP3) 1*P93+3.0*(G(9,4)*SP4-G(3,9)*CP4)*P94+4.0*(G(9,5)*SP5-G(4,9)*CP5)*P 295+5.0*(G(9,6)*SP6-G(5,9)*CP6)*P96+6.0*(G(9,7)*SP7-G(6,9)*CP7)*P97 3+7.0*(G(9,8)*SP8-G(7,9)*CP8)*P98+8.0*(G(9,9)*SP9-G(8,9)*CP9)*P99) IF(NMAX-10) 1,10,10 10 AOR=AOR*AR IF(BM(10)*AOR-ERR) 1,1,1010 1010 SP10=SP2*CP9+CP2*SP9 CP10=CP2*CP9-SP2*SP9 P101=P21*P91-0.25098039*P81 DP101=P21*DP91+DP21*P91-0.25098039*DP81 P102=P21*P92-0.24705882*P82 DP102=P21*DP92+DP21*P92-0.24705882*DP82 P103=P21*P93-0.23529411*P83 DP103=P21*DP93+DP21*P93-0.23529411*DP83 P104=P21*P94-0.21568627*P84 DP104=P21*DP94+DP21*P94-0.21568627*DP84 P105=P21*P95-0.18823529*P85 DP105=P21*DP95+DP21*P95-0.18823529*DP85 P106=P21*P96-0.15294117*P86 DP106=P21*DP96+DP21*P96-0.15294117*DP86 P107=P21*P97-0.10980392*P87 DP107=P21*DP97+DP21*P97-0.10980392*DP87 P108=P21*P98-0.05882352*P88 DP108=P21*DP98+DP21*P98-0.05882352*DP88 P109=P21*P99 DP109=P21*DP99+DP21*P99 P1010=P22*P99 DP1010=9.0*P109 C2=G(10,2)*CP2+G(1,10)*SP2 C3=G(10,3)*CP3+G(2,10)*SP3 C4=G(10,4)*CP4+G(3,10)*SP4 C5=G(10,5)*CP5+G(4,10)*SP5 C6=G(10,6)*CP6+G(5,10)*SP6 C7=G(10,7)*CP7+G(6,10)*SP7 C8=G(10,8)*CP8+G(7,10)*SP8 C9=G(10,9)*CP9+G(8,10)*SP9 C10=G(10,10)*CP10+G(9,10)*SP10 BR=BR-10.0*AOR*(G(10,1)*P101+C2*P102+C3*P103+C4*P104+C5*P105+C6*P1 106+C7*P107+C8*P108+C9*P109+C10*P1010) BT=BT+AOR*(G(10,1)*DP101+C2*DP102+C3*DP103+C4*DP104+C5*DP105+C6*DP 1106+C7*DP107+C8*DP108+C9*DP109+C10*DP1010) BP=BP-AOR*((G(10,2)*SP2-G(1,10)*CP2)*P102+2.0*(G(10,3)*SP3-G(2,10) 1*CP3)*P103+3.0*(G(10,4)*SP4-G(3,10)*CP4)*P104+4.0*(G(10,5)*SP5-G(4 2,10)*CP5)*P105+5.0*(G(10,6)*SP6-G(5,10)*CP6)*P106+6.0*(G(10,7)*SP7 3-G(6,10)*CP7)*P107+7.0*(G(10,8)*SP8-G(7,10)*CP8)*P108+8.0*(G(10,9) 4*SP9-G(8,10)*CP9)*P109+9.0*(G(10,10)*SP10-G(9,10)*CP10)*P1010) 1 BP=BP/P22*1.E-5 BT=BT*1.E-5 BR=BR*1.E-5 B=SQRT (BR*BR+BT*BT+BP*BP) RETURN END SUBROUTINE INVAR(FLAT,FLONG,ALT,ERR,BB,FL) C C REVISED DEC 1965 C NOTE,ERROR IN L IS TYPICALLY LESS THAN 10.*ERR*L (PERCENT) C FLAT=LATITUDE IN DEGREES , FLONG=LONGITUDE IN DEGREES IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V(3,3),B(200),ARC(200),VN(3),VP(3),BEG(200),BEND(200), 1BLOG(200),ECO(200),R1(3),R2(3),R3(3) SAVE C ALT=ALTITUDE=DISTANCE FROM SURFACE OF EARTH IN KILOMETERS V(1,2)=ALT/6371.2 V(2,2)=(90.-FLAT)/57.2957795 V(3,2)=FLONG/57.2957795 ARC(1)=0. ARC(2)=(1.0+V(1,2))*SQRT (ERR)*0.3 DCLT=1.5708-0.2007*COS (V(3,2)+1.239) IF(V(2,2)-DCLT)10,10,11 11 ARC(2)=-ARC(2) 10 CALL START (R1,R2,R3,B,ARC, V) DO 12 I=1,3 VP(I)=V(I,2) 12 VN(I)=V(I,3) CALL LINES (R1,R2,R3,B,ARC,ERR,J,VP,VN) IF(J-200)16,17,17 17 FL=-1.0 GO TO 18 16 JUP=J DO 40 J=1,JUP ARC(J)=ABS (ARC(J)) 40 BLOG(J)=DLOG(B(J)) JEP=JUP-1 DO 21 J=2,JEP ASUM=ARC(J)+ARC(J+1) DX=BLOG(J-1)-BLOG(J) DN=ASUM*ARC(J)*ARC(J+1) BCO =((BLOG(J-1)-BLOG(J+1))*ARC(J)**2-DX*ASUM**2)/DN CCO =(DX*ARC(J+1)-(BLOG(J)-BLOG(J+1))*ARC(J))/DN SA=.75*ARC(J) SC=SA+.25*ASUM DCO=BLOG(J-1)-CCO *SA*SC ECO(J)=BCO +CCO *(SA+SC) BEG(J)=EXP (DCO+ECO(J)*.5*ARC(J)) 21 BEND(J)=EXP (DCO+ECO(J)*.5*(ASUM+ARC(J))) BEG(JUP)=BEND(JEP) BEND(JUP)=B(JUP) ECO(JUP)=(2.0/ARC(JUP))*DLOG(BEND(JUP)/BEG(JUP)) CALL INTEG (ARC,BEG,BEND,B,JEP,ECO,FLINT) 27 CALL CARMEL (B(2),FLINT,FL) 18 BB=B(2) RETURN END SUBROUTINE LINES (R1,R2,R3,B,ARC,ERR,J,VP,VN) C C VERSION OF MAY 1965 IMPLICIT REAL*8 (A-H,O-Z) DIMENSION B(200),ARC(200),R1(3),R2(3),R3(3),VN(3),VP(3),RA(3) SAVE CRE=0.25 IF(ERR-0.15625)74,75,75 74 CRE= (ERR**0.333333333) 75 A3=ARC(3) AAB=ABS (A3) SNA=A3/AAB A1=ARC(1) A2=ARC(2) AO6=A3*A3/6.0 J=3 ILP=1 IS=1 GO TO 87 66 IS=1 J=J+1 AO6=A3*A3/6.0 ARCJ=A1+A2+A3 AD=(ASUM+A1)/AA BD=ASUM/BB CD=A1/CC 36 DO 5 I=1,3 DD=R1(I)/AA-R2(I)/BB+R3(I)/CC GO TO(6,8),IS 6 RT=R1(I)-(AD*R1(I)-BD*R2(I)+CD*R3(I)-DD*ARCJ)*ARCJ RA(I)=R1(I) R1(I)=R2(I) R2(I)=R3(I) R3(I)=RT VP(I)=VN(I) 8 RBAR=(R2(I)+R3(I))/2.-DD*AO6 5 VN(I)=VP(I)+A3*RBAR 87 IF(VN(2))76,77,77 76 VN(2)=-VN(2) 77 IF(VN(2)-3.141592653)78,78,79 79 VN(2)=6.283185307-VN(2) GO TO 77 78 IF(VN(3))80,81,81 80 VN(3)=VN(3)+6.283185307 GO TO 78 81 IF(VN(3)-6.283185307)82,82,83 83 VN(3)=VN(3)-6.283185307 GO TO 81 82 GO TO (9,10),IS 9 SIT=ABS (SIN (VN(2))) PRE1=VN(1) PRE2=PRE1*VN(2) PRE3=PRE1*SIT*VN(3) C**** SSQ=SIT*SIT C**** OER=(6356.912+SSQ*(21.3677+.108*SSQ))/6371.2 C**** AER=VN(1)-OER C**** CALL MAGNET(AER,SIT,VN(3),BR,BT,BP,B(J),VN(2)) CALL NEWMAG(VN ,SIT,VN(3),BR,BT,BP,B(J),VN(2)) R3(1)=BR/B(J) DN=B(J)*VN(1) R3(2)=BT/DN R3(3)=BP/(DN*SIT) ASUM=A3+A2 AA=ASUM*A2 BB=A3*A2 CC=ASUM*A3 IS=2 GO TO 36 10 SIT=ABS (SIN (VN(2))) B(J)=B(J)*((PRE1/VN(1))**3) 59 QRT=.5*ABS (R3(1))/(.1+ABS (R3(2)*VN(1))) X=(ABS (VN(1)-PRE1)+QRT*ABS (VN(1)*VN(2)-PRE2)+ABS (VN(1)*SIT*VN(3 1)-PRE3))/(AAB*ERR*SQRT (1.+QRT*QRT)) GO TO (90,93,90),ILP 93 IF(X-3.3)90,89,89 89 A3=A3*0.2*(8.0+X)/(0.8+X) J=J-1 ILP=3 ASUM=A2+A1 AA=ASUM*A1 BB=A2*A1 CC=ASUM*A2 DO 91 I=1,3 VN(I)=VP(I) R3(I)=R2(I) R2(I)=R1(I) 91 R1(I)=RA(I) GO TO 73 90 IF(J-200)67,60,60 67 A1=A2 IF(B(J)-B(2))49,49,60 49 ILP=2 A2=A3 A3=A3*.2*(8.+X)/(.8+X) AM=(2.-R3(2)*VN(1))*VN(1)*CRE IF(ABS (A3)-AM)84,84,72 72 A3=SNA*AM 84 IF(SNA*R3(1)+.5)85,85,73 85 AM=-.5*SNA*VN(1)/R3(1) IF(ABS (A3)-AM)73,73,86 86 A3=SNA*AM 73 ARC(J+1)=A3 AAB=ABS (A3) GO TO 66 60 RETURN END SUBROUTINE START (R1,R2,R3,B,ARC, V) C VERSION OF DEC 1965 IMPLICIT REAL*8 (A-H,O-Z) DIMENSION B(200),ARC(200),V(3,3),R1(3),R2(3),R3(3) SAVE SIT=ABS (SIN (V(2,2))) AER=V(1,2) SSQ=SIT*SIT OER=(6356.912+SSQ*(21.3677+.108*SSQ))/6371.2 V(1,2)=AER+OER 10 IF(V(3,2))11,12,12 11 V(3,2)=V(3,2)+6.283185307 GO TO 10 C**12 CALL MAGNET(AER,SIT,V(3,2),BR,BT,BP,B(2),V(2,2)) 12 CALL NEWMAG(V(1,2),SIT,V(3,2),BR,BT,BP,B(2),V(2,2)) R2(1)=BR/B(2) DN=B(2)*V(1,2) R2(2)=BT/DN R2(3)=BP/(DN*SIT) IS=0 1 DO 2 I=1,3 2 V(I,1)=V(I,2)-ARC(2)*R2(I) SIT=ABS (SIN (V(2,1))) C** 3 SSQ=SIT*SIT C** OER=(6356.912+SSQ*(21.3677+.108*SSQ))/6371.2 C** AER=V(1,1)-OER C** CALL MAGNET(AER,SIT,V(3,1),BR,BT,BP,B(1),V(2,1)) 3 CALL NEWMAG(V ,SIT,V(3,1),BR,BT,BP,B(1),V(2,1)) IF(B(1)-B(2))4,5,5 4 ARC(2)=-ARC(2) GO TO 1 5 R1(1)=BR/B(1) ARC(3)=ARC(2) DN=B(1)*V(1,1) R1(2)=BT/DN R1(3)=BP/(DN*SIT) DO 6 I=1,3 6 V(I,1)=V(I,2)-ARC(2)*(R1(I)+R2(I))/2. SIT=ABS (SIN (V(2,1))) IS=IS+1 GO TO (3,7),IS 7 DO 8 I=1,3 8 V(I,3)=V(I,2)+ARC(3)*((1.5)*R2(I)-.5*R1(I)) RETURN END SUBROUTINE INTEG (ARC,BEG,BEND,B,JEP,ECO,FI) C C UNCHANGED MAY 1965 IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ARC(200),BEG(200),BEND(200),B(200),ECO(200) SAVE 4 KK=JEP 6 IF(KK-4)14,11,20 11 KK=KK-1 14 A=B(KK-1)/B(2) X2=B(KK)/B(2) X3=B(KK+1)/B(2) ASUM=ARC(KK)+ARC(KK+1) DN=ARC(KK)*ARC(KK+1)*ASUM BB=(-A*ARC(KK+1)*(ARC(KK)+ASUM)+X2*ASUM**2-X3*ARC(KK)**2)/DN C=(A*ARC(KK+1)-X2*ASUM+X3*ARC(KK))/DN FI=1.570796326*(1.-A+BB*BB/(4.*C))/DSQRT (DABS (C)) RETURN 20 T=DSQRT (1.-BEND(2)/B(2)) FI=(2.*T-DLOG((1.+T)/(1.-T)))/ECO(2) IF(B(2)-BEND(KK))21,21,25 25 KK=KK+1 21 T=SQRT (ABS (1.0-BEG(KK)/B(2))) FI=FI-(2.*T-DLOG((1.+T)/(1.-T)))/ECO(KK) KK=KK-1 22 DO 5 I=3,KK ARG1=1.-BEND(I)/B(2) IF(ARG1)26,26,27 26 TE=1.E-5 GO TO 28 27 TE=SQRT (ARG1) 28 ARG1=1.-BEG(I)/B(2) IF(ARG1)29,29,31 31 TB=SQRT (ARG1) GO TO 32 29 TB=1.E-5 32 IF(ABS (ECO(I))-2.E-5) 23,23,24 23 FI=FI+((TE+TB)*(ARC(I)+ARC(I+1)))/4. GO TO 5 24 FI=FI+(2.*(TE-TB)-DLOG((1.+TE)*(1.-TB)/((1.-TE)*(1.+TB))))/ECO(I) 5 CONTINUE 30 RETURN END SUBROUTINE CARMEL (B,XI,VL) C C UNCHANGED MAY 1965 C COMPUTE L IMPLICIT REAL*8 (A-H,O-Z) SAVE IF(XI-1.0E-36)14,14,15 14 VL=(0.311653/B)**(1./3.) RETURN 15 XX=3.0*DLOG(XI) XX=XX+DLOG(B/0.311653) IF(XX+22.)1,1,8 8 IF(XX+3.)2,2,9 9 IF(XX-3.)3,3,10 10 IF(XX-11.7)4,4,11 11 IF(XX-23.)5,5,6 1 GG=.333338*XX+.30062102 GO TO 7 2 GG=((((((((-8.1537735E-14*XX+8.3232531E-13)*XX+1.0066362E-9)*XX+ 18.1048663E-8)*XX+3.2916354E-6)*XX+8.2711096E-5)*XX+1.3714667E-3)* 2XX+.015017245)*XX+.43432642)*XX+.62337691 GO TO 7 3 GG=((((((((2.6047023E-10*XX+2.3028767E-9)*XX-2.1997983E-8)*XX- 15.3977642E-7)*XX-3.3408822E-6)*XX+3.8379917E-5)*XX+1.1784234E-3)* 2XX+1.4492441E-2)*XX+.43352788)*XX+.6228644 GO TO 7 4 GG=((((((((6.3271665E-10*XX-3.958306E-8)*XX+9.9766148E-07)*XX- 11.2531932E-5)*XX+7.9451313E-5)*XX-3.2077032E-4)*XX+2.1680398E-3)* 2XX+1.2817956E-2)*XX+.43510529)*XX+.6222355 GO TO 7 5 GG=(((((2.8212095E-8*XX-3.8049276E-6)*XX+2.170224E-4)*XX-6.7310339 1E-3)*XX+.12038224)*XX-.18461796)*XX+2.0007187 GO TO 7 6 GG=XX-3.0460681 7 VL=(((1.0+EXP (GG))*0.311653)/B)**(1./3.) C END COMPUTE L RETURN END SUBROUTINE SDC1(TIME,E,IFLAG) C C PREPARE ORBITAL 2-LINE ELEMENTS OBTAINED FROM SPACE DEFENSE CENTER C C ELEMENTS ARE PREPARED FOR USE IN TAYLOR SERIES EXPANSIONS. C C ORIGINAL PROGRAM "SDC" ADAPTED IN 81 APRIL BY JDR BOULDING C C INTO A SUBROUTINE FORMAT WHICH:- C 1. SEARCHES LU 33 FOR A SET OF SATELLITE 2-LINE ORBITAL ELEMENTS C WITH AN EPOCH CLOSEST TO INPUT 'TIME'. C 2. PROCESSES THIS SET AND PLACES COEFFICIENTS IN ARRAY 'E' FOR C SUBSEQUENT USE IN A TAYLOR SERIES EXPANSION. C 3. SET 'IFLAG' TO INDICATE ABSENCE OR PRESENCE OF PROBLEMS. C C NOTE: LU 33 MUST CONTAIN SETS OF 2-LINE ELEMENTS IN ORDER OF C INCREASING EPOCH DATE. C 'TIME' MUST BE IN JULIAN DATE FOR SPACE UNITS. C C UPDATE: Converted from CP6 Fortran to Sun Fortran (1992) on June 6 1993 C by Michael G. Fischer (UCalgary) C CHANGES: Decode statements originally used right justified Hollerith C constant notation (ie. 72R1) this was not ansi standard and was changed C to standard A notation, all statements that used the R notation were C replaced with the ASCII equivalent. C Output statements were not ansi std. and were replaced with print * statements C IMPLICIT REAL*8 (B-H,O-Z) SAVE CHARACTER*72 IN1CT,IN2T CHARACTER*72 IN1PT DIMENSION E(27),IN1P(18),IN1C(18),IN2(18),IM(72), 1 TS(27),ITS(54),AS(54) EQUIVALENCE (TS,ITS,AS) EQUIVALENCE(IN1PT,IN1P) EQUIVALENCE(IN1CT,IN1C) EQUIVALENCE(IN2T,IN2) DATA DEG/57.2957795131/,TWOPI/6.28318530717959/ DATA TP/0.0/ C C SPACETRACK EARTH MODEL 1968 (STEM 68R) CONSTANTS DATA RE/6378.145/,F/298.25/ DATA U/398601.2/,TJ2/0.001082549/ C 1000 FORMAT (18A4) c1002 FORMAT (2X,I2,F12) 1002 FORMAT (2X,I2,F12.9) C C C IF FIRST ENTRY GO TO ELEMENT FILE SEARCH ROUTINE IF (TP.EQ.0.0) GO TO 50 T = ABS(TIME-TP) TT = ABS(TIME-TPP) IF (T.LE.TT) GO TO 20 C RESTART SEARCH FROM BEGINNING 10 REWIND 33 TP = 0.0 GO TO 50 C C IF NOT ADVANCING GO TO BOF TEST 20 IF (ABS(TIME-TC).GT.T) GO TO 25 C CONTINUE SEARCH IF EOF HAS NOT BEEN REACHED IF (ABS(TC-TP).GT.0.001) GO TO 30 C ERROR EXIT IF LAST SET OF ELEMENTS CANNOT BE USED IF (2.0*T.GT.TT) GO TO 710 C USE, BUT NO NEED TO PROCESS AGAIN, THE LAST SET OF ELEMENTS WRITE (*,1210) GO TO 700 C C IF NOT AT BOF, PROCESS ELEMENTS 25 IF (TPP.NE.0.0) GO TO 100 C ERROR EXIT IF FIRST SET OF ELEMENTS CANNOT BE USED IF (TIME.GE.TP) GO TO 100 IF (ABS(TIME-TC).LT.2.0*T) GO TO 740 C USE THE FIRST SET OF ELEMENTS WRITE (*,1200) GO TO 100 C C LOOP TO SEARCH FOR PROPER SET OF ELEMENTS C C PUSH TIME STACK 30 TPP = TP TP = TC DO 39, I=1,18 39 IN1P(I) = IN1C(I) C READ LINE 2 OF ELEMENT SET --- PROCESS LATER READ (33,1000) IN2 C C INITIAL ENTRY POINT IN SEARCH LOOP C READ LINE 1 OF ELEMENT SET AND GET EPOCH TIME 50 READ (33,'(a)',END=90) IN1CT C "DECODE" NOT ALLOWED IN MICROSOFT FORTRAN. CHANGED TO INTERNAL READS. read(in1cT(17:32),1002) iy,tc ccc DECODE (16,1002,IN1C(5)) IY,TC C type*,' iy,tc ',iy,tc C CONVERT CURRENT EPOCH YEAR AND DAY TO JDS TC = ((1461*IY)-1)/4 + TC - 21080.0 C LOOP BACK IF CLOSEST EPOCH TIME HAS NOT BEEN FOUND IF (ABS(TIME-TC).LT.ABS(TIME-TP)) GO TO 30 C PROCESS ELEMENTS IN 'IN1P' & 'IN2' IF NOT AT FIRST SET IF (TPP.NE.0.0) GO TO 100 C ERROR EXIT IF FIRST SET OF ELEMENTS CANNOT BE USED IF (ABS(TIME-TC).LT.ABS((TIME-TP)*2.0)) GO TO 740 C WARN OF USE OF FIRST SET OF ELEMENTS WRITE (*,1200) 1200 FORMAT (' FIRST SET OF ELEMENTS USED') C CORRECT ELEMENTS IN 'IN1P' AND 'IN2' GO TO 100 C C END-OF-FILE READ ON ELEMENT FILE LU 33 C 90 WRITE (*,1205) 1205 FORMAT (' EOF REACHED ON ELEMENT FILE') C ERROR EXIT IF LAST ELEMENT SET CANNOT BE USED IF (TIME.LE.TP) GO TO 100 IF (ABS(TIME-TP).GT.ABS((TIME-TPP)/2.0)) GO TO 710 WRITE (*,1210) 1210 FORMAT (' LAST ELEMENT SET USED') C C PROCESS ELEMENTS C DO LINE 1 C "DECODE" NOT ALLOWED IN MICROSOFT FORTRAN. CHANGED TO INTERNAL READS. 100 READ(IN1PT(1:72),1020) (IM(I),I=1,72) c100 DECODE(72,1020,IN1P) (IM(I),I=1,72) C WRITE(*,*) ' IN1PT ' C WRITE(*,'(1X,A72)') IN1PT C WRITE(*,*) ' IM ' C WRITE(*,'(1X,72A1)') (IM(I),I=1,72) 1020 FORMAT(72A1) IF (IM(1).EQ.'1'.AND.IM(53).EQ.' ')GO TO 820 IF (IM(1).EQ.'2'.AND.IM(52).EQ.' ')GO TO 820 WRITE(*,1011) 1011 FORMAT(5X,'FORMAT ERROR') PRINT *,'IM(1)= ' PRINT *, IM(1) PRINT *,'IM(52)= ' PRINT *,IM(52) PRINT *,'IM(53)= ' PRINT *,IM(53) C OUTPUT IM(1),IM(52),IM(53) C CONTINUE IF CONFIRMED THAT LINE 1 IS BEING PROCESSED. 820 IF(IM(1).EQ.'1')GO TO 120 102 WRITE(*,1220) 1220 FORMAT(' SEQUENCE ERROR IN ELEMENT FILE --SEARCH RESTARTED') GO TO 10 C 120 IF(IM(17).EQ.' ')GO TO 130 IM(15)=IM(16) IM(16)=IM(17) 130 IF(IM(16).GE.'0')GO TO 140 IF(IM(16).LT.'J')GO TO 140 IM(15)=IM(15)+1 IF(IM(16).LT.'S')THEN IM(16)=IM(16)-1 GO TO 140 END IF IF(IM(16).EQ.'S')THEN IM(16)=9 GO TO 140 END IF IM(15)=IM(15)+1 IM(16)=IM(16)-3 C "DECODE" NOT ALLOWED IN MICROSOFT FORTRAN. CHANGED TO INTERNAL READS. 140 READ(IN1PT(1:52),1101)ITS(18),ITS(1),IY,TS(4),TS(27),TNDDOT C140 DECODE(52,1101,IN1P)ITS(18),ITS(1),IY,TS(4),TS(27),TNDDOT C TYPE*,' ITS(18),ITS(1),IY,TS(4),TS(27),TNDDOT ',ITS(18),ITS(1), C 1 IY,TS(4),TS(27),TNDDOT c THE VAX REQUIRES THE Fw.d format not just Fw c1101 FORMAT(2(2X,I5),4X,I2,F12,F11,E9.5) 1101 FORMAT(2(2X,I5),4X,I2,F12.8,F11.8,E9.5) C ITS(1)=ITS(1)*100+IAND(IM(15),1ZF)*10+IAND(IM(16),1ZF) ITS(1)=ITS(1)*100+IAND(IM(15),15)*10+IAND(IM(16),15) IY=((1461*IY)-1)/4 TS(4)=TS(4)-21080.D0+DBLE(IY) C C PROCESS LINE 2 C C RESTART SEARCH IF LINE 2 IS NOT BEING PROCESSED C "DECODE" NOT ALLOWED IN MICROSOFT FORTRAN. CHANGED TO INTERNAL READS. READ(IN2T(1:72),1020) (IM(I),I=1,72) C DECODE(72,1020,IN2) (IM(I),I=1,72) C WRITE(*,*) ' IN2T ' C WRITE(*,'(1X,A72)') IN2T C WRITE(*,*) ' IM ' C WRITE(*,'(1X,72A1)') IM IF(IM(1).EQ.'1'.AND.IM(53).EQ.' ')GO TO 821 IF(IM(1).EQ.'2'.AND.IM(52).EQ.' ')GO TO 821 WRITE(*,1011) PRINT *,'IM(1)= ' PRINT *, IM(1) PRINT *,'IM(53)= ' PRINT *,IM(53) PRINT *,'IM(52)= ' PRINT *,IM(52) C OUTPUT IM(1),IM(53),IM(52) 821 IF (IM(1).NE.'2') GO TO 102 C "DECODE" NOT ALLOWED IN MICROSOFT FORTRAN. CHANGED TO INTERNAL READS. READ(IN2T(1:68),1102) I,TS(12),TS(13),TS(11),TS(14),TS(15), 1 TS(21),ITS(5) C DECODE (68,1102,IN2) I,TS(12),TS(13),TS(11),TS(14),TS(15), C 1 TS(21),ITS(5) C WRITE(*,*) ' IN2T AGAIN ' C WRITE(*,'(1X,A72)') IN2T C WRITE(*,*) ' I,TS(12),TS(13),TS(11),TS(14),TS(15), ETC ' C WRITE(*,1102) I,TS(12),TS(13),TS(11),TS(14),TS(15), C 1 TS(21),ITS(5) 1102 FORMAT (2X,I5,2F9.4,F8.7,2F9.4,F12.9,I5) c1102 FORMAT (2X,I5,2F9,F8.7,2F9,F12,I5) C ERROR EXIT IF LINES 1 & 2 CATALOGUE NUMBERS ARE NOT EQUAL IF (I.NE.ITS(18)) GO TO 730 DO 239, I=12,15 239 TS(I) = TS(I)/DEG TS(8) = 1440.0/TS(21) AS(17) = -TS(27)*2880.0/(TS(21)**2) TS(21) = TS(21)*TWOPI TS(27) = TS(27)*TWOPI TNDDOT = TNDDOT*TWOPI C C COMPLETE PROCESSING C C FILL SATELLITE NAME AND CURRENT DATE ITS(2) = ' ' ITS(3) = ' ' ITS(4) = ' ' ITS(6) = 0 TS (5) = RE ITS(11) = ' KM ' AS (12) = F ITS (13) = ' SDC' ITS (14) = 0 T = (((86400.0/TS(21))**2*U)**0.33333333)/RE D = -1.5*TJ2*(1.0-1.5*SIN(TS(12))**2)/(T*T*(1.0-TS(11)**2)**1.5) TS(10) = T*(1.0+(D-D*D)/3.0) TS(16) = -1.33333333*TS(10)*TS(27)/TS(21) TS(17) = -1.33333333*(1.0-TS(11))*TS(27)/TS(21) TS(18) = 0.0 TS(24) = 0.0 T = 1.5*TJ2*TS(21)/(TS(10)*(1.0-TS(11)**2))**2 TS(19) = -T*COS(TS(12)) TS(20) = T*(2.0-2.5*SIN(TS(12))**2) TS(22) = (1.1111111*(TS(27)**2)/TS(21)-TNDDOT)*2.0*TS(10)/TS(21) TS(23) = (0.22222222*TS(27)**2/TS(21)-TNDDOT)* 1 2.0*(1.0-TS(11))/TS(21) T = (1.0+1.33333333*(1.0-TS(11))/(1.0+TS(11)))*TS(27)/TS(21) TS(25) = TS(19)*T TS(26) = TS(20)*T C MOVE DATA TO CALLING PROGRAM STORAGE DO 669, I=1,27 669 E(I) = TS(I) C C EXIT SECTION C C ERROR FREE NORMAL EXIT 700 IFLAG = 0 C WRITE(*,777) (E(I),I=1,27) C777 FORMAT(E13.6) 705 RETURN C C END-OF-FILE ENCOUNTERED ON ELEMENT FILE 710 WRITE (*,1710) 1710 FORMAT (' EOF ON ELEMENT FILE --- UNABLE TO PROCESS') IFLAG = 10 RETURN C C INCONSISTENT CATALOGUE NUMBER 730 WRITE (*,1730) 1730 FORMAT (' ERROR IN CATALOGUE NUMBER --- ELEMENTS NOT PROCESSED') IFLAG = 30 RETURN C C FIRST ELEMENT EPOCH TOO LATE TO BE USED 740 WRITE (*,1740) 1740 FORMAT (' ELEMENT FILE STARTS TOO LATE --- NOT PROCESSED') IFLAG = 40 RETURN END FUNCTION JDS(IDA,IMO,IYR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION M(12) DATA (M(I),I=1,12) / 0,31,59,90,120,151,181,212,243,273,304,334/ JDS=(1461*IYR)/4 + M(IMO)+IDA-21080 IF(IMO.GT.2) GO TO 1 IF((IYR/4)*4.EQ.IYR) JDS=JDS-1 1 RETURN END FUNCTION FDAY(IHR,IMIN,SEC) IMPLICIT REAL*8 (A-H,O-Z) IMIN=ISIGN(IMIN,IHR) SEC=SIGN(SEC,DBLE(IHR)) TM=(60.*DBLE(IHR)+DBLE(IMIN))*60. FDAY=(TM+SEC)/86400. RETURN END SUBROUTINE ELMENT(DT,A0,A1,A2,A) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A0(6),A1(6),A2(6),A(6) DO 10 I=1,6 10 A(I)=(A2(I)*DT+A1(I))*DT+A0(I) DO 20 I=3,6 20 A(I)=ALLOT(A(I)) RETURN END FUNCTION RAG(TM) IMPLICIT REAL*8 (A-H,O-Z) COMMON/CCCRAG/CR1,CR2,CR3,CR4 DATA CR1/.9906780819 /, CR2/100.0021371 / X , CR3/ 1.075231482E-6 /, CR4/1.0027379093/ DOUBLE CPRECISION RG,T JDS=TM FDAY=TM-JDS TWOPI=6.28318530 T=JDS/36525. RG=(CR3*T+CR2)*T+CR1 + CR4*FDAY IRAG=RG RAG=(RG-IRAG)*TWOPI RETURN END SUBROUTINE ELMEFX(A,RAG,X,Y,Z) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(6) A4=A(4)-RAG X3=DSQRT(1.0-A(2)**2) 3 E=XKEP(A(6),A(2),SE,CE) SE=SIN(E) CE=COS(E) SA=SIN(A(5) ) CA=COS(A(5) ) SC=SIN(A4) CC=COS(A4) SB=SIN (A(3)) CB=COS (A(3)) C C COMPUTE POSITION COORDINATES C Q1 = A(1)*(CE-A(2)) Q2 = A(1)*X3*SE V = Q1*CA - Q2*SA W = Q2*CA + Q1*SA Z = CB*W X=CC*V-SC*Z Y=CC*Z + SC*V Z=SB*W RETURN END SUBROUTINE CTNPOLR(X,Y,Z,R,CLON,CLAT) IMPLICIT REAL*8 (A-H,O-Z) DATA TWOPI/6.283185307179/ R=SQRT(X*X+Y*Y+Z*Z) CLON=ATAN2(Y,X) IF(CLON.LT.0) CLON=CLON+TWOPI CLAT=ASIN(Z/R) RETURN END SUBROUTINE GCGD(R,GCLAT,ALT,GDLAT) C CONVERTS FROM GEOCENTRIC TO GEODETIC COORDINATES C INPUT R GEOCENTRIC RADIUS IN EARTH RADII C GCLAT GEOCENTRIC LATITUDE IN RADINNS C OUPUT ALT GEODETIC HEIGHT ABOVE SURFACE EARTH RADII C GDLAT GEODETIC LATITUDE IN RADIANS C ACCURACY 1.E-5 EARTH RADII OR RADIANS C GCTGD IS MORE ACCURATE, BUT ABOUT THREE TIMES SLOWER IMPLICIT REAL*8 (A-H,O-Z) COMMON/CCCGCD/ESQU,F DATA ESQU/.00669454/ R1=ESQU/2. R2=.375*ESQU*ESQU SINPHI=SIN(GCLAT) COSPHI=COS(GCLAT) CSPHSQ=COSPHI*COSPHI DELPRM=ESQU*SINPHI*COSPHI/(1.-ESQU*CSPHSQ ) RE=(6356.7747/6378.160)*((CSPHSQ*R2+R1)*CSPHSQ+1.) DELTA=DELPRM*RE/R GDLAT=GCLAT+DELTA ALT=(R-RE)*(1.-DELTA*DELPRM/2.) RETURN END SUBROUTINE GCTGD(R,AL1,ALTZ, GDLATZ) IMPLICIT REAL*8 (A-H,O-Z) COMMON/CCCGCD/ESQU,F DATA F/.00335289/ C 5 E2=F+F-F*F 6 E4=E2*E2 7 E6=E2*E4 8 E8=E4*E4 10 R1=1./R 11 G1=1.0+E2 12 G2=4.0+3.0*E2 13 G3=35.0*E2 14 G4=E2*R1/32.0 15 G5=G4*E4 16 G6=E2/64.0 17 G7=E4*G6 18 G8=E4*R1 19 A2=G4*((512.0+E2*(128.0+E2*(60.0+G3)))/32.0+G8*(1.0+E2-.375*G2*R1) 1) 20 A4=E4*R1*.0625*(R1*((4.+E2*(1.+G1))+G8*(.9375-R1))-G6*(48.+G3)-1.0 1) 21 A6=G5*(.09375*(4.+5.*E2)+R1*(1.4583333333 *R1*G2-3.*G1)) 22 A8=G6*G5*(R1*(64.-R1*(252.+320.*R1))-5.) 23 AL2=AL1+AL1 24 AL4=AL2+AL2 25 AL6=AL2+AL4 26 AL8=AL4+AL4 AL= AL1+A2* SIN(AL2)+A4* SIN(AL4)+A6* SIN(AL6)+A8* SIN(AL8) 28 A=AL-AL1 CA= COS(A) 30 RCA=R*CA SAL= SIN(AL) 32 SAL2=SAL*SAL 33 ESAL=E2*SAL2 34 X =SQRT(1.-ESAL) 36 ALTZ=RCA-X 37 GDLATZ=AL 38 RETURN END C C CORRECTED TO PROVIDE AN ACCURACY OF 0.01 DEG C FOR PERIOD 1950 TO 1999 C ALGORITHM AND DATA TAKEN FROM THE 1981 ASTRONOMICAL ALMONAC C PAGE C1 AND C20 C SUBROUTINE XYZSUN(TM ,RAG,XSUN,YSUN,ZSUN) IMPLICIT REAL*8 (A-H,O-Z) DATA TMS/8505./, CS1/4.8796/, CS2/0.0172027915/, X CS31/0.03344/,CS32/0.00035/, CS4/6.2302/, CS5/0.01720197/, X CS6/3.495198029E-4/, COSEPS/.9174523215/, SINEPS/ X .3978457463/ T=TM-TMS G=CS4+CS5*T SNLONG=CS1+CS2*T+CS31*SIN(G)+CS32*SIN(2*G) COSLNG=COS(SNLONG) SINLNG=SIN(SNLONG) COSRAG=COS(RAG) SINRAG=SIN(RAG) TEMP1=SINLNG*COSEPS ZSUN=SINLNG*SINEPS XSUN=COSLNG*COSRAG+TEMP1*SINRAG YSUN=TEMP1*COSRAG-COSLNG*SINRAG RETURN END C OLD CHARACTER SET FUNCTION ALLOT(X) IMPLICIT REAL*8 (A-H,O-Z) DATA TWOPI/6.28318530717959/ ALLOT=DMOD(X,TWOPI) IF(X.LT.0) ALLOT=ALLOT+TWOPI RETURN END SUBROUTINE BRWR(A110,DT,A) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION EL(13),S(3),D(23),B(9),G(13),C(10),F(9),A110(6),GP(5), X A11(6),X(28),A(6) COMMON/CCCBRW/ GP DATA GP/1.1467823E+4,5.412400E-4,2.560E-6,.6900E-6, X0.6E-7 / C C COMPUTE INTERMEDIATE QUANTITIES C 5 A2=A110(1)**2 6 A3=A110(1)*A2 F(1) = A110(2) 7 F(2)=A110(2)*A110(2) 8 F(3)=F(2)+F(2) 9 F(4)=F(2)+F(3) 10 F(5)=9.0*F(2) 11 F(6)=4.0+F(4) 12 F(7)=2.0+F(2) 13 B(1)=1.0-F(2) 14 B(2)=SQRT(B(1)) 15 B(3)=B(1)*B(2) 16 B(4)=B(3)*B(3) 17 B(5)=9.0*B(1) 18 B(6)=25.0*B(1) 19 B(7)=126.0*B(1) 20 B(8)=1.0/B(3) 21 B(9)=1.0/(B(1)*B(1)) 22 D(1)=COS(A110(3)) 23 D(2)=D(1)+D(1) 24 D(3)=D(1)*D(1) 25 D(4)=D(3)*D(3) 26 C(1)=GP(2)/A2 27 C(2)=1.5*C(1)*B(9) 28 C(3)=.41666666667E-1*C(2)*C(2) 29 C(4)=.9375*GP(4)/(A2*A2*B(1)*B(4)) 30 C(5)=.8333333333E-1*C(2) 31 C(6)=.6666666667*C(4)/C(2) 32 C(7)=.25*GP(3)/(GP(2)*A110(1)*B(1)) 33 C(8)=.46875*GP(5)/(GP(2)*A3*B(4)) 34 C(9)=.16666666667*C(2) 35 C(10)=C(9)*D(2) 36 D(5)=1.0-D(3) 37 D(6)=3.0*D(3)-1.0 38 D(7)=5.0*D(3)-1.0 39 D(8)=8.0*D(4)/D(7) 40 D(9)=5.0*D(8)/D(7) 41 D(10)=D(8)-D(6) 42 D(11)=C(5)*(5.0*D(10)-4.0+4.0*D(3))-C(6)*D(10) 43 D(12)=.16666666667*C(8)*(3.0*D(10)-2.0) 44 D(13)=.19444444*C(8)*(D(10)+D(10)-D(5)) 45 D(14)=3.0-16.0*D(3)/D(7)+D(9) 46 D(15)=C(8)*D(14) 47 D(16)=F(6)*D(15) 48 D(17)=.33333333333*D(13) 49 D(18)=.12962962963*(D(15)+D(15)-C(8)) 50 D(19)=3.0*D(5) 51 D(20)=D(6)+D(6) 52 D(21)=2.0-D(7) 53 D(22)=C(10)*SQRT(D(5)) 54 D(23)=4.0*C(9)*D(5) 55 F(8)=D(3)*(2.0+F(4)) 56 F(9)=(2.0+5.0*F(2))*D(8)-D(9)*F(3)*D(3) 57 G(1)=SIN(A110(3)) 58 G(2)=G(1)*G(1) 59 G(3)=A110(2)*D(1) 61 G(5)=G(3)/G(1) 60 G(4)=G(3)*D(1) 62 G(6)=G(4)/G(1) 63 G(7)=B(1)*G(1) 64 G(8)=G(7)/A110(2) 65 G(9)=A110(2)*G(1) 66 G(10)=G(3)*G(3) 67 G(11)=B(1)/A110(2) 68 G(12)=.5*G(11) 69 G(13)=G(11)*C(9) C C COMPUTE ''MEAN'' MEAN MOTION C 70 EN0=SQRT(GP(1)/A3) C C COMPUTE COEFFICIENTS OF SECULAR TERMS C 71 S(1)=EN0*(1.0+B(2)*(C(2)*D(6)+C(3)*(16.0*B(2)-15.0+B(6)+(30.0-96.0 1*B(2)-90.0*B(1))*D(3)+(105.0+144.0*B(2)+B(6))*D(4))+C(4)*F(2)*(3.0 2-30.0*D(3)+35.0*D(4)))) 72 S(2)=EN0*(C(2)*D(7)+C(3)*(24.0*B(2)-35.0+B(6)+(90.0-192.0*B(2)-B(7 1))* D(3)+(385.0+360.0*B(2)+45.0*B(1))*D(4))+.33333333*C(4)*(21.0-B 2(5)+ (B(7)-270.0)*D(3)+(385.0-189.0*B(1))*D(4))) 73 S(3)=EN0*(4.0*D(1)*(C(3)*(B(5)-5.0+12.0*B(2)-D(3)*(35.0+36.0*B(2)+ 15.0* B(1)))+.33333333*C(4)*(5.0-3.0*B(1))*(3.0-7.0*D(3)))-C(2)*D( 22)) C C COMPUTE COEFFICIENTS OF LONG PERIOD TERMS C 74 EL(1)=B(3)*D(11) 75 EL(2)=B(2)*G(8)*(-C(7)-(4.0+F(5))*D(12)) 76 EL(3)=G(9)*B(3)*D(13) 77 EL(4)=.5*C(6)*(F(7)-3.0*F(8)+F(9))-.41666667E-1*C(2)*(F(7)-11.0* 1 F(8)+5.0*F(9)) 78 EL(5)=C(7)*(G(1)/F(1)-G(6))+D(12)*((G(8)-G(6))*F(6)+G(9)*(26.0+ 1F(5)))-G(4)*G(1)*D(16) 79 EL(6)=G(9)*(D(18)*G(10)-D(17)*(3.0+F(3)-G(10)/G(2))) 80 EL(7)=G(3)*A110(2)*(C(6)*D(14)-C(5)*(5.0*D(14)-4.0)) 81 EL(8)=G(5)*(C(7)+F(6)*D(12)+G(2)*D(16)) 82 EL(9)=G(5)*F(2)*(-D(17)-D(18)*G(2)) 83 EL(10)=B(1)*F(1)*D(11) 84 EL(11)=G(7)*(C(7)+F(6)*D(12)) 85 EL(12)=-F(2)*G(7)*D(13) 86 EL(13)=-G(3)/G(7) ENTRY BRWR2(A110,DT,A) C C START PROGRAM ********************************* C C C C COMPUTE SECULAR TERMS C A11(6)=ALLOT(A110(6)+S(1)*DT) A11(5)=ALLOT(A110(5)+S(2)*DT) A11(4)=ALLOT(A110(4)+S(3)*DT) A11(3)=A110(3) A11(2)=A110(2) A11(1)=A110(1) C C COMPUTE LONG PERIOD TERMS C X(1)=A11(5)+A11(5) X(2)=X(1)+A11(5) X(3)=SIN(X(1)) X(4)=COS(A11(5)) X(5)=COS(X(2)) AL1 =A11(6)+EL(1)*X(3)+EL(2)*X(4)+EL(3)*X(5) G1 =A11(5)+EL(4)*X(3)+EL(5)*X(4)+EL(6)*X(5) H1 =A11(4)+EL(7)*X(3)+EL(8)*X(4)+EL(9)*X(5) D1E = EL(10)*COS(X(1))+EL(11)*SIN(A11(5))+EL(12)*SIN(X(2)) D1I=EL(13)*D1E AL1=ALLOT(AL1) G1 =ALLOT(G1) H1 =ALLOT(H1) C C COMPUTE SHORT PERIOD TERMS C SA =XKEP(AL1,A110(2),X(6),X(7)) X(8)=1.0-A110(2)*X(7) X(9)=1.0/X(8) X(10)=X(9)*X(9) X(11)=X(9)*X(10) X(12)=X(6)*X(9)*B(2) X(13)=X(9)*(X(7)-A110(2)) X(14)=ATAN2(X(12),X(13)) IF(X(14))370,340,340 370 X(14)=X(14)+6.2831853072 340 X(15)=G1+G1 X(16)=X(15)+X(14) X(17)=X(16)+X(14) X(18)=X(17)+X(14) X(19)=SIN(X(16)) X(20)=SIN(X(18)) X(21)=COS(X(17)) X(22)=D(6)*(X(11)-B(8)) X(23)=D(19)*X(21) X(24)=A110(2)*(3.0*COS(X(16))+COS(X(18))) X(25)=X(9)+B(1)*X(10) X(26)=G(13)*(D(20)*X(12)*(X(25)+1.0)+D(19)*((1.0-X(25))*X(19)+ 1 (X(25)+.33333333)*X(20))) X(27)=6.0*(X(14)-AL1+A110(2)*X(12)) X(28)=3.0*(SIN(X(17))+A110(2)*X(19))+A110(2)*X(20) C COMPUTE OSCULATING ELEMENTS C A(1)=A110(1)*(1.0+C(1)*(X(22)+X(11)*X(23))) A(2) = A110(2)+D1E+G(12)*(C(1)*(X(22)+X(23)*(X(11)-B(9)))-D(23)* 1X(24)) A(3)=A110(3)+D1I+D(22)*(3.0*X(21)+X(24)) UL= AL1-B(2)*X(26) UG= G1+X(26)+C(9)*(D(7)*X(27)+D(21)*X(28)) UH= H1-C(10)*(X(27)-X(28)) A(6)=ALLOT(UL) A(5)=ALLOT(UG) A(4)=ALLOT(UH) RETURN END FUNCTION XKEP(AM,ECC,SE,CE) C IMPLICIT REAL*8 (A-H,O-Z) COMMON/CCCKEP/ ERR DATA ERR/ 1.E-8/ E=AM+ECC*SIN(AM)+ .5*ECC*ECC*SIN(AM+AM) 1 SE = SIN (E) CE= COS (E) DELM=AM-E+ECC*SE E=E+DELM/(1.-ECC*CE) IF ( ABS ( DELM ).GT.ERR) GO TO 1 2 XKEP=E RETURN END SUBROUTINE READ_SCALED(IXPTS,IYR,IDOY,IHR,IMIN,rSEC,FREQX,AP_RANX) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (MAX=50) character*40 oscale CHARACTER*73 HEADER INTEGER AMPX(MAX),AMPO(MAX) DIMENSION AP_RANX(MAX),FREQX(MAX),DELTIMX(MAX),RIOTIMX(MAX) DIMENSION AP_RANO(MAX),FREQO(MAX),DELTIMO(MAX),RIOTIMO(MAX) COMMON/FILE/OSCALE common/times/times(50),iyrs(50),imos(50),idays(50),idoys(50) READ(13,100) IYR,IMO,IDAY,IHR,IMIN,RSEC WRITE(14,100) IYR,IMO,IDAY,IHR,IMIN,RSEC write(7,100) IYR,IMO,IDAY,IHR,IMIN,RSEC c write(8,100) IYR,IMO,IDAY,IHR,IMIN,RSEC c write(60,100) IYR,IMO,IDAY,IHR,IMIN,RSEC c write(9,100) IYR,IMO,IDAY,IHR,IMIN,RSEC c write(63,100) IYR,IMO,IDAY,IHR,IMIN,RSEC c write(10,100) iyr,imo,iday,ihr,imin,rsec c write(62,100) iyr,imo,iday,ihr,imin,rsec JSEC=RSEC IYR=IYR-((IYR/100)*100) jyr=iyr jmo=imo jday=iday jhr=ihr jmin=imin tsec=rsec C CONVERT MONTH AND DAY TO DAY OF YEAR CALL MD2DOY(IYR,IMO,IDAY,IDOY) jdoy=idoy C TYPE*,' DOY = ',IDOY READ(13,101) IXPTS C WRITE(14,101) IXPTS READ(13,102) HEADER WRITE(14,105) C WRITE(14,102) HEADER DO K=1,IXPTS READ(13,*) AP_RANX(K),FREQX(K),DELTIMX(K),RIOTIMX(K),AMPX(K) C CHECK FOR NEGATIVE SCALED FREQ. THIS IS A SIGN THAT FREQ INTERPOLATION C WAS NOT PERFORMED WHEN IONO FILE GENERATED. STOP RUN IF(FREQX(K).LT.0) THEN WRITE(*,*) ' WE HAVE A NEGATIVE FREQ ',FREQX(K) WRITE(*,*) ' STOPPING RUN ' STOP ENDIF WRITE(14,103) AP_RANX(K),FREQX(K),DELTIMX(K),RIOTIMX(K),AMPX(K) C WRITE(*,103) AP_RANX(K),FREQX(K),DELTIMX(K),RIOTIMX(K),AMPX(K) ENDDO c CHECK TO SEE IF INPUT FREQ ARE MONOTONICALLY INCREASING. REMOVE ANY C VIOLATORS. CALL CKFREQ(IXPTS,AP_RANX,FREQX,deltimx,riotimx,ampx) do k=1,ixpts call time_inc(jyr,jmo,jday,jhr,jmin,tsec, 1 iyr,imo,iday,ihr,imin,rsec,riotimx(k)) c type*,' from read_scaled ',iyr,imo,iday,ihr,imin,rsec CALL MD2DOY(IYR,IMO,IDAY,IDOY) jsec=rsec times(k)=ihr+dfloat(imin)/60+rsec/3600. iyrs(k)=iyr imos(k)=imo idays(k)=iday idoys(k)=idoy if(k.eq.1) then iyr1=iyr imo1=imo iday1=iday ihr1=ihr imin1=imin rsec1=rsec jsec=rsec1 endif enddo write(8,100) IYR1,IMO1,IDAY1,IHR1,IMIN1,RSEC1 write(9,100) IYR1,IMO1,IDAY1,IHR1,IMIN1,RSEC1 write(10,100) iyr1,imo1,iday1,ihr1,imin1,rsec1 write(60,100) IYR1,IMO1,IDAY1,IHR1,IMIN1,RSEC1 write(63,100) IYR1,IMO1,IDAY1,IHR1,IMIN1,RSEC1 write(62,100) iyr1,imo1,iday1,ihr1,imin1,rsec1 READ(13,101,END=999) IOPTS C ************************************************** C IF O TRACE NOT SCALED SKIP THIS SECTION. IF(IOPTS.EQ.0) GOTO 999 OPEN(UNIT=15,FILE=OSCALE,STATUS='UNKNOWN',FORM='FORMATTED') c WRITE(15,100) IYR,IMO,IDAY,IHR,IMIN,RSEC c write(61,100) IYR,IMO,IDAY,IHR,IMIN,RSEC write(15,100) jyr,jmo,jday,jhr,jmin,tsec write(61,100) jyr,jmo,jday,jhr,jmin,tsec WRITE(15,106) C WRITE(*,101) IOPTS READ(13,102) HEADER C WRITE(15,102) HEADER DO J=1,IOPTS READ(13,*) AP_RANO(K),FREQO(K),DELTIMO(K),RIOTIMO(K),AMPO(K) WRITE(15,103) AP_RANO(K),FREQO(K),DELTIMO(K),RIOTIMO(K),AMPO(K) WRITE(61,103) AP_RANO(K),FREQO(K),DELTIMO(K),RIOTIMO(K),AMPO(K) ENDDO 999 CONTINUE c reset date and time to original values iyr=jyr idoy=jdoy ihr=ihr1 imin=imin1 jsec=rsec1 rsec=rsec1 C ************************************************************** 100 FORMAT(5I12,F16.7) 101 FORMAT(I8) 102 FORMAT(A) 103 FORMAT(F16.6,F16.7,F16.8,F16.7,I8) 105 FORMAT(' X_SCALED '/' X_SCALE_HGT',12X,'FREQ',10X, 1 'DEL_TIM',10X,'IO_TIM',6X,'AMP') 106 FORMAT(' O_SCALED '/' O_SCALE_HGT',12X,'FREQ',10X, 1 'DEL_TIM',10X,'IO_TIM',6X,'AMP') RETURN END SUBROUTINE MD2DOY(IYR,MONTH,IDAY,IDOY) DIMENSION IDIM(12) DATA IDIM/31,28,31,30,31,30,31,31,30,31,30,31/ IDIM(2)=28 IF(MOD(IYR,4).EQ.0) IDIM(2)=29 IDOY=0 DO 10 I=1,MONTH-1 IDOY=IDOY+IDIM(I) 10 CONTINUE IDOY=IDOY+IDAY RETURN END SUBROUTINE CKFREQ(N,APR,F,delt,riot,amp) REAL*8 F(50),APR(50),BADF(50),BADAPR(50) real*8 delt(50),riot(50),baddelt(50),badriot(50) integer amp(50),badamp(50) c REAL*8 TEMPF(50) C TYPE*,' N GOING INTO CKFREQ ',N C DO JG=1,N c WRITE(*,23) (APR(JG),F(JG),JG=1,N) C ENDDO 23 FORMAT(' ORIGINAL INPUT '/' AP_RANGE FREQ '/(2F15.5)) II=0 50 I=0 DO K=2,N IF(F(K).LE.F(K-1)) THEN I=I+1 II=II+1 BADF(II)=F(K) BADAPR(II)=APR(K) baddelt(ii)=delt(k) badriot(ii)=riot(k) badamp(ii)=amp(k) KK=K GOTO 100 ENDIF ENDDO 100 CONTINUE IF(I.EQ.0) GOTO 99 C REASSIGN VALUES AND RE-CHECK MM=0 DO M=1,KK-1 MM=MM+1 F(MM)=F(M) APR(MM)=APR(M) delt(mm)=delt(m) riot(mm)=riot(m) amp(mm)=amp(m) ENDDO DO M=KK+1,N MM=MM+1 F(MM)=F(M) APR(MM)=APR(M) delt(mm)=delt(m) riot(mm)=riot(m) amp(mm)=amp(m) ENDDO n=mm GOTO 50 99 CONTINUE IF(MM.NE.0) N=MM C DO KB=1,N c WRITE(*,32) (APR(KB),F(KB),KB=1,N) C ENDDO 32 FORMAT(' MODIFIED INPUT '/' AP_RANGE FREQ '/(2F15.5)) cc IF(II.NE.0) THEN cc WRITE(*,34) cc WRITE(*,33) (BADAPR(KK),BADF(KK),KK=1,II) cc ENDIF 33 format(2F12.5) 34 FORMAT(/' FREQ NOT INCREASING: REMOVE ENTRY(S) '/ 1 ' AP_RANGE FREQ') C TYPE*,' N LEAVING CKFREQ ',N RETURN END subroutine time_inc(jyr,jmo,jday,jhr,jmin,tsec, 1 iyr,imo,iday,ihr,imin,rsec,riotim) c routine to take the times of x-trace pts. scaled and adding to the c frame sync time and using these times instead of doing calculations c based on frequency and sweep rate. IMPLICIT REAL*8(A-H,O-Z) integer dmon(12) data dmon/31,28,31,30,31,30,31,31,30,31,30,31/ c type*,' enter yr,mon,day,hh,mm,ss, timeoffset' c accept*,iyr,imo,iday,ihr,imin,rsec,riotim c type*,' riotim ',riotim leapyr=0 if(mod(jyr,4).eq.0) then dmon(2)=29 leapyr=1 endif rsec=tsec+riotim if(rsec.ge.60) then rsec=rsec-60 imin=jmin+1 if(imin.ge.60) then imin=imin-60 ihr=jhr+1 if(ihr.ge.24) then ihr=ihr-24 iday=jday+1 if(iday.gt.dmon(jmo)) then iday=iday-dmon(jmo) imo=jmo+1 if(imo.gt.12) then iyr=jyr+1 imo=imo-12 endif endif endif endif endif c type*,iyr,imo,iday,ihr,imin,rsec return end subroutine get_file character*40 xscale,xtrace,oscale,otrace,neprof,plasf,nhout character*40 inname,neselalt,altselne,wmap character*200 dirname character*16 fid common/file/oscale common/com_icode/ iicode C TYPE*,' ENTER SCALED X TRACE FILENAME ' C READ(5,342) INNAME OPEN(UNIT=50,FILE='idl_tra.txt') READ(50,342) inname 342 format(a) read(50,342) dirname read(50,343) iicode 343 format(i1) write(*,*) iicode CLOSE(UNIT=50) fid=inname(3:18) write(xscale,123) fid 123 format('x_scal_',a,'.dat') write(xtrace,129) fid 129 format('x_fhp_',a,'.dat') write(oscale,124) fid 124 format('o_scal_',a,'.dat') write(otrace,125) fid 125 format('o_tra_',a,'.dat') write(neprof,126) fid 126 format('ne_prof_',a,'.dat') write(plasf,127) fid 127 format('plas_fr_',a,'.dat') write(nhout,128) fid 128 format('nhout_',a,'.dat') write(neselalt,191) fid 191 format('ne_sel_alt_',a,'.dat') write(altselne,192) fid 192 format('alt_sel_ne_',a,'.dat') write(wmap,193) fid 193 format('wmap_',a,'.dat') OPEN(UNIT=20,FILE=NESELALT,STATUS='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=21,FILE=ALTSELNE,STATUS='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=22,FILE=WMAP,STATUS='UNKNOWN',FORM='FORMATTED') C OPEN(UNIT=13,FILE=inname,STATUS='OLD',FORM='FORMATTED') OPEN(UNIT=13,FILE=dirname,STATUS='OLD',FORM='FORMATTED') OPEN(UNIT=14,FILE=xscale,STATUS='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=7,NAME=XTRACE,TYPE='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=8,NAME=OTRACE,TYPE='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=9,NAME=NEPROF,TYPE='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=10,FILE=PLASF,STATUS='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=16,NAME=NHOUT,TYPE='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=4,NAME='Ne1.DAT',TYPE='OLD',FORM='FORMATTED') OPEN(UNIT=3,NAME='NHINPUT.DAT',TYPE='OLD',FORM='FORMATTED') OPEN(UNIT=60,NAME='o_trace.dat',TYPE='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=61,NAME='o_scale.dat',TYPE='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=62,NAME='plas.dat',TYPE='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=63,NAME='ne.dat',TYPE='UNKNOWN',FORM='FORMATTED') OPEN(UNIT=64,NAME='wmap.dat',TYPE='UNKNOWN',FORM='FORMATTED') write(64,*)WMAP c output file for o_scaled will be opened in read_scaled if any c scaled o present c open(unit=15,file=oscale,status='new',form='formatted') c TYPE*,' ' c TYPE*,' O SCALED VALUES WILL BE IN FILE O_SCALED.DAT' c TYPE*,' O TRACED VALUES WILL BE IN FILE O_TRACE.DAT' c TYPE*,' X SCALED VALUES WILL BE IN FILE X_SCALED.DAT' c TYPE*,' X TRACED VALUES WILL BE IN FILE X_TRACE.DAT' c TYPE*,' PLASMA FREQ. VALUES WILL BE IN FILE PLAS_FREQ.DAT' c TYPE*,' NE PROFILE WILL BE IN FILE NE_PROFLIE.DAT' c TYPE*,' ORIGINAL OUTPUT WILL BE IN FILE NHOUTPUT.DAT' c TYPE*,' ' return end