Quote:
Originally Posted by Blue Skies
For solar system objects you need to get a minimum of 3 positions and then some maths is done to compute the orbit. I'm still to learn how that all works as well, but the more observations the better the result.
|
Hi Jacquie,
If you are interested in this topic, have a look at this program below (written in BASIC) that computes preliminary
orbit of Solar system object.
The program is written by
Erich Karkoschka (from Arizona University) back in '80, during a Summer Astronomical camp in Adriatic... I believe it is in public domain since it was published in a summary journal after the occasion.
I tested it (by taking positions of planet from ephemerides, then calculated
orbit, then calculated positions from that
orbit) and it works pretty accurately, however it must be compiled and run in HU-BASIC (that works with 16-digit mantissa) or better otherwise the rounding errors accumulate and the results are not valid, or iteration does not converge at all.
Positional data in a form of Date (string), Time [hours] RA, DEC are at the end of the program.
When finished, the program prompts for saving results in a file.
I found this way convenient when testing....
5 DEFDBL A-H
7 DEFDBL O-Z
10 CLS : PRINT " Program for estimating preliminary orbital "
11 PRINT "elements for members of Solar system from 3 observations": PRINT
12 PRINT "Method developed by Erich Karkoschka, 1987?)": PRINT
14 DIM M$(12): DIM A$(3)
15 M$(1) = "Jan": M$(2) = "Feb": M$(3) = "Mar": M$(4) = "Apr"
16 M$(5) = "May": M$(6) = "Jun": M$(7) = "Jul": M$(8) = "Aug"
17 M$(9) = "Sep": M$(10) = "Oct": M$(11) = "Nov": M$(12) = "Dec"
20 LN = .434294481#: ST = 57.29577951#: DN = .99726957#
22 PI = 3.141592654#: P2 = 2 * PI: EK = 23.44579 / ST
24 SE = SIN(EK): CE = COS(EK)
30 PRINT : REM : INPUT "Name of object ? "; N$
40 DEF FNASN (X#) = ATN(X# / SQR(-X# * X# + 1))
42 DEF FNACS (X#) = -ATN(X# / SQR(-X# * X# + 1)) + PI / 2
45 REM GOTO 130
50 FOR K = 1 TO 3: REM *** TEST***
60 REM ************************
70 READ A$(K), T, RH(K), RM(K), DD(K), DM(K): GOSUB 1940: NEXT
80 REM ************************
85 GOTO 200
90 REM *** INPUT DATA ***
130 HOME: FOR K = 1 TO 3
140 PRINT K; ".": PRINT : INPUT "Datum: (G MJ D) "; A$(K): PRINT
150 IF MID$(A$(K), 5, 1) <> " " THEN 140
160 INPUT "Vrijeme: (UT: Hr,min): "; TH, TM: PRINT : T = TH + TM / 60
170 INPUT "Ra: (h,min,sec )"; RH(K), RM(K), RS(K): PRINT
180 INPUT "Dec: (deg,min,sec )"; DD(K), DM(K), DS(K): PRINT
190 GOSUB 1940: NEXT K
200 CLS : PRINT N$: PRINT : FOR K = 1 TO 3: PRINT : PRINT " "; K; ". ";
210 PRINT A$(K); " (JD="; T(K); ")": PRINT
220 PRINT "Ra= "; RH(K); "hr "; RM(K); "min "; " Dec= "; DD(K); "deg "; DM(K); "min": PRINT
230 NEXT: PRINT "Are the input data O.K. ? (Y/N ?)"
235 INPUT Y$: IF Y$ <> "Y" AND Y$ <> "y" THEN 30
240 O$ = "Eliptic iteration"
250 PRINT : PRINT " Eqx 1950/2000 ?(1/2)": INPUT S$: IF S$ = "1" OR S$ = "2" THEN 270
260 GOTO 250
270 IF S$ = "1" THEN EQ = 1950: GOTO 290
280 IF S$ = "2" THEN EQ = 2000: GOTO 290
290 CH = 0: ON ERROR GOTO 2030
300 REM * Calculation *
310 REM * 1 *
320 FOR K = 1 TO 3
330 B(K) = FNASN(CE * SIN(D(K)) - SE * COS(D(K)) * SIN(R(K)))
340 L(K) = 2 * ATN((SIN(D(K)) * SE + COS(D(K)) * CE * SIN(R(K))) / (COS(B(K)) + COS(D(K)) * COS(R(K))))
350 IF L(K) < 0 THEN L(K) = L(K) + P2
360 IF EQ = 2000 THEN L(K) = L(K) + .7 / ST
370 NEXT
380 REM * 2 *
390 M0 = (T(2) - 2445338.087#) * .985626283# / ST
400 V0 = M0 + (1.915# * SIN(M0) + .02 * SIN(2 * M0)) / ST
410 L0 = V0 - 77.126 / ST
420 R0 = .9997200000000001# / (1 + .016716# * COS(V0))
430 LV = (M0 + V0) / 2 + 12.874 / ST
440 V0 = .01720209895# / R0
450 IF EQ = 1950 THEN L0 = L0 - .7 / ST: LV = LV - .7 / ST
460 REM * 3 *
470 T1 = T(1) - T(2)
480 T3 = T(3) - T(2)
490 REM * 4 *
500 L(1) = L(1) - L(2)
510 L(3) = L(3) - L(2)
520 L0 = L0 - L(2)
530 LV = LV - L(2)
540 REM * 5 *
550 FOR K = 1 TO 3 STEP 2
560 S(K) = 1 / (SIN(B(2)) * TAN(B(K)) + COS(B(2)) * COS(L(K)))
570 X(K) = -S(K) * SIN(L(K))
580 Y(K) = S(K) * (COS(B(2)) * TAN(B(K)) - SIN(B(2)) * COS(L(K)))
590 NEXT
600 REM * 6 *
610 V = SIN(B(2)) / TAN(L0)
620 W = 2 * COS(B(2)) * COS(L0)
630 REM * 7 *
640 F = -(Y(3) - V * X(3)) / (Y(1) - V * X(1))
650 REM * 8 *
660 G = (F * X(1) + X(3)) / (F * T1 * T1 + T3 * T3)
670 REM * 9 *
680 H = (1 / T3 + 1 / ABS(T1)) / (1 + F)
690 REM * 10 *
700 XP = (X(3) - G * T3 * T3) * H
710 YP = (Y(3) - G * V * T3 * T3) * H
720 ZP = (1 / ABS(T1) - F / T3) / (1 + F)
730 REM * 11 *
740 MU = .0002958#
750 H1 = COS(B(2)): H2 = R0 * COS(L0)
760 H3 = R0 * SIN(L0)
770 H4 = SIN(B(2))
780 H5 = ZP * COS(B(2)) - YP * SIN(B(2))
790 H6 = V0 * COS(LV)
800 H7 = V0 * SIN(LV)
810 H8 = ZP * SIN(B(2)) + YP * COS(B(2))
820 IF LEFT$(O$, 1) = "P" THEN GOTO 920
830 REM * Eliptic iteration *
840 G = (2 * G * R0 * R0 * R0) / (MU * SIN(L0))
850 D = 1.5: CLS : PRINT TAB(10); O$
860 FD = ((1 + D * D - W * D) ^ (-1.5) - 1) / D + G
870 FS = -3 * (2 * D - W) * (1 + D * D - W * D) ^ (-2.5) / 2 / D - (1 + D * D - W * D) ^ (-1.5) / D / D + (1 / D / D)
880 D1 = D - FD / FS
890 IF ABS(D - D1) > .0000000000001# THEN D = D1: GOTO 860
900 D = R0 * D
910 GOTO 1030
920 REM * Parabolic iteration *
930 D = 1.2: PRINT TAB(10); : PRINT O$: PRINT TAB(10);
940 U = SQR((D * H1 - H2) * (D * H1 - H2) + H3 * H3 + D * D * H4 * H4)
950 US = (D * (H1 * H1 + H4 * H4) - H1 * H2) / SQR((D * H1 - H2) * (D * H1 - H2) + H3 * H3 + D * D * H4 * H4)
960 FS = US * (D * H5 - H6) * (D * H5 - H6) + U * 2 * (D * H5 * H5 - H5 * H6)
970 FS = FS + US * (-D * XP - H7) * (-D * XP - H7) + U * 2 * (D * XP * XP + XP * H7)
980 FS = FS + US * D * D * H8 * H8 + U * 2 * D * H8 * H8
990 FD = U * (D * H5 - H6) * (D * H5 - H6) + U * (-D * XP - H7) * (-D * XP - H7) + U * D * D * H8 * H8 - 2 * MU
1000 D1 = D - FD / FS
1010 IF ABS(D1 - D) > .0000000000001# THEN D = D1: GOTO 940
1020 REM
1030 REM * 12 *
1040 ON ERROR GOTO 2030
1050 X = D * H1 - H2
1060 Y = -H3
1070 Z = D * H4
1080 R = SQR(X * X + Y * Y + Z * Z)
1090 REM * 13 *
1100 PX = D * H5 - H6
1110 PY = -D * XP - H7
1120 PZ = D * H8
1130 V = SQR(PX * PX + PY * PY + PZ * PZ)
1140 Y1 = R * V * V - 2 * MU
1150 REM * 14 *
1160 XS = Y * PZ - Z * PY
1170 YS = Z * PX - X * PZ
1180 ZS = X * PY - Y * PX
1190 H = SQR(XS * XS + YS * YS + ZS * ZS)
1200 IF LEFT$(O$, 1) = "P" THEN AO = 0: EO = 1: NO = 0: QO = H * H / 2 / MU: GOTO 1290
1210 REM * 15 *
1220 AO = 1 / (2 / R - V * V / MU)
1230 REM * 16 *
1240 NO = SQR(MU / AO / AO / AO)
1250 REM * 17 *
1260 EO = SQR(1 - H * H / MU / AO)
1270 QO = AO * (1 - EO)
1280 REM * 18 *
1290 IO = FNACS(ZS / H)
1300 REM * 19 *
1310 OO = L(2) + SGN(XS) * FNACS(-YS / H / SIN(IO))
1320 IF OO < 0 THEN OO = (OO + P2)
1330 IF OO > P2 THEN OO = (OO - P2)
1340 REM * 20 *
1350 VO = SGN(X * PX + Y * PY + Z * PZ) * FNACS((H * H / MU / R - 1) / EO)
1360 REM * 21 *
1370 WO = SGN(Z) * FNACS((Y * XS - X * YS) / R / H / SIN(IO)) - VO
1380 IF WO < 0 THEN WO = WO + P2
1390 IF WO > P2 THEN WO = WO - P2
1400 IF LEFT$(O$, 1) = "P" THEN T = T(2) - (R * SIN(VO) / 3 / H) * (4 * QO - R * COS(VO)): GOTO 1470
1410 REM * 22 *
1420 EC = 2 * ATN(SQR((1 - EO) / (1 + EO)) * TAN(VO / 2))
1430 REM * 23 *
1440 MO = EC - EO * SIN(EC)
1450 REM * 24 *
1460 T = T(2) - MO / NO
1470 JD = T: GOSUB 1800
1480 REM * Printing resuts *
1490 REM * K *
1500 PRINT TAB(10); : PRINT N$: PRINT "Orbital elements :"
1505 PRINT : PRINT
1510 IF AO <> 0 THEN PRINT " a ="; AO; " AU"
1520 PRINT " q ="; QO; " AU"
1530 PRINT " e ="; EO
1540 IO = IO * ST: PRINT " i ="; IO; " deg"
1550 OO = OO * ST: PRINT " o ="; OO; " deg"
1560 WO = WO * ST: PRINT " w ="; WO; " deg"
1570 PRINT " T ="; JJ; " "; M$(MM); " "; TT; " (JD ="; JD; ")": PRINT
1580 JD = JD - 2433282
1590 PRINT "Do you want to save those data on disk ? (Y/N ?)"
1595 INPUT Y$: IF Y$ <> "Y" THEN 1640
1597 PRINT
1600 'D$ = CHR$(4): PRINT D$; "OPEN"; N$
1610 'PRINT D$; "WRITE"; N$
1620 'PRINT QO: PRINT EO: PRINT IO: PRINT OO: PRINT WO: PRINT JD
1630 'PRINT D$; "CLOSE"; N$
1640 IF LEFT$(O$, 1) = "P" THEN 1670
1650 PRINT
1660 PRINT " n ="; NO * ST; " (deg/24h)"
1670 PRINT : PRINT " Anomalies for T2 ( "; A$(2); " )": PRINT
1680 PRINT " True anomaly An ="; VO * ST; " deg": PRINT
1690 IF LEFT$(O$, 1) = "P" THEN 1710
1700 PRINT " Excentr. anomaly E ="; EC * ST; " deg": PRINT
1710 PRINT " Mean anomaly M ="; MO * ST; " deg"
1712 PRINT : PRINT "More (Y/N) ?"
1717 INPUT Y$: IF Y$ = "Y" THEN RUN
1720 STOP
1730 REM * date -> jul. day *
1740 REM * *
1750 IF MM > 2 THEN C = JJ: D = MM
1760 IF MM <= 2 THEN C = JJ - 1: D = MM + 12
1770 A = INT(C / 100): B = 2 - A + INT(A / 4)
1780 JD = INT(365.25# * C) + INT(30.6001# * (D + 1)) + TT + 1720994.5# + B
1790 RETURN
1800 REM * Jul. day -> date *
1810 JD = JD + .5: Z = INT(JD): F = JD - INT(JD)
1820 IF Z < 2299161 THEN A = Z
1830 IF Z > 2299161 THEN A = INT((Z - 1867216.25#) / 36524.25): A = Z + 1 + A - INT(A / 4)
1840 B = A + 1524
1850 C = INT((B - 122.1) / 365.25)
1860 D = INT(365.25 * C)
1870 E = INT((B - D) / 30.6001)
1880 TT = B - D - INT(30.6001 * E) + F
1890 IF E < 13.5 THEN MM = E - 1
1900 IF E > 13.5 THEN MM = E - 13
1910 IF MM > 2.5 THEN JJ = C - 4716
1920 IF MM < 2.5 THEN JJ = C - 4715
1930 RETURN
1940 REM
1950 REM * Date from A$ into JJ,MM,TT *
1960 JJ = VAL(LEFT$(A$(K), 4)): MM = VAL(MID$(A$(K), 6, 2)): TT = ABS(VAL(MID$(A$(K), 8, 26)))
1970 TT = TT + T / 24
1980 GOSUB 1730
1990 T(K) = JD
2000 R(K) = 15 * (RH(K) + RM(K) / 60 + RS(K) / 3600) / ST
2010 SG = SGN(DD(K)): D(K) = (DD(K) + SG * DM(K) / 60 + SG * DS(K) / 3600) / ST
2020 RETURN
2030 REM * Changes *
2040 IF CH = 1 THEN GOTO 2090
2050 CH = 1: ON ERROR GOTO 10
2060 PRINT "
ORBIT IS NOT ELIPTIC"
2070 O$ = "Parabolic iteration"
2080 GOTO 300
2090 REM * NOT POSSIBLE ! ***
2100 REM
2110 FLASH: PRINT "
ORBIT DOES NOT EXIST !"
2120 STOP
2130 REM
2140 REM ***********************
2150 REM T E S T
2160 REM ***********************
2170 DATA "1983-05-29",10.08,17,57.8,-3,26
2180 DATA "1983-06-05",0.0,17,52.5,-3,21
2190 DATA "1983-06-15",0.0,17,43.1,-3,29