//VISIBLE4 JOB NOTIFY=1212003 /*JOBPARM S=59,L=99,R=5012 /*SERVICE -4 // EXEC FORTVCLG,RC=5012K,RG=1216K //FORT.SYSIN DD * C ----------------------------------------------------------------------00002 C 00001 C 00003 C 00005 C NAME : ALERT 00004 C TYPE : MAIN PROGRAMME 00006 C 00007 C PURPOSE : DETERMINES THE # AND ID OF ALL VISIBLE SATELLITES 00008 C ALONG WITH THE FOUR SATELLITES WHICH YIELD THE 00009 C THE BEST NAV PERFORMANCE AND GDOP VALUES. 00010 C C COMPUTES ALERT FOR DIFFERENT SATELLITES AND C DIFFERENT INPUTS C 00011 C AUTHOR : JUNE 1981 S. MERTIKAS 00012 C 00013 C MODIFIED : FEBR 1983 D. DELIKARAOGLOU 00014 C 00015 C MODIFIED : JUNE 1983 R. SANTERRE C IMPLICIT REAL*8(A-H,O-Z) 00016 LOGICAL*1 ISTAR/'*'/ , IBLANK/'-'/ LOGICAL*1 ICRL(100) C 00017 C THE FOLLOWING ARRAYS HAVE ROW DIMENSION 'MXVS' EQUAL 00018 C TO THE MAXIMUM NUMBER OF ANTICIPATED VISIBLE SATELLITES 00019 C 00020 INTEGER*4 AZBUF(9),ELEBUF(9) ,IVIS(9),IHRVIS(18,100),IHH(25) 00021**2 REAL*4 ELV,AZIM DIMENSION VISAT(9,3) , ELV(9) , AZIM(9) 00022**2 C 00023 C 00024 DIMENSION DTM(7) , REC(16) , EPH(18) , CLK(6) , 00025 # EPHBUF(324) , CLKBUF(108) , IDSAT(18) , Y(5) , 00026 # XSVBUF(18,3) , XRCV(3) , XSV(3) 00027 DIMENSION IBEST(4) , KK(4) , A(9,9) , 00028 # B(9,9) , IAZ(4) , IELEV(4) 00029 DIMENSION IY(18) , MDAY(18) , JSAT(18) , PERIOD(18), # WMOTN(18) , TNODE(18), XLNODE(18) DIMENSION TP(18) , XN(18) , W(18) , WDOT(18) , # EC(18) , AO(18) , OM(18) , OMDOT(18) , # COSI(18) , GAST(18) , SINI(18) DIMENSION XINC(18) , RA(18) , AP(18) , XMA(18) , # PN(18) , WMTN(18) , DAY(18) DIMENSION XLANW(18) , KYAN(18) , KDAYAN(18), KHAN(18) , # XMAN(18) C 00030 C ----------------------------------------------------------------- 00031 C DEFINE I/0 LOGICAL UNITS 00032 C 00033 C ICR - CARD READER 00034**3 C IPR - LINE PRINTER DETAILED OUTPUT 00035**3 C IPR1 - LINE PRINTER SUMMARY OUTPUT 00036**3 C IPR2 IPR3 - OUTPUT FILE UNIT FOR PLOTTING 00037**3 C LEPH - INPUT FILE UNIT FOR SATELLITE EPHEMERIS 00038**3 C LU - OUTPUT FILE UNIT FOR OPTIMUM TRACKING 00039**3 C SCHEDULE IF REQUIRED 00040**3 C 00041**3 C IF ANY OF THE I/O UNITS IS NOT USED SPECIFY IN YOUR JCL CARDS 00042**3 C //GO.FT##F001 DD DUMMY FOR THIS UNIT 00043**3 C 00044**3 C 00045 C ----------------------------------------------------------------- 00046 C 00047 DATA ICR / 5/ , IPR / 6/ , IPR1/16/ , 00048**3 # IPR2/17/ , LEPH/18/ , LU /25/ ,IPR3 /19/ 00049**3 C --------------------------------------------------------------------- C KIND OF INPUT AND KIND OF OUTPUT C C IPARA1=1 : GPS SATELLITES C IPARA1=2 : TRANSIT SATELLITES C IPARA1=3 : LAGEOS SATELLITE(S) C IPARA1=4 : STARLETTE SATELLITE(S) C C IPARA2=1 : NNSS NODAL CROSSING DATA C IPARA2=2 : KEPLERIAN PARAMETERS C IPARA2=3 : GPS THEORETICAL OR TRACKED EPHEMERIS C IPARA2=4 : NASA PREDICTION BULLETIN (GAST USED) C IPARA2=5 : NASA PREDICTION BULLETIN (GAST UNUSED) C C IPARA3=1 : SIMPLE APPROACH (ELEV,AZ) C IPARA3=2 : ELEV,AZ AND DIL. OF PRECISION FOR THE BEST 4 SAT C IPARA3=3 : ELEV,AZ AND DIL. OF PRECISION FOR ALL SATELLITES C --------------------------------------------------------------------- READ(ICR,134) IPARA1,IPARA2,IPARA3 C 00050 C ----------------------------------------------------------------- 00051 C DEFINITION OF SOME BASIC CONSTANTS 00052 C ----------------------------------------------------------------- 00053 C 00054 IGAP=0 PI=DARCOS(-1.D0) 00055 DTR=PI/180.D0 00056 C=299792458.D0 00057 RMEAN=6371.D3 00058 F1=1575.D6 00059 F2=1227.D6 00060 C 00062 ICD=5 00063 SOBS= 4.D0 00064 C 00065 C MXVS - MAXIMUM NUMBER OF VISIBLE SATELLITES 00066 C (SEE NOTE IN DIMENSION STATEMENT ABOVE) 00067 C 00068 MAXSAT = 18 00069**2 MXVS = 9 00070**2 IDA = 9 00071**2 MREC = 0 00072**4 C 00073 C-------------------------------------------------------------------- 00074 C DEFINE REFERENCE ELLIPSOID 00075 C-------------------------------------------------------------------- 00076 C 00077 CALL DATUM(DTM,0,ICR) 00078 BE=DSQRT((1.D0-DTM(2)**2)*DTM(1)**2) 00079 C 00080 C--------------------------------------------------------------------- 00081 C DEFINE LOCATION AND TIME INTERVAL FOR ALERT 00082 C 00083 C LATITUDE IN DEG , MIN , SEC 00084 C LONGITUDE IN DEG , MIN , SEC (+VE EAST) 00085 C HEIGHT IN METRES 00086 C BEGINNING TIME - YR, DAY, HR, MIN, SEC 00087 C END TIME - YR, DAY, HR, MIN, SEC 00088 C TIME INCREMENT (FOR OUTPUT) IN MINUTES 00089 C 00090 C--------------------------------------------------------------------- 00091 C 00092 READ(ICR,1) LATD,LATM,LATS,LOND,LONM,LONS,RHT,ELEVMN 00093 KRHT=IDINT(RHT) KELEV=IDINT(ELEVMN) WRITE(IPR,2) 00094 GO TO (310,320,330,340,350) , IPARA2 310 WRITE(IPR,311) GO TO 360 320 WRITE(IPR,321) GO TO 360 330 WRITE(IPR,331) GO TO 360 340 WRITE(IPR,341) GO TO 360 350 WRITE(IPR,351) 360 CONTINUE GO TO(155,160,165,170) , IPARA1 155 WRITE(IPR,156) GO TO 200 160 WRITE(IPR,161) GO TO 200 165 WRITE(IPR,166) GO TO 200 170 WRITE(IPR,171) 200 CONTINUE WRITE(IPR,8000) KELEV WRITE(IPR,3) LATD,LATM,LATS,LOND,LONM,LONS 00095 WRITE(IPR,7000) KRHT C 00096 C--------------------------------------------------------------------- 00097 C COMPUTE THE Z-COUNT FOR THE : 00098 C BEGINNING DATE (ZCNT1) 00099 C ENDING DATE (ZCNT2) 00100 C--------------------------------------------------------------------- 00101 C 00102 READ(ICR,4) KYEAR,KDAY,KHR,KMIN,KSEC 00103 KKHR=KHR KKYEAR=KYEAR KKDAY=KDAY KKMIN=KMIN KKSEC=KSEC READ(ICR,4) IYEAR,IDAY,IHR,IMIN,ISEC 00104 READ(ICR,7) IDT 00105 IF(IPARA3 .EQ. 1 .AND. IPARA1 .NE.1) GO TO 135 ZCNT1 = ZCOUNT(KYEAR,KDAY,KHR,KMIN,KSEC) 00106 ZZZ=ZCNT1 ZCNT2 = ZCOUNT(IYEAR,IDAY,IHR,IMIN,ISEC) 00107 WRITE(IPR,5) KDAY,KYEAR,KHR,KMIN,KSEC,ZCNT1 00108**2 WRITE(IPR,6) IDAY,IYEAR,IHR,IMIN,ISEC,ZCNT2 00109**2 GO TO 136 135 WRITE(IPR,35) KDAY,KYEAR,KHR,KMIN,KSEC WRITE(IPR,36) IDAY,IYEAR,IHR,IMIN,ISEC 136 CONTINUE C 00110 C_____________________________________________________________________ 00111 C 00112 C HOW MANY SATELLITES YOU HAVE AVAILABLE? 00113 C 00114 C IF SOME ARE AVAILABLE: 00115 C A) DEFINE NSAT=? 00116 C B) DEFINE IDSAT=? 00117 C IF ALL ARE AVAILABLE: 00118 C A) NSAT=18 00119 C B) IDSAT(I)=I ,I=1,18 00120 C C) EPHBUF SHOULD BE FULL 00121 C--------------------------------------------------------------------- 00122 C 00123 READ(ICR,8) NSAT ,(IDSAT(I),I=1,NSAT) 00124 IF (NSAT .NE. MAXSAT) GO TO 9 00125 DO 10 I=1,MAXSAT 00126 IDSAT(I)=I 00127 10 CONTINUE 00128 9 MAXSAT=NSAT 00129 C 00130 C----------------------------------------------------------------- 00131 C READ EPHEMERIS PARAMETERS FOR ALL AVAILABLE SATELLITES 00132 C------------------------------------------------------------------- 00133 C GO TO(43,44,210,131,132) , IPARA2 C 00134 210 CALL READEF(REC,DTM,EPH,EPHBUF,CLK,CLKBUF,MAXSAT,ICR,IPR) **6 C 00136 C 00137 DO 11 I=1,MAXSAT 00138 ISAT=IDSAT(I) 00139 IPSN=IDSAT(I)*6-5 00140 CALL DASET(CLKBUF(IPSN),CLK(1),6) 00141 IPSN=IDSAT(I)*18-17 00142 CALL DASET(EPHBUF(IPSN),EPH(1),18) 00143 11 CONTINUE 00144 GO TO 133 C 43 CALL EPNCD(NSAT,IY,MDAY,JSAT,PERIOD,WMOTN,TNODE,XLNODE,ICR) GO TO 133 C 44 CALL EPKE(NSAT,IY,MDAY,JSAT,TP,XN,W,WDOT,EC,AO,OM,OMDOT, # COSI,GAST,SINI,ICR) GO TO 133 C 131 CALL EPNPBG(NSAT,IDSAT,IY,DAY,XINC,RA,EC,AP,XMA,XN,PN,WMTN,ICR) GO TO 133 C 132 CALL EPNPB(NSAT,IDSAT,IY,DAY,XINC,EC,AP,XMA,XN,PN,WMTN,XLANW, # KYAN,KDAYAN,KHAN,XMAN,ICR) C 133 CONTINUE C IF(IPARA3 .EQ. 1 .AND. IPARA1 .EQ. 1) GO TO 600 IF(IPARA3 .EQ. 1) GO TO 137 IF(IPARA3 .EQ. 2) WRITE(IPR,12) 00145 IF(IPARA3 .EQ. 3) WRITE(IPR,820) WRITE(IPR1,94) 00146**3 GO TO 138 600 WRITE(IPR,601) GO TO 138 137 WRITE(IPR,37) 138 CONTINUE C 00147 C--------------------------------------------------------------------- 00148 C CONVERSION OF RECEIVER LONGITUDE AND LATITUDE 00149 C EXPRESSED IN (DEG, MIN, SEC) INTO RADIANS 00150 C--------------------------------------------------------------------- 00151 C 00152 XLAT= DTR* (LATD+(DFLOAT(LATM)+(DFLOAT(LATS)/60.))/60.) 00153 XLON= DTR* (LOND+(DFLOAT(LONM)+(DFLOAT(LONS)/60.))/60.) 00154 C 00155 C--------------------------------------------------------------------- 00156 C CONVERT RECEIVER GEODETIC CORDINATES (RPHI,RLAT,RHT) 00157 C INTO CARTESIAN ONES (XR ,YR ,ZR ) 00158 C--------------------------------------------------------------------- 00159 C 00160 CALL PLHXYZ(XLAT,XLON,RHT,DTM(3),DTM(4),DTM(5),DTM(1), 00161 # BE,XR,YR,ZR) 00162 XRCV(1)=XR 00163 XRCV(2)=YR 00164 XRCV(3)=ZR 00165 C 00166 C--------------------------------------------------------------------- 00167 C COMPUTE SATELLITE COORDINATES 00168 C--------------------------------------------------------------------- 00169 C 00170 IPP=0 14 J = 0 00171 IP=MOD(KMIN,15) IF(IP .EQ. 0) IPP=IPP+1 IF(IPARA1 .NE. 1) GO TO 139 ZCNT1 = ZCOUNT(KYEAR,KDAY,KHR,KMIN,KSEC) 00172 139 CONTINUE DO 15 K=1,10 00173 IVIS(K) = 0 00174 ELEBUF(K) = 0 00175 AZBUF(K) = 0 00176 ELV(K)=0. AZIM(K)=0. IF(IPARA3 .EQ. 1) GO TO 15 VISAT(K,1)= 0 00177 VISAT(K,2)= 0 00178 VISAT(K,3)= 0 00179 15 CONTINUE 00180 IF(IPARA3 .EQ. 1) GO TO 140 ELVMAX=0.D0 00181 140 CONTINUE DO 17 I=1,NSAT 00182 GO TO(141,142,220,143,144) , IPARA2 220 IPSN=IDSAT(I)*18-17 00183 CALL DASET(EPHBUF(IPSN),EPH(1),18) 00184 Y(1)=EPH(8) 00185 Y(2)=EPH(6) 00186 Y(3)=EPH(11) 00187 Y(4)=EPH(13) 00188 Y(5)=EPH(15) 00189 TIME=ZCNT1*6.-EPH(9) 00190 CALL STXYZ(AMYE,C,DTM,EPH,ICD,IFLAG4,IPR,PHIK,PI,RING, 00191 # RLNG,TIME,XR,YR,ZR,XS,YS,ZS,Y,IPC,IPR) 00192 C GO TO 145 C 141 CALL NNSSNC(IY(I),MDAY(I),TNODE(I),XLNODE(I),PERIOD(I), # WMOTN(I),KYEAR,KDAY,KHR,KMIN,KSEC,PI,XS,YS,ZS) GO TO 145 C 142 CALL NNSSKE(IY(I),MDAY(I),TP(I),XN(I),W(I),WDOT(I),EC(I),AO(I), # OM(I),OMDOT(I),COSI(I),GAST(I),SINI(I),KYEAR,KDAY, # KHR,KMIN,KSEC,DTM(7),PI,XS,YS,ZS,IPR) GO TO 145 C 143 CALL NPBST(IY(I),DAY(I),XINC(I),RA(I),EC(I),AP(I),XMA(I),XN(I), # PN(I),WMTN(I),KYEAR,KDAY,KHR, # KMIN,KSEC,PI,DTM(7),DTM(6),XS,YS,ZS,IPR) GO TO 145 C 144 CALL NPB(IY(I),DAY(I),XINC(I),EC(I),AP(I),XMA(I),XN(I),PN(I), # WMTN(I),XLANW(I),KYAN(I),KDAYAN(I),KHAN(I),XMAN(I),KYEAR, # KDAY,KHR,KMIN,KSEC,PI,DTM(7),DTM(6),XS,YS,ZS,IPR) C 145 CONTINUE C XSV(1)=XS 00193 XSV(2)=YS 00194 XSV(3)=ZS 00195 C 00196 IF(IPARA3 .EQ. 1) GO TO 146 XSVBUF(I,1) = XSV(1) 00197 XSVBUF(I,2) = XSV(2) 00198 XSVBUF(I,3) = XSV(3) 00199 146 CONTINUE C 00200 C--------------------------------------------------------------------- 00201 C CALCULATE SATELLITE ELEVATION ANGLE(ELEV) AND 00202 C AZIMUTH(AZ) IN DEGREES 00203 C------------------------------------------------------------------- 00204 C 00205 CALL DERIV(IPR,DTM(1),BE,XLAT/DTR,XLON/DTR,RHT,XSV,RNG,DSDP, 00206 # DSDL,ELEV,AZ,IER) 00207 C 00208 C---------------------------------------------------------------------- 00209 C DETERMINE NUMBER OF VISIBLE SATELLITES 00210 C--------------------------------------------------------------------- 00211 C 00212 IF(ELEV.GT.ELEVMN) GO TO 20 00213 GO TO 17 00214 C 00215 C-------------------------------------------------------------------- 00216 C IDENTIFY SATELLITE WITH MAXIMUM ELEVATION 00217 C COMPUTE DIRECTION COSINES 00218 C OR COMPILE SATELLITES , ELEVATIONS AND AZIMUTHS ONLY C-------------------------------------------------------------------- 00219 C 00220 20 J = J + 1 00221 IF(IPARA3 .NE. 1) GO TO 148 IF (IPARA3 .EQ. 1 .AND. IPARA1 .EQ. 1) GO TO 300 IF(J .GT. 3) WRITE(IPR,150) ELEBUF(J)=ELEV AZBUF(J)=AZ IVIS(J)=IDSAT(I) GO TO 17 300 IF(J .GT. MXVS) J=J-1 IF(J .GT. MXVS) GO TO 17 ELEBUF(J)=ELEV AZBUF(J)=AZ IVIS(J)=IDSAT(I) ELV(J)=ELEV AZIM(J)=AZ IF(IP .EQ. 0) IHRVIS(I,IPP)=1 GO TO 17 148 CONTINUE IF(J .GT. MXVS) WRITE(IPR,1040) J,MXVS,J 00222**4 IF(J .GT. MXVS) J = J - 1 00223 IF(J .GT. MXVS) GO TO 17 00224 ELEBUF(J) = ELEV 00225 AZBUF(J) = AZ 00226 IVIS(J)=IDSAT(I) 00227 IF(IP .EQ. 0) IHRVIS(I,IPP)=1 IF(ELEV.GT.ELVMAX) GO TO 30 00228 GO TO 32 00229 30 ELVMAX=ELEV 00230 KSAT=J 00231 C 32 VISAT(J,1)= ( XSV(1)-XRCV(1) )/RNG 00232 C VISAT(J,2)= ( XSV(2)-XRCV(2) )/RNG 00233 C VISAT(J,3)= ( XSV(3)-XRCV(3) )/RNG 00234 32 CONTINUE C C DIRECTION COSINES IN THE TANGENTIAL SYSTEM (NORD,EST,HEIGHT) C MODIFICATION 23 MAI 1985 R.SANTERRE. C VISAT(J,1) = DCOS(AZ*DTR) * DCOS(ELEV*DTR) VISAT(J,2) = DSIN(AZ*DTR) * DCOS(ELEV*DTR) VISAT(J,3) = DSIN(ELEV*DTR) C 17 CONTINUE 00235 IF(IPARA3 .NE. 1) GO TO 152 IF(J .EQ. 0) IGAP=IGAP+1 IF(J .NE. 0) IGAP=0 IF(J .EQ. 0 .AND. IGAP .EQ. 1 .AND. IPARA1 .EQ. 1)WRITE(IPR,901) IF(J .EQ. 0 .AND. IGAP .EQ. 1 .AND. IPARA1 .NE. 1)WRITE(IPR,902) IF (J .EQ. 0) GO TO 153 IF(IPARA3 .EQ. 1 .AND. IPARA1 .EQ. 1) GO TO 400 IF(J .GT. 3) GO TO 153 GO TO (471,472,473) ,J 471 WRITE(IPR,481) KHR,KMIN,(IVIS(N),N=1,J),(ELEBUF(N),N=1,J), # (AZBUF(N),N=1,J) GO TO 153 472 WRITE(IPR,482) KHR,KMIN,(IVIS(N),N=1,J),(ELEBUF(N),N=1,J), # (AZBUF(N),N=1,J) GO TO 153 473 WRITE(IPR,483) KHR,KMIN,(IVIS(N),N=1,J),(ELEBUF(N),N=1,J), # (AZBUF(N),N=1,J) GO TO 153 400 IZCNT=IDINT(ZCNT1) WRITE(IPR3,1117) KHR,KMIN,(IVIS(I),I=1,9),(ELV(I),I=1,9), # (AZIM(I),I=1,9) GO TO(451,452,453,454,455,456,457,458,459) ,J 451 WRITE(IPR,461) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J) GO TO 153 452 WRITE(IPR,462) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J) GO TO153 453 WRITE(IPR,463) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J) GO TO153 454 WRITE(IPR,464) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J) GO TO153 455 WRITE(IPR,465) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J) GO TO153 456 WRITE(IPR,466) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J) GO TO 153 457 WRITE(IPR,467) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J) GO TO153 458 WRITE(IPR,468) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J) GO TO 153 459 WRITE(IPR,469) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J) GO TO 153 152 CONTINUE IF(J .EQ. 0) IGAP=IGAP+1 IF(J .NE. 0) IGAP=0 IF(J .EQ. 0 .AND. IGAP .EQ.1) WRITE(IPR,903) IF(J .EQ. 0) GO TO 25 00236**3 WRITE(14,1) J **6 DO 13 L = 1,NSAT **6 IS = IDSAT(L) **6 DO 18 K = 1,J **6 IF(IS .EQ. IVIS(K)) I1 = L **6 IF(IS .EQ. IVIS(K)) GO TO 16 **6 18 CONTINUE **6 GO TO 13 **6 16 WRITE(14,1110) IVIS(K),(XSVBUF(I1,M),M=1,3) **6 1110 FORMAT(I4,3F14.4) **6 13 CONTINUE **6 C 00237 C 00238 IF(IPARA3 .EQ. 3) GO TO 24 C-------------------------------------------------------------------- 00239 C IF IPARA3 EQUAL 2 C C 1. CALCULATION OF THE VOLUME OF THE TETRAHEDRON FORMED BY THE 00240 C FOUR UNIT VECTORS TO EACH SATELLITE FOR ALL COMBINATIONS 00241 C 2. SELECTION OF THE MAXIMUM VOLUME(VOL) 00242 C 3. SELECTION OF THOSE VISIBLE SATELLITES WITH BEST NAVIGATION 00243 C PERFORMANCE (MAX VOLUME) 00244 C-------------------------------------------------------------------- 00245 C 00246 CALL VOLUME(KSAT,VISAT,MXVS,J,KK,VOL) 00247 DO 45 L=1 ,4 00248 IBEST(L) = IVIS( KK(L) ) 00249 IELEV(L) = ELEBUF( KK(L) ) 00250 IAZ(L) = AZBUF( KK(L) ) 00251 45 CONTINUE 00252 IF(J .GE. 4) NVIS = 4 00253**2 IF(J .LT. 4) NVIS = J 00254**3 GO TO (25,21,22,23) , NVIS 00255**3 GO TO 25 00256**2 C 00257**2 C ----------------------------------------------------------------------00258**2 C FOR 3 VISIBLE SATELLITES COMPUTE 00259**2 C 00260**3 C PDOP = SQRT( TRACE COV(X , Y , Z)) 00261**3 C HTDOP = SQRT( TRACE COV(X , Y , CLOCK) ) 00262**3 C HDOP = SQRT( TRACE COV(X , Y)) 00263**3 C 00264**3 C ----------------------------------------------------------------------00265**2 C 00266**2 21 CALL DOP2S(A,B,IDA,VISAT,MXVS,C,HDOP) 00267**3 IZCNT = IDINT(ZCNT1) 00268**3 TOFDAY = ZFRDAY( ZCNT1 ) * 24.D0 00269**3 MREC = MREC + 1 00270**4 WRITE(IPR1,92) MREC,KHR,KMIN,IZCNT,(IVIS(N),N=1,NVIS), 00271**4 # (ELEBUF(N),N=1,NVIS),(AZBUF(N),N=1,NVIS), 00272**3 # HDOP 00273**3 WRITE(IPR ,93) KHR,KMIN,IZCNT,(IVIS(N),N=1,NVIS), 00274**3 # (ELEBUF(N),N=1,NVIS),(AZBUF(N),N=1,NVIS), 00275**3 # HDOP 00276**3 WRITE(IPR2,97) TOFDAY,HDOP,PDOP,HTDOP,GDOP 00277**6 GO TO 25 00278**3 C 00279**3 22 CALL DOP3S(A,B,IDA,VISAT,MXVS,C,HTDOP,PDOP,HDOP) 00280**3 IZCNT = IDINT(ZCNT1) 00281**2 MREC = MREC + 1 00282**4 WRITE(IPR,90) KHR,KMIN,IZCNT,(IVIS(N),N=1,NVIS), 00283**2 # (ELEBUF(N),N=1,NVIS),(AZBUF(N),N=1,NVIS),PDOP,HTDOP 00284**2 TOFDAY = ZFRDAY( ZCNT1 ) * 24.D0 00285**3 WRITE(IPR1,95) MREC,KHR,KMIN,IZCNT,(IVIS(N),N=1,NVIS), 00286**4 # (ELEBUF(N),N=1,NVIS),(AZBUF(N),N=1,NVIS), 00287**3 # HDOP,PDOP,HTDOP 00288**3 WRITE(IPR2,97) TOFDAY,HDOP,PDOP,HTDOP,GDOP 00289**6 GO TO 25 00290**2 C 00291**2 C ----------------------------------------------------------------------00292**2 C 00293**3 C FOR THE FOUR SATELLITES WHICH GIVE THE MAXIMUM 00294**3 C VOLUME OF THE THE FUNDAMENTAL TETRAHEDRON COMPUTE 00295**3 C 00296**3 C GDOP = SQRT( TRACE COV(X , Y , Z , CLOCK) ) 00297**3 C PDOP = SQRT( TRACE COV(X , Y , Z)) 00298**3 C HTDOP = SQRT( TRACE COV(X , Y , CLOCK) ) 00299**3 C HDOP = SQRT( TRACE COV(X , Y)) 00300**3 C 00301**3 C ----------------------------------------------------------------------00302**2 C 00303**2 23 CALL DOP4S(A,B,IDA,VISAT,MXVS,C,KK,GDOP,HTDOP,PDOP,HDOP) 00304**3 MREC = MREC + 1 00305**4 IZCNT = IDINT(ZCNT1) 00306**2 M=J-3 GO TO(514,515,516,517,518,519) ,M 514 WRITE(IPR,524) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # (IBEST(N),N=1,4),GDOP GO TO 26 515 WRITE(IPR,525) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # (IBEST(N),N=1,4),GDOP GO TO 26 516 WRITE(IPR,526) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # (IBEST(N),N=1,4),GDOP GO TO 26 517 WRITE(IPR,527) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # (IBEST(N),N=1,4),GDOP GO TO 26 518 WRITE(IPR,528) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # (IBEST(N),N=1,4),GDOP GO TO 26 519 WRITE(IPR,529) KHR,KMIN,IZCNT,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # (IBEST(N),N=1,4),GDOP 26 CONTINUE TOFDAY = ZFRDAY(ZCNT1) * 24.D0 00310**3 WRITE(IPR1,96) MREC,KHR,KMIN,IZCNT,(IVIS(N),N=1,NVIS), 00311**4 # (ELEBUF(N),N=1,NVIS),(AZBUF(N),N=1,NVIS), 00312**3 # HDOP,PDOP,HTDOP,GDOP 00313**3 WRITE(IPR2,97) TOFDAY,HDOP,PDOP,HTDOP,GDOP 00314**6 WRITE(LU,88) ZCNT1,(IVIS(I),I=1,MXVS),(IBEST(J),J=1,4) **6 GO TO 25 C---------------------------------------------------------------------- C C DILUTION OF PRECISION FOR ALL SATELLITE C C---------------------------------------------------------------------- 24 CALL DOPS(A,B,IDA,VISAT,MXVS,J,C,GDOP,PDOP,HTDOP,HDOP,VDOP) GO TO(801,802,803,804,805,806,807,808,809) ,J 801 WRITE(IPR,811) KHR,KMIN,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # VDOP GO TO 25 802 WRITE(IPR,812) KHR,KMIN,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # HDOP,VDOP GO TO 25 803 WRITE(IPR,813) KHR,KMIN,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # PDOP,HTDOP,HDOP,VDOP GO TO 25 804 WRITE(IPR,814) KHR,KMIN,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # GDOP,PDOP,HTDOP,HDOP,VDOP GO TO 25 805 WRITE(IPR,815) KHR,KMIN,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # GDOP,PDOP,HTDOP,HDOP,VDOP GO TO 25 806 WRITE(IPR,816) KHR,KMIN,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # GDOP,PDOP,HTDOP,HDOP,VDOP GO TO 25 807 WRITE(IPR,817) KHR,KMIN,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # GDOP,PDOP,HTDOP,HDOP,VDOP GO TO 25 808 WRITE(IPR,818) KHR,KMIN,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # GDOP,PDOP,HTDOP,HDOP,VDOP GO TO 25 809 WRITE(IPR,819) KHR,KMIN,(IVIS(N),N=1,J), # (ELEBUF(N),N=1,J),(AZBUF(N),N=1,J), # GDOP,PDOP,HTDOP,HDOP,VDOP 25 CONTINUE 00316**2 GDOP = 0. 00317**3 PDOP = 0. 00318**3 HTDOP = 0. 00319**3 HDOP = 0. 00320**3 VDOP = 0. 153 CONTINUE C 00321**2 C 00322 C------------------------------------------------------------------- 00323 C ITERATION FOR THE SPECIFIED TIME INCREMENT OF (IDT) MINUTES 00324 C------------------------------------------------------------------- 00325 C 00326 IF(KYEAR.GE.IYEAR.AND.KDAY.GE.IDAY.AND.KHR.GE.IHR 00327 # .AND.KMIN.GE.IMIN) GO TO 2221 00328 KMIN = KMIN + IDT 00329 IF(KMIN.LT.60) GO TO 14 00330 KMIN = KMIN - 60 00331 KHR = KHR + 1 00332 IF(KHR.LT.24) GO TO 14 00333 KDAY = KDAY + 1 00334 KHR = KHR - 24 00335 GO TO 14 00336 C----------------------------------------------------------------------- C VISIBILITY SUMMARY FOR GPS SATELLITES C C LIMITATIONS: START TIME MUST BEGIN AT A INTEGER HOUR EX: 12 00 00 C INCREMENT TIME MUST BE (1 OR 3 OR 5 OR 15 OR 30 OR 60)MIN C----------------------------------------------------------------------- 2221 IF(IPARA1 .NE. 1) GO TO 2222 WRITE(IPR,5100) KELEV WRITE(IPR,5110)LATD,LATM,LATS,LOND,LONM,LONS WRITE(IPR,5111) KRHT WRITE(IPR,5120)KKDAY,KKYEAR,KKHR,KKMIN,KKSEC,ZZZ WRITE(IPR,5130)IDAY,IYEAR,IHR,IMIN,ISEC,ZCNT2 DO 5000 J=1,NSAT WRITE(IPR,5002) IDSAT(J) DO 5001 K=1,IPP ICRL(K)=IBLANK IF(IHRVIS(J,K) .EQ. 1) ICRL(K)=ISTAR 5001 CONTINUE WRITE(IPR,5601) (ICRL(K),K=1,IPP) 5601 FORMAT('+',12X,100A1) 5000 CONTINUE NFH=INT(IPP/4)+1 DO 5900 L=1,NFH KKKR=KKHR+L-1 IF(IDT .EQ. 30) KKKR=KKKR*2-KKHR IF(IDT .EQ. 60) KKKR=KKKR*4-KKHR*3 IHH(L)=KKKR IF(IHH(L) .GE. 24) IHH(L)=IHH(L)-24 5900 CONTINUE WRITE(IPR,5600) (IHH(II),II=1,NFH) 5100 FORMAT('1',7X,'VISIBILITY (ELEVATION > ',I2,' DEGREES) OF GPS', # ' SATELLITES',///) 5002 FORMAT('0',5X,'GPS',I3) 5600 FORMAT(10X,25I4) 5110 FORMAT(9X,'LAT=',I3,'DEG',2X,I3,'MIN',2X,I3,'SEC', 0340 # /,(7X,' LONG=',I3,'DEG',2X,I3,'MIN',2X,I3,'SEC')) 0341 5111 FORMAT(11X,'H=',I3,' METRES',/) 5120 FORMAT(8X,'BEGIN ALERT: DAY ',I4,' , ',I5,4X, 0343**2 # 2(I2,':'),I2,' ZCNT1 = ',F7.0,/) 00344**2 5130 FORMAT(8X,' END ALERT: DAY ',I4,' , ',I5,4X, 0345**2 # 2(I2,':'),I2,' ZCNT2 = ',F7.0,/////) 00346**2 C GO TO 2222 C 1 FORMAT(6I4,F6.2,F3.0) 00337 2 FORMAT('1',44X,'YOUR LOCATION AND DESIRED TIME INTERVAL FOR ALERT 00338 # ',/) 00339 3 FORMAT(25X,'LAT=',I3,'DEG',2X,I3,'MIN',2X,I3,'SEC', 00340 # /,(23X,' LONG=',I3,'DEG',2X,I3,'MIN',2X,I3,'SEC')) 00341 7000 FORMAT(27X,'H=',I3,' METRES',/) 4 FORMAT(5I4) 00342 5 FORMAT(14X,'BEGIN ALERT: DAY ',I4,' , ',I5,4X, 00343**2 # 2(I2,':'),I2,' ZCNT1 = ',F7.0,/) 00344**2 6 FORMAT(14X,' END ALERT: DAY ',I4,' , ',I5,4X, 00345**2 # 2(I2,':'),I2,' ZCNT2 = ',F7.0,///) 00346**2 7 FORMAT(I2) 00347 8 FORMAT(19I2) 00348 12 FORMAT(1X,'HR MIN ZCNT',7X,' VISIBLE SATELLITES', 00349**2 # 12X,'ELEVATIONS',22X,'AZIMUTHS',27X,'SELECT',4X,'GDOP', 00350**2 @ /,117X,'(PDOP)',2X,'(HTDOP)',/,117X,'(HDOP)',//) 00351**2 88 FORMAT(F8.0,13I3) **6 90 FORMAT(2I3,I6,2(1X,3I3,18X),3I5,34X,F5.1,3X,F5.1/) 00354**2 92 FORMAT(10X,I4,2I3,I6,2(1X,2I3,8X),2I5,13X,F6.1,18X,/) 00355**5 93 FORMAT(2I3,I6,2(1X,2I3,21X),2I5,39X,F5.1,/) 00356**3 94 FORMAT(15X,'HR MIN ZCNT',3X,'SATELLITES',5X, 00357**5 # 'ELEVATIONS',6X,'AZIMUTHS',11X,' HDOP PDOP HTDOP GDOP',//) 00358**5 95 FORMAT(10X,I4,2I3,I6,2(1X,3I3,5X),3I5,8X,3F6.1,6X,/) 00359**5 96 FORMAT(10X,I4,2I3,I6,2(1X,4I3,2X),4I5,3X,4F6.1,/) 00360**5 97 FORMAT(F10.5,4F7.1) **6 130 FORMAT(5X,G18.9,G18.9) 00362 1000 FORMAT(5X,9G14.8) 00363 1030 FORMAT(//,5X,'SAT #',I3,6(/,5X,'CLK(',I1,')=',E20.10), 00364 # /,18(/,5X,'EPH(',I2,')=',E20.10)) 00365 1040 FORMAT(5X,'*** WARNING : # VISIBLE SATELLITES ',I4, 00366 # ' > # OF MAXIMUM VISIBLE SATELLITES EXPECTED ',I4, 00367 # /,10X,'INCREASE MXVS AND THE ROW DIMENSIONS OF ', 00368**4 # 'THE APPROPRIATE ARRAYS TO AT LEAST',I4,' ***',/) 00369**4 35 FORMAT(14X,'BEGIN ALERT: DAY ',I4,' , ',I5,4X,2(I2,':'),I2,/) 36 FORMAT(14X,' END ALERT: DAY ',I4,' , ',I5,4X,2(I2,':'),I2,///) 37 FORMAT(15X,'HR MIN',6X,'SATELLITE(S)',4X,'ELEVATION(S)', # 6X,'AZIMUTH(S)',//) 134 FORMAT(3I2) 150 FORMAT(5X,'WARNING: VISIBLE SATELLITES GREATER THAN 3') 156 FORMAT(63X,'GPS SATELLITES',/) 161 FORMAT(61X,'TRANSIT SATELLITES',/) 166 FORMAT(61X,'LAGEOS SATELLITE(S)',/) 171 FORMAT(60X,'STARLETTE SATELLITE(S)',/) 8000 FORMAT(53X,'SATELLITE ELEVATION > ',I2,' DEGREES',//) 311 FORMAT(56X,'NNSS NODAL CROSSING DATA',/) 321 FORMAT(45X,'NNSS KEPLERIAN PARAMETERS OR MEAN ORBIT TRACKED',/) 331 FORMAT(62X,'GPS EPHEMERIS',/) 341 FORMAT(50X,'NASA PREDICTION BULLETIN (GAST USED)',/) 351 FORMAT(49X,'NASA PREDICTION BULLETIN (GAST UNUSED)',/) 461 FORMAT(2I3,I6,1X,I3,30X,I3,29X,I5,/) 462 FORMAT(2I3,I6,1X,2I3,27X,2I3,26X,2I5,/) 463 FORMAT(2I3,I6,1X,3I3,24X,3I3,23X,3I5,/) 464 FORMAT(2I3,I6,1X,4I3,21X,4I3,20X,4I5,/) 465 FORMAT(2I3,I6,1X,5I3,18X,5I3,17X,5I5,/) 466 FORMAT(2I3,I6,1X,6I3,15X,6I3,14X,6I5,/) 467 FORMAT(2I3,I6,1X,7I3,12X,7I3,11X,7I5,/) 468 FORMAT(2I3,I6,1X,8I3,9X,8I3,8X,8I5,/) 469 FORMAT(2I3,I6,1X,9I3,6X,9I3,5X,9I5,/) 481 FORMAT(14X,2I3,7X,I3,13X,I3,14X,I5) 482 FORMAT(14X,2I3,7X,2I3,10X,2I3,11X,2I5) 483 FORMAT(14X,2I3,2(7X,3I3),8X,3I5) 524 FORMAT(2I3,I6,1X,4I3,16X,4I3,15X,4I5,25X,4I3,F5.1,/) 525 FORMAT(2I3,I6,1X,5I3,13X,5I3,12X,5I5,20X,4I3,F5.1,/) 526 FORMAT(2I3,I6,1X,6I3,10X,6I3,9X,6I5,15X,4I3,F5.1,/) 527 FORMAT(2I3,I6,1X,7I3,7X,7I3,6X,7I5,10X,4I3,F5.1,/) 528 FORMAT(2I3,I6,1X,8I3,4X,8I3,3X,8I5,5X,4I3,F5.1,/) 529 FORMAT(2I3,I6,2(1X,9I3),9I5,4I3,F5.1,/) 601 FORMAT(1X,'HR MIN ZCNT',6X,' VISIBLE SATELLITES',18X, # 'ELEVATIONS',32X,'AZIMUTHS',//) 811 FORMAT(2I3,1X,I3,25X,I3,24X,I5,60X,F5.1,/) 812 FORMAT(2I3,1X,2I3,22X,2I3,21X,2I5,50X,2F5.1,/) 813 FORMAT(2I3,1X,3I3,19X,3I3,18X,3I5,35X,4F5.1,/) 814 FORMAT(2I3,1X,4I3,16X,4I3,15X,4I5,25X,5F5.1,/) 815 FORMAT(2I3,1X,5I3,13X,5I3,12X,5I5,20X,5F5.1,/) 816 FORMAT(2I3,1X,6I3,10X,6I3,9X,6I5,15X,5F5.1,/) 817 FORMAT(2I3,1X,7I3,7X,7I3,6X,7I5,10X,5F5.1,/) 818 FORMAT(2I3,1X,8I3,4X,8I3,3X,8I5,5X,5F5.1,/) 819 FORMAT(2I3,2(1X,9I3),9I5,5F5.1,/) 820 FORMAT(1X,'HR MIN',6X,' VISIBLE SATELLITES', 00349**2 # 12X,'ELEVATIONS',22X,'AZIMUTHS',24X,'GDOP PDOP HTDP HDOP' 00350**2 @ ,' VDOP') 00351**2 901 FORMAT(1X,122('-'),/) 902 FORMAT(14X,61('-')) 903 FORMAT(1X,129('-'),/) 1117 FORMAT(2I3,9I2,9F6.2,9F8.2) 2222 STOP 00370 END 00371 SUBROUTINE DOPS(A,B,IDA,VISAT,MXVS,J,C,GDOP,PDOP,HTDOP,HDOP,VDOP) C C NAME DOPS C C PURPOSE COMPUTE DILUTION OF PRECISION FACTORS C C VERSION MARCH 1984 - R. SANTERRE C C C PARAMETERS I A , B - WORKING ARRAYS OF DECLARED ROW C IDA .GE. 4 C I VISAT - ARRAY OF DIRECTION COSINES C TO THE VISIBLE SATELLITES C IN THE TANGENTIAL SYSTEM (NORTH,EAST,HGT) C I MXVS - DECLARED ROW DIMENSION OF 'VISAT' C I J - NUMBER OF VISIBLE SATELLITES C I C - SPEED OF LIGHT C I KK - VECTOR OF POSITION INDICATORS WITHIN C 'VISAT' FOR THE FOUR SELECTED SATELLITES C C O GDOP - SQRT( TRACE COV(N,E,H,CLOCK)) C O PDOP - SQRT( TRACE COV(N,E,H)) C O HTDOP - SQRT( TRACE COV(N,E,CLOCK)) C O HDOP - SQRT( TRACE COV(N,E)) C O VDOP - SQRT( TRACE COV(H)) C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(IDA,1) , B(IDA,1) , VISAT(MXVS,1) C C ---------------------------------------------------------------------- C COMPUTE 4-D DILUTION OF PRECISION FACTOR GDOP C ---------------------------------------------------------------------- C IF(J .LT. 4) GO TO 100 MA = J DO 10 I = 1,MA A(I,1) = - VISAT( I , 1 ) A(I,2) = - VISAT( I , 2 ) A(I,3) = - VISAT( I , 3 ) A(I,4) = C 10 CONTINUE C CALL MATMPY(A,A,B,4,MA,4,IDA,IDA,IDA,2) CALL SPIN(B,4,IDA,DET,IEXP) C GDOP = DSQRT( B(1,1) + B(2,2) + B(3,3) + B(4,4) * C*C ) 100 CONTINUE C C ---------------------------------------------------------------------- C COMPUTE 3-D DILUTION OF PRECISION FACTORS C ---------------------------------------------------------------------- C C C COMPUTE PDOP FACTOR FIRST C IF(J .LT. 3) GO TO 200 MA=J DO 15 I = 1,MA A(I,1) = - VISAT( I , 1 ) A(I,2) = - VISAT( I , 2 ) A(I,3) = - VISAT( I , 3 ) 15 CONTINUE C CALL MATMPY(A,A,B, 3,MA, 3,IDA,IDA,IDA,2) CALL SPIN(B, 3,IDA,DET,IEXP) C PDOP = DSQRT( B(1,1) + B(2,2) + B(3,3) ) C C THEN HTDOP FACTOR C DO 20 I = 1,MA A(I,1) = - VISAT( I , 1 ) A(I,2) = - VISAT( I , 2 ) A(I,3) = C 20 CONTINUE C CALL MATMPY(A,A,B, 3,MA, 3,IDA,IDA,IDA,2) CALL SPIN(B, 3,IDA,DET,IEXP) C HTDOP = DSQRT( B(1,1) + B(2,2) + B(3,3) * C*C ) 200 CONTINUE C C C ---------------------------------------------------------------------- C COMPUTE 2-D DILUTION OF PRECISION FACTOR HDOP C ---------------------------------------------------------------------- C IF(J .LT. 2) GO TO 300 MA=J DO 25 I = 1,MA A(I,1) = - VISAT( I , 1 ) A(I,2) = - VISAT( I , 2 ) 25 CONTINUE C CALL MATMPY(A,A,B, 2,MA, 2,IDA,IDA,IDA,2) CALL SPIN(B, 2,IDA,DET,IEXP) C HDOP = DSQRT( B(1,1) + B(2,2) ) 300 CONTINUE C C----------------------------------------------------------------------- C COMPUTE 1-D DILUTION OF PRECISION FACTOR VDOP C----------------------------------------------------------------------- MA=J DO 30 I=1,MA A(I,1) = - VISAT(I,3) 30 CONTINUE CALL MATMPY(A,A,B,1,MA,1,IDA,IDA,IDA,2) VDOP=DSQRT(1.D0/B(1,1)) RETURN END SUBROUTINE DOP4S(A,B,IDA,VISAT,MXVS,C,KK,GDOP,HTDOP,PDOP,HDOP) 00372**3 C 00373**3 C NAME DOP4S 00374**3 C 00375**3 C PURPOSE COMPUTE DILUTION OF PRECISION FACTORS 00376**3 C FROM FOUR SATELLITES 00377**3 C 00378**3 C VERSION FEBRUARY 1983 - D. DELIKARAOGLOU 00379**3 C 00380**3 C THIS VERSION USES THE FOUR SATELLITES WHICH 00381**3 C GIVE THE MAXIMUM VOLUME OF THE FUNDAMENTAL 00382**3 C TETRAHEDRON AND COMPUTES GDOP , HTDOP , 00383**3 C PDOP AND HDOP 00384**3 C 00385**3 C PARAMETERS I A , B - WORKING ARRAYS OF DECLARED ROW 00386**3 C IDA .GE. 4 00387**3 C I VISAT - ARRAY OF DIRECTION COSINES 00388**3 C TO THE VISIBLE SATELLITES 00389**3 C IN THE TANGENTIAL SYSTEM (NORTH,EAST,HGT) C I MXVS - DECLARED ROW DIMENSION OF 'VISAT' 00390**3 C I C - SPEED OF LIGHT 00391**3 C I KK - VECTOR OF POSITION INDICATORS WITHIN 00392**3 C 'VISAT' FOR THE FOUR SELECTED SATELLITES 00393**3 C 00394**3 C O GDOP - SQRT( TRACE COV(N,E,H,CLOCK)) 00395**3 C O HTDOP - SQRT( TRACE COV(N,E,CLOCK)) 00396**3 C O PDOP - SQRT( TRACE COV(N,E,H)) 00397**3 C O HDOP - SQRT( TRACE COV(N,E)) 00398**3 C 00399**3 C 00400**3 IMPLICIT REAL*8(A-H,O-Z) 00401**3 DIMENSION A(IDA,1) , B(IDA,1) , VISAT(MXVS,1) , KK(1) 00402**3 C 00403**3 C ----------------------------------------------------------------------00404**3 C COMPUTE 4-D DILUTION OF PRECISION FACTOR GDOP 00405**3 C ----------------------------------------------------------------------00406**3 C 00407**3 MA = 4 00408**3 DO 10 I = 1,MA 00409**3 A(I,1) = - VISAT( KK(I),1 ) 00410**3 A(I,2) = - VISAT( KK(I),2 ) 00411**3 A(I,3) = - VISAT( KK(I),3 ) 00412**3 A(I,4) = C 00413**3 10 CONTINUE 00414**3 C 00415**3 CALL MATMPY(A,A,B,MA,MA,MA,IDA,IDA,IDA,2) 00416**3 CALL SPIN(B,MA,IDA,DET,IEXP) 00417**3 C 00418**3 GDOP = DSQRT( B(1,1) + B(2,2) + B(3,3) + B(4,4) * C*C ) 00419**3 C 00420**3 C ----------------------------------------------------------------------00421**3 C COMPUTE 3-D DILUTION OF PRECISION FACTORS 00422**3 C ----------------------------------------------------------------------00423**3 C 00424**3 C 00425**3 C COMPUTE PDOP FACTOR FIRST 00426**3 C 00427**3 DO 15 I = 1,MA 00428**3 A(I,1) = - VISAT( KK(I),1 ) 00429**3 A(I,2) = - VISAT( KK(I),2 ) 00430**3 A(I,3) = - VISAT( KK(I),3 ) 00431**3 15 CONTINUE 00432**3 C 00433**3 CALL MATMPY(A,A,B, 3,MA, 3,IDA,IDA,IDA,2) 00434**3 CALL SPIN(B, 3,IDA,DET,IEXP) 00435**3 C 00436**3 PDOP = DSQRT( B(1,1) + B(2,2) + B(3,3) ) 00437**3 C 00438**3 C THEN HTDOP FACTOR 00439**3 C 00440**3 DO 20 I = 1,MA 00441**3 A(I,1) = - VISAT( KK(I),1 ) 00442**3 A(I,2) = - VISAT( KK(I),2 ) 00443**3 A(I,3) = C 00444**3 20 CONTINUE 00445**3 C 00446**3 CALL MATMPY(A,A,B, 3,MA, 3,IDA,IDA,IDA,2) 00447**3 CALL SPIN(B, 3,IDA,DET,IEXP) 00448**3 C 00449**3 HTDOP = DSQRT( B(1,1) + B(2,2) + B(3,3) * C*C ) 00450**3 C 00451**3 C 00452**3 C ----------------------------------------------------------------------00453**3 C COMPUTE 2-D DILUTION OF PRECISION FACTOR HDOP 00454**3 C ----------------------------------------------------------------------00455**3 C 00456**3 DO 25 I = 1,MA 00457**3 A(I,1) = - VISAT( KK(I),1 ) 00458**3 A(I,2) = - VISAT( KK(I),2 ) 00459**3 25 CONTINUE 00460**3 C 00461**3 CALL MATMPY(A,A,B, 2,MA, 2,IDA,IDA,IDA,2) 00462**3 CALL SPIN(B, 2,IDA,DET,IEXP) 00463**3 C 00464**3 HDOP = DSQRT( B(1,1) + B(2,2) ) 00465**3 RETURN END 00467**3 SUBROUTINE DOP3S(A,B,IDA,VISAT,MXVS,C,HTDOP,PDOP,HDOP) 00468**3 C 00469**3 C NAME DOP3S 00470**3 C 00471**3 C PURPOSE COMPUTE DILUTION OF PRECISION FACTORS 00472**3 C FROM THREE SATELLITES 00473**3 C 00474**3 C VERSION FEBRUARY 1983 - D. DELIKARAOGLOU 00475**3 C 00476**3 C PARAMETERS I A , B - WORKING ARRAYS OF DECLARED ROW 00477**3 C IDA .GE. 3 00478**3 C I VISAT - ARRAY OF DIRECTION COSINES 00479**3 C TO THE VISIBLE SATELLITES 00480**3 C IN THE TANGENTIAL SYSTEM(NORTH,EAST,HGT) C I MXVS - DECLARED ROW DIMENSION OF 'VISAT' 00481**3 C I C - SPEED OF LIGHT 00482**3 C 00483**3 C O HTDOP - SQRT( TRACE COV(N,E,CLOCK)) 00484**3 C O PDOP - SQRT( TRACE COV(N,E,H)) 00485**3 C O HDOP - SQRT( TRACE COV(N,E)) 00486**3 C 00487**3 C 00488**3 IMPLICIT REAL*8(A-H,O-Z) 00489**3 DIMENSION A(IDA,1) , B(IDA,1) , VISAT(MXVS,1) 00490**3 C 00491**3 MA = 3 00492**3 DO 10 I = 1,MA 00493**3 A(I,1) = - VISAT( I,1 ) 00494**3 A(I,2) = - VISAT( I,2 ) 00495**3 A(I,3) = C 00496**3 10 CONTINUE 00497**3 C 00498**3 CALL MATMPY(A,A,B,MA,MA,MA,IDA,IDA,IDA,2) 00499**3 CALL SPIN(B,MA,IDA,DET,IEXP) 00500**3 HTDOP = DSQRT( B(1,1) + B(2,2) + B(3,3) * C*C ) 00501**3 C 00502**3 DO 15 I = 1,MA 00503**3 A(I,1) = - VISAT( I,1 ) 00504**3 A(I,2) = - VISAT( I,2 ) 00505**3 A(I,3) = - VISAT( I,3 ) 00506**3 15 CONTINUE 00507**3 C 00508**3 CALL MATMPY(A,A,B,MA,MA,MA,IDA,IDA,IDA,2) 00509**3 CALL SPIN(B,MA,IDA,DET,IEXP) 00510**3 PDOP = DSQRT( B(1,1) + B(2,2) + B(3,3) ) 00511**3 C 00512**3 C ----------------------------------------------------------------------00513**3 C COMPUTE 2-D DILUTION OF PRECISION FACTOR HDOP 00514**3 C ----------------------------------------------------------------------00515**3 C 00516**3 DO 20 I = 1,MA 00517**3 A(I,1) = - VISAT( I,1 ) 00518**3 A(I,2) = - VISAT( I,2 ) 00519**3 20 CONTINUE 00520**3 C 00521**3 CALL MATMPY(A,A,B, 2,MA, 2,IDA,IDA,IDA,2) 00522**3 CALL SPIN(B, 2,IDA,DET,IEXP) 00523**3 C 00524**3 HDOP = DSQRT( B(1,1) + B(2,2) ) 00525**3 RETURN 00526**3 END 00527**3 SUBROUTINE DOP2S(A,B,IDA,VISAT,MXVS,C,HDOP) 00528**3 C 00529**3 C NAME DOP2S 00530**3 C 00531**3 C PURPOSE COMPUTE DILUTION OF PRECISION FACTOR 00532**3 C FROM TWO SATELLITES 00533**3 C 00534**3 C VERSION FEBRUARY 1983 - D. DELIKARAOGLOU 00535**3 C 00536**3 C PARAMETERS I A , B - WORKING ARRAYS OF DECLARED ROW 00537**3 C IDA .GE. 3 00538**3 C I VISAT - ARRAY OF DIRECTION COSINES 00539**3 C TO THE VISIBLE SATELLITES 00540**3 C IN THE TANGENTIAL SYSTEM (NORTH,EST,HGT) C I MXVS - DECLARED ROW DIMENSION OF 'VISAT' 00541**3 C I C - SPEED OF LIGHT 00542**3 C 00543**3 C O HDOP - SQRT( TRACE COV(N,E)) 00544**3 C 00545**3 C 00546**3 IMPLICIT REAL*8(A-H,O-Z) 00547**3 DIMENSION A(IDA,1) , B(IDA,1) , VISAT(MXVS,1) 00548**3 C 00549**3 MA = 2 00550**3 DO 10 I = 1,MA 00551**3 A(I,1) = - VISAT( I,1 ) 00552**3 A(I,2) = - VISAT( I,2 ) 00553**3 10 CONTINUE 00554**3 C 00555**3 CALL MATMPY(A,A,B,MA,MA,MA,IDA,IDA,IDA,2) 00556**3 CALL SPIN(B,MA,IDA,DET,IEXP) 00557**3 HDOP = DSQRT( B(1,1) + B(2,2) ) 00558**3 C 00559**3 RETURN 00560**3 END 00561**3 SUBROUTINE VOLUME(KSAT,VISAT,MXVS,NN,IDSAT,VOL) 00562 C 00563 C NAME VOLUME 00564 C 00565 C PURPOSE TO SELECT THE SET OF FOUR GPS SATELLITES 00566 C WHICH PROVIDE BEST NAVIGATION GEOMETRY 00567 C USING THE MAXIMUM TETRAHEDRON CRITERION 00568 C 00569 C PARAMETERS 00570 C 00571 C KSAT I POSITION OF THE SATELLITE WITH 00572 C MAXIMUM ELEVATION. 00573 C VISAT I MATRIX CONTAINING DIRECTION 00574 C COSINES OF ALL VISIBLE SAT. 00575 C MXVS I MAXIMUM POSSIBLE NUMBER OF VISIBLE 00576 C SATELLITES (EQUAL TO THE DECLARED ROW 00577 C DIMENSION OF MATRIX 'VISAT' IN CALLING 00578 C PROGRAM) 00579 C NN I NUMBER OF VISIBLE SATELLITES 00580 C IDSAT O ID OF VISIBLE SATELLITES WHICH 00581 C PROVIDE BEST NAVIGATION PERFORMANCE. 00582 C VOL O VOLUME OF THE TETRAHEDRON FORMED BY 00583 C UNIT VECTORS TO EACH SATELLITE. 00584 C 00585 IMPLICIT REAL*8(A-H,O-Z) 00586 DIMENSION VISAT(MXVS,3) , A(3), B(3) , C(3) , IDSAT(4) 00587 N1= NN - 1 00588 N2= NN - 2 00589 VOLMX=0.D0 00590 C 00591 DO 1 I=1,N2 00592 IF(I.EQ.KSAT) GO TO 1 00593 K= I + 1 00594 DO 2 J= K ,N1 00595 IF(J.EQ.KSAT) GO TO 2 00596 L= J + 1 00597 DO 3 M= L, NN 00598 IF(M.EQ.KSAT) GO TO 3 00599 C 00600 DO 4 II= 1 , 3 00601 A(II)= VISAT(I,II) - VISAT(KSAT,II) 00602 B(II)= VISAT(J,II) - VISAT(KSAT,II) 00603 C(II)= VISAT(M,II) - VISAT(KSAT,II) 00604 4 CONTINUE 00605 C 00606 VOL= (A(1)*(B(2)*C(3) - B(3)*C(2)) 00607 # - A(2)*(B(1)*C(3) - B(3)*C(1)) 00608 # + A(3)*(B(1)*C(2) - B(2)*C(1))) / 6. 00609 C 00610 IF(DABS(VOL).GT.VOLMX) GO TO 10 00611 GO TO 11 00612 10 VOLMX= DABS(VOL) 00613 IDSAT(1)= KSAT 00614 IDSAT(2)= J 00615 IDSAT(3)= M 00616 IDSAT(4)= I 00617 11 CONTINUE 00618 3 CONTINUE 00619 2 CONTINUE 00620 1 CONTINUE 00621 VOL= VOLMX 00622 RETURN 00623 END 00624 C--------------------------------------------------------------------- SUBROUTINE TIME(IYEAR0,IDAY0,IH0,IM0,IS0,IYEART,IDAYT, # IHT,IMT,IST,ITIME) C--------------------------------------------------------------------- C NAME: TIME C PURPOSE: TO COMPUTE THE NUMBER OF SECONDS BETWEEN TWO C GIVEN INSTANTS C C VERSION : JUNE 1983 , R. SANTERRE C C PARAMETERS: C I IYEAR0: YEAR OF THE INITIAL INSTANT C I IDAY0: DAY OF THE INITIAL INSTANT C I IH0: HOUR OF THE INITIAL INSTANT C I IM0: MINUTE OF THE INITIAL INSTANT C I IS0: SECOND OF THE INITIAL INSTANT C I IYEART: YEAR OF THE FINAL INSTANT C I IDAYT: DAY OF THE FINAL INSTANT C I IHT: HOUR OF THE FINAL INSTANT C I IMT: MINUTE OF THE FINAL INSTANT C I IST: SECOND OF THE FINAL INSTANT C C O ITIME: NUMBER OF SECONDS C C COMMENT: IYEART-IYEAR0 EQUAL -1, 0, OR 1 C--------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C C ----------------------------------------------------------- C TEST IF THE FINAL TIME PRECEDES OR FOLLOWS THE INITIAL TIME C ----------------------------------------------------------- IF(IYEART-IYEAR0) 10,20,30 C 10 IYEAR=IYEART C C ----------------------------------------------------- C TEST IF IYEART IS LEAP AND COMPUTE THE NUMBER OF DAYS C ----------------------------------------------------- IF(IYEAR/4 .NE. IYEAR/4.) GO TO 100 IDD=-366+IDAYT-IDAY0 GO TO 40 100 IDD=-365+IDAYT-IDAY0 GO TO 40 C C -------------------------- C COMPUTE THE NUMBER OF DAYS C -------------------------- 20 IDD=IDAYT-IDAY0 GO TO 40 C 30 IYEAR=IYEAR0 C C ----------------------------------------------------- C TEST IF IYEAR0 IS LEAP AND COMPUTE THE NUMBER OF DAYS C ----------------------------------------------------- IF(IYEAR/4 .NE. IYEAR/4.) GO TO 200 IDD=366+IDAYT-IDAY0 GO TO 40 200 IDD=365+IDAYT-IDAY0 C C ------------------------------------------------------- C COMPUTE THE NUMBER OF SECONDES BETWEEN THE TWO INSTANTS C ------------------------------------------------------- 40 IHH=IDD*24+IHT-IH0 IMM=IHH*60+IMT-IM0 ITIME=IMM*60+IST-IS0 C RETURN END C--------------------------------------------------------------------- SUBROUTINE DAYFRA(DAYF,IDAY,IH,IM,IS) C--------------------------------------------------------------------- C NAME: DAYFRA C PURPOSE: TO TRANSFORM DAY AND FRACTION OF DAY C IN DAY,HOUR,MINUTE AND SECOND C C VERSION : JUNE 1983 , R. SANTERRE C C PARAMETERS: C I DAYF: DAY AND FRACTION OF DAY C C O IDAY: DAY C O IH: HOUR C O IM: MINUTE C O IS: SECONDE C--------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C IDAY=IDINT(DAYF) R=(DAYF-DFLOAT(IDAY))*24.D0 IH=IDINT(R) U=(R-DFLOAT(IH))*60.D0 IM= IDINT(U) T=(U-DFLOAT(IM))*60.D0 IS=IDINT(T) IF((T-IS) .GE. .5) IS=IS+1 C RETURN END C--------------------------------------------------------------------- SUBROUTINE MINFRA(XMINF,IH,IM,IS) C--------------------------------------------------------------------- C NAME: MINFRA C PURPOSE: TO TRANSFORM MINUTE AND FRACTION OF MINUTE C IN HOUR,MINUTE AND SECOND C C VERSION : JUNE 1983 , R. SANTERRE C C PARAMETERS: C I XMINF: MINUTE AND FRACTION OF MINUTE C C O IH: HOUR C O IM: MINUTE C O IS: SECOND C--------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C R=XMINF/60.D0 IH=IDINT(R) C U=(R-DFLOAT(IH))*60.D0 IM=IDINT(U) C T=(U-DFLOAT(IM))*60.D0 IS=IDINT(T) IF((T-IS) .GE. .5) IS=IS+1 C RETURN END C--------------------------------------------------------------------- SUBROUTINE NNSSNC(IY,MDAY,TNODE,XLNODE,PERIOD,WMOTN,KYEAR, # KDAY,KHR,KMIN,KSEC,PI,XS,YS,ZS) C--------------------------------------------------------------------- C NAME: NNSSNC C PURPOSE: TO COMPUTE THE SATELLITE COORDINATES IN THE EARTH FIXED C SYSTEM FROM NNSS NODAL CROSSING DATA C C VERSION : JUNE 1983 , R. SANTERRE C C PARAMETERS: C I IY: INPUT YEAR C I MDAY: INPUT DAY C I TNODE: TIME OF ASCENDING NODE (MINUTES) C I XLNODE: LONGITUDE (EAST) OF ASCENDING NODE (DEGREES) C I PERIOD: NODAL PERIOD (MINUTES) C I WMOTN: WESTWARD MOTION (DEGREES) C I KYEAR: YEAR OF THE ALERT C I KDAY: DAY OF THE ALERT C I KHR: HOUR OF THE ALERT C I KMIN: MINUTE OF THE ALERT C I KSEC: SECOND OF THE ALERT C I PI: VALUE OF PI C C O XS: SATELLITE X COORDINATE (METRES) C O YS: SATELLITE Y CORRDINATE (METRES) C O ZS: SATELLITE Z COORDINATE (METRES) C--------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C C -------------------------------------------- C TRANSFORM TNODE IN HOURS,MINUTES AND SECONDS C -------------------------------------------- CALL MINFRA(TNODE,IH,IM,IS) C C ----------------------------------------------------------------- C COMPUTE THE NUMBER OF SECONDES BETWEEN THE TIME OF ASCENDING NODE C AND THE TIME OF THE ALERT C ----------------------------------------------------------------- CALL TIME(IY,MDAY,IH,IM,IS,KYEAR,KDAY,KHR,KMIN,KSEC,ITIME) C C ----------------------------------------------------------- C COMPUTE LONGITUDE OF ASCENDING NODE (RAD) AT THE ALERT TIME C ----------------------------------------------------------- FF=WMOTN/PERIOD/60.D0*DFLOAT(ITIME) RLANT=(XLNODE-FF)*PI/180.D0 C C --------------------------------------------------- C COMPUTE ARGUMENT OF PERIGEE (RAD) AT THE ALERT TIME C --------------------------------------------------- RPHIT=2.D0*PI/(PERIOD*60.D0)*DFLOAT(ITIME) C C ---------------------------------------------------------- C COMPUTE SATELLITE COORDINATES IN THE ORBITAL PLAN (METRES) C ---------------------------------------------------------- R=7456.D3 XPT=R*DCOS(RPHIT) YPT=R*DSIN(RPHIT) C C ---------------------------------------------------------------- C COMPUTE SATELLITE COORDINATES IN THE EARTH FIXED SYSTEM (METRES) C ---------------------------------------------------------------- XS=XPT*DCOS(RLANT) YS=XPT*DSIN(RLANT) ZS=YPT C RETURN END C--------------------------------------------------------------------- SUBROUTINE NNSSKE(IY,MDAY,TP,XN,W,WDOT,EC,AO,OM,OMDOT,COSI,GAST, # SINI,KYEAR,KDAY,KHR,KMIN,KSEC,WE,PI,XS,YS,ZS, # IPR) C--------------------------------------------------------------------- C NAME: NNSSKE C PURPOSE: TO COMPUTE THE SATELLITE COORDINATES IN THE EARTH FIXED C SYSTEM (WGS-72) FROM NNSS KEPLERIAN PARAMETERS C C VERSION : JUNE 1983 , R. SANTERRE C C PARAMETERS: C I IY: INPUT YEAR C I MDAY: INPUT DAY C I TP: TIME OF PERIGEE (MINUTES) C I XN: MEAN MOTION (DEGREES/MINUTE) C I W: ARGUMENT OF PERIGEE (DEGREES) C I WDOT: ARGUMENT OF PERIGEE DOT (DEGREE/MINUTE) C I EC: ECCENTRICITY C I AO: SEMI-MAJOR AXIS (METRES) C I OM: RIGHT ASCENSION OF ASCENDING NODE (DEGREES) C I OMDOT: RIGHT ASCENSION DOT (DEGREE/MINUTE) C I COSI: COSINE OF INCLINAISON C I GAST: RIGHT ASCENSION OF GREENWICH (DEGREES) C I SINI: SINE OF INCLINAISON C I KYEAR: YEAR OF THE ALERT C I KDAY: DAY OF THE ALERT C I KHR: HOUR OF THE ALERT C I KMIN: MINUTE OF THE ALERT C I KSEC: SECOND OF THE ALERT C I WE: EARTH ROTATION RATE (RAD/SEC) C I PI: VALUE OF PI C I IPR: PRINT CODE FOR LINE CODE C C O XS: SATELLITE X COORDINATE (METRES) C O YS: SATELLITE Y COORDINATE (METRES) C O ZS: SATELLITE Z COORDINATE (METRES) C--------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C C ------------------------------------------------------------- C TRANSFORM THE TIME OF PERIGEE IN HOURS , MINUTES AND SECONDES C ------------------------------------------------------------- CALL MINFRA(TP,IH,IM,IS) C C ---------------------------------------------------------- C COMPUTE THE NUMBER OF SECONDES BETWEEN THE TIME OF PERIGEE C AND THE TIME OF THE ALERT C ---------------------------------------------------------- CALL TIME(IY,MDAY,IH,IM,IS,KYEAR,KDAY,KHR,KMIN,KSEC,ITIME) C C --------------------------------------------------------------- C COMPUTE LONGITUDE OF ASCENDING NODE (RAD) AT THE REFERENCE TIME C --------------------------------------------------------------- RLAN0=(OM-GAST)*PI/180.D0 C C ----------------------------------------------------------- C COMPUTE LONGITUDE OF ASCENDING NODE (RAD) AT THE ALERT TIME C ----------------------------------------------------------- RLANT=RLAN0-(WE-(OMDOT*PI/10800.D0))*DFLOAT(ITIME) C C -------------------------------------------- C COMPUTE MEAN ANOMALY (RAD) AT THE ALERT TIME C -------------------------------------------- RXMAT=XN*PI/10800.D0*DFLOAT(ITIME) C C ---------------------------------------------------- C COMPUTE ECCENTRIC AND TRUE ANOMALY FROM MEAN ANOMALY C ---------------------------------------------------- CALL ANMLY(RXMAT,EC,REAT,RTAT,IFLAG,PI,IPR) C C ---------------------------------------------------- C COMPUTE ARGUMENT OF LATITUDE (RAD) AT THE ALERT TIME C ---------------------------------------------------- RPHIT=RTAT+W*PI/180.D0-DABS(WDOT)*PI/10800.D0*DFLOAT(ITIME) C C ----------------------------- C COMPUTE ORBIT RADIUS (METRES) C ----------------------------- RT=AO*(1.D0-EC*DCOS(REAT))*1000.D0 C C ---------------------------------------------------------- C COMPUTE SATELLITE COORDINATES IN THE ORBITAL PLAN (METRES) C ---------------------------------------------------------- XPT=RT*DCOS(RPHIT) YPT=RT*DSIN(RPHIT) C C ------------------------------------------------------- C COMPUTE SATELLITE COORDINATES IN THE EARTH FIXED SYSTEM C ------------------------------------------------------- XS=XPT*DCOS(RLANT)-YPT*COSI*DSIN(RLANT) YS=XPT*DSIN(RLANT)+YPT*COSI*DCOS(RLANT) ZS=YPT*SINI C RETURN END C--------------------------------------------------------------------- SUBROUTINE NPBST(IY,DAY,XINC,RA,EC,AP,XMA,XN,PN,WMTN, # KYEAR,KDAY,KHR,KMIN,KSEC,PI,WE, # GM,XS,YS,ZS,IPR) C--------------------------------------------------------------------- C NAME: NPBST C PURPOSE: TO COMPUTE THE SATELLITE COORDINATES IN THE EARTH FIXED C SYSTEM (WGS-72) FROM NASA PREDICTION BULLETIN (GAST USED) C C VERSION : JUNE 1983 , R. SANTERRE C C PARAMETERS: C I IY: INPUT YEAR C I DAY: DAY AND FRACTION OF DAY OF THE INPUT C I XINC: INCLINAISON (DEGREES) C I RA: RIGHT ASCENSION OF THE ASCENDING NODE (DEGREES) C I EC: ECCENTRICITY C I AP: ARGUMENT OF PERIGEE (DEGREES) C I XMA: MEAN ANOMALY (DEGREES) C I XN: MEAN MOTION (REVOLUTIONS/DAY) C I PN: NODAL PERIOD (MINUTES) C I WMTN: WESTWARD MOTION (DEGREES) C I KYEAR: YEAR OF THE ALERT C I KDAY: DAY OF THE ALERT C I KHR: HOUR OF THE ALERT C I KMIN: MINUTE OF THE ALERT C I KSEC: SECOND OF THE ALERT C I PI: VALUE OF PI C I WE: EARTH ROTATION RATE (RAD/SEC) C I GM: UNIVERSAL GRAVITATIONAL CONSTANT (M*3/SEC*2) C I IPR: PRINT CODE FOR LINE CODE C C O XS: SATELLITE X COORDINATE (METRES) C O YS: SATELLITE Y COORDINATE (METRES) C O ZS: SATELLITE Z COORDINATE (METRES) C--------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C C ------------------------------------------------------------ C TRANSFORM DAY AND FRACTIONAL PORTION OF DAY IN DAYS , HOURS, C MINUTES AND SECONDS C ------------------------------------------------------------ CALL DAYFRA(DAY,IDAY,IH,IM,IS) C C ---------------------------------------------------------------- C COMPUTE THE NUMBER OF SECONDES BETWEEN THE TIME OF REFERENCE AND C THE TIME OF THE ALERT C ---------------------------------------------------------------- CALL TIME(IY,IDAY,IH,IM,IS,KYEAR,KDAY,KHR,KMIN,KSEC,ITIME) C C ---------------------------------------- C COMPUTE GAST (RAD) AT THE REFERENCE TIME C ---------------------------------------- XMJD0=DJUL(IY-1,12,31.D0) XMJD=XMJD0+DAY RGAST=THETA(XMJD) C C C ------------------------------------- C COMPUTE RIGHT ASCENSION DOT (RAD/SEC) C ------------------------------------- AA=WMTN/PN*PI/10800.D0 RSRAD=WE-AA C C ----------------------------- C COMPUTE MEAN MOTION (RAD/SEC) C ----------------------------- RSXN=XN*2.D0*PI/86400.D0 C C ----------------------------------------- C COMPUTE ARGUMENT OF PERIGEE DOT (RAD/SEC) C ----------------------------------------- BB=2.D0*PI/PN/60.D0 RSAPD=BB-RSXN C C -------------------------------------- C COMPUTE ORBIT SEMI-MAJOR AXIS (METRES) C -------------------------------------- AO=(GM/RSXN**2.D0)**(1.D0/3.D0) C C -------------------------------------------- C COMPUTE MEAN ANOMALY (RAD) AT THE ALERT TIME C -------------------------------------------- RXMAT=XMA*PI/180.D0+RSXN*DFLOAT(ITIME) C C ---------------------------------------------------- C COMPUTE ECCENTRIC AND TRUE ANOMALY FROM MEAN ANOMALY C ---------------------------------------------------- CALL ANMLY(RXMAT,EC,REAT,RTAT,IFLAG,PI,IPR) C C --------------------------------------------------- C COMPUTE ARGUMENT OF PERIGEE (RAD) AT THE ALERT TIME C --------------------------------------------------- RAPT=AP*PI/180.D0+RSAPD*DFLOAT(ITIME) C C ---------------------------------------------------- C COMPUTE ARGUMENT OF LATITUDE (RAD) AT THE ALERT TIME C ---------------------------------------------------- RPHIT=RTAT+RAPT C C ----------------------------- C COMPUTE ORBIT RADIUS (METRES) C ----------------------------- RT=AO*(1.D0-EC*DCOS(REAT)) C C ---------------------------------------------------------- C COMPUTE SATELLITE COORDINATES IN THE ORBITAL PLAN (METRES) C ---------------------------------------------------------- XPT=RT*DCOS(RPHIT) YPT=RT*DSIN(RPHIT) C C ----------------------------------------------------------- C COMPUTE LONGITUDE OF ASCENDING NODE (RAD) AT THE ALERT TIME C ----------------------------------------------------------- CC=(RSRAD-WE)*DFLOAT(ITIME)-RGAST RLANT=RA*PI/180.D0+CC C C ---------------------------------------------------------------- C COMPUTE SATELLITE COORDINATES IN THE EARTH FIXED SYSTEM (METRES) C ---------------------------------------------------------------- XS=XPT*DCOS(RLANT)-YPT*DCOS(XINC*PI/180.D0)*DSIN(RLANT) YS=XPT*DSIN(RLANT)+YPT*DCOS(XINC*PI/180.D0)*DCOS(RLANT) ZS=YPT*DSIN(XINC*PI/180.D0) C RETURN END C--------------------------------------------------------------------- SUBROUTINE NPB(IY,DAY,XINC,EC,AP,XMA,XN,PN,WMTN,XLANW, # KYAN,KDAYAN,KHAN,XMAN,KYEAR,KDAY,KHR,KMIN, # KSEC,PI,WE,GM,XS,YS,ZS,IPR) C--------------------------------------------------------------------- C NAME: NPB C PURPOSE: TO COMPUTE THE SATELLITE COORDINATES IN THE EARTH FIXED C SYSTEM (WGS-72) FROM NASA PREDICTION BULLETIN C C VERSION : JUNE 1983 , R. SANTERRE C C PARAMETERS: C I IY: INPUT YEAR C I DAY: INPUT DAY (DAY AND FRACTION OF DAY) C I XINC: INCLINAISON (DEGREES) C I EC: ECCENTRICITY C I AP: ARGUMENT OF PERIGEE (DEGREES) C I XMA: MEAN ANOMALY (DEGREES) C I XN: MEAN MOTION (REVOLUTIONS/DAY) C I PN: NODAL PERIOD (MINUTES) C I WMTN: WESTWARD MOTION (DEGREES) C I XLANW: LONGITUDE OF ASCENDING NODE (DEGREES WEST) C I KYAN: YEAR OF ASCENDING NODE C I KDAYAN: DAY OF ASCENDING NODE C I KHAN: HOUR OF ASCENDING NODE C I XMAN: MINUTES AND FRACTION OF MINUTE OF ASCENDING NODE C I KYEAR: YEAR OF THE ALERT C I KDAY: DAY OF THE ALERT C I KHR: HOUR OF THE ALERT C I KMIN: MINUTE OF THE ALERT C I KSEC: SECOND OF THE ALERT C I PI: VALUE OF PI C I WE: EARTH ROTATION RATE (RAD/SEC) C I GM: UNIVERSAL GRAVITATIONAL CONSTANT (M*3/SEC*2) C I IPR: PRINT CODE FOR LINE CODE C C O XS: SATELLITE X COORDINATE (METRES) C O YS: SATELLITE Y COORDINATE (METRES) C O ZS: SATELLITE Z COORDINATE (METRES) C--------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) C C --------------------------------------------------------- C TRANSFORM DAY AND FRACTIONAL PORTION OF DAY (TIME INPUTS) C IN HOURS , MINUTES AND SECONDES C --------------------------------------------------------- CALL DAYFRA(DAY,IDAY,IH,IM,IS) C C -------------------------------------------------------- C TRANSFORM MINUTES AND FRACTIONAL PORTION OF MINUTE C (TIME OF ASCENDIND NODE) IN HOURS , MINUTES AND SECONDES C -------------------------------------------------------- CALL MINFRA(XMAN,IHAN,IMAN,ISAN) C C ------------------------------------------------------------- C COMPUTE NUMBER OF SECONDES BETWEEN THE TIME OF ASCENDING NODE C AND THE INPUTS TIME C ------------------------------------------------------------- CALL TIME(KYAN,KDAYAN,KHAN,IMAN,ISAN,IY,IDAY,IH,IM,IS,ITIME1) C C ---------------------------------------------------- C COMPUTE MEAN ANOMALY (DEGREES) AT THE REFERENCE TIME C (TIME OF ASCENDING NODE) C ---------------------------------------------------- XMA=XMA-XN*360.D0/86400.D0*DFLOAT(ITIME1) C C --------------------------------------------------------- C COMPUTE THE NUMBER OF SECONDES BETWEEN THE REFERENCE TIME C (TIME OF ASCENDING NODE) AND THE TIME OF THE ALERT C --------------------------------------------------------- CALL TIME(KYAN,KDAYAN,KHAN,IMAN,ISAN,KYEAR,KDAY,KHR,KMIN,KSEC, # ITIME) C C ------------------------------------- C COMPUTE RIGHT ASCENSION DOT (RAD/SEC) C ------------------------------------- AA=WMTN/PN*PI/10800.D0 RSRAD=WE-AA C C ----------------------------- C COMPUTE MEAN MOTION (RAD/SEC) C ----------------------------- RSXN=XN*2.D0*PI/86400.D0 C C ----------------------------------------- C COMPUTE ARGUMENT OF PERIGEE DOT (RAD/SEC) C ----------------------------------------- BB=2.D0*PI/PN/60.D0 RSAPD=BB-RSXN C C -------------------------------------- C COMPUTE ORBIT SEMI-MAJOR AXIS (METRES) C -------------------------------------- AO=(GM/RSXN**2.D0)**(1.D0/3.D0) C C -------------------------------------------- C COMPUTE MEAN ANOMALY (RAD) AT THE ALERT TIME C -------------------------------------------- RXMAT=XMA*PI/180.D0+RSXN*DFLOAT(ITIME) C C ---------------------------------------------------- C COMPUTE ECCENTRIC AND TRUE ANOMALY FROM MEAN ANOMALY C ---------------------------------------------------- CALL ANMLY(RXMAT,EC,REAT,RTAT,IFLAG,PI,IPR) C C --------------------------------------------------- C COMPUTE ARGUMENT OF PERIGEE (RAD) AT THE ALERT TIME C --------------------------------------------------- RAPT=AP*PI/180.D0+RSAPD*DFLOAT(ITIME-ITIME1) C C ---------------------------------------------------- C COMPUTE ARGUMENT OF LATITUDE (RAD) AT THE ALERT TIME C ---------------------------------------------------- RPHIT=RTAT+RAPT C C ----------------------------- C COMPUTE ORBIT RADIUS (METRES) C ----------------------------- RT=AO*(1.D0-EC*DCOS(REAT)) C C ---------------------------------------------------------- C COMPUTE SATELLITE COORDINATES IN THE ORBITAL PLAN (METRES) C ---------------------------------------------------------- XPT=RT*DCOS(RPHIT) YPT=RT*DSIN(RPHIT) C C ----------------------------------------------------------- C COMPUTE LONGITUDE OF ASCENDING NODE (RAD) AT THE ALERT TIME C ----------------------------------------------------------- CC=(RSRAD-WE)*DFLOAT(ITIME) RLANT=(360.D0-XLANW)*PI/180.D0+CC C C ---------------------------------------------------------------- C COMPUTE SATELLITE COORDINATES IN THE EARTH FIXED SYSTEM (METRES) C ---------------------------------------------------------------- XS=XPT*DCOS(RLANT)-YPT*DCOS(XINC*PI/180.D0)*DSIN(RLANT) YS=XPT*DSIN(RLANT)+YPT*DCOS(XINC*PI/180.D0)*DCOS(RLANT) ZS=YPT*DSIN(XINC*PI/180.D0) C RETURN END C--------------------------------------------------------------------- SUBROUTINE EPNCD(NSAT,IY,MDAY,ISAT,PERIOD,WMOTN,TNODE, # XLNODE,ICR) C C TO READ INPUTS NNSS NODAL CROSSING DATA C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION IY(NSAT) , MDAY(NSAT) , ISAT(NSAT) , # PERIOD(NSAT) , WMOTN(NSAT) , TNODE(NSAT) , # XLNODE(NSAT) C DO 10 I=1,NSAT READ(ICR,20) IY(I) , MDAY(I) , ISAT(I) , # PERIOD(I), WMOTN(I) , TNODE(I), # XLNODE(I) C 20 FORMAT(3I5,4F10.3) 10 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE EPKE(NSAT,IY,MDAY,ISAT,TP,XN,W,WDOT,EC,AO,OM,OMDOT, # COSI,GAST,SINI,ICR) C C TO READ INPUTS NNSS KEPLERIAN PARAMETERS C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION IY(NSAT) , MDAY(NSAT) , ISAT(NSAT) , # TP(NSAT) , XN(NSAT) , W(NSAT) , # WDOT(NSAT) , EC(NSAT) , AO(NSAT) , # OM(NSAT) , OMDOT(NSAT) , COSI(NSAT) , # GAST(NSAT) , SINI(NSAT) C DO 10 I=1,NSAT C READ(ICR,20) IY(I) , MDAY(I) , ISAT(I) , # TP(I) , XN(I) , W(I) , # WDOT(I) , EC(I) , AO(I) , # OM(I) , OMDOT(I) , COSI(I) , # GAST(I) , SINI(I) C 20 FORMAT(3I5,2(F10.4,F10.7),F10.6,/,F10.2,F10.4,F10.7,F10.6, # F10.4,F10.6) C 10 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE EPNPBG(NSAT,IDSAT,IY,DAY,XINC,RA,EC,AP,XMA,XN,PN, # WMTN,ICR) C C TO READ INPUTS NASA PREDICTION BULLETIN (GAST USED) C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION IDSAT(NSAT) , IY(NSAT) , DAY(NSAT) , # XINC(NSAT) , RA(NSAT) , EC(NSAT) , # AP(NSAT) , XMA(NSAT) , XN(NSAT) , # PN(NSAT) , WMTN(NSAT) C DO 10 I=1,NSAT C READ(ICR,20) IDSAT(I) , IY(I) , DAY(I) , # XINC(I) , RA(I) , EC(I) , # AP(I) , XMA(I) , XN(I) , # PN(I) , WMTN(I) C 20 FORMAT(2I5,F13.8,2F10.4,F10.7,F10.4,/,F10.4,F12.8,2F9.3) 10 CONTINUE RETURN END C--------------------------------------------------------------------- SUBROUTINE EPNPB(NSAT,IDSAT,IY,DAY,XINC,EC,AP,XMA,XN,PN,WMTN, # XLANW,KYAN,KDAYAN,KHAN,XMAN,ICR) C C TO READ INPUTS NASA PREDICTION BULLETIN C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION IDSAT(NSAT) , IY(NSAT) , DAY(NSAT) , # XINC(NSAT) , EC(NSAT) , AP(NSAT) , # XMA(NSAT) , XN(NSAT) , PN(NSAT) , # WMTN(NSAT) , XLANW(NSAT) , KYAN(NSAT) , # KDAYAN(NSAT) , KHAN(NSAT) , XMAN(NSAT) C DO 10 I=1,NSAT C READ(ICR,20) IDSAT(I) , IY(I) , DAY(I) , # XINC(I) , EC(I) , AP(I) , # XMA(I) , XN(I) , PN(I) , # WMTN(I) , XLANW(I) , KYAN(I) , # KDAYAN(I) , KHAN(I) , XMAN(I) C 20 FORMAT(2I5,F13.8,F10.4,F10.7,2F10.4,/,F12.8,2F9.3,F8.2,2I5,I3, # F6.2) 10 CONTINUE RETURN END SUBROUTINE DATUM ( DTM , IDTM , LUI ) IMPLICIT REAL * 8 ( A-H , O-Z ) DIMENSION DTM(7) C$ C NAME DATUM C C PURPOSE INITIALIZE DATUM PARAMETERS C C CALLING CALL DATUM ( DTM , IDTM , LUI ) C C PARAMETERS O DTM DATUM PARAMETERS C DTM(1) ELLIPSOID SEMIMAJOR AXIS (M) C DTM(2) ELLIPSOID ECCENTRICITY C DTM(3) X TRANSLATION COMPONENT (M) C DTM(4) Y TRANSLATION COMPONENT (M) C DTM(5) Z TRANSLATION COMPONENT (M) C DTM(6) GRAVITATIONAL CONSTANT (M**3/SEC**2) C DTM(7) EARTH ROTATION RATE (RAD/SEC) C I IDTM DEFAULT SWITCH C 0 = TAKE DEFAULT VALUES C NONZERO = READ USER SPECIFIED VALUES FROM LUI C I LUI INPUT LOGICAL UNIT NUMBER C$ IF(IDTM.NE.0) GO TO 3 C-------------------------------------------------------------- C WGS-72 DATUM PARAMETERS C-------------------------------------------------------------- DTM(1) = 6378135.0D0 DTM(2) = 8.181881123D-2 DTM(3) = 0.D0 DTM(4) = 0.D0 DTM(5) = 0.D0 GO TO 4 C------------------------------------------------------------- C USER DEFINED DATUM C------------------------------------------------------------- 3 READ(LUI,100) (DTM(I),I=1,7) 100 FORMAT(7G12.4) BDVDA = DTM(2) / DTM(1) DTM(2) = DSQRT ( (1.D0 + BDVDA) * (1.D0 - BDVDA) ) C 4 DTM(6) = 398600.8D9 DTM(7) = 7.292115147D-5 RETURN END DOUBLE PRECISION FUNCTION ZCOUNT(IYR,IDAY,IHR,MIN,ISEC) IMPLICIT REAL * 8 ( A - H , O - Z ) C$ C NAME ZCOUNT C C PURPOSE GIVEN A TIME EPOCH, COMPUTE THE Z-COUNT C C CALLING X = ZCOUNT(IYR,IDAY,IHR,MIN,ISEC) C C PARAMETERS I IYR YEAR C I IDAY DAY IN YEAR C I IHR HOUR C I MIN MINUTE C I ISEC SECOND C O ZCOUNT NUMBER OF 6-SECOND INTERVALS C SINCE SATURDAY/SUNDAY MIDNIGHT C$ C REFERENCE YEAR JYEAR=1982 C DAY=0 FOR 1982 IS THURSDAY=5 (I.E. DAY 5 IN THE WEEK) JYEAR = 1982 JDAY0 = 5 C THE TERM (IYR-JYEAR+2)/4 IS TO ACCOUNT FOR LEAP YEARS IWEEK = JDAY0 + (IYR - JYEAR) + (IYR - JYEAR + 2) / 4 IWEEK = MOD( (IWEEK+IDAY) , 7) IF(IWEEK .LE. 0) IWEEK = IWEEK + 7 ZCOUNT = (DBLE(IHR*3600 + MIN*60 + ISEC) # + (IWEEK-1)*86400.D0 ) / 6.D0 RETURN END SUBROUTINE READEF(REC,DTM,EPH,EPHBUF,CLK,CLKBUF,NSAT,LU,IPR) IMPLICIT REAL * 8 ( A-H , O-Z ) DIMENSION REC(1) , DTM(1) , EPH(1) , EPHBUF(1) , CLK(1) , $ CLKBUF(1) DATA PI / 3.141592653589793D0 / C$ C NAME READEF C C PURPOSE READE IS TO READ GPS FILES IN THE EXCHANGE SHELTECH/IBM C FORMAT C C CALLING CALL READEF(REC,DTM,EPH,EPHBUF,CLK,CLKBUF,NSAT,LU,IPR) C C PARAMETERS I REC = RECORD DATA BUFFER C I DTM = DATUM PARAMETERS C DTM(1) = ELLIPSOID SEMIMAJOR AXIS (M) C DTM(2) = ELLIPSOID ECCENTRICITY C DTM(3) = X TRANSLATION COMPONENT (M) C DTM(4) = Y TRANSLATION COMPONENT (M) C DTM(5) = Z TRANSLATION COMPONENT (M) C DTM(6) = GRAVITATIONAL CONSTANT (M**3/SEC**2) C DTM(7) = EARTH ROTATION RATE (RAD/SEC) C I/O CLK = SATELLITE CLOCK PARAMETERS C I LU = INPUT DATA LOGICAL UNIT C I IPR = PRINT CODE FOR LINE PRINTER C O EPH = SATELLITE EPHEMERIS PARAMETERS C EPH(1) = AGE OF EPHEMERIS (SEC) C EPH(2) = RADIUS SIN HARMONIC AMPL (M) C EPH(3) = CORRECT COMPUTED MEAN MOTION (RAD/SEC) C EPH(4) = MEAN ANOMALY AT REF TIME (RAD) C EPH(5) = LATITUDE COS HARMONIC AMPL (RAD) C EPH(6) = ECCENTRICITY C EPH(7) = LATITUDE SIN HARMONIC AMPL (RAD) C EPH(8) = ORBIT SEMIMAJOR AXIS (M) C EPH(9) = EPHEMERIS REF TIME (SEC) C EPH(10)= INCLINATION COS HARMONIC AMPL (RAD) C EPH(11)= RIGHT ASCENSION AT REF TIME (RAD) C EPH(12)= INCLINATION SIN HARMONIC AMPL (RAD) C EPH(13)= INCLINATION AT REF TIME (RAD) C EPH(14)= RADIUS COS HARMONIC AMPL (M) C EPH(15)= ARGUMENT OF PERIGEE C EPH(16)= RATE OF RIGHT ASCENSION (RAD/SEC) C EPH(17)= SATELLITE I.D. C EPH(18)= CORRECTED MEAN MOTION (RAD/SEC) C O REC = RECORD DATA BUFFER C O CLKBUF = BUFFER OF CLOCK PARAMS OF ALL AVAIL SATS C O EPHBUF = BUFFER OF EPHEMERIS PARAMS OF ALL AVAIL SATS C I NSAT = NUMBER OF SATELLITES C I LU = INPUT DATA LOGICAL UNIT C I IPR = PRINT CODE FOR LINE PRINTER C C EXTERNALS NCLOK , DASET , NEPHM C C COMMENTS THE ROUTINE CURRENTLY READS SATELLITE EPHEMERIS DATA C ONLY C$ IPC = 0 N = 3 * NSAT C DO 900 I = 1,N READ ( LU , 1000 , END = 901 ) ( REC(K) , K = 1 , 16 ) 1000 FORMAT ( 4D20.12 ) C------------------------------------------------------------------- C IDENTIFY SATELLITE C------------------------------------------------------------------- C TEMP = REC(1) / 1D4 TEMP1 = DMOD ( TEMP , 1D2 ) ISV = IDINT( TEMP1 ) C------------------------------------------------------------------- C IDENTIFY RECORD TYPE C 0 - MEASUREMENT RECORD C 1 - HEADER INFORMATION C 2 - WEATHER DATA C 3 - EPHEMERIS DATA SUBFRAME 1 C 4 - EPHEMERIS DATA SUBFRAME 2 C 5 - EPHEMERIS DATA SUBFRAME 3 C 6 - MESSAGE SUBFRAME 4 C 7 - ALMANAC SUBFRAME 5 C------------------------------------------------------------------- TEMP = REC(1) / 1D6 ITYPE = IDINT ( TEMP ) IF ( ITYPE .EQ. 0 ) GO TO 10 GO TO (1,2,3,4,5,6,7) , ITYPE 1 CONTINUE 2 CONTINUE GO TO 11 3 CALL NCLOK(IPR,IPC,ISV,REC,CLK) IPSN = ISV*6 - 5 IT = 6 CALL DASET(CLK(1),CLKBUF(IPSN),IT) GO TO 11 C------------------------------------------------------------------- C FOR THIS TEST READ EPHEMERIS ONLY C------------------------------------------------------------------- 4 CONTINUE 5 CONTINUE ISUBFR = ITYPE - 2 CALL NEPHM(ISV,REC,ISUBFR,DTM,PI,EPH) IPSN = ISV*18 - 17 IT = 18 IF ( ITYPE .EQ. 5 ) CALL DASET ( EPH(1) , EPHBUF(IPSN) , IT ) GO TO 11 6 CONTINUE 7 CONTINUE 10 CONTINUE 11 CONTINUE 900 CONTINUE 901 CONTINUE RETURN END SUBROUTINE DASET(DIN,DOUT,LEN) IMPLICIT REAL * 8 ( A-H , O-Z ) INTEGER * 2 ITEMP DIMENSION DIN(1) , DOUT(1) C$ C NAME DASET C C PURPOSE COPY A VECTOR C C CALLING CALL DASET(DIN,DOUT,LEN) C C PARAMETERS I DIN(LEN) C O DOUT(LEN) C I LEN C C$ ITEMP = LEN DO 1 I = 1,ITEMP 1 DOUT(I) = DIN(I) RETURN END SUBROUTINE PLHXYZ(PHI,RLAM,H,XO,YO,ZO,A,B,X,Y,Z) IMPLICIT REAL * 8 ( A-Z ) C$ C NAME PLHXYZ C C PURPOSE COMPUTE CARTESIAN COORDINATES X , Y , Z C GIVEN ELLIPSOIDAL COORDINATES PHI , RLAM , H. C C CALLING CALL PLHXYZ(PHI,RLAM,H,XO,YO,ZO,A,B,X,Y,Z) C C PARAMETERS I PHI LATITUDE (RAD) C I RLAM EAST LONGITUDE (RAD) C I H HEIGHT (M) C I XO,YO,ZO DATUM TRANSLATION COMPONENTS (M) C I A,B ELLIPSOID SEMIMAJOR AND SEMIMINOR AXES (M) C O X,Y,Z CARTESIAN COORDINATES (M) C$ E2 = (A*A - B*B) / (A*A) SP = DSIN(PHI) CP = DCOS(PHI) N = A / DSQRT(1D0 - E2 * SP**2) X = XO + ( N + H ) * CP * DCOS(RLAM) Y = YO + ( N + H ) * CP * DSIN(RLAM) Z = ZO + ( N * (1D0 - E2) + H ) * SP RETURN END SUBROUTINE STXYZ ( AMYE , C , DTM , EPH , ICD , IFLG4 , LUI , $ PHIK , PI , RINC , RLNG , TIME , XR , YR , $ ZR, XS, YS, ZS, Y, IPC, LUO ) IMPLICIT REAL * 8 ( A-H , O-Z ) DIMENSION EPH(18) , DTM(7) , Y(ICD) , TS(3) , TR(3) C$ C NAME STXYZ C C PURPOSE COMPUTE SATELLITE COORDINATES C C CALLING CALL STXYZ ( AMYE , C , DTM , EPH , ICD , C IFLG4 , LUI , PHIK , PI , RINC , C RLNG , TIME , XR , YR , ZR , XS , C YS, ZS, Y, IPC, LUO ) C C PARAMETERS O AMYE ECCENTRIC ANOMALY (RAD) C I C SPEED OF LIGHT (M/SEC) C I DTM DATUM PARAMETERS C DTM(1) ELLIPSOID SEMIMAJOR AXIS (M) C DTM(2) ELLIPSOID ECCENTRICITY C DTM(3) X TRANSLATION COMPONENT (M) C DTM(4) Y TRANSLATION COMPONENT (M) C DTM(5) Z TRANSLATION COMPONENT (M) C DTM(6) GRAVITATIONAL CONSTANT (M**3/SEC**2) C DTM(7) EARTH ROTATION RATE (RAD/SEC) C I EPH SATELLITE EPHEMERIS PARAMETERS C EPH(1) AGE OF EPHEMERIS (SEC) C EPH(2) RADIUS SIN HARMONIC AMPL (M) C EPH(3) CORRECT COMPUTED MEAN MOTION (RAD/SEC) C EPH(4) MEAN ANOMALY AT REF TIME (RAD) C EPH(5) LATITUDE COS HARMONIC AMPL (RAD) C EPH(6) ORBITAL ECCENTRICITY C EPH(7) LATITUDE SIN HARMONIC AMPL (RAD) C EPH(8) ELLIPSE SEMIMAJOR AXIS (M) C EPH(9) EPHEMERIS REF TIME (SEC) C EPH(10)INCLINATION COS HARMONIC AMPL (RAD) C EPH(11) RIGHT ASCENSION AT REF TIME (RAD) C EPH(12) INCLINATION SIN HARMONIC AMPL (RAD) C EPH(13) INCLINATION AT REF TIME (RAD) C EPH(14) RADIUS COS HARMONIC AMPL (M) C EPH(15) ARGUMENT OF PERIGEE (RAD) C EPH(16) RATE OF RIGHT ASCENSION (RAD/SEC) C EPH(17) SATELLITE I.D. C EPH(18) CORRECTED MEAN MOTION (RAD/SEC) C I ICD COLUMN DIMENSION OF D MATRIX C O IFLG4 FLAG INDICATES PROCESSING C CONVERGENCE PROBLEMS C 0 ==> CONTINUE PROCESSING C 1 ==> ABORT PROCESSING C I LUI INPUT LOGICAL UNIT C O PHIK ARGUMENT OF LATITUDE (RAD) C I PI OBVIOUS C O RINC INCLINATION (RADIANS) C O RLNG LONGITUDE OF ASCENDING NODE (RAD) C I TIME GPS TIME MEASURED FROM REF EPOCH (SEC) C I XR RECEIVER WGS-72 X COORDINATE (M) C I YR RECEIVER WGS-72 Y COORDINATE (M) C I ZR RECEIVER WGS-72 Z COORDINATE (M) C O XS SATELLITE WGS-72 X COORDINATE (M) C O YS SATELLITE WGS-72 Y COORDINATE (M) C O ZS SATELLITE WGS-72 Z COORDINATE (M) C I Y VECTOR OF BIAS ESTIMATES C Y(2) ORBIT ELLIPSE SEMIMAJOR AXIS (M)OF C Y(3) ORBIT ECCENTRICITY C Y(4) ASC NODE RIGHT ASC AT REF TIME (RAD) C Y(5) ORBIT INCLINATION AT REF TIME (RAD) C Y(6) ARGUEMENT OF PERIGEE OF ELLIPSE (RAD) C Y(1) FOR DOPPLER = FREQUENCY OFFSET (HZ) C FOR RANGE = TIME OFFSET UNKNOWN (SEC) C I IPC PRINT CODE: C 1 = NORMAL C 2 = EXTENDED C 3 = SUPER EXTENDED C 4 = EVEN LONGER|| C I LUO LISTING LU C C EXTERNALS ANMLY , RANGE C$ IFLG4 = 0 C----------------------------------------------------------------------- C MEAN ANOMALY C----------------------------------------------------------------------- ANYM = EPH(4) + EPH(18) * TIME IF ( IPC .LT. 5 ) GO TO 100 WRITE(LUO,1000) WRITE(LUO,1001) ( I, EPH(I), I=1,18 ) WRITE(LUO,1002) ( I, Y(I), I=1,7 ) WRITE(LUO,1102)(I,DTM(I), I = 1,7) WRITE(LUO,1003) TIME C WRITE(LUO,1004) C WRITE(LUO,1005) ANYM, Y(1), AMYE, AMYT 100 CONTINUE C----------------------------------------------------------------------- C COMPUTE ECCENTRIC AND TRUE ANOMALIES C----------------------------------------------------------------------- CALL ANMLY ( ANYM , Y(2) , AMYE , AMYT , IFLG2 , PI , LUI ) C----------------------------------------------------------------------- C COMPUTE ARGUMENT OF LATITUDE C----------------------------------------------------------------------- PHIK = AMYT + Y(5) IF ( IPC .LT. 5 ) GO TO 200 WRITE(LUO,1010) WRITE(LUO,1005) ANYM, Y(2), AMYE, AMYT WRITE(LUO,1011) PHIK 200 CONTINUE C----------------------------------------------------------------------- C COMPUTE 2ND HARMONIC PERTURBATION CORRECTION C----------------------------------------------------------------------- CPK = DCOS ( 2.D0 * PHIK ) SPK = DSIN ( 2.D0 * PHIK ) DUK = EPH(7) * SPK + EPH(5) * CPK DRK = EPH(14) * CPK + EPH(2) * SPK DIK = EPH(10) * CPK + EPH(12) * SPK C----------------------------------------------------------------------- C CORRECT ARGUMENT OF LATITUDE C----------------------------------------------------------------------- PHIK = PHIK + DUK C----------------------------------------------------------------------- C CORRECT RADIUS C----------------------------------------------------------------------- RK = Y(1) * ( 1.D0 - Y(2) * DCOS(AMYE) ) + DRK C----------------------------------------------------------------------- C CORRECT INCLINATION C----------------------------------------------------------------------- RINC = Y(4) + DIK C----------------------------------------------------------------------- C POSITIONS IN ORBITAL PLANE C----------------------------------------------------------------------- XK = RK * DCOS(PHIK) YK = RK * DSIN(PHIK) C----------------------------------------------------------------------- C BEGIN ITERATIONS TO PROPERLY CORRECT COORDS FOR PROPAGATION TIME C----------------------------------------------------------------------- CINC = DCOS(RINC) SINC = DSIN(RINC) PROP= 0.D0 IF ( IPC .LT. 5 ) GO TO 300 WRITE(LUO,1020) CPK, SPK, DUK, DRK, DIK WRITE(LUO,1011) PHIK WRITE(LUO,1021) RK WRITE(LUO,1022) RINC, XK, YK WRITE(LUO,1023) CINC, SINC 300 CONTINUE DO 1 I=1,10 C----------------------------------------------------------------------- C SOLVE FOR CORRECTED LONGITUDE C----------------------------------------------------------------------- RLNG = Y(3) + ( EPH(16) - DTM(7) ) * TIME + DTM(7) * ((-1.D0) * & PROP - EPH(9) ) C----------------------------------------------------------------------- C COMPUTE EARTH FRAME COORDINATES OF SATELLITE C----------------------------------------------------------------------- CLNG = DCOS(RLNG) SLNG = DSIN(RLNG) XS = XK * CLNG - YK * CINC * SLNG YS = XK * SLNG + YK * CINC * CLNG ZS = YK * SINC C----------------------------------------------------------------------- C COMPUTE RANGE TO SATELLITE IN SECONDS C----------------------------------------------------------------------- TS(1) = XS TS(2) = YS TS(3) = ZS TR(1) = XR TR(2) = YR TR(3) = ZR C = 299792458.D0 RNG = RANGE ( TS , TR ) / C IF ( IPC .LT. 5 ) GO TO 400 WRITE(LUO,1030) RLNG, CLNG, SLNG WRITE(LUO,1031) XS, YS, ZS WRITE(LUO,1032) XR, YR, ZR WRITE(LUO,1033) RNG , C 400 CONTINUE IF ( DABS( PROP - RNG ) .LT. 1.D-12) GO TO 2 PROP = RNG 1 CONTINUE IF ( I .LT. 10 ) GO TO 2 IFLG4 = 1 2 CONTINUE IF ( IPC .EQ. 5 ) WRITE(LUO,1040) C----------------------------------------------------------------------- C FORMAT STATEMENTS C----------------------------------------------------------------------- 1000 FORMAT( '0', 'ENTERING STCRD (FORM)' ) 1001 FORMAT( ' EPH', I3, ' = ', F30.15 ) 1102 FORMAT( ' DTM', I3, ' = ', F30.15 ) 1002 FORMAT( ' Y ', I3, ' = ', F30.15 ) 1003 FORMAT( ' TIME = ', F30.15 ) 1004 FORMAT( ' BEFORE ANMLY ' ) 1005 FORMAT( ' ANYM, Y(1), AMYE, AMYT ', / 2F30.15 / 2F30.15 ) 1010 FORMAT( ' AFTER ANMLY ' ) 1011 FORMAT( ' PHIK = ', F30.15 ) 1020 FORMAT( ' CPK, SPK, DUK, DRK, DIK ' / 5D15.11 ) 1021 FORMAT( ' RK = ', D20.10 ) 1022 FORMAT( ' RINC, XK, YK ', 3D20.10 ) 1023 FORMAT( ' CINC, SINC ', 2F30.15 ) 1030 FORMAT( ' RLNG, CLNG, SLNG ', 3D20.10 ) 1031 FORMAT( ' XS, YS, ZS ', 3F20.5 ) 1032 FORMAT( ' XR, YR, ZR ', 3F20.3 ) 1033 FORMAT( ' RNG = ', F20.15 , ' C = ' , F20.5 ) 1040 FORMAT( ' LEAVING STCRD (FORM) ', // ) RETURN END SUBROUTINE DERIV(LU,AE,BE,XLAT,XLON,HGT,XS, $ S,DSDP,DSDL,ELEV,AZ,IER) IMPLICIT REAL * 8 ( A-H , O-Z ) DIMENSION ANG(3) , DRDL(3) , DRDP(3) , DX(3) , NAXIS(3) , $ ROT(3,3) , XR(3) , XS(1) , XT(3) DATA NAXIS /3,2,-2/, $ PI /3.14159265358979D0/ C$ C NAME DERIV C C PURPOSE COMPUTE SLANT RANGE(S) TO SATELLITE, ITS C DERIVATIVES WITH RESPECT TO LAT (DSDP), LONG (DSDL), C AND ELEVATION ANGLE (ELEV) C C CALLING CALL DERIV (LU,AE,BE,XLAT,XLON,HGT,XS, C S,DSDP,DSDL,ELEV,AZ,IER) C C PARAMETERS I LU LISTING LU C I AE EARTH SEMIMAJOR AXIS (M) C I BE EARTH SEMIMINOR AXIS (M) OR RECIP FLATNG C I XLAT RECEIVER LATITUDE (DEG) C I XLON RECEIVER LONGITUDE (DEG) C I HGT RECEIVER ELLIPSOID HEIGHT (M) C I XS SATELLITE LOCAL TERRESTRIAL COORDS (M) C O S RECEIVER TO SATELLITE SLANT RANGE (M) C O DSDP DERIV OF S WITH RESPECT TO LAT (M/DEG) C O DSDL DERIV OF S WITH RESPECT TO LONG (M/DEG) C O ELEV SATELLITE ELEVATION ABOVE HORIZON (DEG) C O AZ SATELLITE AZIMUTH FROM NORTH (DEG) C O IER 0 SUCCESSFUL RETURN C C EXTERNALS ROTREF C$ IER = 0 C---------------------------------------------------------------- C COMPUTE RECEIVER LOCAL TERRESTRIAL COORDS XR C---------------------------------------------------------------- CP = DCOS(XLAT*PI/180.D0) SP = DSIN(XLAT*PI/180.D0) CL = DCOS(XLON*PI/180.D0) SL = DSIN(XLON*PI/180.D0) BOA = (BE / AE)**2 IF(BE .LT. 6.D3) BOA = (1. - 1./BE)**2 RN = AE / DSQRT(CP**2 + BOA * SP**2) RM = RN * BOA * (RN/AE)**2 XR(1) = (RN+HGT) * CP * CL XR(2) = (RN+HGT) * CP * SL XR(3) = (RN*BOA + HGT) * SP C---------------------------------------------------------------- C COMPUTE DERIVS OF XR WRT TO LAT AND LONG C---------------------------------------------------------------- DRDP(1) = -(RM+HGT) * SP * CL * PI / 180.D0 DRDP(2) = -(RM+HGT) * SP * SL * PI / 180.D0 DRDP(3) = (RM+HGT) * CP * PI / 180.D0 DRDL(1) = -XR(2) * PI / 180.D0 DRDL(2) = XR(1) * PI / 180.D0 DRDL(3) = 0.D0 C---------------------------------------------------------------- C COMPUTE SLANT RANGE AND ITS DERIVATIVES C---------------------------------------------------------------- DO 10 I = 1,3 10 DX(I) = XS(I) - XR(I) S = DSQRT(DX(1)**2 + DX(2)**2 + DX(3)**2) DSDP = -(DX(1)*DRDP(1) + DX(2)*DRDP(2) + DX(3)*DRDP(3)) / S DSDL = -(DX(1)*DRDL(1) + DX(2)*DRDL(2) + DX(3)*DRDL(3)) / S C---------------------------------------------------------------- C COMPUTE ELEVATION AND AZIMUTH ANGLE C---------------------------------------------------------------- ANG(1) =(XLON - 180.)*PI/180.D0 ANG(2) =(XLAT - 90.)*PI/180.D0 ANG(3) = 0.D0 IT = 3 CALL ROTREF(IT,NAXIS,ANG,ROT) DO 20 I = 1,3 XT(I) = 0.D0 DO 20 J = 1,3 20 XT(I) = XT(I) + ROT(I,J) * DX(J) ELEV = DATAN2(XT(3),DSQRT(XT(1)**2 + XT(2)**2)) * 180.D0/PI AZ = DATAN2(XT(2),XT(1)) * 180.D0/PI RETURN END SUBROUTINE MATMPY(M1,M2,M3,L,M,N,JL,JM,JN,ICODE) REAL*8 M1(JL,1) , M2(JM,1) , M3(JN,1) C$ C NAME MATMPY C C PURPOSE TO COMPUTE THE PRODUCT OF TWO MATRICES IN ANY C ALLOWABLE TRANSPOSE COMBINATION AS FOLLOWS: C C CALLING CALL MATMPY(M1,M2,M3,L,M,N,JL,JM,JN,ICODE) C C OPTION ICODE PRODUCT M3 C C 1 M1*M2 C 2 (M1)T*M2 C 3 M1*(M2)T C 4 (M1)T*(M2)T C C C (L,M),(M,N),(L,N) ARE THE DIMENSIONS OF THE C PRE- , POST- AND PRODUCT-MATRICES C RESPECTIVELY C JL,JM,JN ARE CORRESPONDING DECLARED ROW C DIMENSIONS AT THE CALLINE PROGRAM C$ DO 11 I = 1,L DO 11 J = 1,N M3(I,J) = 0.D0 DO 11 K = 1,M GO TO(1,2,3,4),ICODE 1 CONTINUE C M3 = M1 * M2 M3(I,J) = M3(I,J) + M1(I,K)*M2(K,J) GO TO 11 2 CONTINUE C M3 = M1 TRANSPOSE * M2 M3(I,J) = M3(I,J) + M1(K,I)*M2(K,J) GO TO 11 3 CONTINUE C M3 = M1 * M2 TRANSPOSE M3(I,J) = M3(I,J) + M1(I,K)*M2(J,K) GO TO 11 4 CONTINUE C M3 = M1 TRANSPOSE * M2 TRANSPOSE M3(I,J) = M3(I,J) + M1(K,I)*M2(J,K) 11 CONTINUE RETURN END SUBROUTINE ANMLY(AMYM,ORBEC,AMYE,AMYT,IFLAG,PI,LU) IMPLICIT REAL*8(A-H,O-Z) C$ C NAME ANMLY C C PURPOSE CONVERT ECCENTRICITY AND MEAN ANOMALY TO ECCENTRIC AND C TRUE ANOMALY C C CALLING CALL ANMLY(AMYM,ORBEC,AMYE,AMYT,IFLAG,PI,LU) C C PARAMETERS I AMYM = MEAN ANOMALY (RAD) C I ORBEC= ORBITAL ECCENTRICITY C O AMYE = ECCENTRIC ANOMALY (RAD) C O AMYT = TRUE ANOMALY (RAD) C O IFLAG= 0 OK C 1 NONCONVERGENCE ON KEPLERS EQUATION FOR ECC ANOM C I PI C I LU = DIAGNOSTIC PRINTOUT LU C$ C IFLAG = 0 PI2 = PI * 2.0D0 ORBECS = ORBEC * ORBEC AMYM = DMOD( AMYM,PI2 ) C----------------------------------------------------------------- C COMPUTE FIRST APPROXIMATION TO ECCENTRIC ANOMALY C----------------------------------------------------------------- S = DSIN(AMYM) A = ORBEC * S B = A * ORBEC * DCOS(AMYM) C = ORBECS * A * (1.0D0 - 1.5D0 * S * S) AMYE = AMYM + A + B + C C ---------------------------------------------------------------------- C ITERATE ON ECCENTRIC ANOMALY C ---------------------------------------------------------------------- I = 0 10 I = I + 1 E1 = AMYE AMYE = E1+(AMYM-E1+ORBEC * DSIN(E1))/(1.0D0 - ORBEC * DCOS(E1)) IF ( I .LE. 8 ) GO TO 30 IFLAG = 1 GO TO 40 30 IF ( DABS(AMYE-E1) .GT. 1.0D-10 ) GO TO 10 C--------------------------------------------------------------------- C ITERATIONS COMPLETE C--------------------------------------------------------------------- 40 CONTINUE SE = DSIN(AMYE) DIFF = DABS(AMYM-(AMYE-ORBEC*SE)) IF ( DIFF .LT. 1.0D-10 ) GO TO 60 IFLAG = 1 C WRITE(LU,50) DIFF, I, AMYM, AMYE C 50 FORMAT(" /ANML2: ECCENTRICITY MISCLOSURE="D20.10/ C # " ITERATIONS="I4/" AMYM="D20.10" AMYE="D20.10) 60 CONTINUE C------------------------------------------------------------------- C COMPUTE TRUE ANOMALY C------------------------------------------------------------------- SF = DSQRT ( 1.0D0 - ORBECS ) * SE CF = DCOS ( AMYE ) - ORBEC AMYT = DATAN2 ( SF , CF ) RETURN END SUBROUTINE ROTREF(NUM,NAXIS,ANGLE,ROT) IMPLICIT REAL * 8 ( A-H , O-Z ) DIMENSION ROT(3,3) , R1(3,3) , R2(3) , $ ANGLE(NUM) , NAXIS(NUM) C$ C NAME ROTREF C C PURPOSE COMPUTE PRODUCT MATRIX OF SEQUENCE OF C ROTATIONS AND REFLECTIONS C C CALLING CALL ROTREF(NUM,NAXIS,ANGLE,ROT) C C PARAMETERS I NUM = NUMBER OF ROTATIONS AND REFLECTIONS IN SEQUENCE C I NAXIS = SEQUENCE OF ROTATION AND REFLECTION AXES C FOR ROTATIONS USE 1, 2, OR 3 C FOR REFLECTIONS USE -1, -2, OR -3 C I ANGLE = SEQUENCE OF ROTATION ANGLES IN RADIANS C IGNORE (ASSUME ZERO) FOR REFLECTIONS C O ROT = 3X3 PRODUCT MATRIX C$ EPS = 1D-15 C-------------------------------------------------------------------- C INITIALIZE (SET 'ROT' = IDENTITY MATRIX) C-------------------------------------------------------------------- DO 10 I=1,3 DO 10 J=1,3 ROT(I,J)=0.D0 IF(I.EQ.J) ROT(I,J)=1.D0 10 CONTINUE C-------------------------------------------------------------------- C PROCESS SEQUENCE OF ROTATIONS AND REFLECTIONS C-------------------------------------------------------------------- DO 100 N=1,NUM C-------------------------------------------------------------------- C DEFINE 3 AXES FOR CURRENT ROTATION OR REFLECTION C-------------------------------------------------------------------- N1=IABS(NAXIS(N)) N2=N1+1 IF(N2.GT.3) N2=N2-3 N3=N2+1 IF(N3.GT.3) N3=N3-3 C------------------------------------------------------------------- C DEFINE DIAGONAL ELEMENTS C------------------------------------------------------------------- R1(N1,N1) = 1.D0 IF(NAXIS(N).LT.0.D0) R1(N1,N1) = -1.D0 R1(N2,N2) = DCOS(ANGLE(N)) IF(NAXIS(N).LT.0.D0) R1(N2,N2) = 1.D0 R1(N3,N3) = R1(N2,N2) C------------------------------------------------------------------ C DEFINE NONZERO OFF-DIAGONAL ELEMENTS C------------------------------------------------------------------ R1(N2,N3) = DSIN(ANGLE(N)) IF(NAXIS(N).LT.0.D0) R1(N2,N3) = 0.D0 R1(N3,N2) = -R1(N2,N3) C----------------------------------------------------------------- C DEFINE ZERO OFF-DIAGONAL ELEMENTS C------------------------------------------------------------------ R1(N1,N2) = 0.D0 R1(N1,N3) = 0.D0 R1(N2,N1) = 0.D0 R1(N3,N1) = 0.D0 C------------------------------------------------------------------ C FORM PRODUCT (SET 'ROT' = 'R1' * 'ROT') C------------------------------------------------------------------ DO 100 J=1,3 DO 30 I=1,3 R2(I) = 0.D0 DO 30 K=1,3 30 R2(I) = R2(I) + R1(I,K)*ROT(K,J) C DO 100 I=1,3 ROT(I,J) = R2(I) IF(DABS(ROT(I,J)).LT.EPS) ROT(I,J) = 0.D0 100 CONTINUE RETURN END SUBROUTINE NCLOK ( LUO, IPC, ISAT, REC, CLK ) IMPLICIT REAL * 8 ( A-H , O-Z ) DIMENSION CLK(*) , REC(*) C$ C NAME NCLOK C C PURPOSE EXTRACT CLOCK VALUES FROM INPUT RECORD C C CALLING CALL NCLOK (LUO,IPC,ISAT,REC,CLK) C C PARAMETERS I LUO LISTING LU C I IPC PRINT CONTROL C I ISAT SATELLITE ID C I REC INPUT RECORD C O CLK OUTPUT ARRAY C C EXTERNAL DASET C$ C----------------------------------------------------------------------- C EXTRACT THE CLOCK VALUES C----------------------------------------------------------------------- CALL DASET ( REC(12), CLK, 5 ) CLK(6) = FLOAT(ISAT) C----------------------------------------------------------------------- C PRINT CLOCK VALUES IF REQUIRED C----------------------------------------------------------------------- IF ( IPC .GE. 3 ) WRITE(LUO,10100) (REC(J),J=3,11) RETURN C----------------------------------------------------------------------- C FORMAT STATEMENTS C----------------------------------------------------------------------- 10100 FORMAT( / ' CLOCK VALUES: ALPHA 0-3, BETA 0-3, TGD' / 1 2(4D20.10 /) , D20.10 ) END SUBROUTINE NEPHM ( ISAT, REC, J, DTM, PI, EPH ) IMPLICIT REAL * 8 ( A-H , O-Z ) DIMENSION EPH(*) , DTM(*) , REC(*) C$ C NAME NEPHM C C PURPOSE EXTRACT EPHEMERIS VALUES FROM INPUT RECORD C C CALLING CALL NEPHM(ISAT,REC,J,DTM,PI,EPH) C C PARAMETERS I ISAT SATELLITE ID. C I REC INPUT RECORD C I J SUBFRAME ID. C I DTM DATUM PARAMETERS C I PI OBVIOUS C O EPH OUTPUT ARRAY C C EXTERNALS DASET, DFLOT C$ C----------------------------------------------------------------------- C SEE IF FIRST OR SECOND HALF C----------------------------------------------------------------------- IF ( J .EQ. 3 ) GO TO 200 C----------------------------------------------------------------------- C FIRST HALF - GET VALUES C----------------------------------------------------------------------- IT = 9 CALL DASET ( REC(3), EPH , IT) C----------------------------------------------------------------------- C CONVERT TO BASE UNITS C----------------------------------------------------------------------- EPH(3) = EPH(3) * PI EPH(4) = EPH(4) * PI EPH(8) = EPH(8) * EPH(8) C----------------------------------------------------------------------- C COMPUTE CORRECTED MEAN POSITION C----------------------------------------------------------------------- EPH(18) = DSQRT(DTM(6)) / EPH(8)**1.5D0 + EPH(3) GO TO 300 C----------------------------------------------------------------------- C SECOND HALF - GET DATA AND C CONVERT TO BASE UNITS C----------------------------------------------------------------------- 200 CONTINUE IT = 7 CALL DASET ( REC(3), EPH(10), IT) EPH(11) = EPH(11) * PI EPH(13) = EPH(13) * PI EPH(15) = EPH(15) * PI EPH(16) = EPH(16) * PI EPH(17) = FLOAT(ISAT) 300 RETURN END DOUBLE PRECISION FUNCTION RANGE(XSV,XRV) IMPLICIT REAL * 8 ( A-H , O-Z ) DIMENSION XRV(3) , XSV(3) C$ C NAME RANGE C C PURPOSE COMPUTE STATION-TO-SATELLITE RANGE, GIVEN C STATION AND SATELLITE CARTESIAN COORDINATES C C CALLING X = RANGE(XSV,XRV) C C PARAMETERS I XRV(3) STATION CARTESIAN COORDINATES C I XSV(3) SATELLITE CARTESIAN COORDINATES C O RANGE STATION-TO-SATELLITE RANGE C$ R = 0D0 DO 10 I = 1 , 3 10 R = ( XSV(I) - XRV(I) ) ** 2 + R RANGE = DSQRT(R) RETURN END DOUBLE PRECISION FUNCTION ZFRDAY(ZCNT) C C PURPOSE CONVERT 'ZCNT' TIME INTO FRACTION OF DAY C IMPLICIT REAL*8(A-H,O-Z) C TSEC = ZCNT * 6. TIDD = TSEC / 86400.D0 ITIDD = IDINT(TIDD) ZFRDAY = TIDD - DBLE(ITIDD) RETURN END SUBROUTINE SPIN(Q,N,MM,DET,IDEXP) IMPLICIT REAL * 8 ( A-H , O-Z ) REAL*8 Q(MM,MM) C$ C NAME SPIN C C PURPOSE MATRIX INVERSION ROUTINE FOR SYMMETRIC POSITIVE-DEFINITE C MATRICES. THE MATRIX INVERTED IS THE UPPER N BY N C PORTION OF THE MATRIX Q (DIMENSIONED MM BY MM C IN CALLING ROUTINE). C C CALLING CALL SPIN (Q,N,MM,DET,IDEXP) C C PARAMETERS I/O Q = MATRIX TO BE INVERTED (DIMENSIONED MM BY MM) C ON OUTPUT, UPPER LEFT N BY N SUBMATRIX CONTAINS C INVERSE OF INPUT UPPER LEFT N BY N SUBMATRIX. C I N = DIMENSION OF THE ACTUAL PART (UPPER LEFT CORNER) C OF Q TO BE INVERTED (N .LE. MM) C I MM = DIMENSIONED SIZE OF Q IN THE CALLING ROUTINE. C O DET = NON-EXPONENT PORTION OF THE DETERMINANT OF THE C INPUT N BY N (UPPER LEFT CORNER) MATRIX. C O IDEXP = EXPONENT (OF 10) PART OF THE DETERMINANT. THUS C DETERMINANT = DET * 10 ** IDEXP . C C$ DET = 0.0D0 DO 4 J = 1,N DO 4 I = 1,J IF (I.EQ.1) GO TO 2 M = I - 1 SUM = 0.0D0 DO 1 K = 1,M 1 SUM = SUM + Q(K,I) * Q(K,J) Q(I,J) = Q(I,J) - SUM 2 IF (I.EQ.J) GO TO 3 Q(I,J) = Q(I,J) / Q(I,I) GO TO 4 3 CONTINUE DET = DET + DLOG10(Q(I,I)) Q(I,I) = DSQRT(Q(I,I)) 4 CONTINUE IDEXP = DET RPART = DET - IDEXP APART = DABS(RPART) IF (APART.LT.1.0D-20) DET = 1.D0 IF (APART.LT.1.0D-20) GO TO 10 DET = 10.0D0 ** RPART 10 CONTINUE DO 7 J = 1,N DO 7 I = 1,J IF (I.LT.J) GO TO 5 Q(J,J) = 1.0D0 / Q(J,J) GO TO 7 5 SUM = 0.0D0 M = J-1 DO 6 K = I,M 6 SUM = SUM - Q(I,K) * Q(K,J) Q(I,J) = SUM / Q(J,J) 7 CONTINUE DO 9 J = 1,N DO 9 I = 1,J SUM = 0.0D0 DO 8 K = J,N 8 SUM = SUM + Q(I,K) * Q(J,K) Q(I,J) = SUM IF (I.EQ.J) GO TO 9 Q(J,I) = SUM 9 CONTINUE RETURN END FUNCTION DJUL(J1,M1,T) C* CC NAME : DJUL (FUNCTION) CC CC XMJD= DJUL(J1,M1,T) CC CC PURPOSE : COMPUTES THE MODIFIED JULIAN DATE (MJD) CC FROM YEAR,MONTH AND DAY CC MJD = JULIAN DATE - 2400000.5 CC CC PARAMETERS : J1 : YEAR (E.G. 1984) I*4 CC M1 : MONTH(E.G. 2 FOR FEBRUARY) I*4 CC T : DAY OF MONTH R*8 CC CC AUTHOR : DJUL IS A MEMBER OF ASTLIB CC ASTRONOMICAL INSTITUTE , UNIVERSITY OF BERNE CC SWITZERLAND CC REAL*8 DJUL,T J=J1 M=M1 IF(M.GT.2)GO TO 1 J=J-1 M=M+12 1 I=J/100 K=2-I+I/4 DJUL=(365.25D0*J-DMOD(365.25D0*J,1.D0))-679006.D0 DJUL=DJUL+AINT(30.6001*(M+1))+T+K RETURN END FUNCTION THETA(T) C* CC NAME : THETA (FUNCTION) CC CC TH=THETA(T) CC CC PURPOSE : COMPUTATION OF MEAN SIDERAL TIME GREENWICH CC AT EPOCH T (MODIFIED JULIAN DATE) CC SEE "EXPLANATORY SUPPLEMENT TO THE EPHEMERIS" CC EDITION 1961 CC CC PARAMETER : T : EPOCH IN MODIFIED JULIAN DATE CC T AND THETA: REAL*8 CC CC AUTHOR : THETA IS A MEMBER OF ASTLIB CC ASTRONOMICAL INSTITUTE, UNIVERSITY OF BERNE CC SWITZERLAND CC REAL*8 T,TH,THETA,PI,ZPI PI=3.141592653589793D0 ZPI=2*PI TH=T-15019.5D0 THETA=279.6909833D0+360.9856473353D0*TH+.2902D-12*TH**2 THETA=DMOD(THETA/180*PI,ZPI) RETURN END //LKED.SYSLIB DD DSN=UNB1.FORTLIB,DISP=SHR // DD // DD // DD // DD DSN=A.M12129.SELIBOBJ,DISP=SHR //GO.FT14F001 DD DUMMY //GO.FT16F001 DD DUMMY //GO.FT17F001 DD DUMMY //GO.FT19F001 DD DUMMY //*GO.FT19F001 DD DSN=A.M1212B.QUE.ALERT,UNIT=SYSDA,DISP=(NEW,CATLG), //* SPACE=(TRK,(3,2),RLSE),DCB=(RECFM=FB,LRECL=150,BLKSIZE=15000) //GO.FT25F001 DD DUMMY //GO.SYSIN DD * 010401 45 26 34 283 44 42 50.0015. 1983 213 21 10 00 1983 214 2 10 00 05 5 4 1983 192.49547046 63.1997 167.5589 .0039040 214.2885 145.4408 2.00575775 717.980 180.002 6 1983 193.34519674 63.7477 45.8370 .0024403 119.4508 240.8805 2.00564376 717.958 179.997 8 1983 192.41335369 63.2219 167.5884 .0031790 335.1053 24.6996 2.00563908 717.943 179.993 5 1983 186.05356987 63.7280 168.3037 .0040514 227.8061 131.8517 2.00561920 717.967 179.998 9 1983 191.40597304 63.5613 46.0623 .0092091 73.5887 287.5003 2.00570750 717.977 180.001 //