PROGRAM TEST C C PURPOSE: ILLUSTRATE USAGE OF SUBROUTINES HEIGHTDRIVER, HEIGHT, AND C ASSOCIATED SUBROUTINES C THIS PROGRAM DUPLICATES THE FIRST 10 ENTRIES C (TIME,LAT,LON,HEIGHT) C IN THE GTRACK FILE (GRNDTK). THIS COULD BE USED C AS A FIRST TEST TO MAKE SURE THE CODE WORKS C C GERARD L.H. KRUIZINGA 10/20/93 C C MODIFIED BY: GERARD L.H. KRUIZINGA 05/08/94 C MODIFICATION TO TAKE INTO ACCOUNT PRECISION C OTHER THAN CRAY 64-BIT MACHINES C MODIFIED BY: MICHAEL J. GABOR 3/20/95 C MODIFICATION TO TAKE BURN TIMES INTO ACCOUNT C IMPLICIT NONE DOUBLE PRECISION EPOCHORB,EPOCHREC DOUBLE PRECISION DEPOCH,DEPOCH1,EPOCH,TORB(4) DOUBLE PRECISION TIME,DELTT,TBOFF(1000),TBON(1000) DOUBLE PRECISION Y(4,3),HT,QLAT,QLON DOUBLE PRECISION SECONDS,AE,RFLAT,DELT INTEGER ISATID,IDATE(5) INTEGER INITFLG,IUNITG,IUNITO,I INTEGER IFLAG,IUNITM,NUMBRN CHARACTER*10 JOBTIM,JOBDAT,NAMESAT,VERSION LOGICAL MFLAG COMMON /HTSTUFF/ EPOCHORB,DEPOCH,DEPOCH1,TORB,Y,INITFLG COMMON /GTRINFO/ DELT,ISATID,IDATE,AE,RFLAT,EPOCH,SECONDS, . VERSION,JOBDAT,JOBTIM,NAMESAT DATA MFLAG/.TRUE./ DELTT=60.D0 IUNITG = 11 IUNITO = 12 IUNITM = 13 C INITIALIZE NUMBRN BEFORE THE FIRST CALL OF HEIGHTDRIVER SUBROUTINE. C SUBROUTINE UPDATES IT TO ACTUAL VALUE. MAKE SURE THAT INITIALLY C IT IS LARGE NUMBRN=1000 OPEN(IUNITO,FILE='OUTPUT') OPEN(IUNITG,FILE='GRNDTK',STATUS='OLD') OPEN(IUNITM,FILE='MANEUVERS',STATUS='OLD') REWIND(IUNITM) CALL CHECK_HEADER(IUNITG) REWIND (IUNITG) C C SOME COMMENTS ON CALLING SUBROUTINE HEIGHTDRIVER C C VARIABLE "TIME" IS IN SECONDS PAST EPOCHREC (IN MJD) C VARIABLE "HT" IS IN METERS C VARIABLES "QLAT" AND "QLON" ARE IN DEGREES (QLAT IS GEODETIC LATITUDE) C VARIABLE "DELTT" SPECIFIES REGION OF PERMITTED TIME AROUND BURN IN C SECONDS. TOTAL REGION IS DELTT ON EITHER SIDE OF BURN. C ARRAY "TBON" AND "TBOFF" SPECIFY THE START AND STOP TIMES OF THE C BURNS AND ARE AVAILABLE FROM THE MANEUVERS TABLE C EPOCHREC= EPOCH DO 10 I = 110530,110570 C DO 10 I = 1,20 TIME= DFLOAT(I-1)*10.0D0 CALL HEIGHTDRIVER (IUNITM,DELTT,MFLAG,NUMBRN,TBON,TBOFF, & TIME,EPOCHREC,IUNITG,HT,QLAT,QLON,IFLAG,IUNITO) WRITE(IUNITO,1000) TIME,QLAT,QLON,HT 10 CONTINUE CLOSE(IUNITG) CLOSE(IUNITO) CLOSE(IUNITM) 1000 FORMAT (F15.6,2F15.8,F15.4) STOP END SUBROUTINE HEIGHTDRIVER (IUNITM,DELTT,MFLAG,NUMBRN,TBON,TBOFF, & TIME,EPOCHREC,IUNITG,HT,QLAT,QLON,IFLAG,IUNITO) ************************************************************************ * MODIFIED BY: MICHAEL J. GABOR 3/20/95 * MODIFICATION TO TAKE BURN TIMES INTO ACCOUNT * THIS IS A MODIFICATION TO THE HEIGHT SUBROUTINE (DETAILS UNDER * HEIGHT) THAT ALLOWS FOR DATA TO BE IGNORED THAT IS WITHIN * A SPECIFIED NUMBER OF SECONDS OF A BURN. DATA IS ARBITRARILY * REMOVED DURING THE ACTUAL BURN. TO INCLUDE DATA THAT IS DURING * A BURN, SPECIFY 0.D0 SECONDS FOR DELTT. * VARIABLES SPECIFICALLY USED IN HEIGHTDRIVER: * ARRAY "TBON" AND "TBOFF" SPECIFY THE START AND STOP TIMES OF THE * BURNS AND ARE AVAILABLE FROM THE MANEUVERS TABLE * VARIABLE "MFLAG" DETERMINES IF THIS IS THE FIRS PASS THROUGH THE * LIST OF BURNS * VARIABLE "DELTT" SPECIFIES REGION OF PERMITTED TIME AROUND BURN IN * SECONDS. TOTAL REGION IS DELTT ON EITHER SIDE OF BURN. ************************************************************************ IMPLICIT NONE DOUBLE PRECISION DELTT,TBON(1000),TBOFF(1000),TIME,EPOCHREC,HT DOUBLE PRECISION QLAT,QLON INTEGER IUNITM,IUNITG,NUMBRN,I,IFLAG,IUNITO LOGICAL MFLAG,STATUS C DATA STATUS/.TRUE./ SAVE STATUS = .TRUE. QLAT = -999.0D0 QLON = -999.0D0 HT = -999.0D0 DO 100 I=1,NUMBRN IF (MFLAG.EQ..TRUE.) THEN READ(IUNITM,5,END=200) TBON(I),TBOFF(I) TBON(I)=(TBON(I)-2400000.5D0-EPOCHREC)*86400.D0 TBOFF(I)=(TBOFF(I)-2400000.5D0-EPOCHREC)*86400.D0 NUMBRN=I END IF IF (DABS(TIME-TBOFF(I)).LT.DELTT.OR.DABS(TIME-TBON(I)) & .LT.DELTT) THEN STATUS=.FALSE. END IF IF (TIME.GT.TBON(I).AND.TIME.LT.TBOFF(I).AND.DELTT.GT.0.D0) & THEN STATUS=.FALSE. END IF 100 CONTINUE 200 MFLAG=.FALSE. IF (STATUS.EQ..TRUE.) THEN CALL HEIGHT (TIME,EPOCHREC,IUNITG,HT,QLAT,QLON,IFLAG) ELSE END IF 5 FORMAT(2X,F18.10,4X,F18.10) 1000 FORMAT (F15.6,2F15.8,F15.4) END SUBROUTINE HEIGHT (TREC,EPOCHREC,UNIT,HT,LAT,LONG,IFLAG) ************************************************************************ * CODED BY: DON CHAMBERS -- 2/4/92 -- VERSION 1.3 * * MODIFIED BY: GERARD L.H. KRUIZINGA 4/1/92 * INCLUSION OF CORRECT SPLINE INTERPOLATION ACROSS * THE PRIME (GREENWHICH MERIDIAN) * MODIFIED BY: GERARD L.H. KRUIZINGA 5/8/94 * MODIFICATIONS TO TAKE INTO ACCOUNT PRECISION * DIFFERENCES BETWEEN CRAY AND SUN WORKSTATIONS * * THIS SUBROUTINE READS THROUGH A UTOPIA GTRACK FILE CONTAINING HEIGHTS * ABOVE THE REFERENCE ELLIPSOID SPACED AT 10 SECOND INTERVALS FROM THE * EPOCH OF THE CYCLE. * * THIS SUBROUTINE INTERPOLATES (USING A CUBIC SPLINE) TO FIND THE ORBIT * HEIGHT AT TREC SO THAT THE ORBIT HEIGHT CAN BE COMPARED TO THE GDR * HEIGHT. * * * VARIABLES: * INPUT-- * TREC - TIME OF RECORD (SECONDS PAST EPOCHREC) * EPOCHREC - EPOCH OF RECORD (IN MJD) * UNIT - UNIT NUMBER FOR ORBIT HEIGHTS FILE * OUTPUT-- * HT - ORBITAL HEIGHT ABOVE REFERENCE ELLIPSOID (M) * LAT - GEODETIC LATITUDE OF SUB-SATELLITE POINT (IN DEGREES) * LONG - GEOCENTRIC LONGITUDE OF SUB-SATELLITE POINT (IN DEGREES) * IFLAG - ERROR FLAG * IFLAG = 0 : SUCCESFUL INTERPOLATION * IFLAG =-1 : REQUESTED TIME BEFORE EPOCH * IFLAG = 1 : REQUESTED TIME AFTER LAST RECORD IN FILE * OTHERS-- * EPOCHORB - EPOCH OF ORBIT CYCLE (IN JD) * TORB - TIME OF ORBIT PAST EPOCHORB (IN SECONDS) * HTERR - ERROR MESSAGES * * MODIFICATION-- * HEIGHT SUBROUTINE CAN READ GTRACK FILES WITH AND WITHOUT * GTRACK HEADERS. SUBROUTINE WILL CHECK IF CORRECT TIME SPACING * ON APPEARS ON GTRACK HEADER * * MODIFICATION * INTERPOLATION WILL TAKE INTO ACCOUNT LEAP SECONDS * **************************************************************************** IMPLICIT NONE DOUBLE PRECISION TREC,EPOCHORB,EPOCHREC DOUBLE PRECISION DEPOCH,DEPOCH1,TORB(4),TOLD,EPOCH,DTNEW DOUBLE PRECISION LAT,LONG DOUBLE PRECISION SECONDS,AE,RFLAT,DELT DOUBLE PRECISION Y(4,3),HT DOUBLE PRECISION DT,DY2,DY3,DY4,A,B,C,VAL INTEGER UNIT INTEGER ISATID,IDATE(5) INTEGER INITFLG,I,J,IFLAG CHARACTER*80 RECORD,HTERR CHARACTER*10 JOBTIM,JOBDAT,NAMESAT,VERSION CHARACTER*5 TEXT LOGICAL FOUND,FLAG COMMON /HTSTUFF/ EPOCHORB,DEPOCH,DEPOCH1,TORB,Y,INITFLG COMMON /GTRINFO/ DELT,ISATID,IDATE,AE,RFLAT,EPOCH,SECONDS, . VERSION,JOBDAT,JOBTIM,NAMESAT DATA NAMESAT /'NOTDET'/ DATA TOLD /999999999999999.0/ FOUND = .FALSE. HTERR = 'ERROR IN READING ORBIT HEIGHT FILE' IFLAG = 0 DO WHILE (.NOT. FOUND) * * READ IN FIRST 4 POINTS THE FIRST TIME THE SUBROUTINE IS CALLED * IF (INITFLG .EQ. 0) THEN TOLD = 99999999999999.0 CALL CHECK_HEADER(UNIT) IF (NAMESAT .EQ. 'NOTDET') THEN READ (UNIT,1000,END=99,ERR=199) RECORD READ(RECORD,1200,ERR=299) TEXT,EPOCHORB DEPOCH = (EPOCHORB - EPOCHREC)*86400.0D0 DEPOCH1 = DEPOCH ELSE IF (DELT .NE. 10.0) THEN WRITE(*,*) ' INCORRECT TIME INTERVAL FOR HEIGHT' WRITE(*,*) ' INTERVAL : ',DELT,' SHOULD BE 10.0 SEC' STOP ENDIF DEPOCH = (EPOCH - EPOCHREC)*86400.0D0 DEPOCH1 = DEPOCH ENDIF DO 10 J =1,4 9 READ (UNIT,1000,END=99,ERR=199) RECORD READ (RECORD,1100)TORB(J),Y(J,1),Y(J,2),Y(J,3) TORB(J) = TORB(J) + DEPOCH IF (TOLD .EQ. TORB(J)) GOTO 9 TOLD = TORB(J) 10 CONTINUE INITFLG = 1 WRITE(*,*) ' ' WRITE(*,*) 'INITIALIZATION INFORMATION FROM HEIGHT' WRITE(*,*) ' ' IF (NAMESAT .EQ. 'NOTDET') THEN WRITE(*,*) 'EPOCH READ FROM GTRACK FILE :',EPOCHORB ELSE WRITE(*,*) 'EPOCH READ FROM GTRACK FILE :',EPOCH ENDIF WRITE(*,*) 'EPOCH FOR TIME TAG :',EPOCHREC WRITE(*,*) 'RELATIVE EPOCH :',DEPOCH WRITE(*,*) ' ' END IF IF (TREC.LT.DEPOCH1) THEN HTERR = 'RECORD TIME IS BEFORE ORBIT TIME' IFLAG=-1 GO TO 199 END IF * * SEE IF T LIES IN THE FOUR POINTS; IF SO, INTERPOLATE TO FIND * HT, LAT, AND LONG; IF NOT, SEARCH THROUGH NEXT FOUR POINTS * IF (TREC.GE.TORB(1) .AND. TREC.LT.TORB(4)) THEN DO 30 I =1,3 CALL CHECK_LONGITUDE DY2 = Y(2,I) - Y(1,I) DY3 = Y(3,I) - Y(1,I) DY4 = Y(4,I) - Y(1,I) DT = TREC - TORB(1) CALL LEAPSEC(TREC,TORB,EPOCHREC,DY2,DY3,DY4, & A,B,C,FLAG,DTNEW) IF (FLAG) THEN DT=DTNEW GOTO 29 ENDIF A = (DY4 + 3.0*DY2 - 3.0*DY3)/6000.0 B = (4.0*DY3 - 5.0*DY2 - DY4)/200.0 C = (DY4 - 4.5*DY3 + 9.0*DY2)/30.0 29 VAL = A*DT*DT*DT + B*DT*DT + C*DT + Y(1,I) IF (I.EQ.1) LAT = VAL IF (I.EQ.2) LONG = VAL IF (I.EQ.2 .AND. VAL .GT. 360.0) LONG = VAL - 360.0 IF (I.EQ.3) HT = VAL 30 CONTINUE IF (Y(4,2) .GT. 360.0) Y(4,2)=Y(4,2)-360.0 FOUND = .TRUE. ELSE TORB(1) = TORB(4) Y(1,1) = Y(4,1) Y(1,2) = Y(4,2) Y(1,3) = Y(4,3) DO 35 J = 2,4 721 READ (UNIT,1000,END=99,ERR=199) RECORD IF (RECORD(1:5) .EQ. 'EPOCH') THEN READ(RECORD,1200) TEXT,EPOCHORB DEPOCH = (EPOCHORB - EPOCHREC)*86400.0 * * SKIP FIRST RECORD AFTER EACH EPOCH CARDS PAST THE FIRST ONE * READ (UNIT,1000,END=99,ERR=199) RECORD READ (UNIT,1000,END=99,ERR=199) RECORD END IF READ (RECORD,1100)TORB(J),Y(J,1),Y(J,2),Y(J,3) TORB(J) = TORB(J) + DEPOCH IF (TOLD .EQ. TORB(J)) GOTO 721 TOLD = TORB(J) 35 CONTINUE END IF END DO IF (.NOT. FOUND) THEN HTERR = 'RECORD TIME IS BIGGER THAN LARGEST ORBIT TIME' IFLAG=1 GO TO 199 END IF RETURN C 99 PRINT *,'END OF ORBIT HEIGHT FILE REACHED' C PRINT *,'LAST REQUESTED TIME :',TREC 99 CONTINUE IFLAG=1 RETURN C 199 PRINT *,'ERROR OCCURRED IN SUBROUTINE HEIGHT' C PRINT *,HTERR 199 CONTINUE RETURN 299 PRINT *,'ERROR READING FIRST LINE OF GTRACK FILE; NO EPOCH HEADER' STOP 1000 FORMAT (A80) 1100 FORMAT (F15.6,2F15.8,F15.4) 1150 FORMAT (F15.7,2F15.8,F15.4) 1200 FORMAT (A5,1X,F25.15) END SUBROUTINE CHECK_LONGITUDE * * SUBROUTINE TO CHECK IF ALL 4 LONGITUDES DO NOT CONTAIN A * DISCONTINUITY CAUSED BY CROSSING THE PRIME MERIDIAN * CODED BY : GERARD L.H. KRUIZINGA 4/2/92 IMPLICIT NONE DOUBLE PRECISION EPOCHORB,DEPOCH,DEPOCH1,TORB(4) DOUBLE PRECISION Y(4,3) INTEGER INITFLG,I COMMON /HTSTUFF/ EPOCHORB,DEPOCH,DEPOCH1,TORB(4),Y(4,3),INITFLG IF (ABS(Y(1,2)-Y(4,2)).GT. 100.0) THEN DO 10 I=1,4 IF ( Y(I,2) .LT. 100) Y(I,2)=Y(I,2)+360.0 10 CONTINUE ENDIF RETURN END C********************************************************************** SUBROUTINE KALDAY (ITIME, STIME, UTC) C C...PURPOSE: TO CONVERT BETWEEN CALENDAR DATE AND JULIAN DATE C C...PROGRAMMED BY JAMES D MCMILLAN - UNIV. OF TEXAS - 7/23/1973 C JULDAY ADDED BY BRIAN CUTHBERTSON - UNIV. OF TEXAS - 9/29/1979 C C...KALDAY INPUT (JULDAY OUTPUT): C UTC THE DOUBLE PRECISION JULIAN DATE C C...KALDAY OUTPUT (JULDAY INPUT): C ITIME INTEGER ARRAY WITH MONTH, DAY, YEAR, HOUR AND MINUTE C STIME SINGLE PRECISION SECONDS C DIMENSION ITIME(5) DOUBLE PRECISION UTC C C...COMPUTE THE HOUR, MINUTE, AND SECOND. JD=UTC XTIME=AMOD(SNGL(UTC-JD)*24.+12.,24.) ITIME(4)=XTIME XTIME=(XTIME-ITIME(4))*60.0 ITIME(5)=XTIME STIME=(XTIME-ITIME(5))*60.0 C C...ROUND 60 SECONDS TO NEXT MINUTE IF (STIME.LT.59.99999999999) GO TO 10 STIME=0.0 ITIME(5)=ITIME(5)+1 IF (ITIME(5).NE.60) GO TO 10 ITIME(5)=0 ITIME(4)=ITIME(4)+1 IF (ITIME(4).NE.24) GO TO 10 ITIME(4)=0 C C...COMPUTE THE YEAR, MONTH, AND DAY. 10 IF (ITIME(4).LT.12) JD=JD+1 LX=JD+68569 NX=4*LX/146097 LX=LX-(146097*NX+3)/4 ITIME(3)=4000*(LX+1)/1461001 LX=LX-1461*ITIME(3)/4+31 ITIME(1)=80*LX/2447 ITIME(2)=LX-2447*ITIME(1)/80 LX=ITIME(1)/11 ITIME(1)=ITIME(1)+2-12*LX ITIME(3)=100*(NX-49)+ITIME(3)+LX RETURN C ENTRY JULDAY (ITIME, STIME, UTC) C JD=ITIME(2)-32075+1461*(ITIME(3)+4800+(ITIME(1)-14)/12)/4+367*(ITI 1ME(1)-2-(ITIME(1)-14)/12*12)/12-3*((ITIME(3)+4900+(ITIME(1)-14)/12 2)/100)/4 UTC=DBLE(FLOAT(JD))-0.5D+00+DBLE(FLOAT(ITIME(4)))/24.0D+00+DBLE(FL 1OAT(ITIME(5)))/1440.0D+00+DBLE(STIME)/86400.0D+00 RETURN C C.....END KALDAY/JULDAY END C********************************************************************** SUBROUTINE GET_UTC(KYR,KDAY,KSEC,KMSEC,TGDR) C C THIS SUBROUTINE COMPUTES GEOSAT GDR TIME (SEC PAST 01/01/1985) C C CODED BY GERARD L.H. KRUIZINGA JULY 199OA2 C CENTER FOR SPACE RESEARCH UNIVERSITY OF TEXAS AT AUSTIN C C INPUT : KYR YEAR C KDAY DAY IN THE YEAR C KSEC SECONDS PAST MIDNIGHT OF KDAY C KMSEC MICRO SECONDS C OUTPUT : TGDR SECONDS PAST 01/01/1985 C DOUBLE PRECISION TGDR,UTC,JAN1,GEOSAT_EPOCH DIMENSION ITIME(5) DATA ITEST /0/ DATA GEOSAT_EPOCH /0.0D0/ IF (ITEST .EQ. 0) THEN ITIME(1) = 1 ITIME(2) = 1 ITIME(3) = 1985 ITIME(4) = 0 ITIME(5) = 0 STIME = 0.0 CALL JULDAY (ITIME,STIME,UTC) GEOSAT_EPOCH = UTC ITEST = 1 ENDIF ITIME(1) = 1 ITIME(2) = 1 ITIME(3) = 1900 + KYR ITIME(4) = 0 ITIME(5) = 0 STIME = 0.0 CALL JULDAY (ITIME, STIME, JAN1) TGDR = (JAN1 - GEOSAT_EPOCH)*86400.0 + DBLE(KDAY-1)*86400.0 & + DBLE(KSEC) + FLOAT(KMSEC)/1E6 RETURN END SUBROUTINE CHECK_HEADER(IUNIT) C C PURPOSE: CHECK FOR THE EXISTANCE OF AN GTRACK HEADER C AND IF SO READ THE INFO C C CODED BY: GERARD L.H. KRUIZINGA 05/07/93 C CENTER FOR SPACE RESEARCH UNIVERSITY OF TEXAS AT AUSTIN C IMPLICIT NONE DOUBLE PRECISION EPOCH DOUBLE PRECISION SECONDS,AE,RFLAT,DELT INTEGER ISATID,IDATE(5),IUNIT CHARACTER*10 JOBTIM,JOBDAT,NAMESAT,VERSION CHARACTER*60 CARD COMMON /GTRINFO/ DELT,ISATID,IDATE,AE,RFLAT,EPOCH,SECONDS, . VERSION,JOBDAT,JOBTIM,NAMESAT READ(IUNIT,390) CARD READ(IUNIT,390) CARD READ(IUNIT,390) CARD READ(IUNIT,390) CARD PRINT*,'CARD = ',CARD IF (CARD(1:7) .EQ. 'VERSION') THEN IF (CARD(1:8) .EQ. 'VERSION:') THEN READ(CARD,400) VERSION,DELT READ(IUNIT,410) JOBDAT,JOBTIM,ISATID,NAMESAT READ(IUNIT,390) CARD READ(IUNIT,420) AE READ(IUNIT,430) RFLAT READ(IUNIT,390) CARD READ(IUNIT,440) EPOCH,IDATE,SECONDS READ(IUNIT,390) CARD READ(IUNIT,190) CARD READ(IUNIT,390) CARD ELSE READ(CARD,200) VERSION,DELT READ(IUNIT,210) JOBDAT,JOBTIM,ISATID,NAMESAT READ(IUNIT,190) CARD READ(IUNIT,220) AE READ(IUNIT,230) RFLAT READ(IUNIT,190) CARD READ(IUNIT,240) EPOCH,IDATE,SECONDS READ(IUNIT,190) CARD READ(IUNIT,190) CARD READ(IUNIT,190) CARD ENDIF ELSE REWIND(IUNIT) ENDIF WRITE(*,*) DELT,ISATID,AE,RFLAT,EPOCH,IDATE,SECONDS, . VERSION,JOBDAT,JOBTIM,NAMESAT 190 FORMAT (A79) 200 FORMAT (10X,A10,30X,F9.1) 210 FORMAT (9X,2A10,23X,I7,1X,A10) 220 FORMAT (49X,F10.1) 230 FORMAT (49X,F10.5) 240 FORMAT (8X,F20.12,12X,I3,1X,I2,1X,I2,I4,1X,I2,1X,F10.7) 390 FORMAT (A60) 400 FORMAT (9X,A10,22X,F9.1) 410 FORMAT (8X,2A10,15X,I7,A10) 420 FORMAT (49X,F11.2) 430 FORMAT (49X,F11.6) 440 FORMAT (8X,F20.12,6X,I3,1X,I2,1X,I2,I3,1X,I2,1X,F10.7) RETURN END C***************************************************************** SUBROUTINE CRAMER(A,B,X) C C PURPOSE: SOLVE 3X3 SYSTEM A.X = B USING CRAMER'S RULE C C CODED BY: GERARD L.H. KRUIZINGA 06/15/1993 C CENTER FOR SPACE RESEARCH UNIVERSITY OF TEXAS AT AUSTIN C C INPUT: C C A 3X3 MATRIX C B 3X1 VECTOR C C OUTPUT: C C X 3X1 SOLUTION X = A-1.B C IMPLICIT NONE DOUBLE PRECISION A(3,3),HELP(3,3),X(3),B(3) DOUBLE PRECISION DETA,DETX,DETY,DETZ INTEGER I,J C... CHECK IF SOLUTION EXITS CALL DETERMINANT(A,DETA) IF (DETA .EQ. 0.0) THEN WRITE(*,*) ' SOLUTION DOES N0T EXIST BECAUSE ' WRITE(*,*) ' DETERMINANT OF : ' CALL WRITEMAT(A,3) WRITE(*,*) ' IS ZERO' RETURN ENDIF DO 20 I=1,3 DO 10 J=1,3 HELP(I,J)=A(I,J) IF (J .EQ. 1) HELP(I,J)=B(I) 10 CONTINUE 20 CONTINUE CALL DETERMINANT(HELP,DETX) X(1) = DETX/DETA DO 40 I=1,3 DO 30 J=1,3 HELP(I,J)=A(I,J) IF (J .EQ. 2) HELP(I,J)=B(I) 30 CONTINUE 40 CONTINUE CALL DETERMINANT(HELP,DETY) X(2) = DETY/DETA DO 60 I=1,3 DO 50 J=1,3 HELP(I,J)=A(I,J) IF (J .EQ. 3) HELP(I,J)=B(I) 50 CONTINUE 60 CONTINUE CALL DETERMINANT(HELP,DETZ) X(3) = DETZ/DETA RETURN END C***************************************************************** SUBROUTINE DETERMINANT(A,DET) C C PURPOSE: DETERMINE DETERMINANT OF A 3X3 MATRIX C C CODED BY: GERARD L.H. KRUIZINGA 06/15/1993 C CENTER FOR SPACE RESEARCH UNIVERSITY OF TEXAS AT AUSTIN C C INPUT: C C A 3X3 MATRIX C C OUTPUT: C C DET DETERMINANT OF MATRIX A C IMPLICIT NONE DOUBLE PRECISION A(3,3),DET DET = A(1,1)*A(2,2)*A(3,3) + . A(1,2)*A(2,3)*A(3,1) + . A(1,3)*A(2,1)*A(3,2) - . A(1,3)*A(2,2)*A(3,1) - . A(1,1)*A(2,3)*A(3,2) - . A(1,2)*A(2,1)*A(3,3) RETURN END C***************************************************************** SUBROUTINE WRITEMAT (A,N) C C PURPOSE: WRITE OUT A 3X3 MATRIX A C C CODED BY: GERARD L.H. KRUIZINGA 06/15/1993 C CENTER FOR SPACE RESEARCH UNIVERSITY OF TEXAS AT AUSTIN C IMPLICIT NONE DOUBLE PRECISION A(N,N) INTEGER N,IROW,JCOL CHARACTER*1 NUMB(0:9) CHARACTER*80 FORMAT DATA NUMB /'0','1','2','3','4','5','6','7','8','9'/ IF (N .GT. 9) THEN WRITE(*,*) ' SUBROUTINE WRITEMAT CAN ONLY PRINT ARRAYS' WRITE(*,*) ' WITH DIMENSION LESS OR EQUAL THAN 9' WRITE(*,*) ' CURRENT DIMENSION : ',N RETURN ENDIF FORMAT(1:1)='(' FORMAT(2:2)=NUMB(N) FORMAT(3:7)='F9.3)' WRITE(6,*) ' ' DO 10 IROW=1,N WRITE(6,FORMAT)(A(IROW,JCOL),JCOL=1,N) 10 CONTINUE WRITE(6,*) ' ' RETURN END C***************************************************************** SUBROUTINE LEAPSEC(TREC,TORB,EPOCHREC,DY2,DY3,DY4, & A,B,C,FLAG,DTNEW) C C PURPOSE: DETERMINE IF INPUT TIME IS CLOSE TO A TIME C WHEN A LEAPSECOND WAS INSERTED C C CODED BY: GERARD L.H. KRUIZINGA 06/16/1993 C CENTER FOR SPACE RESEARCH C C NOTE: LEAP SECONDS UP TO DATE UNTIL JULY 1993 C IMPLICIT NONE INTEGER NLEAP PARAMETER (NLEAP=19) DOUBLE PRECISION TREC,TORB,EPOCHREC,MJDLIM(NLEAP),TIME DOUBLE PRECISION EPSLEAP,TIME1,TIME2,TPREV,DTNEW DOUBLE PRECISION TRECMJD DOUBLE PRECISION A,B,C,DY2,DY3,DY4,DT1,DT2,DT3 DOUBLE PRECISION X(3),Y(3),MAT(3,3) INTEGER INDEX LOGICAL FLAG,FIRST DIMENSION TORB(4) DATA MJDLIM /41317.0E+0, 41499.0E+0, 41683.0E+0, . 42048.0E+0, 42413.0E+0, 42778.0E+0, 43144.0E+0, . 43509.0E+0, 43874.0E+0, 44239.0E+0, 44786.0E+0, . 45151.0E+0, 45516.0E+0, 46247.0E+0, 47161.0E+0, . 47892.0E+0, 48257.0E+0, 48804.0E+0, 49169.0E+0 . / DATA EPSLEAP /8640.0D0/ DATA INDEX /0/ DATA TPREV /0.0D0/ DATA FIRST /.TRUE./ C DETERMINE THE LEAP-SECOND ARRAY INDEX WHICH IS AHEAD FOR C THE INPUT TIME WHEN LEAPSEC IS CALLED FOR THE FIRST TIME. C THIS MEANS THAT THE FIRST TIME OF THE INTERPOLATION WINDOW C MUST BE .GE. TIME OF LEAP + 10.0 SEC IN UTC IF (FIRST .OR. (TPREV - TREC) .GT. 0) THEN TIME1 = TORB(1)/86400.0D0 + EPOCHREC 5 IF (INDEX .GE. NLEAP) GOTO 10 INDEX = INDEX + 1 IF (TIME1 .GE. MJDLIM(INDEX) + 1.0/EPSLEAP) GOTO 5 10 CONTINUE FIRST = .FALSE. ENDIF FLAG = .FALSE. TPREV = TREC TIME2 = TORB(4)/86400.0D0 + EPOCHREC C CHECK IF LAST TIME OF INTERPOLATION WINDOW IS .LE. THAN C TIME OF LEAP-SECOND IF (TIME2 .LE. MJDLIM(INDEX)) RETURN C CHECK IF FIRST TIME OF INTERPOLATION WINDOW IS .GE. THAN C TIME OF LEAP-SECOND + 10.0 IN UTC TIME1 = TORB(1)/86400.0D0 + EPOCHREC IF (TIME1 .GE. MJDLIM(INDEX) + 1.0/EPSLEAP) THEN INDEX = INDEX + 1 IF (INDEX .GT. NLEAP) INDEX = NLEAP RETURN ENDIF C INTERPOLATION AT THIS POINT CONTAINS A 11 SECOND INTERVAL C NOTIFY USER OF OCCURENCE OF THE LEAP SECOND WRITE(*,*) ' ' WRITE(*,*) ' LEAP SECOND ENCOUNTERED !!!!' WRITE(*,*) ' TIME OF LEAP-SECOND (MJD) : ',MJDLIM(INDEX) WRITE(*,*) ' FIRST TIME OF INTERVAL(MJD): ',TIME1 WRITE(*,*) ' LAST TIME OF INTERVAL(MJD) : ',TIME2 WRITE(*,*) ' PLEASE VERIFY RESULTS !!!!!' WRITE(*,*) ' ' C DETERMINE IN WHICH INTERVAL THE 11 SECONDS OCCURS TRECMJD = TREC/86400.0D0 + EPOCHREC C INTERVAL BETWEEN TORB(1),TORB(2) TIME = TORB(2)/86400.0D0 + EPOCHREC IF ( MJDLIM(INDEX) .LE. TIME) THEN DT1 = 11.0 DT2 = 21.0 DT3 = 31.0 GOTO 20 ENDIF C INTERVAL BETWEEN TORB(2),TORB(3) TIME = TORB(3)/86400.0D0 + EPOCHREC IF ( MJDLIM(INDEX) .LE. TIME) THEN DT1 = 10.0 DT2 = 21.0 DT3 = 31.0 GOTO 20 ENDIF C INTERVAL BETWEEN TORB(3),TORB(4) DT1 = 10.0 DT2 = 20.0 DT3 = 31.0 GOTO 20 20 CONTINUE IF (TRECMJD .GT. MJDLIM(INDEX)) THEN DTNEW = TREC - TORB(1) + 1.0 ELSE DTNEW = TREC - TORB(1) ENDIF MAT(1,1) = DT1*DT1*DT1 MAT(2,1) = DT2*DT2*DT2 MAT(3,1) = DT3*DT3*DT3 MAT(1,2) = DT1*DT1 MAT(2,2) = DT2*DT2 MAT(3,2) = DT3*DT3 MAT(1,3) = DT1 MAT(2,3) = DT2 MAT(3,3) = DT3 Y(1)=DY2 Y(2)=DY3 Y(3)=DY4 CALL CRAMER(MAT,Y,X) A = X(1) B = X(2) C = X(3) FLAG = .TRUE. RETURN END