DEFINT I-M
DEFDBL A-H, O-Z
OPTION BASE 1
DIM SHARED ROT(3, 3), XO(3), XN(3), PI2
DECLARE SUB R3 (TH)
DECLARE SUB R2 (TH)
DECLARE SUB R1 (TH)
DECLARE SUB M1 (ROT)
DECLARE SUB M2 (ROT)
DECLARE SUB M3 (ROT)
DECLARE SUB MATMULT (XO, XN, ROT)
DECLARE SUB XHAT (XLON, XLAT, XO)
DECLARE SUB CLR (XO, XN)
DECLARE FUNCTION PLG (X, Y)
PI = ATN(1#) * 4#
DEGRAD = PI / 180!
PI2 = PI / 2#
CLS
'
' Input angles:
'    PHI = site latitude in DEGREES
'    DEC = declination of star in DEGREES
'    RA = right ascension of star in DEGREES (e.g., 12h = 180 deg)
'    ST = local sidereal time (actual or mean depending on which
'         equinox has been chosen for reference) in DEGREES
'
PHI = 45#: PHIR = PHI * DEGRAD
DEC = 0#: DECR = DEC * DEGRAD
RA = 180#: RAR = RA * DEGRAD
ST = 0#: STR = ST * DEGRAD
'
' First rotate the RA and DEC (Q system) coordinates into HA and DEC
' (A system) by the formula X(HA,DEC)=M2*R3(ST)*X(RA,DEC)
'
CALL XHAT(RAR, PI2 - DECR, XO)
CALL R3(STR)
CALL MATMULT(XO, XN, ROT)
CALL CLR(XO, XN)
CALL M2(ROT)
CALL MATMULT(XO, XN, ROT)
CALL CLR(XO, XN)
'
' Vector XO now contains the components of the star position in the A system
'
' Now we convert from the A system to the H (horizon) system by the formula
'    X(AZ,ALT)=R3(180)*R2(90-PHI)*X(HA,DEC)
'
ANG = PI2 - PHIR
CALL R2(ANG)
CALL MATMULT(XO, XN, ROT)
CALL CLR(XO, XN)
CALL R3(PI)
CALL MATMULT(XO, XN, ROT)
CALL CLR(XO, XN)
'
' Vector XO contains the components of the star position in the H system
'
' Now we need to get the altitude and azimuth angles
'
X = XO(1)
Y = XO(2)
IF (X <> 0# OR Y <> 0#) THEN
  ARG = XO(3) / (SQR(XO(1) * XO(1) + XO(2) * XO(2)))
  ALTR = ATN(ARG)
  ALT = ALTR / DEGRAD
  PRINT "Altitude angle is "; : PRINT USING "###.###"; ALT; : PRINT "degrees."
ELSE
  IF (XO(3) > 0) THEN
    PRINT "Object is at the zenith."
  ELSE
    PRINT "Object is at the nadir."
  END IF
END IF
AZ = PLG(X, Y) / DEGRAD
IF (AZ <= 360#) THEN
  PRINT "Azimuth angle is "; : PRINT USING "###.###"; AZ; : PRINT "degrees."
ELSE
  PRINT "Azimuth is undefined."
END IF
SUB CLR (XO, XN)
'
' This subroutine replaces vector XO with vector XN
'
FOR I = 1 TO 3
    XO(I) = XN(I)
    TEST = ABS(XO(I))
    IF (TEST < .000000000000001#) THEN
       XO(I) = 0#
    END IF
NEXT I
END SUB
SUB M1 (ROT)
ROT(1, 1) = -1#
ROT(1, 2) = 0#
ROT(1, 3) = 0#
ROT(2, 1) = 0#
ROT(2, 2) = 1#
ROT(2, 3) = 0#
ROT(3, 1) = 0#
ROT(3, 2) = 0#
ROT(3, 3) = 1#
END SUB
SUB M2 (ROT)
ROT(1, 1) = 1#
ROT(1, 2) = 0#
ROT(1, 3) = 0#
ROT(2, 1) = 0#
ROT(2, 2) = -1#
ROT(2, 3) = 0#
ROT(3, 1) = 0#
ROT(3, 2) = 0#
ROT(3, 3) = 1#
END SUB
SUB M3 (ROT)
ROT(1, 1) = 1#
ROT(1, 2) = 0#
ROT(1, 3) = 0#
ROT(2, 1) = 0#
ROT(2, 2) = 1#
ROT(2, 3) = 0#
ROT(3, 1) = 0#
ROT(3, 2) = 0#
ROT(3, 3) = -1#
END SUB
SUB MATMULT (XO, XN, ROT)
'
' This subroutine computes the product of a vector and a rotation matrix
'
XN(1) = ROT(1, 1) * XO(1) + ROT(1, 2) * XO(2) + ROT(1, 3) * XO(3)
XN(2) = ROT(2, 1) * XO(1) + ROT(2, 2) * XO(2) + ROT(2, 3) * XO(3)
XN(3) = ROT(3, 1) * XO(1) + ROT(3, 2) * XO(2) + ROT(3, 3) * XO(3)
END SUB
FUNCTION PLG (X, Y)
  SGNAX = SGN(ABS(X))
  SGNAY = SGN(ABS(Y))
  SGNX = SGN(X)
  SGNY = SGN(Y)
  C = PI2 * (2# - (1# + SGNX) * (SGNY - SGNAY + 1#))
  IF (X <> 0#) THEN
    PLG = SGNAX * ATN(Y / X) + C
  ELSE
    IF (Y <> 0) THEN
       PLG = C
    ELSE
'      PLG not defined for X and Y both zero.
       PLG = 9999#
    END IF
  END IF
END FUNCTION
SUB R1 (TH)
'DIM ROT(3, 3)
CTH = COS(TH)
ST = SIN(TH)
ROT(1, 1) = 1!
ROT(1, 2) = 0!
ROT(1, 3) = 0#
ROT(2, 1) = 0#
ROT(2, 2) = CTH
ROT(2, 3) = STH
ROT(3, 1) = 0#
ROT(3, 2) = -STH
ROT(3, 3) = CTH
END SUB
SUB R2 (TH)
'DIM ROT(3, 3)
CTH = COS(TH)
STH = SIN(TH)
ROT(1, 1) = CTH
ROT(1, 2) = 0#
ROT(1, 3) = -STH
ROT(2, 1) = 0#
ROT(2, 2) = 1#
ROT(2, 3) = 0#
ROT(3, 1) = STH
ROT(3, 2) = 0#
ROT(3, 3) = CTH
END SUB
SUB R3 (TH)
CTH = COS(TH)
STH = SIN(TH)
ROT(1, 1) = CTH
ROT(1, 2) = STH
ROT(1, 3) = 0#
ROT(2, 1) = -STH
ROT(2, 2) = CTH
ROT(2, 3) = 0#
ROT(3, 1) = 0#
ROT(3, 2) = 0#
ROT(3, 3) = 1#
END SUB
SUB XHAT (XLON, XLAT, XO)
'
' This subroutine computes the unit vector given the longitude and
' co-latitude angles (in radians)
'
XO(1) = SIN(XLAT) * COS(XLON)
XO(2) = SIN(XLAT) * SIN(XLON)
XO(3) = COS(XLAT)
END SUB