'    EMEPOL.BAS by G3SEK

'    *** MOON TRACKING PROGRAM ***  from "Amateur Radio Software"
'    by John Morris, GM4ANB (c) RSGB (available from ARRL).
'    For full explanation of astronomical aspects, read the book!
'    Microsoft BASIC version with mods by G3SEK and G4PMK,
'    September 1983 onwards

'    Optionally uses standard HOME.DAT file: an ASCII file containing -
'      latitude (decimal degrees N), longitude (decimal degrees WEST),
'       location (string, in quotes), callsign (string, in quotes)

'    EMEPOL calculates polarization from home station to selected
'    DX locations.


    DECLARE FUNCTION HHMM$ (hour%, MINUTES%)
    DECLARE FUNCTION MONTH$ (M%)
    DECLARE FUNCTION ATAN2 (SI, CO)
    DECLARE FUNCTION FRAC (X)

    COMMON SHARED PI, PT, TP

    ' Standard printer controls
    ESC$ = CHR$(27)
    FF$ = CHR$(12)
    ' Insert any special printer initiation codes here.
    LPINIT$ = ""
    ' Printer width: 100 is useful, but 80 is OK
    LWIDTH% = 100
   
    PI = 4 * ATN(1)
    TP = PI * 2
    PT = PI / 2
    RD = 180 / PI
    DR = PI / 180
    C1 = .9174640600000014#
    S1 = .397818675#
   
    DEF FNB (X) = TP * (X / TP - INT(X / TP))
 
    LOCATION$ = ""
    CALL$ = ""
 
    DIM SKY(24, 6)
    GOSUB FillSKY

    DX% = 21  ' 1 + number of DX locations for polarization-shift
    DIM DX$(DX%), DXN(DX%), DXE(DX%)
    GOSUB FillDX
  
    KEY OFF
    CLS
    LOCATE 2, 27
    COLOR 14
    PRINT " Moon Tracker by G3SEK "
    COLOR 7, 0

    GOSUB Where
   
    OPEN "SCRN:" FOR OUTPUT AS #3

WhatDay:
   
    YR% = VAL(RIGHT$(DATE$, 4))
    DY% = VAL(MID$(DATE$, 4, 2))
    MN% = VAL(LEFT$(DATE$, 2))
    LOCATE 6, 1
    PRINT SPACE$(75)
    LOCATE 6, 1
    PRINT "Day of month (1-31, Enter = today) : ";
    INPUT "", T$
    IF T$ = "" THEN
       LOCATE 6, 38
       PRINT SPACE$(20)
       LOCATE 6, 38
       PRINT USING "## & ####"; DY%; MONTH$(MN%); YR%
       GOTO SaveOld
    END IF
    DY% = INT(VAL(T$))
    IF (DY% <= 0 OR DY% > 31) THEN BEEP: GOTO WhatDay
    LOCATE 6, 38
    PRINT SPACE$(20)
    LOCATE 6, 38
    PRINT USING "##"; DY%

WhatMonth:
    LOCATE 8, 1
    PRINT SPACE$(75)
    LOCATE 8, 1
    PRINT USING "        Month number (Enter = ##)  : "; MN%;
    INPUT "", T$
    IF T$ <> "" THEN
       MN1% = INT(VAL(T$))
       IF (MN% <= 0 OR MN% > 12) THEN BEEP: GOTO WhatMonth
       MN% = MN1%
    END IF
    LOCATE 8, 38
    PRINT SPACE$(20)
    LOCATE 8, 38
    PRINT USING "##"; MN%
    OMN% = MN%
   
Whatyear:
    LOCATE 10, 1
    PRINT SPACE$(75)
    LOCATE 10, 1
    PRINT USING "Year               ( Enter = ####) : "; YR%;
    INPUT "", T$
    IF T$ <> "" THEN
       YR1 = VAL(T$)
       IF YR1 <= 0 THEN BEEP: GOTO Whatyear
       YR% = INT(YR1)
       IF YR% < 100 THEN YR% = YR% + 1900
    END IF
    LOCATE 10, 38
    PRINT SPACE$(20)
    LOCATE 10, 38
    PRINT USING "####"; YR%

