Pole of satellites' tracks, for video rec. etc.

Bjorn Gimle (Bjorn_Gimle@lector.kth.se)
Wed, 5 Jul 95 20:10:12 +0100

I have written a program to demonstrate the possibilities of aligning
a telescope mount's polar axis, so a satellite can be tracked almost
exclusively along the RA axis.

This was suggested to simplify camcorder shots of Mir, but can also be
useful to record flashing periods/patterns of other (bright) objects.
If you can fix your binoculars or telescope to a similar mount, it can
offer a safer way of searching for faint and/or long-period flashers.

The trigonometric formula presented by David Moore, about June 25, is
correct, and probably the most compact one possible, but it is limited
to use two prediction points. 

I prefer using vector algebra methods, because most problems are solved
with elementary methods (the "inner" or "dot", and the "outer" or
"vector" products, which use only * and +, and square roots) and the
straightforward conversion from angles, and back after computation.

My program may seem long, but most of it is reading output from a 
prediction program (currently QuickSat 2.10, but the program identifies
the subroutines that depend on that format), the elementary operations
mentioned above, and printing the results.

If only one point is predicted on a pass, and it is identified as the
culmination, or two points are predicted, an approximate pole, 90 
degrees from the predictions, is computed (like David's result)

For three points, one (of the two opposite) point at equal distance
from the three points is computed.

For more than three points, a least squares solution is obtained.

The computed pole's RA and Decl. at the time of the first predicted
track point, is given. If you set your instrument NN minutes in
advance, subtract that many minutes from the RA (in HHMM format).

For each point predicted, the program prints the "latitude" or 
"declination" (90-distance) of the point from the computed pole.
This will show you how much declination adjustment you will need, and
even where the track will be located. With Mir and several 1400-
2100 km satellites, the declination was between -12 and +5, and the
variation +- 0.2 degrees or less. For an excentric orbit like Cosmos 
382, on a pass from 3200 to 5200 km, the declination was -17, and
the variation +- 0.5 degrees.

If you have MS QBasic ( MS-DOS 5.0 or later ), split the text below 
into files pole.ctl, pole.qou, pole.log, and pole.bas, 
and use this command: QBASIC /RUN pole.bas
Answer pole.ctl and xxxx.log to the prompts, and compare pole.log 
to xxxx.log.

If you need an .exe-cutable file (by QuickBasic 4.0) instead, I will 
send that to Bart for publication in the archive.

-------------------- Start of pole.ctl ----------------------------
 1995 05 Year, month number
 24 25 Start date, end date
  5.3  3.6 Start time, end time
   47.2238  -18.2284       44.    Someplace
  0 UT  24 correction and time zone name and 12/24 for UT to CDT
 2000 Epoch of predicted RA, Dec
  7.9 Magnitude limit
  8 Altitude cut-off value
 1.2 The search/step parameter value
 T f True means accept only the most recent elements for each object
 t True means ignore shadow test
 T  2 True means generate multiple prediction points, how many each way
 f True means output distance values in miles
 T True means generate a blank line before each object's prediction.
 M V                   B E b N p        N           v V p       f x E N      
v V p  Non-blank selects a class of objects
 A   D      Output format 
C:\TEMP\SATPROGN\qs57.mag
pole.qou
C:\TLE\cs950527.tle
EOF                                     End of input file list
-------------------- Start of pole.qou ----------------------------
  47.224 -18.228   44.    Someplace             2000  7.9  8 F T T F T
  
***  1995 May  25  Thu morning  *** Times are AM UT  ***  1934  152

 H  M  S  Tim Al Azi C Dir  Mag Dys F  Hgt Shd  Rng  EW Phs  R A  Dec

16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0  
 3 58 45   .0 44 148 C 271 -1.2   0 2  408 408  573 2.0 102 2252   5.0      

16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0  
 3 57 57   .0 35 187   299 -1.1   0 2  408 408  672 1.9  70 2056  -7.7      
 3 59 32   .0 36 107   241   .1   0 2  408 408  663 1.3 134  057  15.5      

16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0  
 3 57 57   .0 35 187   299 -1.1   0 2  408 408  672 1.9  70 2056  -7.7      
 3 58 45   .0 44 148 C 271 -1.2   0 2  408 408  573 2.0 102 2252   5.0      
 3 59 32   .0 36 107   241   .1   0 2  408 408  663 1.3 134  057  15.5      

16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0  
 3 57  9   .0 23 206   311  -.5   0 2  407 407  899 1.3  50 1940 -15.7      
 3 57 57   .0 35 187   299 -1.1   0 2  408 408  672 1.9  70 2056  -7.7      
 3 58 45   .0 44 148 C 271 -1.2   0 2  408 408  573 2.0 102 2252   5.0      
 4  0 20   .0 24  87   229  2.1   0 2  408 408  886  .7 155  222  19.0 

16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0  
 3 57  9   .0 23 206   311  -.5   0 2  407 407  899 1.3  50 1940 -15.7      
 3 57 57   .0 35 187   299 -1.1   0 2  408 408  672 1.9  70 2056  -7.7      
 3 58 45   .0 44 148 C 271 -1.2   0 2  408 408  573 2.0 102 2252   5.0      
 3 59 32   .0 36 107   241   .1   0 2  408 408  663 1.3 134  057  15.5      
 4  0 20   .0 24  87   229  2.1   0 2  408 408  886  .7 155  222  19.0 
-------------------- Start of pole.log ----------------------------
  47.224 -18.228   44.    Someplace             2000  7.9  8 F T T F T
***  1995 May  25  Thu morning  *** Times are AM UT  ***  1934  152
 H  M  S  Tim Al Azi C Dir  Mag Dys F  Hgt Shd  Rng  EW Phs  R A  Dec
16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0  
Consider the following QuickSat prediction:
  47.224 -18.228   44.    Someplace             2000  7.9  8 F T T F T
***  1995 May  25  Thu morning  *** Times are AM UT  ***  1934  152
 H  M  S  Tim Al Azi C Dir  Mag Dys F  Hgt Shd  Rng  EW Phs  R A  Dec
16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0  
16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0     1546  72.9 = Track
pole.
 -0.138294 -0.208933  0.813420 236.50  72.88
 3 58 45   .0 44 148 C 271 -1.2   0 2  408 408  573 2.0 102 2252   5.0  -0.0
---------------------------------------------------------------------------
16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0     1611  67.1 = Track
pole.
 -0.159588 -0.309323  0.825777 242.71  67.14
 3 57 57   .0 35 187   299 -1.1   0 2  408 408  672 1.9  70 2056  -7.7  -0.0
 3 59 32   .0 36 107   241   .1   0 2  408 408  663 1.3 134  057  15.5  -0.0
---------------------------------------------------------------------------
16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0     1525  65.8 = Track
pole.
 -0.281718 -0.351409  1.000000 231.28  65.75
 3 57 57   .0 35 187   299 -1.1   0 2  408 408  672 1.9  70 2056  -7.7  -4.0
 3 58 45   .0 44 148 C 271 -1.2   0 2  408 408  573 2.0 102 2252   5.0  -4.0
 3 59 32   .0 36 107   241   .1   0 2  408 408  663 1.3 134  057  15.5  -4.0
---------------------------------------------------------------------------
16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0     1519  65.6 = Track
pole.
 -0.293334 -0.345506  1.000000 229.67  65.62
 3 57  9   .0 23 206   311  -.5   0 2  407 407  899 1.3  50 1940 -15.7  -4.6
 3 57 57   .0 35 187   299 -1.1   0 2  408 408  672 1.9  70 2056  -7.7  -4.6
 3 58 45   .0 44 148 C 271 -1.2   0 2  408 408  573 2.0 102 2252   5.0  -4.6
 4  0 20   .0 24  87   229  2.1   0 2  408 408  886  .7 155  222  19.0  -4.6
---------------------------------------------------------------------------
16609 Mir Complex    86 17 A 32.7 4.2     294 1.2 V-1.0     1518  65.6 = Track
pole.
 -0.294336 -0.345857  1.000000 229.60  65.57
 3 57  9   .0 23 206   311  -.5   0 2  407 407  899 1.3  50 1940 -15.7  -4.6
 3 57 57   .0 35 187   299 -1.1   0 2  408 408  672 1.9  70 2056  -7.7  -4.6
 3 58 45   .0 44 148 C 271 -1.2   0 2  408 408  573 2.0 102 2252   5.0  -4.7
 3 59 32   .0 36 107   241   .1   0 2  408 408  663 1.3 134  057  15.5  -4.5
 4  0 20   .0 24  87   229  2.1   0 2  408 408  886  .7 155  222  19.0  -4.7
---------------------------------------------------------------------------
'------------------- Start of pole.bas ----------------------------
DEFSTR A-H, O-Q
TYPE VectorPos
        H AS INTEGER
        M AS INTEGER
        S AS INTEGER
        RA AS DOUBLE
        Dec AS DOUBLE
        Pred AS STRING * 70
        X AS DOUBLE
        Y AS DOUBLE
        Z AS DOUBLE
END TYPE
DECLARE FUNCTION acos! (c!)
DECLARE FUNCTION atan2# (Y AS DOUBLE, X AS DOUBLE)
DECLARE FUNCTION DotProduct! (p AS VectorPos, Q AS VectorPos)
DECLARE SUB FindZenith (a AS VectorPos, SiteLat!, R AS VectorPos)
DECLARE FUNCTION GetPredLine$ (QSline$)
DECLARE SUB GetQSdata (QSline$, SatPos() AS VectorPos, Lines)
DECLARE SUB NormVect (p AS VectorPos)
DECLARE SUB PrintPredDecl (Object, SatPos() AS VectorPos, Lines!)
DECLARE SUB RAdecData (p() AS VectorPos, i!)
DECLARE SUB SetupForQuickSat (Format, Lat, Lon)
DECLARE SUB SolveXYexact (a AS VectorPos, B AS VectorPos, c AS VectorPos, R AS
VectorPos)
DECLARE SUB SolveXYLSQ (p() AS VectorPos, Lines!)
DECLARE SUB VectorProduct (a AS VectorPos, B AS VectorPos, X AS VectorPos)
DECLARE SUB XYZdata (p() AS VectorPos, i!)
PRINT
PRINT "Reads prediction program output, and computes apparent pole of motion
for"
PRINT "each pass with at least two points, and 'declination' of each
pass/point. "
PRINT "Currently supports QuickSat 2.10 'A' format output (Az.in col.18-20,"
PRINT "RA in col.61-64). Easily adjusted to other prediction formats."
PRINT
' SetupForQuicksat, FindZenith, GetPredLine, and GetQSdata depend, more or
' less, on the prediction program format. For other programs, add your own
' subroutines, and calls to them, and comment out the old calls.
' Or, modify these routines to recognize and check for alternative formats
' by setting unique codes for PredFileFormat and PredType.
' The crucial points are that (for PredType="P") 'Lines' count the no. of
' prediction points for each pass; the .H .M .S .RA .Dec components of
' SatPos() are set before XYZdata is called; and that the processing of the
' points for one pass ( PredType="N" and Lines>0 ) is untouched.

CONST MaxLines = 20
DIM SatPos(MaxLines) AS VectorPos
SatPos(2).Pred = "-"    ' for debugging

SetupForQuickSat PredFileFormat, SiteLat, SiteLong

INPUT "Log calculations to file : ", PolFile
IF PolFile = "" THEN PolFile = "NUL"
OPEN "A", #2, PolFile

FOR ever! = 0 TO 1 STEP 0
    IF EOF(1) THEN
        PredType = "N"
    ELSE
        PredType = GetPredLine(QSline)
    END IF
    SELECT CASE PredType
    CASE "H"
        PRINT QSline; " Pole dec."
    CASE "P"
        'PRINT QSline
        Lines = 1 + Lines
        IF Lines > MaxLines THEN
            PRINT "Too many prediction lines for "; ObjectName
            PRINT "ignoring : "; LEFT$(QSline, 50); "......"
            Lines = Lines - 1
        ELSE
            GetQSdata QSline, SatPos(), Lines
            XYZdata SatPos(), Lines
        END IF
    CASE "N"
        SELECT CASE Lines
        CASE 1
            Culm = MID$(SatPos(1).Pred, 22, 1)
            IF Culm = "C" THEN
                PRINT "--- Only one prediction.  This is the culmination, an
approximate ---"
                PRINT "--- pole has been computed from the prediction and
SiteLat.       ---"
                FindZenith SatPos(1), SiteLat, SatPos(0)        ' (0)=zenith
                XYZdata SatPos(), 0
                VectorProduct SatPos(0), SatPos(1), SatPos(2)   ' (2)=90 deg
off
                RAdecData SatPos(), 2
                VectorProduct SatPos(2), SatPos(1), SatPos(0)   ' (0)="pole"
                IF SatPos(0).Z < 0 THEN
                    SatPos(0).X = -SatPos(0).X
                    SatPos(0).Y = -SatPos(0).Y
                    SatPos(0).Z = -SatPos(0).Z
                    SatPos(0).Dec = -SatPos(0).Dec
                    SatPos(0).RA = SatPos(0).RA - 180
                    IF SatPos(0).RA < 0 THEN SatPos(0).RA = 360 + SatPos(0).RA
                END IF
                PrintPredDecl ObjectName, SatPos(), 1
            ELSE
                PRINT ObjectName: PRINT #2, ObjectName
                PRINT "--- Only one prediction. If this is the culmination, an
approximate ---"
                PRINT "--- pole could be computed from the prediction and
SiteLat.         ---"
                PRINT SatPos(1).Pred, STRING$(75, "-")
                PRINT #2, SatPos(1).Pred, STRING$(75, "-")
            END IF
        CASE 2
            VectorProduct SatPos(1), SatPos(2), SatPos(0)
            PrintPredDecl ObjectName, SatPos(), 2
        CASE 3  ' This case can removed, if the next one is > 2 !
                ' SolveXYexact can then be erased.
            SolveXYexact SatPos(1), SatPos(2), SatPos(3), SatPos(0)
            PrintPredDecl ObjectName, SatPos(), 3
        CASE IS > 3
            SolveXYLSQ SatPos(), Lines
            PrintPredDecl ObjectName, SatPos(), Lines
        CASE ELSE                           ' Lines = 0
            PRINT QSline: PRINT #2, QSline
        END SELECT                          ' Lines
        ObjectName = LEFT$(QSline + SPACE$(80), 64)
        IF EOF(1) THEN EXIT FOR
        IF Lines > 0 AND Continuous <> "N" THEN
            PRINT
            PRINT "Press 'n' for non-stop printing, any other key for one pass
at a time : ";
            Continuous = UCASE$(INPUT$(1))
            PRINT
        END IF
        Lines = 0
    CASE ELSE                               ' PredType = ""
        '
    END SELECT                              ' PredType
NEXT ever!

CLOSE

'
-----------------------------------------------------------------------------
FUNCTION acos! (c!)
    S! = SQR(1 - c! * c!)
    IF c! = 0 THEN
        acos! = 3.14159265# / 2
    ELSE
        IF c! > 0 THEN
            acos! = ATN(S! / c!)
        ELSE
            acos! = 3.14159265# + ATN(S! / c!)
        END IF
    END IF
END FUNCTION

'
-----------------------------------------------------------------------------
FUNCTION atan2# (Y AS DOUBLE, X AS DOUBLE)
    v90# = 3.14159265357989# / 2#
    IF X = 0 THEN
        a# = v90#
        IF Y < 0 THEN a# = -a#
    ELSE
        a# = ATN(Y# / X#)
        IF X < 0 THEN a# = a# + 2# * v90#
    END IF
    IF a# < 0 THEN
        atan2# = a# + v90# * 4#
    ELSE
        atan2# = a#
    END IF
END FUNCTION

'
-----------------------------------------------------------------------------
FUNCTION DotProduct! (p AS VectorPos, Q AS VectorPos)
    DotProduct! = p.X * Q.X + p.Y * Q.Y + p.Z * Q.Z
END FUNCTION

SUB FindZenith (a AS VectorPos, SiteLat, R AS VectorPos)
' Find siderial time from Az and Elev in prediction line, and SiteLat
' If not in prediction line, ask for it.
    SHARED PredFileFormat
    Rad# = 3.1415926535# / 180
    Lat! = SiteLat * Rad#
    R.Dec = SiteLat
    IF PredFileFormat <> " A" THEN
        INPUT "Enter siderial time,converted to degrees : ", R.RA
    ELSE
        az! = VAL(MID$(a.Pred, 18, 3)) * Rad#
        el! = VAL(MID$(a.Pred, 14, 3)) * Rad#
        ha! = SIN(el!) * COS(Lat!) - COS(el!) * SIN(Lat!) * COS(az!)
        ha! = acos!(ha! / COS(a.Dec * Rad#)) / Rad#
        R.RA = a.RA - ha!
        IF az! > 3.1415926535# THEN R.RA = a.RA + ha!
    END IF
    R.H = a.H: R.M = a.M: R.S = a.S
END SUB

'.
-----------------------------------------------------------------------------
FUNCTION GetPredLine (QSline)
' Read one line from the prediction file #1 (QuickSat format assumed here).
' Return "P" for a prediction point line, "N" for a (possible) object name.
        LINE INPUT #1, QSline
        GetPredLine = "N"
        SELECT CASE LEN(RTRIM$(QSline))
        CASE 0:
            GetPredLine = " "
        CASE 70:
            GetPredLine = "N"
            IF MID$(QSline$, 69, 1) = "." THEN
                GetPredLine = "P"
            ELSEIF LEFT$(QSline, 9) = " H  M  S " THEN
                GetPredLine = "H"
            END IF
        CASE ELSE
            '
        END SELECT
END FUNCTION

'
-----------------------------------------------------------------------------
SUB GetQSdata (QSline, p() AS VectorPos, Lines)
' Get HMS, RA and Dec from Quicksat 2.10 'A' format line.

'3 57  9   .0 23 206   311  -.5   0 2  407 407  899 1.3  50 1940 -15.7
'          ^12   ^18                                        ^61  ^66
p(Lines).H = VAL(MID$(QSline, 1, 2))
p(Lines).M = VAL(MID$(QSline, 4, 2))
p(Lines).S = VAL(MID$(QSline, 7, 2))
p(Lines).Pred = QSline
' Add code for a/p in col.10 ?
p(Lines).RA = VAL(MID$(QSline, 61, 2)) * 15 + VAL(MID$(QSline, 63, 2)) / 4
p(Lines).Dec = VAL(MID$(QSline, 66, 5))

END SUB

'
-----------------------------------------------------------------------------
SUB NormVect (p AS VectorPos)
    R# = SQR(p.X * p.X + p.Y * p.Y + p.Z * p.Z)
    p.X = p.X / R#
    p.Y = p.Y / R#
    p.Z = p.Z / R#
END SUB

'
-----------------------------------------------------------------------------
SUB PrintPredDecl (ObjectName, SatPos() AS VectorPos, Lines)
        RAdecData SatPos(), 0
        Rad# = 180# / 3.14159265357989#
        RAh = INT(SatPos(0).RA / 15)
        RAm = 4 * (SatPos(0).RA - 15 * RAh)
        WHILE RAm > 59.5: RAm = RAm - 60: RAh = RAh + 1: WEND
        PRINT LEFT$(ObjectName, 60);
        PRINT USING "##"; RAh; RAm;
        PRINT USING " ###.# &"; SatPos(0).Dec; "= Pole."
        PRINT #2, LEFT$(ObjectName, 60);
        PRINT #2, USING "##"; RAh; RAm;
        PRINT #2, USING " ###.# &"; SatPos(0).Dec; "= Track pole."
        PRINT #2, USING "###.######"; SatPos(0).X; SatPos(0).Y; SatPos(0).Z;
        PRINT #2, USING "####.##"; SatPos(0).RA; SatPos(0).Dec
        FOR i = 1 TO Lines
            NormVect SatPos(0)
            Pdec! = 90 - Rad# * acos!(DotProduct(SatPos(0), SatPos(i)))
            'PRINT #2, USING "###.######"; SatPos(i).X; SatPos(i).Y;
SatPos(i).Z;
            'PRINT #2, USING "####.##"; SatPos(i).RA; SatPos(i).Dec
            PRINT #2, SatPos(i).Pred;
            PRINT #2, USING "####.#"; Pdec!
            PRINT SatPos(i).Pred;
            PRINT USING "####.#"; Pdec!
        NEXT i
        PRINT STRING$(75, "-")
        PRINT #2, STRING$(75, "-")
END SUB

'
-----------------------------------------------------------------------------
SUB RAdecData (SatPos() AS VectorPos, i)
    Rad# = 3.14159265357989# / 180#
    R# = SQR(SatPos(i).X * SatPos(i).X + SatPos(i).Y * SatPos(i).Y)
    IF R# > 0 THEN
        SatPos(i).Dec = ATN(SatPos(i).Z / R#) / Rad#
    ELSE
        SatPos(i).Dec = 90
    END IF
    SatPos(i).RA = atan2#(SatPos(i).Y, SatPos(i).X) / Rad#
END SUB

'
-----------------------------------------------------------------------------
SUB SetupForQuickSat (PredFileFormat, SiteLat, SiteLong)

' Open the prediction file as file #1. Determine PredFileFormat
' and SiteLat and SiteLon (normally not needed)

INPUT "Name of QuickSat control file (quicksat.ctl) : ", CtlFile
IF CtlFile = "" THEN CtlFile = "quicksat.ctl"   ' "H:\QUICKSKY\pole.qfg"  
OPEN "I", #1, CtlFile
FOR i = 1 TO 4: LINE INPUT #1, Cline: NEXT i
SiteLat = VAL(LEFT$(Cline, 10))
SiteLon = VAL(MID$(Cline, 11, 10))
FOR i = 1 TO 2: LINE INPUT #1, Cline: NEXT i
IF UCASE$(LEFT$(Cline, 2)) = " N" THEN
        PRINT "Older Quicksat formats not tested ?"
        LINE INPUT #1, Cline
END IF
IF UCASE$(LEFT$(Cline, 2)) = " Y" THEN
        PRINT "Quicksat Radio formats not supported"
        STOP
END IF
FOR i = 1 TO 10: LINE INPUT #1, Cline: NEXT i
PredFileFormat = UCASE$(LEFT$(Cline, 2))
IF PredFileFormat <> " A" THEN
        PRINT "Only Quicksat 'A' format supported (line 16 in "; CtlFile; ")"
        STOP
END IF
FOR i = 1 TO 2: LINE INPUT #1, Cline: NEXT i
CLOSE
QouFile = LEFT$(Cline, INSTR(Cline + " ", " ") - 1)
PRINT "Quicksat output file = "; QouFile
'3 57  9   .0 23 206   311  -.5   0 2  407 407  899 1.3  50 1940 -15.7
'          ^12   ^18                                        ^61  ^66
OPEN "I", #1, QouFile

END SUB '.
-----------------------------------------------------------------------------

DEFDBL X-Z
SUB SolveXYexact (a AS VectorPos, B AS VectorPos, c AS VectorPos, R AS
VectorPos)
    z2 = B.Z - a.Z
    z3 = c.Z - a.Z
    x2 = a.X - B.X
    x3 = a.X - c.X
    y2 = a.Y - B.Y
    y3 = a.Y - c.Y
    det# = x2 * y3 - y2 * x3
    IF det# = 0 THEN
        PRINT "--- Solution with Zpole==1 assumed impossible. Try with X or
Y==1. ---"
    ELSE
        R.X = (z2 * y3 - z3 * y2) / det#
        R.Y = (x2 * z3 - x3 * z2) / det#
        R.Z = 1
    END IF
END SUB '
-----------------------------------------------------------------------------

SUB SolveXYLSQ (p() AS VectorPos, Lines)
DIM a AS VectorPos, B AS VectorPos
    a = p(1)
    FOR i = 2 TO Lines
        B = p(i)
        x0 = a.X - B.X
        y0 = a.Y - B.Y
        x2 = x2 + x0 * x0
        y2 = y2 + y0 * y0
        xy = xy + x0 * y0
        xz = xz + x0 * (B.Z - a.Z)
        yz = yz + y0 * (B.Z - a.Z)
    NEXT i
    det# = x2 * y2 - xy * xy
    IF det# = 0 THEN
        PRINT "--- Solution with Zpole==1 assumed impossible. Try with X or
Y==1. ---"
    ELSE
        p(0).X = (xz * y2 - yz * xy) / det#
        p(0).Y = (x2 * yz - xy * xz) / det#
        p(0).Z = 1
    END IF
END SUB '
-----------------------------------------------------------------------------

DEFSNG X-Z
SUB VectorProduct (a AS VectorPos, B AS VectorPos, X AS VectorPos)
    X.X = a.Y * B.Z - a.Z * B.Y
    X.Y = a.Z * B.X - a.X * B.Z
    X.Z = a.X * B.Y - a.Y * B.X
    R# = SQR(X.X * X.X + X.Y * X.Y)
    Rad# = 180# / 3.14159265357989#
    X.Dec = Rad# * ATN(X.Z / R#)
    X.RA = Rad# * atan2#(X.Y, X.X)
    X.Pred = " "    ' To guard against garbage like hex 00 !
END SUB '
-----------------------------------------------------------------------------

SUB XYZdata (p() AS VectorPos, i)
    Rad# = 3.14159265357989# / 180#
' Adjust RA for siderial time between first and current point, since the
' telescope's pole is fixed in Alt/Az, not RA/dec !
    tdiff = (p(i).H - p(1).H) * 15 + (p(i).M - p(1).M) / 4
    tdiff = 366.24 / 365.24 * (tdiff + (p(i).S - p(1).S) / 240)
    R# = COS(p(i).Dec * Rad#)
    p(i).Z = SIN(p(i).Dec * Rad#)
    RArad# = (p(i).RA - tdiff) * Rad#
    p(i).Y = R# * SIN(RArad#)
    p(i).X = R# * COS(RArad#)
END SUB
'------------------- End of pole.bas ----------------------------



--