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 ---------------------------- --