//GPSTATS JOB /*JOBPARM M=4,R=2048,L=45 /*SERVICE -4 //STEP1 EXEC FORTVCLG,RC=2048K,RL=1024K,RG=1024K, // PARM.FORT='NOMAP,NOXREF' //FORT.SYSIN DD * IMPLICIT REAL*8(A-H,O-Z) REAL*4 EPHRT DIMENSION DTM(7), OBS(13,4), EPH(20,4), CLK(5,4), # DAT(4), ICOUN(12), IRECD(12,38), DEPH(20,4), # EPHRT(12,38), DCLK(5,4), TIMETR(4) C DIMENSION ISAT(4), TIM(4) DIMENSION TIME(25), DLTT(4), XRV(3) DIMENSION XSV(3,4), XSVD(3,4) DIMENSION AA(4,4), DD(4), DL(4) DIMENSION IWA1(4), IWA2(4) DIMENSION XMED(4), OMED(4), OBSU(4) C DIMENSION PSEUDO(4) , DX(3) , BIAS(4) , DELTAX(4) DIMENSION DXX(2000) , DYY(2000) , DZZ(2000), BU(2000) DIMENSION XUU(3) , DL1(2000) DIMENSION DL2(2000), DL3(2000), DL4(2000) DIMENSION G(2000) DIMENSION GDOPR(2000), XDOPR(2000), YDOPR(2000), ZDOPR(2000) C C--------------------------------------------------------------------- C CONSTANT VALUES C---------------------------------------------------------------------- IREAD = 5 IPRINT = 6 IDISK1 = 11 IDISK2 = 12 IPC = 1 MSAT = 12 MEPH = 38 LREC = 3000 NREC = 20 PI = 4.D0*DATAN(1.D0) F1 = 1.57542D+9 F2 = 1.2276D+9 C = 299792458.D0 C------------------------------------------- IDLAT = 45 IMLAT = 23 SLAT = 55.840D0 IDLOG =-75 IMLOG = 55 SLOG = 20.269D0 HGT = 45.64 RCHGT = 1.00 T = 273.D0 + 11.D0 P = 980.D0 E = 50.D0 C--------------------------------------------------------------------- C OPEN DISK FILES C--------------------------------------------------------------------- OPEN(IDISK1,ERR=3010,STATUS='OLD',ACCESS='SEQUENTIAL', # IOSTAT=IERR) OPEN(IDISK2,ERR=3020,STATUS='OLD',ACCESS='DIRECT', # RECL=3000,IOSTAT=IERR,FORM='FORMATTED') C------------------------------------------------------------------------ C DEFINE DATUM PARAMETERS C------------------------------------------------------------------------ CALL DATUM(DTM,0,IREAD) C----------------------------------------------------------------------- C DEFINE RECEIVER'S APPROXIMATE COORDS (METRES) C----------------------------------------------------------------------- DLAT = DMSRA(IDLAT,IMLAT,SLAT) DLOG = DMSRA(IDLOG,IMLOG,SLOG) A = DTM(1) B = A * DSQRT(1.D0 - DTM(2) * DTM(2)) CALL PLXYZ(DLAT,DLOG,HGT,0.D0,0.D0,0.D0,A,B,XRV(1),XRV(2),XRV(3)) WRITE(IPRINT,55) (XRV(I),I=1,3) 55 FORMAT(//,2X,'RECEIVER COORDINATES',G30.10,/) C---------------------------------------------------------------------- C EPHEM. TABLES FOR EASY AND FAST ACCESS C---------------------------------------------------------------------- C CALL ETABL(IPRINT,IDISK2,LREC,NREC,ICOUN,EPHRT,IRECD, # MSAT,MEPH) C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$4 C READ OBSERVATION RECORDS FROM K=2 TO K=2000 C------------------------------------------------------------------------ C ICHECK = 0 C DO 1 K=2,2000 READ(IDISK1,600,ERR=2010,IOSTAT=IERROR) (DAT(I),I=1,4), # ((OBS(I,J),J=1,4),I=1,13) C------------------------------------------------------------------------ C DEFINE SATELLITE IDENTITIES IN OBSRV. RECORD C------------------------------------------------------------------------ DO 2 J=1,4 ISAT(J) = OBS(1,J) IF(ISAT(J).EQ.0) GO TO 1 2 CONTINUE C---------------------------------------------------------------------- C COMPUTE TIME OF RECEPTION (TIMERC) C---------------------------------------------------------------------- TIMERC = RECEPT(DAT(1),DAT(2),DAT(3)) C--------------------------------------------------------------------- C COMPUTE TIME OF TRANSMISSION (IONOSPHERIC EFFECT) C-------------------------------------------------------------------- C CALL TRANSM(OBS,F1,F2,ISAT,IPRINT,TIMETR) CALL SET(TIM,4,TIMERC) C-------------------------------------------------------------------- C FIND BEST AVAILABLE BRODCAST EPHEMERIDES C-------------------------------------------------------------------- CALL FINDE(IPRINT,IDISK2,ISAT,TIM ,ICOUN,EPHRT,IRECD, # MSAT,MEPH,DTM,EPH,CLK) C------------------------------------------------------------------- C CORRECT CLOCK COEFFICIENTS FOR RELATIVISTIC EFFECTS C------------------------------------------------------------------ DO 7 J=1,4 IF(ISAT(J).EQ.0) GO TO 7 CALL CLKAN(CLK(1,J), EPH(1,J), PI, IPRINT, IFLG3) 7 CONTINUE C---------------------------------------------------------------------- C COMPUTE SATELITE COORDINATES AND SV-TIME OFFSET(DLTT) C---------------------------------------------------------------------- CALL SATOR(IPRINT,IPC,ISAT,TIMETR,EPH,CLK,DTM,XRV,XSV,XSVD, # DLTT) WRITE(IPRINT,10) K 10 FORMAT(5X,'NUMBER OF ITERATIONS=',I8,3X,'SATELLITE COORDINATES') C----------------------------------------------------------------------- C COMPUTE PSEUDO-RANGES C----------------------------------------------------------------------- C CALL PSRAN(IPRINT,ISAT,TIMETR,TIMERC,PSEUDO) C----------------------------------------------------------------------- C CORRECT FOR TROPOSHERIC DELAYS C----------------------------------------------------------------------- DO 11 I=1,4 IF(ISAT(I).EQ.0) GO TO 11 CALL ELEVT(XSV(1,I),XRV,DOT,CROSS,ELEV) C----------------------------------------------- C CHECK FOR MIN ELEVATION ANGLE C---------------------------------------------- CHECK = ELEV * 180.0 / PI IF(CHECK .LT. 5.0) GO TO 1 CALL HOPF(T,P,E,ELEV,HGT,RCHGT,PI,DELTA) PSEUDO(I) = PSEUDO(I) - DELTA 11 CONTINUE C------------------------------------------------ C COMPUTE NAVIGATIONAL SOLUTION C------------------------------------------------ ICHECK = ICHECK + 1 DO 12 I=1,3 XUU(I) = XRV(I) 12 CONTINUE CALL ADJUST(IPRINT,XSV,PSEUDO,XUU,BIASU,DELTAX,GDOP, # XDOP,YDOP,ZDOP,AA) C GDOPR(ICHECK) = GDOP XDOPR(ICHECK) = XDOP YDOPR(ICHECK) = YDOP ZDOPR(ICHECK) = ZDOP C 5 FORMAT(7X,'THE ICHECK VALUE=',I8,4X,'GDOP=',F10.3,2X,'K=',I6) C------------------------------------------------- C COMPUTE DIFFERENCES (DX,DY,DZ) FROM MEDIAN C------------------------------------------------- XMED(1) = XRV(1) + 16.834D0 - 8.045D0 XMED(2) = XRV(2) - 4.81D0 + 7.208D0 XMED(3) = XRV(3) + 17.281D0 - 7.794D0 BMED = 3993590.D0 +3.3888 C DXX(ICHECK) = XUU(1) - XMED(1) DYY(ICHECK) = XUU(2) - XMED(2) DZZ(ICHECK) = XUU(3) - XMED(3) BU( ICHECK) = BIASU C DO 555 JJ=1,4 DD(1) = DXX(ICHECK) DD(2) = DYY(ICHECK) DD(3) = DZZ(ICHECK) DD(4) = (BU(ICHECK) - BMED) 555 CONTINUE C DO 556 MD=1,4 DL(MD)= PSEUDO(MD) -(RANGE(XSV(1,MD),XMED)- BMED) 556 CONTINUE DL1(ICHECK) = DL(1) DL2(ICHECK) = DL(2) DL3(ICHECK) = DL(3) DL4(ICHECK) = DL(4) 1 CONTINUE C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C---------------------------------------------------------------------- C WRITE DEVIATIONS FROM RECEIVER'S COORDINATES C-------------------------------------------------------------------- WRITE(IPRINT,141) DO 14 I= 1, ICHECK WRITE(17,139) I, GDOPR(I), XDOPR(I), YDOPR(I), ZDOPR(I) C DO 14 I=1, ICHECK C WRITE(22 ,142) I, DXX(I), DYY(I), DZZ(I), BU(I),DL1(I), C # DL2(I), DL3(I), DL4(I) C WRITE(22,142) I, DXX(I), DYY(I), DZZ(I), BU(I),DL1(I), C # DL2(I), DL3(I), DL4(I) 14 CONTINUE C GO TO 2020 2010 WRITE(IPRINT,9040) IERROR GO TO 2020 3010 WRITE(IPRINT,9050) IERR GO TO 2020 3020 WRITE(IPRINT,9060) IERR GO TO 2020 C-------------------------------------------------------------------- C FORMAT STATEMENTS C------------------------------------------------------------------- 100 FORMAT(2X,'OBSERVATION RECORD=',I6) 111 FORMAT(120A) 139 FORMAT(2X,I8,4(F20.2,2X)) 140 FORMAT(2X,I8,2X,3(F10.3,2X),F20.2,F20.2,1X,3(F12.2)) 142 FORMAT(I8,1X,3(F10.3,2X),F20.4,F20.4,1X,3(F14.2)) 141 FORMAT('1',2X,'------I----------DX-----------DY-----------DZ--- #-----BIAS(METRES)----') 500 FORMAT(2X,'ICOUNT=',I8,'I=',I3) 501 FORMAT(2X,12(F10.2)) 502 FORMAT(2X,12(I5)) 600 FORMAT(4G30.15) 9040 FORMAT(' ','ERROR',I4,3X,'WHILE READING FILE OF OBSERVATIONS') 9050 FORMAT(' ','ERROR',I4,3X,'WHILE OPENINING FILE OF OBSERVATIONS') 9060 FORMAT(' ','ERROR',I4,3X,'WHILE OPENINING FILE OF EPHEMERIDES') 9070 FORMAT(' ','ERROR',I4,3X,'WHILE CLOSING FILE OF OBSERVATIONS') 9080 FORMAT(' ','ERROR',I4,3X,'WHILE CLOSING FILE OF EPHEMERIDES') C------------------------------------------------- C CLOSE FILES C------------------------------------------------- 2020 CLOSE(IDISK1,ERR=4010,STATUS='KEEP',IOSTAT=IERROR) CLOSE(IDISK2,ERR=4020,STATUS='KEEP',IOSTAT=IERROR) GO TO 2222 4010 WRITE(IPRINT,9070) IERROR GO TO 2222 4020 WRITE(IPRINT,9080) IERROR 2222 CONTINUE STOP END C SUBROUTINE ADJMED(IPRINT,XSV,XMED,BMED,OMED) C C FUNCTION : COMPUTES OBSERVATIONS BASED ON A GIVEN POSITION C C HISTORY : VERSION DATE PROGRAMMER C 1 MARCH, 86 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I LOG. UNIT OF PRINTER C XSV I ARRAY OF SATELLITE COORDINATES C XMED I VECTOR OF RECEIVER COORDINATES C BMED I BIAS OF RECEIVER (METRE) C OMED O VECTOR OF OBSERVATIONS C C EXTERNALS : RANGE , MMULD C------------------------------------------------------------- C IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION A(4,4) , DELTAX(4), U(4), V(4), W(4) DIMENSION XSV(3,4) , XMED(3) DIMENSION DOT(4) , RMED(4) DIMENSION OMED(4) DIMENSION C(4) DIMENSION ERROR(3) ITER = 0 C C--------------------------------------------------------- C COMPUTE DIRECTIONAL COSINES C--------------------------------------------------------- 1 CONTINUE ITER = ITER + 1 DO 2 I = 1,4 DIST = RANGE(XSV(1,I), XMED) U(I) = ( XSV(1,I) - XMED(1) ) / DIST V(I) = ( XSV(2,I) - XMED(2) ) / DIST W(I) = ( XSV(3,I) - XMED(3) ) / DIST 2 CONTINUE C-------------------------------------------------------- C GENERATE THE A DESIGN MATRIX C-------------------------------------------------------- DO 3 I=1,4 A(I,1) = U(I) A(I,2) = V(I) A(I,3) = W(I) A(I,4) = 1.D0 3 CONTINUE C------------------------------------------------------ C COMPUTE DOT PRODUCTS AND GENERATE OBSERVATION VECTOR C------------------------------------------------------ DO 4 I=1,4 DOT(I) = U(I) * XSV(1,I) + # V(I) * XSV(2,I) + # W(I) * XSV(3,I) 4 CONTINUE RMED(1) = XMED(1) RMED(2) = XMED(2) RMED(3) = XMED(3) RMED(4) = BMED C------------------------------------------------------ C COMPUTE OBSERVATION CORRESPONDING TO MEDIAN C------------------------------------------------------ CALL MMULD(C,4,A,4,RMED,4,4,4,1) C DO 5 I=1,4 OMED(I) = C(I) -DOT(I) 5 CONTINUE 999 RETURN END C SUBROUTINE ADJUST(IPRINT,XSV,PSEUDO,XU,BIASU,DELTAX,GDOP, # XDOP,YDOP,ZDOP,A) C C FUNCTION : COMPUTES NAVIGATION SOLUTION GIVEN 4 RANGES C C HISTORY : VERSION DATE PROGRAMMER C 1 MARCH, 86 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I LOG. UNIT OF PRINTER C XSV I ARRAY OF SATELLITE COORDINATES C BIAS I SV- TIME OFFSET (METRES) C XU I/O VECTOR OF APPR. RECEIVER COORDINATES C BIASU O TIME OFFSET OF RECEIVER (METRE) C A O A DESIGN MATRIX C DELTAX O VECTOR OF CORRECTIONS (METRE) C U O VECTOR OF DIRECTIONAL COSINES(1) C V O VECTOR OF DIRECTIONAL COSINES(2) C W O VECTOR OF DIRECTIONAL COSINES(3) C C I SPEED OF LIGHT C PL I COVARIANCE MATRIX OF OBSERVATIONS C C EXTERNALS : MINVD, RANGE , MMULD, MTXMD C------------------------------------------------------------- C IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION A(4,4) , DELTAX(4), U(4), V(4), W(4) DIMENSION XSV(3,4) , XU(3) DIMENSION DOT(4) , BIAS(4) , PSEUDO(4) DIMENSION OBSERV(4) , IW1(4) , IW2(4) DIMENSION B(4,4) , IWB1(4) , IWB2(4) DIMENSION ERROR(3) BIASU = 41514654.D0 ITER = 0 C C--------------------------------------------------------- C COMPUTE DIRECTIONAL COSINES C--------------------------------------------------------- 1 CONTINUE ITER = ITER + 1 DO 2 I = 1,4 DIST = RANGE(XSV(1,I), XU) U(I) = ( XSV(1,I) - XU(1) ) / DIST V(I) = ( XSV(2,I) - XU(2) ) / DIST W(I) = ( XSV(3,I) - XU(3) ) / DIST 2 CONTINUE C-------------------------------------------------------- C GENERATE THE A DESIGN MATRIX C-------------------------------------------------------- DO 3 I=1,4 A(I,1) = U(I) A(I,2) = V(I) A(I,3) = W(I) A(I,4) = 1.D0 3 CONTINUE C------------------------------------------------------- C COMPUTE GDOP C------------------------------------------------------- CALL MTXMD(B,4,A,4,4,4) CALL MINVD(B,4,4,DETB,IWB1,IWB2) GDOP = DSQRT( B(1,1) + B(2,2) + B(3,3) + B(4,4) ) XDOP = DSQRT( B(1,1) ) YDOP = DSQRT( B(2,2) ) ZDOP = DSQRT( B(3,3) ) C------------------------------------------------------ C COMPUTE DOT PRODUCTS AND GENERATE OBSERVATION VECTOR C------------------------------------------------------ DO 4 I=1,4 DOT(I) = U(I) * XSV(1,I) + # V(I) * XSV(2,I) + # W(I) * XSV(3,I) 4 CONTINUE DO 5 I=1,4 OBSERV(I) = DOT(I) - PSEUDO(I) 5 CONTINUE C------------------------------------------------------ C COMPUTE SOLUTION DELTAX = (A-INVERSE) * L C------------------------------------------------------ CALL MINVD(A,4,4,DETA,IW1,IW2) C IF(DABS(DETA).LT.001D0) WRITE(IPRINT,200) DETA CALL MMULD(DELTAX,4,A,4,OBSERV,4,4,4,1) C C WRITE(IPRINT,90) ITER C WRITE(IPRINT,100) (DELTAX(I),I=1,4) C DO 6 I=1,3 ERROR(I) = XU(I) - DELTAX(I) XU(I) = DELTAX(I) 6 CONTINUE C DO 7 I=1,3 IF(DABS(ERROR(I)).GT.0.001D0) GO TO 1 7 CONTINUE C BIASU = DELTAX(4) 90 FORMAT(2X,'-------------------I T E R A T I O N#',I8,'FROM ADJU---- #-------------------') 100 FORMAT(2X,'DELTAX(I)=',G30.15) 200 FORMAT(2X,'A-MATRIX PROBABLY ILL-CONDITIONED WITH DETA=',G30.7) RETURN END C SUBROUTINE HOPF(T,P,E,ELEV,HGT,RCHGT,PI,DELTA) IMPLICIT REAL * 8 (A-H,O-Z) C FUNCTION : ORRECT RANGES FOR TROPOSHERIC REFRACTION C C HISTORY : VERSION DATE PROGRAMMER C 1 1980 U.N.B. C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C T I ATMOSPH. TEMPERATURE (KELVIN DEGREES) C P I ATMOSPHERIC PRESSURE (MBARS) C E I WATER VAPOUR PRESSURE (MBARS) C ELEV I ELEVATION ANGLE (RADIANS) C HGT I HEIGHT ABOVE REFER. ELLIPSOID(M) C RCHGT I RECEIVER HGT ABOVE SURVEY MARK(M) C PI I 3.14159..... C DELTA O TROPOSHERIC DELAY (METRES) C C EXTERNALS : NONE C------------------------------------------------------------- C HGTT = HGT + RCHGT DEG = 180D0 / PI DEND = DSQRT ( ( ELEV * DEG ) ** 2 + 6.25D0 ) / DEG DENW = DSQRT ( ( ELEV * DEG ) ** 2 + 2.25D0 ) / DEG HW = 11000D0 HD = 148.72D0 * T - 488.3552D0 RKD = 1.552D-05 * ( P / T ) * ( HD - HGTT) RKW = 7.46512D-02 * ( E / T ** 2 ) * ( HW - HGTT) DELTA= ( RKD / DSIN ( DEND ) ) + RKW / DSIN( DENW ) RETURN END C SUBROUTINE ELEVT(XSV,XSTAT,DOT,CROSS,ELEV) C------------------------------------------------------------ C FUNCTION : COMPUTES THE ELEVATION ANGLE OF A SATELLITE C C HISTORY : VERSION DATE PROGRAMMER C 1 MARCH 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C XSV I SATELLITE COORDINATES (METRE) C XSTAT I STATION (RECEIVER) COORDINATES(M) C DOT O DOT PRODUCT OF (XSTAT) WITH (DX) C CROSS O CROSS PRODUCT (XSTAT) WITH (DX) C ELEV O ELEVATION ANGLE (RADIANS) C C EXTERNALS : NONE C------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) DIMENSION XSV(3) , XSTAT(3), DX(3) PI = 4.D0 * DATAN(1.D0) C----------------------------------------- C COMPUTE RANGE VECTOR DX C----------------------------------------- DOT = 0.D0 DO 1 I=1,3 DX(I) = XSV(I) - XSTAT(I) 1 CONTINUE C----------------------------------------- C COMPUTE DOT PRODUCT C----------------------------------------- DO 2 I=1,3 DOT = XSTAT(I) * DX(I) + DOT 2 CONTINUE C---------------------------------------- C COMPUTE CROSS PRODUCT OF XSTAT AND DX C---------------------------------------- A1 = XSTAT(1) B1 = XSTAT(2) C1 = XSTAT(3) C A2 = DX(1) B2 = DX(2) C2 = DX(3) C CROSS = (B1 * C2 - B2 *C1) **2 + # (C1 * A2 - C2 *A1) **2 + # (A1 * B2 - A2 *B1) **2 CROSS = DSQRT (CROSS) C----------------------------------------- C COMPUTE ELEVATION ANGLE C----------------------------------------- ELEV = DATAN( - DOT / CROSS ) IF(ELEV.LT.0.D0) ELEV = DABS(ELEV) RETURN END SUBROUTINE PSRAN(IPRINT,ISAT,TIMETR,TIMERX,PSEUDO) C------------------------------------------------------------ C FUNCTION : COMPUTES PSEUDORANGES FROM TIME READINGS C C HISTORY : VERSION DATE PROGRAMMER C 1 JULY 1984 AL. KLEUSBERG C 2 MARCH 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I L. UNIT OF PRINTER. C C I SPEED OF LIGHT C ISAT I SATELLITES' IDENTITIES C TIMETR I TIME OF SV SIGNAL TRANSMISSION C TIMERX I TIME OF SIGNAL RECEPTION OF RECEIVER C PSEUDO O PSEUDO-RANGE (METRES) C C EXTERNALS : NONE C------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) DIMENSION ISAT(4), TIMETR(4), PSEUDO(4) C = 299792458.D0 C DO 100 I=1,4 IF(ISAT(I).EQ.0) GO TO 100 PSEUDO(I) = (TIMERX - TIMETR(I)) * C 100 CONTINUE RETURN END C FUNCTION DMSRA(ID,IM,S) IMPLICIT REAL*8(A-H,O-Z) C RHO = 45.D0 / DATAN(1.D0) D = DBLE(ID) SIG = DSIGN(1.D0,D) C DMSRA = (DABS(D) + DBLE(IM) / 60.D0 + S/3600.D0) * SIG / RHO RETURN END C SUBROUTINE PLXYZ(PHI,RLAM,H,X0,Y0,Z0,A,B,X,Y,Z) IMPLICIT REAL*8(A-H,O-Z) C E2 = (A*A -B*B) / (A*A) SP = DSIN(PHI) CP = DCOS(PHI) N = A / DSQRT(1D0 - E2*SP**2) X = X0 + (N + H) * CP * DCOS(RLAM) Y = Y0 + (N + H) * CP * DSIN(RLAM) Z = Z0 + (N * (1D0 -E2) + H ) * SP RETURN END C SUBROUTINE KLORB(EPH,XSV,XSVD,C,XRV,DTM,TIM,ITIM,IPRINT,IVEL) C------------------------------------------------------------ C FUNCTION : COMPUTES SATELLITE COORDINATES AND VELOCITY C IN EARTH-FIXED, EARTH-CENTERED SYSTEM DEFINED C AT THE EPOCH OF SIGNAL TRANSMISSION C C HISTORY : VERSION DATE PROGRAMMER C 1 1981 U.N.B.(DELIKARAOGLOU) C 2 JUNE 1984 AL. KLEUSBERG C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I L. UNIT OF PRINTER. C EPH I SATELLITE EPHEMERIS C XSV O SATELLITE COORDINATES (METRES) C XSVD O SATELLITE VELOCITIES (M/SEC) C C I SPEED OF LIGHT C XRV I RECEIVER'S COORDINATES (METRES) C DTM I DATUM PARAMETRES C TIM I TIME SINCE EPH REFERENCE TIME C ITIM I TIME FLAG = 0 TIME OF SIGNAL TRANS C TIME FLAG > 0 TIME OF SIGNAL RECEPT C IVEL I COMPUTATION CODE C = 0 COMPUTE POSITION COORDINATES C > 0 COMPUTE COORDINATES & VELOCITIES C EXTERNALS : ANMLY, RANGE C------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) DIMENSION EPH(20), XSVD(3), XSV(3), DTM(7) , XRV(3) DIMENSION XR(3) C A = EPH(8) ECC = EPH(6) RINC= EPH(13) CALL ANMLY(EPH,DTM,TIM,IPRINT,AMYM,AMYE,AMYT,PHIK0) RADIUS = A * (1.D0 - ECC * DCOS(AMYE)) C------------------------------------------------- C TIME DERIV OF RADIUS AND TRUE ANOMALY C------------------------------------------------- RDOT = A*ECC*DSIN(AMYE)*EPH(18)/(1.D0-ECC*DCOS(AMYE)) VDOT = DSQRT(1-ECC*ECC)*EPH(18)/(1.D0-ECC*DCOS(AMYE))**2 C------------------------------------------------ C COMPUTE PERTUBATION FROM EPHEMERIS CORRECTION PARAMETERS C------------------------------------------------ CPK = DCOS(2.D0 * PHIK0) SPK = DSIN(2.D0 * PHIK0) C DUK = EPH(5) * CPK + EPH(7) * SPK DRK = EPH(14) * CPK + EPH(2) * SPK DIK = EPH(10) * CPK + EPH(12) * SPK C------------------------------------------------- C CORRECT ARGUMENT OF LATITUDE, RADIUS, INCLINATION C------------------------------------------------- RADIUS = RADIUS + DRK RINC = RINC + DIK + EPH(19)* TIM ITER = 0 PROP = 0.D0 10 PHIK = PHIK0 + DUK - PROP*VDOT COSU = DCOS(PHIK) SINU = DSIN(PHIK) COSI = DCOS(RINC) SINI = DSIN(RINC) C RLNG = EPH(11) + (EPH(16) - DTM(7) )*(TIM - PROP) # - DTM(7) * EPH(9) COSW = DCOS(RLNG) SINW = DSIN(RLNG) C---------------------------------------------- C COORDINATES IN ORBITAL PLANE C---------------------------------------------- X = RADIUS*COSU Y = RADIUS*SINU C---------------------------------------------- C COORDINATES IN EARTH FIXED SYSTEM C---------------------------------------------- XSV(1) = X*COSW - Y*COSI*SINW XSV(2) = X*SINW + Y*COSI*COSW XSV(3) = Y*SINI C---------------------------------------------- C EXIT IF TIME IS TIME OF SIGNAL TRANSMISSION C---------------------------------------------- IF(ITIM.EQ.0) GO TO 30 DR = PROP*DTM(7) SINDR = DSIN(DR) COSDR = DCOS(DR) XR(1) = XRV(1)*COSDR + XRV(2)*SINDR XR(2) = XRV(2)*COSDR - XRV(1)*SINDR XR(3) = XRV(3) RNG = RANGE (XSV,XR) TPROP = RNG / C IF(DABS(PROP - TPROP) .LT. 1.D-12) GO TO 30 IF(ITER .GE. 10 ) GO TO 20 PROP = TPROP ITER = ITER + 1 GO TO 10 20 CONTINUE WRITE(IPRINT,1001) DABS(PROP - TPROP), ITER 1001 FORMAT(2X,'KLORB:ABS(PROP-TPROP) =',G20.8,'AFTER',I4, # 'ITERATIONS') 30 CONTINUE IF(IVEL .EQ. 0) RETURN C------------------------------------------------------- C TIME OF ORBITAL ELEMENTS C------------------------------------------------------- DUDOT = 2.D0*VDOT*(-EPH(5) * SPK + EPH(7) * CPK) DRDOT = 2.D0*VDOT*(-EPH(14) * SPK + EPH(2) * CPK) DIDOT = 2.D0*VDOT*(-EPH(10) * SPK + EPH(12) * CPK) DUDOT = DUDOT + VDOT DRDOT = DRDOT + RDOT DIDOT = DIDOT + EPH(19) OMDOT = EPH(16) -DTM(7) C------------------------------------------------------- C DERIVATIVES IN ORBITAL PLANE C------------------------------------------------------- XDOT = DRDOT*COSU - RADIUS*SINU*DUDOT YDOT = DRDOT*SINU + RADIUS*COSU*DUDOT C------------------------------------------------------ C VELOCITIES IN EARTH FIXED SYSTEM C------------------------------------------------------ FAC1 = XDOT - Y*COSI*OMDOT FAC2 = -X*OMDOT - YDOT*COSI + Y*SINI*DIDOT XSVD(1) = COSW*FAC1 + SINW*FAC2 XSVD(2) = SINW*FAC1 - COSW*FAC2 XSVD(3) = YDOT*SINI + Y*COSI*DIDOT RETURN END C SUBROUTINE ANMLY(EPH,DTM,TIME,IPRINT,AMYM,AMYE,AMYT,PHIK) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION EPH(20), DTM(7) PI = 3.141592653589793D0 ITER = 0 ECC = EPH(6) C----------------------------------------------------------------- C SOLVE KEPLER'S LAW FOR THE ECCENTRIC ANOMALY C----------------------------------------------------------------- AMYM = EPH(4) + EPH(18) * TIME AMYM = DMOD(AMYM, 2.D0* PI) AMYE0= AMYM + ECC * DSIN(AMYM) + ECC*2 * DSIN(2.D0*AMYM)/2.D0 5 CONTINUE AMYM1 = AMYE0 - ECC * DSIN(AMYE0) DE =(AMYM - AMYM1) / (1.D0 - ECC*DCOS(AMYE0)) AMYE0 = AMYE0 + DE ITER = ITER + 1 IF(ITER .LT. 10 .AND. DABS(DE) .GT. 1.D-10) GO TO 5 IF(DABS(DE) .GT. 1.D-10) WRITE(IPRINT,1000) DABS(DE) C AMYE = AMYE0 C------------------------------------------------------------- C COMPUTE TRUE ANOMALY C------------------------------------------------------------ COSV = DCOS(AMYE) - ECC SINV = DSQRT(1.D0 - ECC**2) * DSIN(AMYE) AMYT = DATAN2(SINV,COSV) C PHIK = AMYT + EPH(15) 1000 FORMAT(//,5X,'NO CONVERGENCE IN ANMLY AFTER 10', # 'ITERATIONS DABS(DE) =',G20.8) RETURN END C SUBROUTINE CKCOR(CLK ,EPH, IOBS, NMOBS, NSTN, TIME, PI, # IPRINT,IFLG2, ICKOP, OBS, C, IPC, IDPRG,NTRVL,DLTSI) C------------------------------------------------------------ C FUNCTION : CONVERT SV-TIME TO GPS-TIME, C AND APPLY CLOCK CORRECTION PARAMETERS C TRANSMITTED FROM SATELLITE C C HISTORY : VERSION DATE PROGRAMMER C 1 1981 U.N.B. C 2 JUNE 1984 AL. KLEUSBERG C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I L. UNIT OF PRINTER. C IPC I PRINTING CODE(=1) C EPH I SATELLITE EPHEMERIS C CLK I CORRECTED CLOCK PARAMETERS C IOBS I NUMBER OF OBSERVATIONS/ SEGMENT C NMOBS I NUMBER OF OBSERVATIONS/ SATELLITE C NSTN I NUMBER OF STATIONS SIMULT. OBSERVING C ICKOP I CLOCK MODEL OPTION C -1 --> IMPROVED MODEL C 0 --> STANDARD MODEL C C I SPEED OF LIGHT C IDPRG I =0 ; WHEN OBS=DOPPLER C =1 ; WHEN OBS=PSEUDO RANGES C =2 ; WHEN OBS=P-CODE READINGS C NTRVL I DOPPLER INTERVAL C TIME O VECTOR OF TIMES SINCE EPHEMERIS EPOCH C CORRESPONDING TO EITHER RANGES OR TO C THE BEGINNING OF DOPPLER INTERVAL C IFLG2 O CONVERGENCE INDICATOR FROM ANML2 C = 0 WHEN PROPER CONVERGENCE C = 1 WHEN IMPROPER CONVERGENCE C OBS I/O OBSERV TO BE CORRECTED (DOPPLER, C PSEUDO RANGE, CODE) C DLTSI O SV TIME OFFSET C EXTERNALS : ANML2 C------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) DIMENSION TIME(IOBS,NSTN), CLK(5), EPH(20), OBS(1) LOGICAL RESET, FIRST IFLG2 = 0 FIRST = .TRUE. C------------------------------------------------------------- C IMPROVED CLOCK MODEL C------------------------------------------------------------- DO 600 I=1,NMOBS RESET = .FALSE. PTIME = TIME(I,1) C------------------------------------------------------------- C CORRECT FOR END OF WEEK CROSSOVERS C------------------------------------------------------------- 5 CONTINUE DLTR = 0.0D0 TKE = PTIME - EPH(9) IF( TKE. LE. 3024.D2 ) GO TO 10 TKE = TKE -6048.D2 GO TO 20 10 CONTINUE IF( TKE .LT. -3024.D2 ) TKE = TKE + 6048.D2 20 CONTINUE TKC = PTIME - CLK(2) IF( TKC .LE. 3024.D2 ) GO TO 30 TKC = TKC - 6048.D2 GO TO 40 30 CONTINUE IF( TKC. LT. -3024.D2 ) TKC = TKC + 6048.D2 40 CONTINUE IF( ICKOP. EQ. 0) GO TO 50 C--------------------------------------------------------- C COMPUTE MEAN ANOMALY C---------------------------------------------------------- ANYM = EPH(4) + EPH(18) * TKE C--------------------------------------------------------- C COMPUTE ECCENTRIC ANOMALY C--------------------------------------------------------- CALL ANML2 ( ANYM, EPH(6), AMYE, AMYT, IFLG, PI, IPRINT) IF( IFLG. EQ. 1 ) IFLG2 = 1 C--------------------------------------------------------- C COMPUTE IMPROVED CORRECTION TERM C--------------------------------------------------------- RK = -4.443D-10 * EPH(6) * DSQRT(EPH(8)) DLTR= RK * DSIN(AMYE) 50 CONTINUE C-------------------------------------------------------- C COMPUTE CORRECTED TIME FROM EPHEMERIS EPOCH C-------------------------------------------------------- ODLTSI = DLTSI DLTSI = CLK(5) + CLK(4) * TKC + CLK(3) * TKC * TKC + DLTR C-------------------------------------------------------- C SEE IF DOPPLER OR RANGE C-------------------------------------------------------- IF( IDPRG. GT. 0 ) GO TO 500 J = I - 1 IF( .NOT. RESET ) TIME (I,1) = TKE - DLTSI IF(FIRST .AND. .NOT. RESET ) GO TO 100 C------------------------------------------------------- C APPLY DLTSI TO DOPPLER C------------------------------------------------------- IF( RESET ) J = J + 1 TC = ( DLTSI - ODLTSI )/ (FLOAT(NTRVL) + DLTSI -ODLTSI) OBS(J) = OBS(J) * TC + OBS(J) IF (IPC .GE. 2 ) WRITE(IPRINT,10060) J, OBS(J), DLTSI, ODLTSI, # TC C------------------------------------------------------ C SEE IF WE JUST RESET AFTER BREAK C------------------------------------------------------ 100 CONTINUE IF(RESET) GO TO 600 FIRST = .FALSE. C------------------------------------------------------ C CHECK FOR LAST OBSERVATION C------------------------------------------------------ IF( I .EQ. NMOBS) GO TO 200 C------------------------------------------------------ C SEE IF BREAK COMING UP C------------------------------------------------------ IF( (TIME(I+1,1)- PTIME) .LE. FLOAT(NTRVL) ) GO TO 600 C------------------------------------------------------ C BREAK RESET PARAMETERS C----------------------------------------------------- 200 CONTINUE PTIME = PTIME + FLOAT(NTRVL) RESET = .TRUE. FIRST = .TRUE. GO TO 5 C----------------------------------------------------- C CORRECT RANGES AND P-CODE READINGS C----------------------------------------------------- 500 CONTINUE TIME(I,1) = TKE - DLTSI IF( IDPRG. EQ. 1) OBS(I) = OBS(I) + DLTSI * C IF( IDPRG. GT. 1) OBS(I) = OBS(I) - DLTSI IF( IPC. GE. 2) WRITE(IPRINT,10500) I, OBS(I), DLTSI 600 CONTINUE RETURN C----------------------------------------------------- C FORMAT STATEMENTS C----------------------------------------------------- 10060 FORMAT(I6,4D20.10) 10500 FORMAT(I6,D20.10, D20.10) END C SUBROUTINE SATOR(IPRINT,IPC,ISAT,TIME,EPH,CLK,DTM,XRV,XSV,XSVD, # DLTT) C------------------------------------------------------------ C FUNCTION : CONVERT SV-TIME TO GPS-TIME, C COMPUTE SATELLITE COORDINATES AND C VELOCITIES IN EARTH-FIXED EARTH-CENTERED SYSTEM C C HISTORY : VERSION DATE PROGRAMMER C 1 JUNE 1984 AL. KLEUSBERG C 2 FEB 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I L. UNIT OF PRINTER. C IPC I PRINTING CODE C ISAT I SATELLITE IDENTITIES C TIME I/O SV-TIMES OF SIGNAL TRANSMISSION C EPH I SATELLITE EPHEMERIS C CLK I CORRECTED CLOCK PARAMETERS C DTM I DATUM PARAMETERS C XRV I RECEIVER COORDINATES C XSV O SATELLITES' COORDINATES C XSVD O SATELLITES' VELOCITIES C DLTT O SV-TIME OFFSETS. C EXTERNALS : CKCOR , KLORB, DMOVE C------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) C DIMENSION ISAT(4), TIME(4), DLTT(4) DIMENSION XRV(3) , EPH(20,4), CLK(5,4) DIMENSION XSV(3,4) , XSVD(3,4) DIMENSION DTM(7) , TIM(4) PI = 4.D0 * DATAN(1.D0) C DATA MONE, NULL, IONE, ITWO / -1, 0, 1, 2/ DATA C /299792458.D0/ DATA ZERO /0.0D0 / C DO 130 I=1,4 IF(ISAT(I).NE.0) GO TO 140 130 CONTINUE RETURN 140 CALL DMOVE (TIME,4,TIM,1,4) C DO 100 I=1,4 IF(ISAT(I) .EQ. 0) GO TO 110 C------------------------------------------------ C CONVERT SV-TIME TO GPS-TIME C------------------------------------------------ CALL CKCOR(CLK(1,I),EPH(1,I),IONE,IONE,IONE,TIM(I),PI,IPRINT, # IFLG3,MONE,TIME(I),C,IPC,ITWO,NULL,DLTT(I)) C------------------------------------------------ C COMPUTE SATELLITE COORDINATES AND VELOCITIES C------------------------------------------------ CALL KLORB(EPH(1,I),XSV(1,I),XSVD(1,I),C,XRV,DTM,TIM(I),NULL, # IPRINT,NULL) GO TO 100 110 CALL SET(XSV(1,I) ,3,ZERO) CALL SET(XSVD(1,I),3,ZERO) C 100 CONTINUE RETURN END C SUBROUTINE CLKAN(CLK, EPH, PI, IPRINT, IFLG3) C------------------------------------------------------------ C FUNCTION : CORRECT CLOCK COEFFICIENTS FOR C RELATIVISTIC EFFECTS. C C HISTORY : VERSION DATE PROGRAMMER C 1 1981 U.N.B.? C 2 JUNE 1984 AL. KLEUSBERG C 3 FEB 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I L. UNIT OF PRINTER. C CLK(1) I/O AGE OF CLOCK DATA (SEC) C CLK(2) I/O CLOCK REFERENCE TIME (SEC) C CLK(3) I/O A2 COEFFICIENT (SEC/SEC**2) C CLK(4) I/O A1 COEFFICIENT (SEC/SEC) C CLK(5) I/O A0 COEFFICIENT (SEC) C PI I GREEK PI=3.14159... C EPH(1) I AGE OF EPHEMERIS C EPH(2) I RADIUS SIN HARMONIC AMPL.(METRE) C EPH(3) I CORRECT COMPUTED MEAN MOTION(RAD/SEC) C EPH(4) I MEAN ANOMALY AT REF TIME (RAD) C EPH(5) I ARG LATITUDE COS HARMONIC AMPL(RAD) C EPH(6) I ECCENTRICITY C EPH(7) I ARG LATITUDE SIN HARMONIC AMPL(RAD) C EPH(8) I SEMI-MAJOR AXIS(METRE) C EPH(9) I REFERENCE TIME= T0(SEC) C EPH(10) I INCLINATION COS HARM AMPL (RAD) C EPH(11) I RIGHT ASCENSION AT REF TIME (RAD) C EPH(12) I INCLINATION SIN HARM AMPL (RAD) C EPH(13) I INCLINATION AT REF TIME(RAD) C EPH(14) I RADIUS COS HARM AMPL (METRE) C EPH(15) I ARGUMENT OF PERIGEE(RAD) C EPH(16) I RATE OF RIGHT ASCENSION (RAD/SEC) C EPH(17) I AODE C EPH(18) I CORRECTED MEAN MOTION(RAD/SEC) C EPH(19) I RATE OF INCLINATION(RAD/SEC) C EPH(20) I SATELLITE ID C IFLAG3 O MISCOSURE FLAG FROM ANML2 C C EXTERNALS : ANML2 C------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) DIMENSION CLK(5), EPH(20) C------------------------------------------------------------- C CORRECT COEFFICIENTS FOR RELATIVISTIC EFFECTS C COMPUTE ECCENTRIC ANOMALY AT REFERENCE TIME C------------------------------------------------------------- CALL ANML2(EPH(4),EPH(6),AMYE,AMYT,IFLG3,PI,IPRINT) ASQRT = DSQRT (EPH(8)) SINE = DSIN (AMYE) RK = -4.443D-10 * EPH(6) * ASQRT C------------------------------------------------------------ C A0 COEFFICIENT C------------------------------------------------------------ CLK(5) = CLK(5) - RK * SINE C------------------------------------------------------------ C A1 COEFFICIENT C------------------------------------------------------------ COSE = DCOS (AMYE) DENOM = (1.D0 - EPH(6) * COSE) CLK(4)= CLK(4) - RK * EPH(18) * COSE / DENOM C------------------------------------------------------------ C A2 COEFFICIENT C------------------------------------------------------------ RMM2 = EPH(18) * EPH(18) CLK(3) = CLK(3) + RK * RMM2 * SINE / (2.D0 * DENOM**2) 1 CONTINUE RETURN END C SUBROUTINE ANML2(AMYM,ORBEC,AMYE,AMYT,IFLAG,PI,IPRINT) C------------------------------------------------------------ C FUNCTION : CONVERT ECCENTRICITY AND MEAN ANOMALY TO C ECCENTRIC AND TRUE ANOMALY. C C HISTORY : VERSION DATE PROGRAMMER C 1 1981 U.N.B.? C 2 JUNE 1984 AL. KLEUSBERG C 3 FEB 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I LOG. UNIT OF PRINTER C AMYM I MEAN ANOMALY (RAD) C PI I GREEK PI=3.14159... C IPRINT I L. UNIT OF PRINTER. C ORBEC I ORBITAL ECCENTRICITY C AMYE O ECCENTRIC ANOMALY (RAD) C AMYT O TRUE ANOMALY (RAD) C IFLAG O =0 IF EVERYTHING OK C =1 IF NO CONVERGENCE ON KEPLERS C EQUATION FOR ECC. ANOMALY C C EXTERNALS : NONE C------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) 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.D0-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 COMPLETED C----------------------------------------------------------- 40 CONTINUE SE = DSIN (AMYE) DIFF = DABS(AMYM-(AMYE-ORBEC*SE)) IF( DIFF. LT. 1.0D-10 ) GO TO 60 IFLAG = 1 WRITE(IPRINT,50) DIFF, I, AMYM, AMYE 50 FORMAT(2X,'ANML2: ECCENTRICITY MISCOSURE=',D20.10,'ITERATION=', # I4,'AMYM=',D20.10,'MAYE=',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 C SUBROUTINE ETABL(IPRINT,LUEPH,LREC,NREC,ICOUN,EPHRT,IRECD, # MSAT,MEPH) C------------------------------------------------------------ C FUNCTION : SCAN THE EPHEMERIS DISC FILE AND CREATE TABLES C ICOUN(MSAT),EPHRT(MSAT,MEPH) AND IRECD(MSAT,MEPH) C FOR EASY AND FAST ACCESS TO EPHEMERIS OF A C CERTAIN EPOCH. C C HISTORY : VERSION DATE PROGRAMMER C 1 JUNE 1984 AL. KLEUSBERG C 2 FEB 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I LOG. UNIT OF PRINTER C LUEPH I LOG. UNIT OF EPHEM. FILE C LREC I LENGTH RECORD EPHEM. FILE C NREC I NUMBER OF RECORDS IN EPHM.FILE C MSAT I ROW DIMENSION OF ICOUN, C EPHRT, IRECD. C MEPH I COLUMN DIMENSION OF EPHRT,IRECD C ICOUN O COUNTER VECTOR ICOUN(MSAT) C CONTAINING THE # OF ENTRIES IN C EPHRT AND IRECD. C IRECD O ARRAY CONTAINING RECORD # OF C EPHEMERIS DISK FILE,CORRESPONDING C TO EPHRT. C EPHRT O ARRAY CONTAINING EPHEMERIS C REFERENCE TIMES (TOE). C C EXTERNALS : NONE C------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) REAL*4 EPHRT DIMENSION ICOUN(MSAT), EPHRT(MSAT,MEPH), IRECD(MSAT,MEPH), # DEPH(20,4) , DCLK(5,4) DO 100 K=1,NREC IREC = K READ(LUEPH,678,IOSTAT=IERROR,ERR=300,REC=IREC) # ((DEPH(I,J),J=1,4),I=1,20) 678 FORMAT(100G30.15) DO 200 J=1,4 IDSAT = DEPH(20,J) IF(IDSAT.EQ.0) GO TO 200 IF(IDSAT.GT.MSAT) GO TO 190 ICOUN(IDSAT) = ICOUN(IDSAT) + 1 IF(ICOUN(IDSAT).GT.MEPH) GO TO 310 C EPHRT(IDSAT,ICOUN(IDSAT)) = SNGL(DEPH(9,J)) IRECD(IDSAT,ICOUN(IDSAT)) = IREC GO TO 200 190 WRITE(IPRINT,1020) IDSAT 200 CONTINUE 100 CONTINUE GO TO 320 300 WRITE(IPRINT,1000) IERROR GO TO 320 310 WRITE(IPRINT,1010) 320 CONTINUE 600 FORMAT(4G30.15) 800 FORMAT(2X,'ICOUN=',I8,3X,'ID OF SATELLITE=',I5) 700 FORMAT(2X,'OBSERVATION RECORD=',I6) 1000 FORMAT(2X,'ETABL-ERROR',I5,'WHILE READING EPHEMERIS FILE') 1010 FORMAT(2X,'ETABL-ERROR EHHRT AND IRECD ISUFFICIENT DIMEN') 1020 FORMAT(2X,'ETABL-ERROR SUSPICIOUS SAT ID',I5,'NOT INCLUDED') RETURN END C SUBROUTINE FINDE(IPRINT,LUEPH,ISAT,TIM ,ICOUN,EPHRT,IRECD, # MSAT,MEPH,DTM,EPH,CLK) C------------------------------------------------------------ C FUNCTION : EXTRACT EPHEMERIS AND CLOCK PARAMETERS FOR A C CERTAIN EPOCH OF SIGNAL TRANSMISSION, USING THE C EPHEMERIS REFERENCE TABLE C C HISTORY : VERSION DATE PROGRAMMER C 1 JUNE 1984 AL. KLEUSBERG C 2 FEB 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I LOG. UNIT OF PRINTER C LUEPH I LOG. UNIT OF EPHEM. FILE C ISAT I/O SAT ID'S ARE SET TO ZERO IF NO C EPHEMERIS ARE READ C TIM I EPOCH EPHEMERIS ARE NEEDED C MSAT I ROW DIMENSION OF ICOUN, C EPHRT, IRECD. C MEPH I COLUMN DIMENSION OF EPHRT,IRECD C ICOUN I COUNTER VECTOR ICOUN(MSAT) C CONTAINING THE # OF ENTRIES IN C EPHRT AND IRECD. C IRECD I ARRAY CONTAINING RECORD # OF C EPHEMERIS DISK FILE,CORRESPONDING C TO EPHRT. C EPHRT I ARRAY CONTAINING EPHEMERIS C REFERENCE TIMES (TOE). C EPH O EPHEMERIS FOR 'ISAT' AND 'TIM' C CLK O CLOCK PARAM. FOR 'ISAT' AND 'TIM' C C EXTERNALS : NONE C------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) REAL*4 EPHRT DIMENSION ICOUN(MSAT), EPHRT(MSAT,MEPH), IRECD(MSAT,MEPH) DIMENSION DEPH(20,4) , EPH(20,4) DIMENSION DCLK(5,4), CLK(5,4) DIMENSION ISAT(4), DTM(7), TIM(4) DO 100 I=1,4 IF(ISAT(I).EQ.0) GO TO 100 IF((ISAT(I).LT.0).OR.(ISAT(I).GT.MSAT)) GO TO 155 TIME = TIM(I) IREC = 0 TDIF = 10800.D0 IDSAT = ISAT(I) MAXSAT = ICOUN(IDSAT) DO 110 IC=1,MAXSAT TDUF = DABS(DBLE(EPHRT(IDSAT,IC))- TIME) IF(TDUF.GT.TDIF) GO TO 110 TDIF = TDUF IREC = IRECD(IDSAT,IC) 110 CONTINUE IF(IREC.NE.0) GO TO 120 WRITE(IPRINT,1000) IDSAT GO TO 140 120 CONTINUE READ(LUEPH,9,IOSTAT=IERROR,ERR=130,REC=IREC) ((DEPH(L,M),M=1,4), # L=1,20), ((DCLK(L,M),M=1,4),L=1,5) 9 FORMAT(100(G30.15)) DO 150 J=1,4 ISA = DEPH(20,J) IF(ISA.EQ.IDSAT) GO TO 160 150 CONTINUE WRITE(IPRINT,1000) IDSAT GO TO 140 155 WRITE(IPRINT,1020) ISAT(I) GO TO 140 160 CONTINUE DO 666 II=1,20 EPH(II,I) = DEPH(II,J) 666 CONTINUE C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ WRITE(6,*) IDSAT, IREC C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ A = EPH(8,I) EPH(18,I) = DSQRT(DTM(6)/A/A/A) + EPH(3,I) DO 777 II=1,5 CLK(II,I) = DCLK(II,J) 777 CONTINUE GO TO 100 130 CONTINUE WRITE(IPRINT,1010) IERROR, IDSAT 140 CONTINUE ISAT(I) = 0 DO 888 II=1,20 EPH(II,I) = 0.D0 888 CONTINUE DO 999 II=1,5 CLK(II,I) = 0.D0 999 CONTINUE 100 CONTINUE RETURN 1000 FORMAT(2X,'FINDE: EPHEMERIS OF SAT #',I3,'OLDER THAN 3 HOURS #OR NOT IN RECORD') 1010 FORMAT(2X,'FINDE: ERROR',I5,'WHILE READING EPH OF SAT #',I3) 1020 FORMAT(2X,'FINDE: SUSPICIOUS SAT ID#',I5,'NO EPHEMERIS') END C FUNCTION RECEPT(DAT1,DAT2,DAT3) C COMPUTES TIME OF RECEPTION IMPLICIT REAL*8(A-H,O-Z) RECEPT = (DAT2-DAT1)*0.02D0 + DAT3*0.02D0 RETURN END C SUBROUTINE TRANSM(OBS,F1,F2,ISAT,IPRINT,TIMETR) C------------------------------------------------------------ C FUNCTION : COMPUTES TIME OF TRANSMISSION C C HISTORY : VERSION DATE PROGRAMMER C 1 FEB 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IPRINT I LOG. UNIT OF PRINTER C F1 I L1 GPS TRANSMISSION FREQUENCY C F2 I L2 GPS TRANSMISSION FREQUENCY C ISAT I SAT IDENTITIES C TIMETR O VECTOR OF TRANSMISSION TIMES C C EXTERNALS : TION C------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) DIMENSION OBS(13,4), TIMETR(4), ISAT(4) DO 100 I=1,4 IF(ISAT(I). EQ. 0) GO TO 100 CODE1 = OBS(2,I) + OBS(3,I) + OBS(4,I) CODE2 = OBS(5,I) + OBS(6,I) + OBS(7,I) CALL TION(CODE1,CODE2,F1,F2) TIMETR(I) = CODE1 100 CONTINUE RETURN END C SUBROUTINE TION(T1,T2,F1,F2) C------------------------------------------------------------ C FUNCTION : COMPUTES IONOSHERIC CORRECTION FOR P-CODE C C HISTORY : VERSION DATE PROGRAMMER C 1 1980 U.N.B.?????? C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C T1 I/O P-CODE MEASUREMENT ON L1 C T2 I P-CODE MEASUREMENT ON L2 C F1 I L1 GPS TRANSMISSION FREQUENCY C F2 I L2 GPS TRANSMISSION FREQUENCY C EXTERNALS : NONE C------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) DT = T1 - T2 F1SQ = F1 ** 2 F2SQ = F2 ** 2 DT = DT * F2SQ / (F2SQ - F1SQ) T1 = T1 + DT RETURN END C SUBROUTINE DATUM(DTM,IDTM,IREAD) IMPLICIT REAL*8(A-H,O-Z) DIMENSION DTM(7) 1 CONTINUE 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 3 CONTINUE C----------------------------------------- C USER DEFINED DATUM C----------------------------------------- READ(IREAD,100) (DTM(I),I=1,7) 100 FORMAT(7G12.4) BDVDA = DTM(2) / DTM(1) DTM(2) = DSQRT( (1.D0 + BDVDA) * (1.D0 - BDVDA) ) 4 CONTINUE DTM(6) = 398600.8D9 DTM(7) = 7.292115147D-5 RETURN END C FUNCTION RANGE(XSV,XRV) C------------------------------------------------- C COMPUTES THE STATION TO SATELLITE RANGE C------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) DIMENSION XRV(3), XSV(3) R = 0D0 DO 10 I=1,3 10 R = ( XSV(I) - XRV(I) )** 2 + R RANGE = DSQRT(R) RETURN END C SUBROUTINE SET(ARR,N,CONST) C------------------------------------------------------------ C FUNCTION : SETS ARRAY TO A CONSTANT C C HISTORY : VERSION DATE PROGRAMMER C 1 JUNE 1984 AL. KLEUSBERG C 2 FEB 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C IARR I/O VECTOR C N I DIMENSION OF VECTOR C CONST I CONSTANT C C EXTERNALS : NONE C------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) DIMENSION ARR(100) DO 10 I=1,N ARR(I) = CONST 10 CONTINUE RETURN END C SUBROUTINE DMOVE(ARRAY1,NROW,ARRAY2,IFIRST,ILAST) C------------------------------------------------------------ C FUNCTION : MOVE ARRAY2 TO ARRAY1 C C HISTORY : VERSION DATE PROGRAMMER C 1 FEB 1986 S. MERTIKAS C C PARAMETERS: NAME I/O DESCRIPTION C ---------------------------------- C ARRAY2 O ARRAY TO BE FILLED C ARRAY1 I ARRAY GINEN C NROW I ROW DIMENSION C IFIRST I FIRST ELEMENT TO BE FILLED C ILAST I LAST ELEMENT TO BE FILLED C C EXTERNALS : NONE C------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) DIMENSION ARRAY1(100),ARRAY2(100) IF(ILAST.LT.NROW) WRITE(6,1000) DO 1 I=IFIRST,ILAST ARRAY2(I) = ARRAY1(I) 1 CONTINUE 1000 FORMAT(2X,'DMOVE: ERROR ENCOUNTERED IN DIMENSION') RETURN END C //* //LKED.USERLIB DD DSN=A.M12129.SELIBOBJ,DISP=SHR //GO.FT22F001 DD DSN=A.M12128.MERTIKAS.FINAL, // UNIT=P3350,DISP=(NEW,CATLG,DELETE), // DCB=(RECFM=FB,LRECL=128,BLKSIZE=17920), // SPACE=(TRK,(6,1,1),RLSE) //GO.FT17F001 DD DSN=A.M12128.MERTIKAS.GEDOP, // UNIT=P3350,DISP=(NEW,CATLG,DELETE), // DCB=(RECFM=FB,LRECL=98,BLKSIZE=17640), // SPACE=(TRK,(6,1,1),RLSE) //GO.FT11F001 DD DSN=A.M12128.MERTIKAS.OBSRV, // UNIT=P3350,DISP=OLD //GO.FT12F001 DD DSN=A.M12128.MERTIKAS.EPHEM, // UNIT=P3350,DISP=OLD //GO.SYSIN DD * 310 20 35 27.8 17785.2 25775.9 17307.6