DOUBLE PRECISION TRISE, TCA, TSET, $ TNODE, TSTART, TSTOP, $ TNEXT, EQXYDA, DMOD DIMENSION ISAT(10), PER(10), WMTN(10), $ EQXYDA(10), EQXMIN(10), EQXLON(10), $ TRISE(10), TCA(10), TSET(10), $ TNODE(10), IWEEK(2,7), MONTH(2,12), $ IECA(10), ICADIR(10), IDIR(10) DIMENSION ITCLOS(5), ITSET(5), IP(5), $ ISIDE(10), ITIME(5), ITRISE(5), $ JYEAR(10) EQUIVALENCE (IP(1),LUCON), $ (IP(2),LULIST), $ (IP(3),IONLIN), $ (IP(4),NSAT) DATA $ IWEEK /2HMO,2HN ,2HTU,2HE ,2HWE,2HD ,2HTH,2HU , $ 2HFR,2HI ,2HSA,2HT ,2HSU,2HN /, $ MONTH /2HJA,2HN ,2HFE,2HB ,2HMA,2HR ,2HAP,2HR , $ 2HMA,2HY ,2HJU,2HN ,2HJU,2HL ,2HAU,2HG , $ 2HSE,2HP ,2HOC,2HT ,2HNO,2HV ,2HDE,2HC / C$ C NAME ALERT C C TYPE PROGRAM C C PURPOSE COMPUTE TRANSIT SATELLITE ALERTS FOR UP TO 10 SATELLITES C C AUTHOR/DATE D.WELLS 800130 (DEFAULT DATA 800124) C C EXTERNALS RMPAR,EXEC,FMOD,FLOAT,UPDAT,NEXT,MOD C C CALLING RU,ALERT,LUCON,LULIST,IONLIN,NSAT C C PARAMETERS LUCON = CONSOLE LU C LULIST = LISTING LU (DEFAULT = LUCON) C IONLIN = 0 FOR OFFLINE ALERTS (DEFAULT) C 1 FOR ONLINE ALERTS C NSAT = NUMBER OF SATELLITES (DEFAULT = 6) C C LANGUAGE FORTRAN C C COMMENTS PROGRAM FIRST CHECKS TO SEE IF FILE "NODAL" C EXISTS, AND IF SO, WHETHER IT CONTAINS UP-TO-DATE C NODAL CROSSING DATA (PLACED THERE BY PROGRAM SATFX). C IF NOT, DEFAULT DATA FROM DATA STATEMENT IS USED. C C ONLINE MODE COMPUTES NEXT 10 ALERTS FOR CURRENT POSITION C C OFFLINE MODE LETS OPERATOR SELECT POSITION, START AND C STOP TIMES. C C PROGRAM "ACTIV" SCHEDULES ALERT IN ONLINE MODE ONCE C PER HOUR. C C C**********PROGRAM ALERT INPUT********** C C COLUMNS FORMAT NAME DESCRIPTION C _________________________________________ C C CARD 1 C 1-3 I3 MAXEL MAX. ELEV. USED TO EDIT ALERTS C 4-6 I3 MINEL MIN. ELEV. USED TO EDIT ALERTS C 7-9 I3 IOPTEL OPTIMAL ELEV. FOR ALERTS C C CARD 2 C C 1-4 I4 LUCON INPUT LOGICAL UNIT NO.(DEFAULT=5) C 5-8 I4 LULIST OUTPUT LOGICAL UNIT NO.(DEFAULT=6) C 9-12 I4 IONLIN =0 FOR OFFLINE ALERTS C =1 FOR ONLINE ALERTS C 13-16 I4 NSAT NUMBER OF SATELLITES C C C CARD 3 C C 1-10 I10 IYEAR YEAR FOR WHICH ALERTS REQUIRED C 11-20 F10.4 TSTAX FIRST DAY FOR WHICH ALERTS REQUIRED C 21-30 F10.4 TENDX LAST DAY FOR WHICH ALERTS REQUIRED C (DEFAULT=TSTAX) C 31-40 F10.4 OBSLAX LATITUDE OF OBSERVING STATION C 41-50 F10.4 OBSLOX LONGITUDE OF OBSERVING STATION C C C CARD 4 + NOTE: ONE CARD OF THIS TYPE IS NEEDED FOR EACH SATELLITE C FOR WHICH ALERTS ARE TO BE COMPUTED C C 1-10 F10.4 SATX SATELLITE NUMBER C 11-20 F10.4 PERX SATELLITE ORBITAL PERIOD(MINUTES) C 21-30 F10.4 WMNTX SATELLITE WESTWARD MOTION(DEGREES) C 31-40 F10.4 EQXDAX DATE OF OBSERVATION (YEAR . DAY) C 41-50 F10.4 EQXMIX TIME OF EQUATOR CROSSING(MINUTES) C 51-60 F10.4 EQXLOX LONGITUDE OF EQUATOR CROSSING C (DEGREES (EAST)) C C C*********************************************************************** C C$ C GET EDIT CRITERIA READ(5,*) MAXEL,MINEL,IOPTEL C GET RUN PARAMETERS READ(5,*)(IP(I),I=1,4) IF(LUCON.EQ.0)LUCON=5 IF(LULIST.EQ.0)LULIST=6 C READ FILE NODAL FOR INPUT CALL NODAL(LUCON,LULIST,IONLIN,NSAT, $ ISAT,PER,WMTN,EQXYDA,EQXMIN,EQXLON, $ OBSLAT,OBSLON,IYEAR,TSTART,TSTOP) C FOR OFFLINE MODE CONFIRM INPUT WITH OPERATOR IF(IONLIN .NE. 0) GO TO 5 CALL OPERA(LUCON,LULIST,NSAT, $ ISAT,PER,WMTN,EQXYDA,EQXMIN,EQXLON, $ OBSLAT,OBSLON,IYEAR,TSTART,TSTOP) 5 CONTINUE OBSLON = FMOD(OBSLON,360.,N) NALERT = 0 LBLDAY = TSTART - 1D0 C COMPUTE TNODE IN DAYS DO 10 I = 1,NSAT JYEAR(I) = EQXYDA(I) EQXDAY = DMOD(EQXYDA(I),1D0) * 1000D0 CALL XYEAR(IYEAR,JYEAR(I),EQXDAY) 10 TNODE(I) = EQXDAY + EQXMIN(I) / 1440D0 C FILL ALERT BUFFER DO 20 I = 1,NSAT NORB = (TSTART - TNODE(I)) * 1440D0 / PER(I) - 1 CALL UPDAT(LULIST,PER(I),WMTN(I),TNODE(I),EQXLON(I),NORB) 15 CALL NEXT(LULIST,OBSLAT,OBSLON,MAXEL,MINEL, $ PER(I),WMTN(I),TNODE(I),EQXLON(I), $ TRISE(I),TCA(I),TSET(I), $ IECA(I),ICADIR(I),IDIR(I),ISIDE(I)) IF(TCA(I) .LT. TSTART) GO TO 15 20 CONTINUE C FIND EARLIEST ALERT IN BUFFER 30 TNEXT = TCA(1) J = 1 DO 40 I = 2,NSAT IF(TCA(I) .GT. TNEXT) GO TO 40 TNEXT = TCA(I) J = I 40 CONTINUE C QUIT AFTER TSTOP IF(TCA(J) .GT. TSTOP) GO TO 999 C LABEL EACH DAYS ALERTS IF(IDINT(TCA(J)) .LE. LBLDAY) GO TO 45 LBLDAY = LBLDAY + 1 ND = LBLDAY IY = IYEAR CALL DATE(IY,ND,MNUM,MDAY,IW) OUTLON = OBSLON IF(OUTLON .GT. 180.) OUTLON = OUTLON - 360. C COMPUTE RISE/CA/SET HR/MIN 45 CALL TIME(TRISE(J),ITRISE) CALL TIME(TCA(J),ITCLOS) CALL TIME(TSET(J),ITSET) C WRITE LATEST ALERT WRITE(LULIST,1002) ISAT(J),IECA(J),ND,(ITRISE(7-I),I=3,4), $ (ITSET(7-I),I=3,4) C REFILL EMPTY BUFFER SLOT CALL NEXT(LULIST,OBSLAT,OBSLON,MAXEL,MINEL, $ PER(J),WMTN(J),TNODE(J),EQXLON(J), $ TRISE(J),TCA(J),TSET(J), $ IECA(J),ICADIR(J),IDIR(J),ISIDE(J)) C CHECK FOR MAXIMUM OF 10 ALERTS ONLINE NALERT = NALERT + 1 IF(IONLIN .NE. 0 .AND. NALERT .GE. 10) GO TO 999 GO TO 30 1002 FORMAT(1X,I2,1X,I2,1X,I3,2(1X,I2,1X,I2)) 999 STOP END SUBROUTINE NODAL(LUCON,LULIST,IONLIN,NSAT, $ ISAT,PER,WMTN,EQXYDA,EQXMIN,EQXLON, $ OBSLAT,OBSLON,IYEAR,TSTART,TSTOP) DOUBLE PRECISION TSTART,TSTOP,DAY,TSTAX,EQXYDA,EQXDAX DIMENSION ISAT(10), PER(10), WMTN(10), $ EQXYDA(10), EQXMIN(10), EQXLON(10), $ NPROG(3), NFILE(3), IDCB(144), $ IBUF(30), ITIME(6), IBLANK(30) DATA NPROG /2HAL,2HER,2HT /, $ NFILE /2HNO,2HDA,2HL /, $ IBLANK /30*2H / C$ C NAME NODAL C C TYPE SUBROUTINE C C PURPOSE READ FILE "NODAL" IF IT EXISTS C C CALLING CALL NODAL(LUCON,LULIST,IONLIN,NSAT, C ISAT,PER,WMTN,EQXYDA,EQXMIN,EQXLON, C OBSLAT,OBSLON,TSTART,TSTOP) C$ C CHECK WHETHER FILE NODAL EXISTS READ(5,*) IYEAR,TSTAX,TENDX,OBSLAX,OBSLOX CALL TIME(TSTAX,ITIME) OBSLAT = OBSLAX OBSLON = OBSLOX C READ REST OF FILE (LATEST NODAL PARAMETERS) DO 30 I=1,NSAT READ(5,*) SATX,PERX,WMTNX,EQXDAX,EQXMIX,EQXLOX J=I ISAT(J)=SATX PER(J) = PERX WMTN(J) = WMTNX EQXYDA(J) = EQXDAX EQXMIN(J) = EQXMIX EQXLON(J) = EQXLOX 30 CONTINUE C OFFLINE- COMPUTE 5 DAYS ALERTS FROM LAST FIX TIME ITIME(1) = 0 ITIME(2) = 0 ITIME(3) = 0 TSTART = DAY(ITIME) TSTOP = TSTART + 7D0 / 24D0 IF(IONLIN.EQ.0) TSTOP = TENDX + 1.D0 IF(TENDX.EQ.0.) TSTOP = TSTART + 1.D0 RETURN END SUBROUTINE OPERA(LUCON,LULIST,NSAT, $ ISAT,PER,WMTN,EQXYDA,EQXMIN,EQXLON, $ OBSLAT,OBSLON,IYEAR,TSTART,TSTOP) DOUBLE PRECISION TSTART, TSTOP, EQXYDA DIMENSION ISAT(10),PER(10), WMTN(10), $ EQXYDA(10),EQXMIN(10),EQXLON(10) C$ C NAME OPERA C C TYPE SUBROUTINE C C PURPOSE INPUT-CHECKING DIALOGUE WITH OPERATOR C C CALLING CALL OPERA(LUCON,LULIST,NSAT, C OBSLAT,OBSLON,TSTART,TSTOP) C$ C ECHO INPUT ON LULIST IDAY = TSTART IHR = (TSTART - IDAY) * 24D0 JDAY = TSTOP JHR = (TSTOP - JDAY) * 24D0 RETURN END SUBROUTINE NEXT(LU,OBSLAT,OBSLON,MAXEL,MINEL, $ PER,WMTN,TNODE,EQXLON, $ TRISE,TCA,TSET, $ IECA,ICADIR,IDIR,ISIDE) DOUBLE PRECISION TNODE,TRISE,TCA,TSET DIMENSION INS(2),IEW(2) DATA PI/3.14159/, $ REARTH /6365./, $ GRAV /398600.5/, $ INS /2HNS,2HSN/, $ IEW /2HW ,2HE / C$ C NAME NEXT C C TYPE SUBROUTINE C C PURPOSE COMPUTE ALERT PARAMETERS FOR NEXT PASS OF ONE SATELLITE C C AUTHOR/DATE D.WELLS 790517 C C EXTERNALS COS,SIN,ATAN,ATAN2 C C CALLING CALL NEXT(LU,OBSLAT,OBSLON, C PER,WMTN,TNODE,EQXLON, C TRISE,TCA,TSET,IECA,ICADIR,IDIR,ISIDE) C C PARAMETERS LU = LISTING LU C OBSLAT = OBSERVER'S LATITUDE (DEGREES NORTH) C OBSLON = OBSERVER'S LONGITUDE (DEGREES EAST) C MAXEL = MAXIMUM ACCEPTABLE PASS ELEVATION C MINEL = MINIMUM ACCEPTABLE PASS ELEVATION C PER = SATELLITE ORBITAL PERIOD (MINUTES) C WMTN = SATELLITE WESTWARD MOTION (DEGREES) C TNODE = NODAL CROSSING TIME (DAYS) C EQXLON = NODAL CROSSING LONGITUDE (DEGREES EAST) C TRISE = RISE TIME (DAYS) C TCA = CLOSEST APPROACH TIME (DAYS) C TSET = SET TIME (DAYS) C IECA = ELEVATION AT TCA (DEGREES) C ICADIR = SATELLITE DIRECTION AT TCA (DEGREES) C IDIR = SATELLITE DIRECTION (1 = NS, 2 = SN) C ISIDE = SATELLITE SIDE (1 = WEST, 2 = EAST) C C LANGUAGE FORTRAN C$ C USE KEPLER'S THIRD LAW TO GET ORBIT RADIUS FROM PERIOD RORBIT = ((PER * 60. / (2. * PI))**2 * GRAV) ** (1./3.) COSAMX = REARTH / RORBIT SOBLAT = SIN(OBSLAT*PI/180.) COBLAT = COS(OBSLAT*PI/180.) TOBLAT = SOBLAT / COBLAT C COMPUTE SATELLITE POSITION AT TCA (SUBLAT,SUBLON) 10 IVIS = 0 DLONG = FMOD(EQXLON-OBSLON,360.,N) SDLON = SIN(DLONG*PI/180.) CDLON = COS(DLONG*PI/180.) TDLON = SDLON / CDLON SUBLAT = OBSLAT IF(ABS(DLONG) .GE. 90. .AND. ABS(DLONG) .LE. 270.) $ SUBLAT = 180. - OBSLAT DO 20 ITER = 1,10 SUB2 = SUBLAT SUBLON = FMOD(DLONG - WMTN * SUBLAT / 360., 360.,N) SLON = SIN(SUBLON*PI/180.) CLON = COS(SUBLON*PI/180.) SUBLAT = ATAN2(TOBLAT + WMTN * SLON / 360.,CLON)*180. / PI IF(ABS(SUBLAT - SUB2) .LT. 0.1) GO TO 25 20 CONTINUE GO TO 40 C IS PASS VISIBLE AT TCA 25 SLAT = SIN(SUBLAT*PI/180.) CLAT = COS(SUBLAT*PI/180.) CALF = CLAT * CLON * COBLAT + SLAT * SOBLAT IF(CALF .LT. COSAMX) GO TO 40 C VISIBLE PASS FOUND IVIS = 1 TCA = TNODE + (PER / 1440D0) * (SUBLAT / 360D0) DURAT = (PER / 1440D0) * (ARCOS(COSAMX / CALF) / 360.) TRISE = TCA - DURAT TSET = TCA + DURAT SALF = SQRT(1. - CALF **2) IECA = 90 AZCA= 0. IF( ABS(SALF) .LT. 1.E-5) GO TO 35 IECA = ATAN((CALF - COSAMX) / SALF) * 180. / PI IF (IECA .LT. MINEL .OR. IECA .GT. MAXEL) THEN IVIS = 0 GO TO 40 END IF AZCA = ARCOS((SLAT - SOBLAT * CALF) / (COBLAT * SALF)) 35 AZCA = (AZCA+ 45.*(1.-SIGN(1.,OBSLAT)))*SIGN(1.,SDLON/CDLON) ICADIR = AZCA - 90. * SIGN(1.,SDLON) IF(ICADIR .LT. 0) ICADIR = ICADIR + 360 IDIR = (IFIX(CDLON/ABS(CDLON)) + 3) / 2 ISIDE = (IFIX(TDLON/ABS(TDLON)) + 3) / 2 IDIR = INS(IDIR) ISIDE = IEW(ISIDE) C UPDATE TNODE,EQXLON 40 CALL UPDAT(LU,PER,WMTN,TNODE,EQXLON,1) C IF NO ALERT TRY NEXT PASS IF(IVIS .EQ. 0) GO TO 10 RETURN END SUBROUTINE UPDAT(LU,PER,WMTN,TNODE,EQXLON,NORB) DOUBLE PRECISION TNODE C$ C NAME UPDATE C C TYPE SUBROUTINE C C PURPOSE UPDATE TNODE,EQXLON BY NORB ORBITS C C AUTHOR/DATE D.WELLS 790517 C C CALLING CALL UPDAT(LU,PER,WMTN,TNODE,EQXLON,NORB) C C LANGUAGE FORTRAN C$ TNODE = TNODE + NORB * PER / 1440D0 EQXLON = FMOD(EQXLON - NORB * WMTN,360.,N) RETURN END FUNCTION ARCOS(COSA) DATA PI/3.14159/ C$ C NAME ARCOS C C TYPE FUNCTION C C PURPOSE RETURNS ARC COSINE IN DEGREES C C CALLING A = ARCOS(COSA) C$ ARCOS = 90. IF(COSA .EQ. 0.) RETURN ARCOS = 0. IF(ABS(COSA) .GE. 1.) RETURN ARCOS = ATAN(SQRT(1. - COSA**2) / COSA) * 180. / PI RETURN END FUNCTION FMOD(VALUE,CYCLE,N) C$ C NAME FMOD C C TYPE FUNCTION C C PURPOSE GENERALIZED MOD FUNCTION C C CALLING X = FMOD(VALUE,CYCLE,N) C C PARAMETERS VALUE = INPUT VALUE C CYCLE = INPUT MODULUS RANGE C N = OUTPUT NUMBER OF CYCLES TO NORMALIZE VALUE C FMOD = OUTPUT NORMALIZED VALUE, C IN THE RANGE (0,CYCLE) IF CYCLE POSITIVE C IN THE RANGE (CYCLE,0) IF CYCLE NEGATIVE C$ N = VALUE / CYCLE FMOD = VALUE - N * CYCLE IF(FMOD * CYCLE .GE. 0.) RETURN FMOD = FMOD + CYCLE N = N - 1 RETURN END FUNCTION JERR(LUCON,NPROG,NFILE,ID,IER) DIMENSION NPROG(3),NFILE(3) C$ C NAME JERR C C TYPE FUNCTION C C PURPOSE PRINT FMP ERROR MESSAGES C C CALLING I = JERR(LUCON,NPROG,NFILE,ID,IER) C C PARAMETERS LUCON = CONSOLE LU C NPROG(3) = NAME OF CALLING PROGRAM C NFILE(3) = NAME OF DISC FILE WITH FMP ERROR C ID = FMP OPERATION RESULTING IN ERROR C IER = ERROR CODE C$ JERR = IER IF(IER .GE. 0) RETURN WRITE(LUCON,1001) NPROG,NFILE,IER,ID RETURN 1001 FORMAT(1X,3A2,' - FILE ',3A2,' HAS ERROR ',I6,' ON ',A2) END SUBROUTINE XYEAR(IYEAR,JYEAR,XDAY) C$ C NAME XYEAR C C TYPE SUBROUTINE C C PURPOSE GIVEN XDAY REFERRED TO JYEAR, REFER IT TO IYEAR C C CALLING CALL XYEAR(IYEAR,JYEAR,XDAY) C$ IF(IYEAR .EQ. JYEAR) RETURN NYEAR = IABS(IYEAR - JYEAR) KYEAR = JYEAR IF(IYEAR .LT. JYEAR) KYEAR = IYEAR IF(IYEAR .LT. JYEAR) XDAY = -XDAY DO 10 I = 1,NYEAR XDAY = XDAY + 365. IF(MOD(KYEAR + I - 1,4) .EQ. 0) XDAY = XDAY + 1. 10 CONTINUE RETURN END SUBROUTINE DATE(IY,ND,M,ID,IW) DIMENSION MONTH(12) DATA MONTH /31,28,31,30,31,30,31,31,30,31,30,31/ C$ C NAME DATE C C TYPE SUBROUTINE C C PURPOSE GIVEN YEAR AND DAY OF YEAR, C FIND MONTH, DAY OF MONTH, DAY OF WEEK C$ C CHECK LEAP YEAR 10 MONTH(2) = 28 IF(MOD(IY,4) .EQ. 0) MONTH(2) = 29 LEAP = 0 IF(MOD(IY,4) .EQ. 0) LEAP = 1 C MODIFY IY AND ND IF ND PAST END OF YEAR IF(ND .LE. 365 + LEAP) GO TO 20 IY = IY + 1 ND = ND - 365 - LEAP GO TO 10 C MODIFY IY AND ND IF ND BEFORE START OF YEAR 20 IF(ND .GT. 0) GO TO 30 IY = IY - 1 ND = ND + 365 IF(MOD(IY,4) .EQ. 0) ND = ND + 1 GO TO 10 C SUBTRACT MONTHS UNTIL ND NEGATIVE. C THEN NUMBER OF SUBTRACTIONS = NUMBER OF CURRENT MONTH 30 ND1 = ND DO 40 I = 1,12 IM = I ND1 = ND1 - MONTH(IM) IF(ND1 .LE. 0) GO TO 50 40 CONTINUE 50 M = IM C COMPUTE DAY OF WEEK. DAY ZERO FOR 1971 IS THU (IW=4) ID = ND1 + MONTH(IM) JY = 1971 IW = 4 + (IY - JY) IF(IY .GT. JY) IW = IW + (IY - JY + 2) / 4 IF(IY .LT. JY) IW = IW + (IY - JY - 1) / 4 IW = IW + ND IW = MOD(IW-1,7) + 1 RETURN END SUBROUTINE TIME(DAYS,ITIME) DOUBLE PRECISION DAYS,DBLE,SCALE,XDAY DIMENSION ITIME(5),SCALE(4) DATA SCALE/100D0,2*60D0,24D0/ C C NAME TIME C C TYPE SUBROUTINE C C PURPOSE CONVERTS TIME IN DAYS TO RTE TIME VECTOR C C AUTHOR/DATE D. WELLS 780605 C C EXTERNALS DBLE,FLOAT,IDINT C C CALLING CALL TIME(DAYS,ITIME) C C PARAMETERS DAYS = DOUBLE PRECISION TIME IN DAYS C ITIME(5) = EQUIVALENT RTE TIME VECTOR C C LANGUAGE FORTRAN C XDAY = DAYS DO 5 I=1,4 J=5-I ITIME(J+1) = IDINT(XDAY) XDAY=(XDAY-DBLE(FLOAT(ITIME(J+1)))) * SCALE(J) 5 CONTINUE ITIME(1)=IDINT(XDAY) RETURN END DOUBLE PRECISION FUNCTION DAY(ITIME) DOUBLE PRECISION DBLE,SCALE,XDAY DIMENSION ITIME(5),SCALE(4) DATA SCALE/100D0,2*60D0,24D0/ C C NAME DAY C C TYPE DOUBLE PRECISION FUNCTION C C PURPOSE CONVERT RTE TIME VECTOR TO TIME IN DAYS C C AUTHOR/DATE D. WELLS 780605 C C EXTERNALS DBLE,FLOAT C C CALLING X = DAY(ITIME) C C PARAMETERS ITIME(5) = RTE TIME VECTOR C DAY = EQUIVALENT TIME IN DAYS C C LANGUAGE FORTRAN C C XDAY = 0.1D0 DO 5 I=1,4 5 XDAY=(XDAY+DBLE(FLOAT(ITIME(I)))) / SCALE(I) XDAY=XDAY+DBLE(FLOAT(ITIME(5))) DAY = XDAY RETURN END