SaveOld:      ' Save starting values for repeat
    ODY% = DY%
    OMN% = MN%
    OYR% = YR%
    GOSUB InitGST

WhatStartTime:
    LOCATE 12, 1
    PRINT SPACE$(75)
    LOCATE 12, 1
    T2$ = LEFT$(TIME$, 2)
    PRINT "  GMT start (hh00,  Enter = "; T2$; "00)  :   00";
    LOCATE 12, 38
    INPUT "", T$
    T$ = LEFT$(T$, 2)

    IF (T$ = "0") OR (T$ = "00") THEN
       TSTART% = 0
    ELSEIF T$ = "" THEN
       TSTART% = VAL(T2$)
    ELSE
       TSTART% = INT(VAL(T$))
       IF (TSTART% <= 0) OR (TSTART% > 23) THEN
          BEEP
          GOTO WhatStartTime
       END IF
    END IF

    LOCATE 12, 38
    PRINT SPACE$(20)
    LOCATE 12, 38
    PRINT HHMM$(TSTART%, 0)

    TSTART% = TSTART% * 100
    OTSTART% = TSTART%

WhatInterval:
    LOCATE 14, 1
    PRINT SPACE$(75)
    LOCATE 14, 1
    PRINT "  Interval  (minutes, Enter = 60)  : ";
    INPUT "", T$
    IF T$ <> "" THEN
       DT = INT(VAL(T$))
       IF (DT < 5) OR (DT > 60) THEN BEEP: GOTO WhatInterval
       DT = 5 * INT(DT / 5)
    ELSE
       DT = 60
    END IF
    LOCATE 14, 38
    PRINT SPACE$(20)
    LOCATE 14, 40
    PRINT USING "##"; DT
    ODT% = DT

WhatMoonsets:
    LOCATE 16, 1
    PRINT SPACE$(75)
    LOCATE 16, 1
    PRINT "  Number of moonsets  (Enter = 1)  : ";
    INPUT "", T$
    IF T$ <> "" THEN
       SS% = VAL(T$)
       IF (SS% < 1) OR (SS% > 5) THEN BEEP: GOTO WhatMoonsets
    ELSE
       SS% = 1
    END IF
    LOCATE 16, 38
    PRINT SPACE$(20)
    LOCATE 16, 40
    PRINT SS%
    OSS% = SS%

    SLEEP (1)
    CLS


'*************** MAIN CALCULATION LOOP **********************************


StartCalculation:
   
    IF PRINTING% THEN
       COLOR 0, 3
       LOCATE 24, 1
       PRINT " PRINTING... ";
       COLOR 7, 0
    END IF
        
    MOONSETS% = 0
    MOONUP% = 0
    JUSTSET% = 0
    T1 = TSTART% / 100
    GM = (60 * INT(T1) + FRAC(T1) * 100) / 60
    DT = DT / 60
   
OncePerDay:
    GOSUB DayOfWeek
    GOSUB DaynumToDate
    MN$ = MONTH$(MN%)

    GOSUB NewScreen
    PRINT #3, USING "## & ####"; DY%; MONTH$(MN%); YR%
    PRINT #3,
    PRINT #3, "  GMT"; TAB(15); "Azimuth"; TAB(25); "Elevation"
    PRINT #3,

