'===========================================================================
' Subject: CELESTIAL COORD TO ALT-AZIMUTH     Date: 06-29-99 (19:37)       
'  Author: Jim Phillips                       Code: QB, QBasic, PDS        
'  Origin: jphillip@northernnet.com         Packet: ALGOR.ABC
'===========================================================================
' CEL2ALAZ - Converts celestial coordinates to alt-azimuth.
'            From Sky & Telescope magazine, May and June 1984.
'            Much modified and enhanced by Jim Phillips.
  
COLOR 4            'Display in red, for night vision.
Pi = 3.14159265#
DegToRad = Pi / 180

PRINT
PRINT "     Note: To accept any default value, press Enter."
PRINT

StartProgram:
GOSUB GetLatLongTZone
GOSUB GetDate

DO
  GOSUB GregorianToJulian
  GOSUB GetLocalTime
  GOSUB AdjustForLocalTime
  GOSUB JulianToGMST
    TimeRads = S * 15 * DegToRad
  GOSUB GetRAAndDec
  GOSUB CalcAltAndAz
GetChoice:
  PRINT "        Find another   Restart program   Quit program"
  INPUT "Choose:  F              R                 Q          "; ChoiceStr$
  SELECT CASE UCASE$(ChoiceStr$)
    CASE "F"
      'Continue with DO loop....
    CASE "R"
      GOTO StartProgram
    CASE "Q"
      END
    CASE ELSE          'Input was not valid....
      GOTO GetChoice
  END SELECT
LOOP

AdjustForLocalTime:
  TimeAdjust = (LocalTime + TZoneAdjust) / 24     'Fraction of a day.
    JulFrac = JulFrac + TimeAdjust
    IF JulFrac >= 1 THEN
      JulFrac = JulFrac - 1
      JulDay = JulDay + 1
    END IF
RETURN

CalcAltAndAz:
  ' From Sky & Telescope, June 1984, p. 558.  Converts Right Ascension/
  ' Declination to Altitude/Azimuth.
  T5 = TimeRads - RARads + LongRads
  S1 = SIN(LatRads) * SIN(DecRads)
  S1 = S1 + COS(LatRads) * COS(DecRads) * COS(T5)
  C1 = 1 - S1 * S1
  IF C1 > 0 THEN C1 = SQR(C1)
  IF C1 > 0 THEN
    AltRads = ATN(S1 / C1)
  ELSE
    AltRads = SGN(S1) * Pi / 2
  END IF
  C2 = COS(LatRads) * SIN(DecRads)
  C2 = C2 - SIN(LatRads) * COS(DecRads) * COS(T5)
  S2 = -COS(DecRads) * SIN(T5)
  IF C2 = 0 THEN
    AzRads = SGN(S2) * Pi / 2
  ELSE
    AzRads = ATN(S2 / C2)
    IF C2 < 0 THEN AzRads = AzRads + Pi
  END IF
  IF AzRads < 0 THEN AzRads = AzRads + 2 * Pi
  PRINT "ALTITUDE: ";
  PRINT USING "###.#"; AltRads / DegToRad;
  PRINT "     AZIMUTH: ";
  PRINT USING "###.#"; AzRads / DegToRad
  PRINT
RETURN

GetRAAndDec:
  INPUT "Right Ascension HH MM SS "; RAStr$
    A$ = LEFT$(RAStr$, 2)          'Keep consistent with other entry format.
    A2 = VAL(MID$(RAStr$, 4, 2))
    A3 = VAL(RIGHT$(RAStr$, 2))
    GOSUB SexaToDeci
    RARads = DeciEquiv * 15 * DegToRad
  INPUT "Declination (-)DD MM SS "; DecStr$
    MinusFudge = 0
      IF LEFT$(DecStr$, 1) = "-" THEN MinusFudge = 1   'Negative declination.
    A$ = LEFT$(DecStr$, 2 + MinusFudge)
    A2 = VAL(MID$(DecStr$, 4 + MinusFudge, 2))
    A3 = VAL(RIGHT$(DecStr$, 2))
    GOSUB SexaToDeci
    DecRads = DeciEquiv * DegToRad
RETURN

GetDate:
  DateStr$ = DATE$            'Get date from computer to offer as a default.
    Year$ = RIGHT$(DateStr$, 4)   'Work with strings to make it simpler for
    Month$ = LEFT$(DateStr$, 2)   'user to accept default or change date.
    Day$ = MID$(DateStr$, 4, 2)
    PRINT "Date (-)YYYY MM DD ( default is "; Year$; " "; Month$; " "; Day$;
      INPUT " ) "; NewDateStr$
      MinusFudge = 0
        IF LEFT$(NewDateStr$, 1) = "-" THEN MinusFudge = 1   'Negative year.
      IF NewDateStr$ <> "" THEN
        Year$ = LEFT$(NewDateStr$, 4 + MinusFudge)
        Month$ = MID$(NewDateStr$, 6 + MinusFudge, 2)
        Day$ = RIGHT$(NewDateStr$, 2)
      END IF
  Year = VAL(Year$)
  Month = VAL(Month$)
  Day = VAL(Day$)
