' AMSAT ORBITAL PREDICTION PROGRAM de W3IWI - May 1980
' COPYRIGHT 1980 by Dr. Thomas A. Clark, W3IWI
'                   6388 Guilford Road
'                   Clarksville, MD. 21029

' REVISED & MODIFIED FOR IBM-PC by R. D. Welch, W0SL - May 1983
'                                  908 Dutch Mill Dr.
'                                  Manchester, Mo. 63011

' ENHANCED AND DEBUGGED BY Ing. H.F.STRASSER, OE1HSI - JAN. 1985
'                          A 1238
'                          VIENNA, AUSTRIA

' Modified to use the NASA/NORAD two-line element set format (as distributed
' by Dr. Kelso of the Air Force Institute of Technology) and revised for
' QBASIC by Antonio Querubin, AH6BW - Dec. 1991
'    Internet:  tony@mpg.phys.hawaii.edu
'    AMPRnet:  ah6bw@uhm.ampr.org
'    PBBS:  ah6bw@ah6bw.hi.usa.oc
'    BITnet:  querubin@uhunix

' Permission granted for non-commercial use providing
' credit is given to the author(s), AMSAT and ORBIT Magazine.

DEFDBL A-Z

DECLARE FUNCTION gets$ ()
DECLARE SUB getkey (k$)
DECLARE SUB timezone (d0$, t0$, d1$, t1$, tzcor!)
DECLARE FUNCTION juliantodate$ (day%, year%)
DECLARE FUNCTION datetojulian% (d$)
DECLARE FUNCTION trim$ (i$)
DECLARE FUNCTION readtle% (tlefile$, name$(), epochyear%(), epochjulianday#(), inclination#(), TLERAAN#(), eccentricity#(), TLEArgOfPerigee#(), TLEMA#(), MeanMotion#(), TLERevolution&())
DECLARE SUB readgrnd ()
DECLARE SUB plotloc (longitude!, latitude!, xcoord%, ycoord%)
DECLARE FUNCTION arctan# (y#, x#)

DEF fnyear% = VAL(RIGHT$(d$, 2))
DEF fnhour% = VAL(LEFT$(t$, 2))
DEF fnmonth% = VAL(LEFT$(d$, 2))
DEF fnday% = VAL(MID$(d$, 4, 2))
DEF fnminute% = VAL(MID$(t$, 4, 2))
DEF fnsecond% = VAL(RIGHT$(t$, 2))
DEF fnt$ (d%) = MID$(STR$(d% + 100), 3)
DEF fnleap% (year%) = ((year% MOD 4) = 0)
DEF fnradians (degrees) = degrees / 180# * pi#
DEF fndegrees (radians) = radians / pi# * 180#
DEF fnarcsin (sine) = ATN(sine / SQR(1 - sine * sine))
' Greenwich Mean Sidereal Time
DEF fngmst# (YR%) = (99.6367# - .2387# * (YR% - 1989) + .9856# * INT((YR% - 1989) / 4#)) / 360!

CLEAR

maxsats% = 17
DIM name$(maxsats%)
DIM epochyear%(maxsats%), TLERevolution&(maxsats%), satl(maxsats%, 2)
DIM epochjulianday#(maxsats%), inclination#(maxsats%), TLERAAN#(maxsats%), eccentricity#(maxsats%), TLEArgOfPerigee#(maxsats%)
DIM TLEMA#(maxsats%), MeanMotion#(maxsats%)
DIM BeaconFrequency&(maxsats%), a0(maxsats%)
DIM sat%(7), pkt%(7), kep%(7) ' Image arrays (9 x 5 pixels)

DIM cc(3, 2) ' Used in orbit determination routines
   
' ****** NUMERIC CONSTANTS ******
pi# = 3.141592653589793# ' Value of PI
r0# = 6378.16#           ' Earth's radius
f# = 1# / 298.16#        ' 1/Earth flattening coef.
g0# = 75369793000000#    ' GM of Earth in (orbits/day)^2/km^3
g1# = 1.0027379093#      ' Sidereal/Solar time rate ratio
c# = 299792.5#           ' Speed of light (km/sec).

' String constants
tlefile$ = "USAT.TLE"   ' Default NASA/NORAD element file

' YOU WILL NEED THE FILES USATCGA.MAP, USAT.TLE AND USATGRND.DAT
' TO RUN THIS PROGRAM
' **** SATAUS Menuprogramm de OE1HSI VERSION 1.5  26.JAN. 1985
 
' ****** ESTABLISH SATELLITE ELEMENT MATRIX ******
PRINT "Reading satellite elements from "; tlefile$
numsats% = readtle%(tlefile$, name$(), epochyear%(), epochjulianday#(), inclination#(), TLERAAN#(), eccentricity#(), TLEArgOfPerigee#(), TLEMA#(), MeanMotion#(), TLERevolution&())

' Initialize beacon matrix and check age of satellite elements.
FOR i% = 1 TO numsats%
  BeaconFrequency&(i%) = 0
  IF epochyear%(i%) < fnyear% - 1 THEN
    PRINT "Warning:  Elements for satellite "; trim(name$(i%)); " are very old."
    PRINT "You should update your "; tlefile$; " file."
  END IF
NEXT

' Initialize the ground station data
CALL readgrnd

' ****** DRAW AND STORE SATELLITE SYMBOLS ******
CLS
SCREEN 2
LINE (4, 0)-(4, 4)
LINE (0, 2)-(8, 2)
GET (0, 0)-(8, 4), sat%
PUT (0, 0), sat%
CLS
LINE (4, 1)-(4, 3)
LINE (3, 2)-(5, 2)
GET (0, 0)-(8, 4), pkt%
PUT (0, 0), pkt%

' ***** Initializations done. *****

50
' ****** DERIVED CONSTANTS ******
' These values depend on the ground station data.
' C$=GROUND STATION CALL+LOCATION STRING
c$ = trim$(gsid$) + " " + trim$(gsloc$)
pi2# = 2 * pi#
SinGSLat# = SIN(fnradians(gslat!))    ' sin/cos of ground station latitude
CosGSLat# = COS(fnradians(gslat!))
SinGSLong# = SIN(fnradians(-gslong!))  ' sin/cos of ground station longitude
CosGSLong# = COS(fnradians(gslong!))
R9 = r0# * (1 - (f# / 2) + (f# / 2) * COS(2 * fnradians(gslat!))) + gsheight% / 1000
L8 = ATN((1 - f#) ^ 2 * SinGSLat# / CosGSLat#)
GSzcoord = R9 * SIN(L8)
GSxcoord = R9 * COS(L8) * CosGSLong#
GSycoord = R9 * COS(L8) * SinGSLong#

80
CLS
SCREEN 0
KEY(9) OFF
KEY(10) OFF
PRINT "                        USAT90.BAS - Version 1989-91"
PRINT
PRINT "======================================================================"
PRINT
PRINT "            SELECT ONE OF THE FOLLOWING OPTIONS:"
PRINT
PRINT " (P) ORBITAL PREDICTION PROGRAM"
PRINT
PRINT " (R) REALTIME TRACKING AND HIGH RESOLUTION SCREEN"
PRINT
PRINT " (S) Display ELEMENTS OF SATELLITES"
PRINT
PRINT " (G) CHANGE OR ENTER GROUNDSTATION DATA"
PRINT
PRINT " (D) RETURN TO DOS"
PRINT
PRINT
PRINT "ENTER SELECTION (P,R,C,G,D)--> "
CALL getkey(k$)
ON INSTR("PRSGD", k$) GOTO 4820, 2090, 240, 1620, 9000
BEEP
GOTO 50

240
' ****** SATFILE.BAS - VERSION 1.0, ISSUE 1.0 - HSIMODIF.1/25/85 ******
CLS
DO
  PRINT "Elements of the following satellites are available:"
  PRINT
  FOR j% = 1 TO numsats%
    PRINT name$(j%),
  NEXT
  PRINT
  ' ****** FIND SATELLITE RECORD ******
  DO
    PRINT : INPUT "Which satellite"; v1$: IF v1$ = "" THEN 80
    v1$ = UCASE$(trim$(v1$))
    FOR j% = 1 TO numsats%
      IF UCASE$(trim$(name$(j%))) = v1$ THEN
        CLS
        ' ****** DISPLAY SATELLITE RECORD ******
        PRINT "Satellite       = "; name$(j%)
        PRINT "Epoch year      = "; epochyear%(j%)
        PRINT "Epoch day       = "; epochjulianday#(j%)
        PRINT "Inclination     = "; inclination#(j%)
        PRINT "R.A.A.N.        = "; TLERAAN#(j%)
        PRINT "Eccentricity    = "; eccentricity#(j%)
        PRINT "Arg. of perigee = "; TLEArgOfPerigee#(j%)
        PRINT "Mean anomaly    = "; TLEMA#(j%)
        PRINT "Mean motion     = "; MeanMotion#(j%)
        PRINT "Epoch orbit no. = "; TLERevolution&(j%)
        PRINT "Beacon freq.    = "
        PRINT
        EXIT DO
      END IF
    NEXT
    PRINT "Satellite not found."
  LOOP
LOOP

1620
' ******* Groundsation data change v.1.0 OE1HSI  jan.-1985**********
CLS
PRINT "CURRENT GROUND STATION DATA"
PRINT
GOSUB 2080
DO
  PRINT "Do you want to CHANGE this DATA ? (Y/N)"
  CALL getkey(k$)
  ON INSTR("YN", k$) GOTO 1780, 2060
  BEEP
LOOP
1710
PRINT
PRINT
GOSUB 2080
DO
  PRINT "Do you want to make further changes ? (Y/N) "
  CALL getkey(k$)
  ON INSTR("YN", k$) GOTO 1780, 2050
  BEEP
LOOP
1780
PRINT
PRINT "ENTER NEW DATA OR <RETURN> FOR UNCHANGED DATA":
PRINT
PRINT "Ground Station Identifier                    : ";
LINE INPUT u$
IF trim$(u$) <> "" THEN gsid$ = trim$(u$)
PRINT "Location of station                          : ";
LINE INPUT u$
IF trim$(u$) <> "" THEN gsloc$ = trim$(u$)
PRINT
PRINT "Longitude WEST of Greenwich (max +360) or East of Greenw. entered as -0 to -180"
PRINT
INPUT "Enter (with decimals)                 : ", u$
IF trim$(u$) <> "" THEN gslong! = VAL(u$)
IF gslong! < 0 THEN gslong! = 360 + gslong!
PRINT
PRINT "LATITUDE NORTH of Equator + (max 90) SOUTH of Equator - (max 90)"
PRINT
INPUT "ENTER (With decimals)                 : ", u$
IF trim$(u$) <> "" THEN gslat! = VAL(u$)
PRINT
INPUT "Station height above sea level in meters     : ", u$
IF trim$(u$) <> "" THEN gsheight% = VAL(u$)
PRINT
PRINT "Correction from local computer time to UTC"
PRINT
INPUT "Enter (in hours with decimals)        : ", u$
IF trim$(u$) <> "" THEN tzcor! = VAL(u$)
GOTO 1710
2050
PRINT
PRINT "DATA SAVED AS USATGRND.DAT"
OPEN "USATGRND.DAT" FOR OUTPUT AS #1
PRINT #1, gsid$
PRINT #1, gsloc$
PRINT #1, gslong!, gslat!, gsheight%, tzcor!
CLOSE #1
GOTO 50
2060
PRINT
PRINT "DATA NOT CHANGED"
GOTO 50

2080
  PRINT "Ground Station Identifier            : "; gsid$
  PRINT "Location name                        : "; gsloc$
  PRINT USING "West longitude (degrees)             : +###.##"; gslong!
  PRINT USING "Latitude (degrees)                   :  +##.##"; gslat!
  PRINT USING "Height above Mean Sea Level (meters) :   #####"; gsheight%
  PRINT USING "Local time correction to UTC (hours) :  +##.##"; tzcor!
  PRINT
RETURN
'**** END PROGRAM GROUNDSTATION DATA CHANGE/STORAGE OE1HSI  JAN. 1985 ****

2090
' ****** ORBITS2 - VERSION 1.0, ISSUE 1.2 -11/1/83 ******
KEY OFF
SCREEN 2, 0
CLS

' Initialize D$ and T$ from date$ and time$
CALL timezone(DATE$, TIME$, d$, t$, tzcor!)
' ****** SET UP KEY TRAPPING ******
ON KEY(9) GOSUB 4760
KEY(9) STOP
ON KEY(10) GOSUB 4790
KEY(10) OFF
flg9 = 0
flg10 = 0
GOSUB 4290
   
' ****** ORBIT DETERMINATION LOOP STARTS HERE ******
q$ = ""
DO UNTIL q$ = CHR$(27)
  FOR j% = 1 TO numsats%
    q$ = UCASE$(INKEY$)
    IF q$ = CHR$(27) THEN EXIT FOR
    CALL timezone(DATE$, TIME$, d$, t$, tzcor!)
    t = INT(30.55 * (fnmonth% + 2)) - 2 * (INT(.1 * (fnmonth% + 7))) - 91
    IF fnmonth% > 2 AND fnleap%(fnyear%) THEN t = t + 1
    IF epochyear%(j%) = fnyear% - 1 THEN
      t = t + 365
      IF fnleap%(epochyear%(j%)) THEN t = t + 1
    END IF
    t = t + fnday% + fnhour% / 24# + fnminute% / 1440# + fnsecond% / 86400#
    GOSUB 6780
    GOSUB 6910
    IF flg9 THEN CALL plotloc(longitude!, latitude!, plotx%, ploty%): GOSUB 3890
    GOSUB 4040
  NEXT
LOOP
GOTO 50
   
' ****** PUT SATELLITE ON SCREEN ROUTINE ******
3890
  GET (plotx% - 4, ploty% - 2)-(plotx% + 4, ploty% + 2), kep%
  PUT (plotx% - 4, ploty% - 2), sat%, PRESET
  PUT (plotx% - 4, ploty% - 2), sat%
  PUT (plotx% - 4, ploty% - 2), kep%, PSET
  PUT (plotx% - 4, ploty% - 2), sat%
  IF FLG0 <> 0 THEN
    PUT (satl(j%, 1), satl(j%, 2)), sat%
    PUT (satl(j%, 1), satl(j%, 2)), pkt%, OR
  END IF
  satl(j%, 1) = plotx% - 4
  satl(j%, 2) = ploty% - 2
  IF j% = numsats% THEN FLG0 = 1
RETURN
 
' ****** PRINT SATELLITE DETAILS ROUTINE ******
4040
  KEY(9) ON
  KEY(10) ON
4050
  KEY(10) STOP
  KEY(9) STOP
  IF FLGK = 1 THEN GOSUB 4270
  IF flg9 = 0 THEN
    LOCATE 2, 49
    PRINT d$
    LOCATE 2, 63
    PRINT t$
    IF elevation! >= 0 THEN
      IF elevation! > 0 AND elevation! < 1 THEN BEEP
    END IF
    IF j% + 6 > 23 GOTO 4230
    LOCATE j% + 6, 15
    PRINT SPACE$(50)
    LOCATE j% + 6, 8 - (LEN(trim$(name$(j%))) / 2)
    PRINT trim$(name$(j%))
    LOCATE j% + 6, 15
  ELSE
    IF flg10 = 1 THEN GOSUB 4540
    IF flg10 <> 0 THEN
      IF UCASE$(trim$(name$(j%))) <> u$ THEN RETURN
      LOCATE 25, 69
      PRINT "SELECTED";
    END IF
    LOCATE 25, 1
    PRINT SPACE$(68);
    LOCATE 25, 8 - (LEN(trim$(name$(j%))) / 2)
    PRINT trim$(name$(j%));
    LOCATE 25, 15
  END IF
  PRINT USING "###   ###  #####    #####  ###.#   ###.#    ######"; azimuth!; elevation!; range#; (r - r0#); latitude!; longitude!; revolution&;
4230
  IF flg9 <> 0 THEN LOCATE 20, 48: PRINT t$;
RETURN
   
' ****** SET UP SCREEN DISPLAY ROUTINE ******
4270
  CLS
  FLG0 = 0
  FLGK = 0
4290
  IF flg9 = 1 THEN
    SCREEN 2, 0
    DEF SEG = &HB800
    BLOAD "USATCGA.MAP", 0
    DEF SEG
    CALL plotloc(gslong!, gslat!, plotx%, ploty%)
    CIRCLE (plotx%, ploty%), 2
    LOCATE , , 1, 12, 13
    L1 = 22
    L2 = 23
    L3 = 24
    LOCATE 20, 3
    PRINT "Data for Groundstation             At Time=           UTC  On: "; d$
    LOCATE 20, 26
    PRINT gsid$
  ELSE
    ' FLG9=0
    SCREEN 0, 1
    LOCATE 1, 40 - LEN(c$) / 2, 0
    PRINT c$
    LOCATE 2, 5
    PRINT "Real-time satellite tracking coordinates on            at          UTC"
    LOCATE 25, 3
    PRINT "F9 TOGGLES THE GRAPH/TABLE     F10 TO SELECT SINGLE SAT IN GRAPH   ESC TO END";
    L1 = 4
    L2 = 5
    L3 = 6
  END IF
  LOCATE L1, 3
  PRINT "  NAME OR   AZ    EL   RANGE   HEIGHT   LAT.   LONG.    ORBIT"
  LOCATE L2, 3
  PRINT "DESIGNATOR  DEG   DEG    KM      KM     DEG     DEG       NO."
  LOCATE L3, 3
  PRINT "----------  ---   ---  -----   ------  -----   -----    ------";
RETURN
   
' ****** SELECT SINGLE SATELLITE ROUTINE ******
4540
  LOCATE 25, 1
  PRINT SPACE$(79);
  LOCATE 25, 1
  INPUT ; "WHICH SATELLITE? (CR FOR ALL)"; i$
  i$ = UCASE$(trim$(i$))
  IF i$ = "" THEN
    flg10 = 0
  ELSE
    u$ = i$
    flg10 = 2
  END IF
  LOCATE 25, 1
  PRINT SPACE$(79);
RETURN

4760 ' toggle FLG9
  flg9 = -(flg9 = 0)
  FLGK = 1
RETURN 4050

4790 flg10 = flg9
RETURN 4050

4820
  '****** ORBIT2 - VERSION 2.0, ISSUE 1.0/HSI - 17/01/85 *****
  KEY OFF
  SCREEN 0
  CLS
  ' ****** HOUSEKEEPING ITEMS ******
  C8$ = CHR$(10) + CHR$(10) + CHR$(10) + CHR$(10)
  C9$ = CHR$(12) + CHR$(7)

  DO
    PRINT "Enter start date (mm-dd-yy):  ";
    u$ = gets$
    year% = VAL(RIGHT$(u$, 2))
    month% = VAL(LEFT$(u$, 2))
    day% = VAL(MID$(u$, 4, 2))
    IF month% >= 1 AND month% <= 12 AND day% >= 1 AND day% <= 31 AND year% >= 91 AND MID$(u$, 3, 1) = "-" AND MID$(u$, 6, 1) = "-" THEN EXIT DO
  LOOP
  YY = year% + 1900
  t$ = fnt$(year%) + "/" + fnt$(month%) + "/" + fnt$(day%) + " at "
  D8 = day% + INT(30.55 * (month% + 2)) - 2 * (INT(.1 * (month% + 7))) - 91
  IF month% > 2 THEN IF fnleap%(year%) THEN D8 = D8 + 1
  PRINT "    Day #"; D8
  PRINT
  PRINT "Enter start time (hh:mm):  ";
  DO
    u$ = gets$
    h% = VAL(LEFT$(u$, 2))
    m% = VAL(RIGHT$(u$, 2))
    IF h% >= 0 AND h% <= 23 AND m% >= 0 AND m% <= 59 AND MID$(u$, 3, 1) = ":" THEN EXIT DO
    PRINT "Valid format is:  hh:mm"
    PRINT "Re-enter:  ";
  LOOP
  T7 = D8 + h% / 24# + m% / 1440#
  t$ = t$ + fnt$(h%) + fnt$(m%) + ":00 H"
  PRINT "Enter duration in hours and minutes (hh:mm):  ";
  DO
    u$ = gets$
    h% = VAL(LEFT$(u$, 2))
    m% = VAL(RIGHT$(u$, 2))
    IF h% >= 0 AND m% >= 0 AND m% <= 59 AND MID$(u$, 3, 1) = ":" THEN EXIT DO
    PRINT "Valid format is:  hh:mm"
    PRINT "Re-enter:  ";
  LOOP
  T8 = T7 + h% / 24# + m% / 1440#
  PRINT USING "    From ###.####### to ###.#######"; T7; T8
  DO
    PRINT "Enter time increment (minutes):  ";
    INPUT m%
    IF m% > 0 THEN EXIT DO
  LOOP
  T9 = m% / 1440#
  PRINT
  INPUT "MINIMUM ELEVATION ? (DEFAULT 0) Deg. = ", E8
  DO
    PRINT
    INPUT "Output to Printer (P) or Screen (S) ?-->", P$
    P$ = UCASE$(P$)
    IF (P$ = "P" OR P$ = "S") THEN
      EXIT DO
    ELSE BEEP
    END IF
  LOOP
  IF P$ = "P" THEN outfile$ = "prn:" ELSE outfile$ = "scrn:"
  CLOSE #2
  OPEN outfile$ FOR OUTPUT AS #2
  IF outfile$ = "prn:" THEN
    PRINT "Put the printer on-line and align the page,"
    PRINT "then press any key when ready."
    DO
      getkey (u$)
      IF u$ <> "" THEN EXIT DO
    LOOP
  END IF
 
' ****** SELECT SATELLITE FROM MENU ******
  CLS
  PRINT "SATELLITE SELECTION MENU"
  PRINT
  FOR j% = 1 TO numsats%
    PRINT name$(j%),
  NEXT
  PRINT
  DO
    PRINT
    INPUT "SELECT satellite >", y3$
    y3$ = UCASE$(trim$(y3$))
    FOR j% = 1 TO numsats%
      IF UCASE$(trim$(name$(j%))) = y3$ THEN EXIT DO
    NEXT
    PRINT "Satellite not found."
  LOOP
  PRINT
  PRINT "Doppler calculated for frequ. = "; BeaconFrequency&(j%); " MHz"
  INPUT " Change frequency to: (0 for default) ", d
  IF d <> 0 THEN BeaconFrequency&(j%) = d
  PRINT #2,
  PRINT #2, "Orbital ELEMENTS for "; name$(j%)
  PRINT #2,
  PRINT #2, "Reference epoch = "; epochyear%(j%); " +"; epochjulianday#(j%)
  PRINT #2, "Starting  epoch = "; year%; " +"; T7; " = "; t$
  PRINT #2,
  PRINT #2, "Parameter"; TAB(20); "Reference"; TAB(40); "Starting"
  t = T7
  IF epochyear%(j%) = year% - 1 THEN
    t = t + 365
    T8 = T8 + 365
    IF fnleap%(epochyear%(j%)) THEN t = t + 1: T8 = T8 + 1
  END IF
  ' Calculate start time parameters
  GOSUB 6780
  PRINT #2, "Orbit Number "; TAB(20); TLERevolution&(j%); TAB(40); revolution&
  PRINT #2, "Mean Anomaly "; TAB(20); TLEMA#(j%); TAB(40); fndegrees(MARadians#)
  PRINT #2, "Inclination  "; TAB(20); inclination#(j%)
  PRINT #2, "Eccentricity "; TAB(20); eccentricity#(j%)
  PRINT #2, "Mean Motion  "; TAB(20); MeanMotion#(j%)
  PRINT #2, "S.M.A.,km    "; TAB(20); a0(j%)
  PRINT #2, "Arg. Perigee "; TAB(20); TLEArgOfPerigee#(j%); TAB(40); ArgOfPerigee#
  PRINT #2, "R. A. A. N.  "; TAB(20); TLERAAN#(j%); TAB(40); RAAN#
  PRINT #2, "Freq.,MHz    "; TAB(20); BeaconFrequency&(j%)
  OldRevolution& = -1
  OldJulianDay% = -1
  
  '****** COMPUTATION LOOP ******
6400
  t = t + T9
  JulianDay% = INT(t)
  IF JulianDay% <> OldJulianDay% THEN
    OldJulianDay% = -1
    OldRevolution& = -1
  END IF
  GOSUB 6780
  GOSUB 6910
  IF elevation! >= E8 THEN
    IF JulianDay% <> OldJulianDay% OR revolution& <> OldRevolution& THEN
      IF JulianDay% <> OldJulianDay% THEN
        GOSUB 6690
        OldJulianDay% = JulianDay%
        PRINT #2, " U.T.C.    AZ    EL  DOPPLER   RANGE   HEIGHT  LAT.  LONG.  PHASE"
        PRINT #2, "HHMM:SS   deg   deg    Hz       km       km    deg    deg   <256>"
      END IF
      PRINT #2, TAB(21); "- - - ORBIT #"; revolution&; "- - -"
    END IF
    OldRevolution& = revolution&
    T4 = t - JulianDay%
    S4 = INT(T4 * 86400)
    H4 = INT(S4 / 3600 + .000001)
    M4 = INT((S4 - H4 * 3600) / 60 + .000001)
    S4 = S4 - 3600 * H4 - 60 * M4
    doppler# = -BeaconFrequency&(j%) * 1000000 * R8 / c#
    PRINT #2, fnt$(H4) + fnt$(M4) + ":" + fnt$(S4);
    PRINT #2, USING "   ###   ###  #####"; azimuth!; elevation!; doppler#;
    PRINT #2, USING "    #####    #####  ###.#  ###.#  ###"; range#; (r - r0#); latitude!; longitude!; phase%
  END IF
  IF t < T8 THEN GOTO 6400
  PRINT #2, C9$
  PRINT "Do YOU have another INQUIRY  ?  (Y/N) "
  PRINT
  PRINT "Else you return to the MAIN MENU !"
  PRINT
  DO
    CALL getkey(k$)
    ON INSTR("YN", k$) GOTO 4820, 6670
    BEEP
  LOOP
6670
  CLOSE #2
GOTO 80
    
'****** PAGE HEADER SUBROUTINE ******
6690
  PRINT #2, C9$; c$; "   Lat.="; gslat!; "  W.Long.="; gslong!; "  Ht.="; gsheight%;
  pagenum% = pagenum% + 1
  PRINT #2, TAB(70); "Page # "; pagenum%
  PRINT #2, TAB(15); " - - - Minimum Elevation = "; E8; "Deg. - - -"
  PRINT #2,
  PRINT #2, TAB(14); "- - - DAY #"; JulianDay%; "- - -  "; juliantodate$(JulianDay%, INT(YY)); "  - - -"
  PRINT #2,
RETURN
 
'****** ORBIT DETERMINATION AND UTILITY ROUTINES ******

6780
  a0(j%) = ((g0# / (MeanMotion#(j%) * MeanMotion#(j%))) ^ (1 / 3))
  e2 = 1 - eccentricity#(j%) * eccentricity#(j%)
  E1 = SQR(e2)
  TLEOrbitPos# = TLEMA#(j%) / 360 + TLERevolution&(j%)
  K2 = 9.95 * ((r0# / a0(j%)) ^ 3.5) / (e2 * e2)
  S1 = SIN(fnradians(inclination#(j%)))
  C1 = COS(fnradians(inclination#(j%)))
  RAAN# = TLERAAN#(j%) - (t - epochjulianday#(j%)) * K2 * C1
  S0 = SIN(fnradians(RAAN#))
  C0 = COS(fnradians(RAAN#))
  ArgOfPerigee# = TLEArgOfPerigee#(j%) + (t - epochjulianday#(j%)) * K2 * (2.5 * (C1 * C1) - .5)
  S2 = SIN(fnradians(ArgOfPerigee#))
  C2 = COS(fnradians(ArgOfPerigee#))
  cc(1, 1) = (C2 * C0) - (S2 * S0 * C1)
  cc(1, 2) = -(S2 * C0) - (C2 * S0 * C1)
  cc(2, 1) = (C2 * S0) + (S2 * C0 * C1)
  cc(2, 2) = -(S2 * S0) + (C2 * C0 * C1)
  cc(3, 1) = (S2 * S1)
  cc(3, 2) = (C2 * S1)
  OrbitPos# = TLEOrbitPos# + MeanMotion#(j%) * (t - epochjulianday#(j%))
  revolution& = INT(OrbitPos#)
  phase% = INT((OrbitPos# - revolution&) * 256)
  MARadians# = (OrbitPos# - revolution&) * pi2#
RETURN

6910
  E = MARadians# + eccentricity#(j%) * SIN(MARadians#) + .5 * (eccentricity#(j%) * eccentricity#(j%)) * SIN(2 * MARadians#)
  DO
    S3 = SIN(E)
    C3 = COS(E)
    R3 = 1 - eccentricity#(j%) * C3
    M1 = E - eccentricity#(j%) * S3
    M5 = M1 - MARadians#
    IF ABS(M5) < .000001 THEN EXIT DO ELSE E = E - M5 / R3
  LOOP
  X0 = a0(j%) * (C3 - eccentricity#(j%))
  y0 = a0(j%) * E1 * S3
  r = a0(j%) * R3
  X1 = X0 * cc(1, 1) + y0 * cc(1, 2)
  y1 = X0 * cc(2, 1) + y0 * cc(2, 2)
  Z1 = X0 * cc(3, 1) + y0 * cc(3, 2)
  G7 = t * g1# + fngmst#(1900 + epochyear%(j%))
  G7 = (G7 - INT(G7)) * pi2#
  S7 = -SIN(G7)
  C7 = COS(G7)
  x = (X1 * C7) - (y1 * S7)
  y = (X1 * S7) + (y1 * C7)
  z = Z1
  X5 = (x - GSxcoord)
  Y5 = (y - GSycoord)
  Z5 = (z - GSzcoord)
  range# = SQR(X5 * X5 + Y5 * Y5 + Z5 * Z5)
  ' Speed to/from ground station
  IF OldT# <> t THEN R8 = ((range# - OldRange#) / (t - OldT#)) / 86400 ELSE R8 = -900000000
  OldRange# = range#
  OldT# = t
  ' Translate to ground station coordinate system
  z8 = (X5 * CosGSLong# * CosGSLat#) + (Y5 * SinGSLong# * CosGSLat#) + (Z5 * SinGSLat#)
  x8 = -(X5 * CosGSLong# * SinGSLat#) - (Y5 * SinGSLong# * SinGSLat#) + (Z5 * CosGSLat#)
  y8 = (Y5 * CosGSLong#) - (X5 * SinGSLong#)
  elevation! = fndegrees(fnarcsin(z8 / range#))
  azimuth! = fndegrees(arctan(y8, x8))
  longitude! = 360 - fndegrees(arctan(y, x))
  latitude! = fndegrees(fnarcsin(z / r))
RETURN

9000 END

FUNCTION arctan (y, x)
' Modified arctangent function.  Returns a value between 0 and 2 PI.
SHARED pi#, pi2#

IF x > 0 AND y >= 0 THEN
  arctan = ATN(y / x)
ELSEIF x > 0 AND y < 0 THEN
  arctan = ATN(y / x) + pi2#
ELSEIF x < 0 THEN
  arctan = ATN(y / x) + pi#
ELSEIF y >= 0 THEN
  arctan = pi# / 2
ELSE arctan = 1.5# * pi#
END IF
END FUNCTION

DEFSNG A-Z
' Converts a date string to the julian day of the year
FUNCTION datetojulian% (d$)

day% = VAL(MID$(d$, 4, 2))

SELECT CASE LEFT$(d$, 2)
CASE "01"
CASE "02"
  day% = day% + 31
CASE "03"
  days% = day% + 59
CASE "04"
  day% = day% + 90
CASE "05"
  day% = day% + 120
CASE "06"
  day% = day% + 151
CASE "07"
  day% = day% + 181
CASE "08"
  day% = day% + 212
CASE "09"
  day% = day% + 243
CASE "10"
  day% = day% + 273
CASE "11"
  day% = day% + 304
CASE "12"
  day% = day% + 334
END SELECT

IF fnleap%(VAL(RIGHT$(d$, 4))) AND LEFT$(d$, 2) > "02" THEN
  datetojulian% = day% + 1
ELSE
  datetojulian% = day%
END IF

END FUNCTION

SUB getkey (k$)

k$ = ""
WHILE k$ = ""
        k$ = UCASE$(INKEY$)
WEND

END SUB

' Gets an input line, trims leading and trailing blanks
' then converts to upper-case before returning.
FUNCTION gets$
LINE INPUT i$
gets$ = UCASE$(trim$(i$))
END FUNCTION

' Converts a julian day and year to a date string.
FUNCTION juliantodate$ (day%, year%)

IF day% = 60 THEN
  IF fnleap%(year%) THEN tempdate$ = "02-29" ELSE tempdate$ = "03-01"
  GOTO appendyear
ELSE IF fnleap%(year%) AND day% > 60 THEN day% = day% - 1
END IF

SELECT CASE day%
CASE 1 TO 31
  tempdate$ = "01-" + MID$(STR$(day% + 100), 3)
CASE 32 TO 59
  tempdate$ = "02-" + MID$(STR$(day% - 31 + 100), 3)
CASE 60 TO 90
  tempdate$ = "03-" + MID$(STR$(day% - 59 + 100), 3)
CASE 91 TO 120
  tempdate$ = "04-" + MID$(STR$(day% - 90 + 100), 3)
CASE 121 TO 151
  tempdate$ = "05-" + MID$(STR$(day% - 120 + 100), 3)
CASE 152 TO 181
  tempdate$ = "06-" + MID$(STR$(day% - 151 + 100), 3)
CASE 182 TO 212
  tempdate$ = "07-" + MID$(STR$(day% - 181 + 100), 3)
CASE 213 TO 243
  tempdate$ = "08-" + MID$(STR$(day% - 212 + 100), 3)
CASE 244 TO 273
  tempdate$ = "09-" + MID$(STR$(day% - 243 + 100), 3)
CASE 274 TO 304
  tempdate$ = "10-" + MID$(STR$(day% - 273 + 100), 3)
CASE 305 TO 334
  tempdate$ = "11-" + MID$(STR$(day% - 304 + 100), 3)
CASE 335 TO 365
  tempdate$ = "12-" + MID$(STR$(day% - 334 + 100), 3)
END SELECT

appendyear: juliantodate$ = tempdate$ + "-" + MID$(STR$(year%), 2)

END FUNCTION

SUB plotloc (longitude!, latitude!, xcoord%, ycoord%)
' Returns the screen x and y coordinates for a given map location

ycoord% = CINT(.7111 * (90 - latitude!) + 3)
IF longitude! >= 0 AND longitude! <= 270 THEN
  xcoord% = CINT(477 - longitude! * 1.7444)
ELSE xcoord% = CINT(1105 - longitude! * 1.7444)
END IF

END SUB

SUB readgrnd
SHARED gsid$, gsloc$, gslong!, gslat!, gsheight%, tzcor!

' Format of USATGRND.DAT:
'   gsid$    = Station Identifier
'   gsloc$   = Station location (name)
'   gslong!   = Longitude in degrees
'   gslat!    = Latitude in degrees
'   gsheight% = Height above sea level in meters
'   tzcor!    = Time zone correction to UTC
OPEN "USATGRND.DAT" FOR INPUT AS #1
LINE INPUT #1, gsid$
LINE INPUT #1, gsloc$
INPUT #1, gslong!, gslat!, gsheight%, tzcor!
CLOSE #1

END SUB

FUNCTION readtle% (tlefile$, name$(), epochyear%(), epochjulianday#(), inclination#(), TLERAAN#(), eccentricity#(), TLEArgOfPerigee#(), TLEMA#(), MeanMotion#(), TLERevolution&())
' ****** ESTABLISH SATELLITE ELEMENT MATRIX ******
' tlefile$ conforms to the NASA/NORAD three-line element set format.

OPEN tlefile$ FOR INPUT AS #1

FOR i% = 1 TO UBOUND(name$)
  IF EOF(1) THEN EXIT FOR
  LINE INPUT #1, name$(i%)
  IF trim$(name$(i%)) = "" THEN EXIT FOR
  LINE INPUT #1, y3$
  LINE INPUT #1, t0$
  epochyear%(i%) = VAL(MID$(y3$, 19, 2))
  epochjulianday#(i%) = VAL(MID$(y3$, 21, 12))
  inclination#(i%) = VAL(MID$(t0$, 9, 8))
  TLERAAN#(i%) = VAL(MID$(t0$, 18, 8))
  eccentricity#(i%) = VAL("." + MID$(t0$, 27, 7))
  TLEArgOfPerigee#(i%) = VAL(MID$(t0$, 35, 8))
  TLEMA#(i%) = VAL(MID$(t0$, 44, 8))
  MeanMotion#(i%) = VAL(MID$(t0$, 53, 11))
  TLERevolution&(i%) = VAL(MID$(t0$, 64, 5))
NEXT

CLOSE #1
readtle% = i% - 1
           
END FUNCTION

' Converts an input date and time string from one time zone to another.
SUB timezone (d0$, t0$, d1$, t1$, tzcor!)

JulianDay% = datetojulian(d0$)

year% = VAL(RIGHT$(d0$, 4))

hour! = VAL(LEFT$(t0$, 2)) + INT(tzcor!) + (VAL(MID$(t0$, 4, 2)) + tzcor! - INT(tzcor!)) / 60
SELECT CASE hour!
CASE IS < 0
  hour! = hour! + 24
  JulianDay% = JulianDay% - 1
CASE IS >= 24
  hour! = hour! - 24
  JulianDay% = JulianDay% + 1
END SELECT
t1$ = fnt$(INT(hour!)) + ":" + fnt$(INT(60 * (hour! - INT(hour!)))) + MID$(t0$, 6)

SELECT CASE JulianDay%
CASE 0
  year% = year% - 1
  IF fnleap%(year%) THEN JulianDay% = 366 ELSE JulianDay% = 365
CASE IS > 365
  IF NOT fnleap%(year%) THEN JulianDay% = 1: year% = year% + 1
END SELECT

d1$ = juliantodate$(JulianDay%, year%)
END SUB

' Trims leading and trailing spaces from a string.
FUNCTION trim$ (i$)
          
trim$ = RTRIM$(LTRIM$(i$))

END FUNCTION