OncePerInterval:
    TM = GM
    GOSUB HoursMinutes
    GOSUB GMTtoGST
    GOSUB MoonPosition
    GOSUB EclipticToEquatorial
    GOSUB EqToEL
    GOSUB MoonParallax

    GOSUB CheckIfSet

    IF MOONUP% THEN  ' finish off this calculation
       GOSUB EqToAZ
       EL = CE
      
       GOSUB NewScreen
       PRINT #3, HHMM$(HR%, MT%); TAB(16);
       PRINT #3, USING "###.#"; RD * AZ;
       PRINT #3, TAB(28);
       PRINT #3, USING "##.#"; RD * EL
   
      ' do this just before every full hour,
      ' and also just before every half-hour if DT < 0.5h
      IF (MT% >= (30 - 60 * DT) AND MT% < 30) OR MT% >= 60 * (1 - DT) THEN
          GOSUB MoonDoppler
          'NB Pol_shift alters variables relevant to home-stn Doppler
          GOSUB PolarizationShift
          PRINT #3, USING "Echos +#### Hz"; DV
       END IF

    ELSE  ' moon not up
   
       IF JUSTSET% THEN GOSUB Lastline

       IF MOONSETS% >= SS% THEN           ' tidy-up and exit via Menu
          CLOSE #3
          GOSUB Menu
       END IF

    END IF

    ' Next time step
    GM = GM + DT
    IF GM < 24 THEN
       GOTO OncePerInterval
    ELSE
       DO
          GM = GM - 24
          DN% = DN% + 1
       LOOP UNTIL GM < 24   ' In case of time steps > 1 day
       GOTO OncePerDay
    END IF

Cutshort:
    GOSUB Menu

' ************************ SUBROUTINES **********************************
   

InitGST:
     T1 = YR% - 1
     DE = INT(365.25 * (T1 - 1980)) - INT(T1 / 100) + INT(T1 / 400) + 381
     T1 = (DE + 29218.5) / 36525!
     T1 = 6.6460656# + T1 * (2400.051262# + T1 * 2.581E-05)
     SE = T1 - 24 * (YR% - 1900)
'    Date to Day Number
     T1 = YR%
     T2 = MN%
     IF T2 <= 2.5 THEN
        T1 = T1 - 1
        T2 = T2 + 12
     END IF
     DN% = INT(365.25 * (T1 - 1980)) - INT(T1 / 100) + INT(T1 / 400) - 16
     DN% = DN% + DY% + 30 * T2 + INT(.6 * T2 - .3)
     RETURN

                      
DayOfWeek:
     T1 = INT(DN% - 7 * INT(DN% / 7) + .5)
     D$ = MID$("   Mon  TuesWednes Thurs   Fri Satur   Sun", 6 * T1 + 1, 6)
     D$ = D$ + "day"
     RETURN
                               

DaynumToDate:
     T2 = INT(DN% - 39410!)
     T1 = INT((T2 + 32044.8) / 36524.3)
     T1 = T1 + T2 - INT(T1 / 4) + 1486
     YR% = INT((T1 - 122.1) / 365.25)
     T1 = T1 - INT(365.25 * YR%)
     MN% = INT(T1 / 30.6001)
     DY% = T1 - INT(30.6001 * MN%)
     YR% = YR% + 2084
     MN% = MN% - 1
     IF MN% > 12 THEN
        YR% = YR% + 1
        MN% = MN% - 12
     END IF
     RETURN


HoursMinutes:
     MT% = INT(TM * 60 + .5)
     HR% = INT(MT% / 60)
     MT% = MT% - HR% * 60
     RETURN


GMTtoGST:
     T1 = (SE + .0657098 * (DN% - DE) + GM * 1.00274) / 24
     GS = 24 * FRAC(T1)
     RETURN


