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') ERR=0.01 L=0 RE=6371.2 PI=3.14159 DR=PI/180. 665 WRITE (*,666) ' 1 FOR ALOUETTE 1' WRITE(*,666) ' 2 FOR ALOUETTE 2' WRITE(*,666) ' 3 FOR ISIS 1' WRITE(*,666) ' 4 FOR ISIS 2' 666 FORMAT (1X,'Please enter the ',A) READ (*,10) ISAT C IF(ISAT.EQ.1) WRITE(75,*) ' ISIS 1' C IF(ISAT.EQ.2) WRITE(75,*) ' ISIS 2' IF (ISAT.GE.5 .OR. ISAT.EQ.0) GOTO 665 IF (ISAT.EQ.1) THEN ! ALOUETTE 1 OPEN(UNIT=32,FILE='AL1.',STATUS='OLD') C OPEN(UNIT=32,FILE='AL1B.',STATUS='OLD') ELSEIF (ISAT.EQ.2) THEN ! ALOUETTE 2 OPEN(UNIT=32,FILE='AL2.',STATUS='OLD') ELSEIF (ISAT.EQ.3) THEN ! ISIS 1 ISAT = 5 OPEN (UNIT=32,FILE='ISIS1.',STATUS='OLD') OPEN (UNIT=33,FILE='ISIS1N.',STATUS='OLD') ELSEIF (ISAT.EQ.4) THEN ! ISIS 2 ISAT = 6 OPEN (UNIT=32,FILE='ISIS2.',STATUS='OLD') OPEN (UNIT=33,FILE='ISIS2N.',STATUS='OLD') 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) WRITE (*,666) 'year LAST 2 DIGITS (ZERO TO END)' READ(*,10,END=200) LYEAR IF(LYEAR.EQ.00) GOTO 999 WRITE (*,666) 'day' READ(*,10,END=200) LDAY 10 FORMAT(I3) WRITE (*,666) 'hour' READ(*,10) LHOUR WRITE (*,666) 'min' READ(*,10) MIN WRITE (*,666) 'sec' READ(*,10) LSEC 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) WRITE(*,20)LYEAR,LDAY,LHOUR,MIN,LSEC,JTIME,DLAT,DLONG,HGT, 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 20 FORMAT(I3,I4,I3,I2,1H/,I2,I5,F7.2,F8.2,F6.0,I5,F7.2,F8.2,F7.2,F6.3 C,2I4,1X,A3,F8.2) c f5.2 GO TO 1 999 CONTINUE 200 CLOSE(UNIT=32) CLOSE(UNIT=33) STOP 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 FIELDG(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 FIELDG (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 FIELDG(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 A1(6),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