RETURN

GetLocalTime:
  TimeStr$ = TIME$           'Get time from computer to offer as a default.
    THrs$ = LEFT$(TimeStr$, 2)     'Work with strings to make it simpler for
    TMins$ = MID$(TimeStr$, 4, 2)  'user to accept default or change time.
    TSecs$ = RIGHT$(TimeStr$, 2)
    PRINT "Local 24-hour-clock time HH MM SS ( default is "; THrs$; " ";
      PRINT TMins$; " "; TSecs$; " ) ";
    INPUT NewTimeStr$
      IF NewTimeStr$ <> "" THEN
        THrs$ = LEFT$(NewTimeStr$, 2)
        TMins$ = MID$(NewTimeStr$, 4, 2)
        TSecs$ = RIGHT$(NewTimeStr$, 2)
      END IF
  A$ = THrs$      'Leave as a string, to be compatible with SexaToDeci Sub.
  A2 = VAL(TMins$)
  A3 = VAL(TSecs$)
  GOSUB SexaToDeci
  LocalTime = DeciEquiv   'Local time as a decimal.
RETURN

GregorianToJulian:
  ' From Sky & Telescope, May 1984, p. 455.  Converts a Gregorian calendar
  ' date to a Julian date.
  D1 = INT(Day)
  JulFrac = Day - D1 - .5
  JulDay = -INT(7 * (INT((Month + 9) / 12) + Year) / 4)
  S = SGN(Month - 9)
  A = ABS(Month - 9)
  J1 = INT(Year + S * INT(A / 7))
  J1 = -INT((INT(J1 / 100) + 1) * 3 / 4)
  JulDay = JulDay + INT(275 * Month / 9) + D1 + J1
  JulDay = JulDay + 1721027 + 2 + 367 * Year
  IF JulFrac < 0 THEN
    JulFrac = JulFrac + 1
    JulDay = JulDay - 1
  END IF
  PRINT "Julian Date: "; JulDay; JulFrac
RETURN

JulianToGMST:
  ' From Sky & Telescope, June 1984, p. 558.  Converts a Julian date to
  ' Greenwich Mean Sidereal Time.
  D = JulDay - 2451545
  T = D / 36525
  T1 = INT(T)
  J0 = T1 * 36525 + 2451545
  T2 = (JulDay - J0 + .5) / 36525
  S = 24110.54841# + 184.812866# * T1
  S = S + 8640184.812866# * T2
  S = S + .093104 * T * T
  S = S - .0000062# * T * T * T
  S = S / 86400
  S = S - INT(S)
  S = 24 * (S + (JulFrac - .5) * 1.002737909#)
  IF S < 0 THEN S = S + 24
  IF S > 24 THEN S = S - 24
  Hrs = INT(S)
  MinsSecs = 60 * (S - Hrs)
  Mins = INT(MinsSecs)
  Secs = 60 * (MinsSecs - Mins)
  PRINT "Greenwich Mean Sidereal Time: "; Hrs; Mins; Secs
RETURN

GetLatLongTZone:
  ' Latitude and Longitude default values below can be edited to whatever
  ' values you want (typically, your local coordinates).
  Latitude = 47.1                                     'Observer's latitude.
    PRINT "Latitude ( default is "; Latitude; ") ";
    INPUT LatitudeStr$
    IF LatitudeStr$ <> "" THEN       'User entered a different latitude.
      Latitude = VAL(LatitudeStr$)
    END IF
    LatRads = Latitude * DegToRad
  Longitude = -93.5                                   'Observer's longitude.
    PRINT "Longitude ( default is "; Longitude; ") ";
    INPUT LongitudeStr$
    IF LongitudeStr$ <> "" THEN        'User entered a different longitude.
      Longitude = VAL(LongitudeStr$)
    END IF
    LongRads = Longitude * DegToRad
  PRINT "Time zone adjustment (Eastern Standard Time 5, Central 6, etc."
    INPUT "  1 hour less if Daylight Savings.) "; TZoneAdjust
RETURN

SexaToDeci:
  ' From Sky & Telescope, June 1984, p. 558.  Converts a sexagesimal
  ' number (Hours/Degrees, Minutes, Seconds) to a decimal.  The
  ' Hours/Degrees value A$ is handled as a string in order to deal with
  ' Degree values -1 < value < 0.  (e.g. Basic will treat -00 35 20 as
  ' 00 35 20.  The method below works around this.)
  A1 = ABS(VAL(A$))
  DeciEquiv = A1 + A2 / 60 + A3 / 3600
  IF LEFT$(A$, 1) = "-" THEN DeciEquiv = -DeciEquiv
RETURN