MoonPosition:
     D = DN% + GM / 24
     EW = FNB(1.134193# + D * .2299715060000012#)
     MM = FNB(1.319238# + D * .228027135#)
     T1 = FNB(6.217512000000011# + D * .01720196977#)
     T2 = 2 * FNB(2.550677# + D * .212768711#)
     T3 = FNB(4.7652214# + D * .230895723#)
     EW = EW + .01148 * SIN(T2) + .10976 * SIN(MM) - .022235 * SIN(MM - T2)
     EW = EW - .003246 * SIN(T1) + .003735 * SIN(2 * MM) - .0019897 * SIN(2 * T3)
     EW = EW - .0010297 * SIN(2 * MM - T2) - .0009948 * SIN(MM + T1 - T2)
     EN = T3 + .011507 * SIN(T2) + .10873924# * SIN(MM) - .0222006 * SIN(MM - T2)
     EN = .0897797 * SIN(EN) - .002548 * SIN(T3 - T2)
     RETURN


EclipticToEquatorial:
     SI = C1 * SIN(EN) + S1 * COS(EN) * SIN(EW)
     CO = SQR(1 - SI * SI)
     DC = ATAN2(SI, CO)
     RA = ATAN2((SIN(EW) * C1 - TAN(EN) * S1), COS(EW))
     IF RA < 0 THEN RA = RA + TP
     RETURN


EqToEL:
     T1 = GS / 24 - RA / TP
     GH = TP * FRAC(T1)
     LH = GH + E
     SI = COS(LH) * COS(DC) * COS(N) + SIN(DC) * SIN(N)
     CO = SQR(1 - SI * SI)
     EL = ATAN2(SI, CO)
     RETURN


MoonParallax:
     RO = .996986 / (1 + .0549 * COS(MM + .10976 * SIN(MM)))
     CE = EL - RO * .0166 * COS(EL)
     RETURN


CheckIfSet:
     IF CE <= 0 THEN
        IF MOONUP% THEN             ' moon has just set
           JUSTSET% = -1
           MOONSETS% = MOONSETS% + 1
           MOONUP% = 0
        ELSE
           JUSTSET% = 0
           RETURN
        END IF
     ELSE                           ' moon has just risen, or is already up
        MOONUP% = -1
        RETURN
     END IF
     RETURN


EqToAZ:
     AZ = ATAN2((-SIN(LH) * COS(DC) * COS(N)), (SIN(DC) - SIN(N) * SIN(EL)))
     IF AZ < 0 THEN AZ = AZ + TP
     RETURN


MoonDoppler:
     T2 = .10976
     T1 = MM + T2 * SIN(MM)
     DMDT = .01255 * RO * RO * SIN(T1) * (1 + T2 * COS(MM))
     DMDT = DMDT * 4449             ' moon's radial velocity, m/s
     T1 = 6378                      ' earth's mean radius
     T2 = 384401                    ' mean orbital radius
     FREQ = 432

' Old Doppler calculation
     T3 = T1 * T2 * (COS(DC) * COS(N) * SIN(LH)) / SQR(T2 * T2 - 2 * T2 * T1 * SIN(EL))
'    NB  Minor bug in original code - there ought to be a 2 here ^
'        but it makes only about 2Hz difference to results!
     DV = DMDT + T3 * .0753125
     DV = -DV * 2 * FREQ / 299.8:   ' Doppler (Hz); *2 is for echo
     RETURN
' New Doppler calculation   -  see notes for derivation
     T1 = 1000 * T1                           ' Earth's mean radius, m
     T2 = T2 * RO * 1000                      ' Corrected orbital radius, m
     R = SQR(T2 * T2 + 2 * T1 * T2 * SIN(EL)) ' Slant range, observer to
                                              ' moon centre, m

     ' Correct the moon's centre-centre radial velocity for slant range
     DRDT1 = DMDT * (T2 + (T1 * SIN(EL))) / R ' Slant radial velocity, m/s
    
     ' Rate of change of hour angle, radian/sec - from book
     DLHDT = 7.53125E-05 * (SIN(LH) * COS(DC) * COS(N))

     ' Calculate rate of change of declination, assuming sinusoidal monthly
     ' variation up to current max value (obtain by running program).
     DCMAX = 20 * DR                          ' Present max. declination, rad
     DDCDTMAX = .0000716 * DR                 ' Present max. dDC/dt, rad/sec
     T3 = DC / DCMAX
     IF (T3 < 1) THEN
        T3 = SIN(T3)
        DDCDT = DDCDTMAX * SQR(1 - T3 * T3)
        T3 = DDCDT * (COS(DC) * SIN(N) - COS(LH) * SIN(DC) * COS(N))
     ELSE
        T3 = 0   ' Zero rate of change of DC, if DC is near DCMAX
     END IF

     ' Add terms due to change of hour angle and change of declination
     DSINELDT = DLHDT + T3
     ' ... to give observer's relative motion term DRDT2, m/s
    
     DRDT2 = T1 * T2 * DSINELDT / R
    
     DV = DRDT1 + DRDT2                 ' Total relative velocity, m/s

     DV = -DV * 2 * FREQ / 299.8        ' Echo Doppler

'     IF DV < 0 THEN STOP                ' For testing only...

     RETURN
    

PolarizationShift:
      COLOR 9
      PRINT #3, "Pol rotator: ";
      GOSUB Polarization
      HOMEPOL = PO
      EHOME = E
      NHOME = N
      J% = 0
      FOR I% = 1 TO DX%
         IF DX$(I%) = "" THEN EXIT FOR
         E = DXE(I%)
         N = DXN(I%)
         GOSUB EqToEL
         GOSUB MoonParallax
         IF CE > 0 THEN
            GOSUB EqToAZ
            GOSUB Polarization
            OFFSET = PO - HOMEPOL
            OFFSET = (RD * OFFSET) MOD 180
            IF ABS(OFFSET) > 90 THEN OFFSET = OFFSET - (180 * SGN(OFFSET))
            PRINT #3, USING "&+###ø   "; DX$(I%); OFFSET;
            J% = J% + 1
            IF (J% = 6) THEN ' start next line
                  PRINT #3,
                  PRINT #3, SPACE$(13);
                  J% = 0
            END IF
         END IF
      NEXT I%
      E = EHOME
      N = NHOME
      PRINT #3,
      COLOR 7
      RETURN
  
Polarization:
      PO = ATAN2((SIN(LH) * COS(DC) * COS(N)), (SIN(N) - SIN(DC) * SIN(EL)))
      RETURN

FillDX:
      RESTORE DXlocations

      FOR I% = 1 TO DX%
         READ DX$(I%)
         IF DX$(I%) = "" THEN EXIT FOR
         READ DXN(I%), DXE(I%)
         DXN(I%) = DXN(I%) * DR
         DXE(I%) = DXE(I%) * DR
      NEXT I%

DXlocations:
' Degrees      N      E      in order E to W.  Max number = DX%
DATA "VK3",  -37.82, 144.97
DATA "JA ",   35.67, 139.75
DATA "9M2",    3.0,  101.7
DATA "4X ",   31.75,  35.25
DATA "UB5",   50.43,  30.52
' DATA "A22",  -24.75,  25.9
DATA "ZS6",  -34.0,   18.5
' DATA "EA9",   35.88,  -5.32
' DATA "EA8",   28,    -15.5
DATA "PY ",   -16,   -48
DATA "W1 ",   41.75, -71.70
DATA "HP3",    8.95, -79.5
DATA "W4 ",   30.43, -84.32
DATA "W5 ",   36.17, -85.83
DATA "W0 ",   41.58, -93.58
DATA "W7 ",  39.17, -119.77
DATA "W6 ",   38.53, -121.50
DATA "KH6",   21.32, -157.83
'Last entry is a trap
DATA "", 0, 0
     
      RETURN


Lastline:  ' Print data after moonset
    IF PRINTING% THEN PRINT #3,
    COLOR 2
    PRINT #3, USING "RA  +###.#ø  Dec +##.#ø      "; RD * RA; RD * DC;
    PRINT #3, USING "Signals +##.# dB      "; -17.37 * LOG(RO);
    GOSUB SkyTemperature
    PRINT #3, USING "Sky temp #### K"; SKYTEMP
    COLOR 7
    IF PRINTING% THEN PRINT #3, FF$;
    RETURN


SkyTemperature:
      GHA% = INT(1 + RA * RD / 15)
      DEC% = INT(DC * RD / 10) + 4
      '   400MHz temperatures corrected to 432MHz by (400/432)^2.4 = 1/1.2
      SKYTEMP = SKY(GHA%, DEC%) / 1.2
      RETURN


FillSKY:
      RESTORE Skydata
      FOR I% = 1 TO 24
         FOR J% = 6 TO 1 STEP -1
            READ SKY(I%, J%)
         NEXT J%
      NEXT I%
      '   400MHz sky temperatures, from R E Taylor, Proc IEEE Apr 1973 p469
      '   I=GHA; J=Dec 30to20,20to10,...,-20to-30; ie grid is 15x10 degrees.
Skydata:
      DATA 23,23,23,25,25,22, 25,25,25,25,25,20, 30,30,27,24,20,19
      DATA 30,30,24,22,18,18
      DATA 35,33,25,23,19,16, 43,41,34,30,24,17, 35,38,35,35,32,25
      DATA 25,25,28,30,34,32
      DATA 17,17,17,19,22,28, 15,16,15,18,20,23, 17,17,17,17,19,22
      DATA 18,18,19,20,23,23
      DATA 20,24,26,25,25,25, 24,30,25,25,26,30, 33,35,28,28,32,34
      DATA 38,42,35,32,36,43
      DATA 42,55,44,40,48,80, 38,50,70,90,110,200, 60,80,120,160,150,180
      DATA 70,75,100,90,70,50
      DATA 55,40,40,34,29,29, 45,27,22,28,25,22, 27,27,20,20,21,20
      DATA 23,23,21,23,25,22
      RETURN


Where:   ' Observer's location input
      LOCATE 4, 10
      PRINT "Press <Enter> to use HOME.DAT file.    Press N for new QTH."
      GOSUB Keypress
      IF T$ <> "N" THEN
         ON ERROR GOTO NoHome
         OPEN "HOME.DAT" FOR INPUT AS #1
         INPUT #1, LA, LO, LOCATION$, CALL$
         CLOSE
         ON ERROR GOTO 0
         N = LA * DR
         E = -LO * DR
      ELSE
         PRINT
         INPUT "        Latitude (whole degrees)   : ", T1
         INPUT "                 (decimal minutes) : ", T2
         PRINT "      North or South?  (Press N/S) > ";
         DO
            GOSUB Keypress
         LOOP UNTIL T$ = "N" OR T$ = "S"
         PRINT T$
         PRINT
         N = DR * (T1 + T2 / 60)
         IF T$ = "S" THEN N = -N
         INPUT "       Longitude (whole degrees)   : ", T1
         INPUT "                 (decimal minutes) : ", T2
         PRINT "        East or West?  (Press E/W) > ";
         DO
            GOSUB Keypress
         LOOP UNTIL T$ = "E" OR T$ = "W"
         PRINT T$
         E = DR * (T1 + T2 / 60)
         IF T$ = "W" THEN E = -E
      END IF
      RETURN


NoHome:
      BEEP
      PRINT
      PRINT "Can't find HOME.DAT on this directory"
      END

      
Menu:    ' for repeat calculations - returns to various labels in MAIN
    LOCATE 24, 1
    PRINT SPACE$(30);
    LOCATE 24, 1
    PRINT "   Repeat ("; : COLOR 11: PRINT "Y"; : COLOR 7: PRINT "/n)? : ";
    GOSUB Keypress
    IF T$ = "N" THEN
       CLS
       END
    ELSE
       LOCATE 24, 20
       PRINT "Y";
    END IF
 
    PRINTING% = 0
    LOCATE 25, 1
    PRINT "Printout (y/"; : COLOR 11: PRINT "N"; : COLOR 7: PRINT ")?  : ";
    GOSUB Keypress
    IF T$ = "Y" THEN
       PRINTING% = -1
       OPEN "LPT1:" FOR OUTPUT AS #3
       WIDTH #3, LWIDTH%
       PRINT #3, LPINIT$;
    ELSE
       OPEN "SCRN:" FOR OUTPUT AS #3
    END IF

    CLS
    PRINT
    PRINT TAB(10); "REPEAT CALCULATION"
    PRINT
    PRINT "Old parameters were:"
    PRINT "Start date         ";
    PRINT USING "## & ####"; ODY%; MONTH$(OMN%); OYR%
    PRINT "Start time         "; HHMM$(OTSTART% / 100, 0)
    PRINT "Interval (mins)   "; ODT%
    PRINT "Moonsets          "; OSS%
    PRINT
    PRINT "Press   "; : COLOR 11: PRINT "R"; : COLOR 7
      PRINT "         "; : COLOR 11: PRINT "R"; : COLOR 7
      PRINT "epeat old calculation": PRINT
    PRINT "        "; : COLOR 11: PRINT "D"; : COLOR 7
      PRINT "         new start "; : COLOR 11: PRINT "D"; : COLOR 7
      PRINT "ate, including-": PRINT
    PRINT "        "; : COLOR 11: PRINT "T"; : COLOR 7
      PRINT "         new start "; : COLOR 11: PRINT "T"; : COLOR 7
      PRINT "ime, including-": PRINT
    PRINT "        "; : COLOR 11: PRINT "I"; : COLOR 7
      PRINT "         new "; : COLOR 11: PRINT "I"; : COLOR 7
      PRINT "nterval,   including-": PRINT
    PRINT "        "; : COLOR 11: PRINT "N"; : COLOR 7
      PRINT "         new "; : COLOR 11: PRINT "N"; : COLOR 7
      PRINT "umber of moonsets": PRINT
    PRINT "        "; : COLOR 11: PRINT "Q"; : COLOR 7
      PRINT "         "; : COLOR 11: PRINT "Q"; : COLOR 7
      PRINT "uit"
 
    DO
       GOSUB Keypress
       SELECT CASE T$

          CASE "R"   ' Repeat old calculation
             CLS
             DY% = ODY%
             MN% = OMN%
             YR% = OYR%
             GOSUB InitGST
             TSTART% = OTSTART%
             DT = ODT%
             SS% = OSS%
             RETURN StartCalculation

          CASE "D"   ' New start date etc
             CLS
             RETURN WhatDay

          CASE "T"   ' New start time etc
             CLS
             DY% = ODY%
             MN% = OMN%
             YR% = OYR%
             GOSUB InitGST
             RETURN WhatStartTime

          CASE "I"   ' New time interval etc
             CLS
             DY% = ODY%
             MN% = OMN%
             YR% = OYR%
             GOSUB InitGST
             TSTART% = OTSTART%
             RETURN WhatInterval

          CASE "N"   ' New number of moonsets
             CLS
             DY% = ODY%
             MN% = OMN%
             YR% = OYR%
             GOSUB InitGST
             TSTART% = OTSTART%
             DT = ODT%
             RETURN WhatMoonsets

          CASE "Q"   ' Quit
             CLS
             END

          CASE ELSE  ' Not a menu key
             BEEP

       END SELECT

    LOOP ' until valid keypress


Keypress:
      DO
        T$ = UCASE$(INPUT$(1))
      LOOP UNTIL T$ <> ""
      RETURN


NewScreen:
      IF PRINTING% OR (CSRLIN < 21) THEN RETURN
      LOCATE 25, 1
      COLOR 0, 3
      PRINT " Press Enter or Esc ";
      COLOR 7, 0
      T$ = INPUT$(1)
      IF T$ = CHR$(27) THEN
         CLOSE #3
         LOCATE 25, 1
         PRINT SPACE$(25);
         RETURN Cutshort
      END IF
      CLS
      RETURN


' ************************ FUNCTIONS ***********************************
' Filed separately - press F2

FUNCTION ATAN2 (SI, CO)
    
     T1 = ABS(SI)
     T2 = ABS(CO)
     IF T1 > T2 THEN
        TH = PT - ATN(T2 / T1)
     ELSE
        TH = ATN(T1 / T2)
     END IF
     IF CO < 0 THEN TH = PI - TH
     IF SI < 0 THEN TH = -TH
     ATAN2 = TH

END FUNCTION

FUNCTION FRAC (X)
' Fractional part of X
   FRAC = X - INT(X)
END FUNCTION

FUNCTION HHMM$ (HOURS%, MINUTES%)
  
   ' Returns 24-hour clock time with leading zeros

   T% = 10000 + (100 * HOURS%) + MINUTES%
   HHMM$ = RIGHT$(STR$(T%), 4)

END FUNCTION

    FUNCTION MONTH$ (M%)
       MONTH$ = MID$("JanFebMarAprMayJunJulAugSepOctNovDec", 3 * M% - 2, 3)
    END FUNCTION