REM GE Copeland REM Dept of Physics REM Old Dominion University REM Norfolk VA REM 3/21/2001 REM 10 REM *** INPUT NUMERICAL CONSTANTS *** 20 REM *** OB = OBLIQUITY 32 OB = .40914 33 RD = 57.2957795# PI = 4! * ATN(1!) 34 TPI = 2! * PI 35 TP = TPI 40 REM *** INPUT ORBITAL ELEMENTS--ANGLES IN RADIANS *** 50 REM E = ECCENTRICITY; O1 = NODE LONGITUDE; O2 = PERIHELION LONGITUDE 60 REM IN = INCLINATION; A = SEMI-MAJOR AXIS; P = PERIOD IN YEARS 70 REM PD = JULIAN DAY OF PERIHELION; PH = ADDITIONAL DAY FRAC. AT PERI. 80 E = .967276 82 O1 = 1.01482798# 84 O2 = 1.95211743# 90 IN = 2.83160961# 92 A = 17.941104# 94 P = 75.993 100 PD = 2446470.5# 102 PH = .45174 110 GOSUB 6000 120 N = TP / P / 365.25 REM MEAN MOTION 121 PRINT CHR$(12); PRINT " THIS PROGRAM CALCULATES INFORMATION ON " 122 PRINT " FOR COMET HALLEY FOR ANY LOCATION ON EARTH. THE" 123 PRINT " OPERATION IS SELF EXPLANATORY." PRINT 124 PRINT " SOME DEFINITIONS: " PRINT " X,Y,Z ARE COORDINATES IN EQUATORIAL" 125 PRINT " SYSTEM IN A.U. R IS DISTANCE FROM SUN." 126 PRINT " DELTA IS DISTANCE FROM EARTH." PRINT 127 PRINT " INPUT YOUR LATITUDE IN DECIMAL DEGREES" 128 INPUT " NEGATIVE IF IN SOUTHERN HEMISPHERE"; LA 129 LA = LA / RD 130 INPUT " DATE OF INTEREST? (MM,DD,YYYY) "; MM, DD, YY S = 1 GOSUB 5000 140 INPUT " TIME OF INTEREST? (HH,MM)"; HH, M1 150 JH = (HH + M1 / 60) / 24 160 M = ((JD - PD) + (JH - PH)) * N 165 GOSUB 7500 REM FIND POSITION OF SUN 170 KY = 0 180 GOSUB 7000 REM SOLVE KEPLER'S EQUATION 190 GOSUB 8000 REM FIND X,Y,Z OF COMET 195 REM SC = SUN - COMET DISTANCE 200 SC = SQR(XC * XC + YC * YC + ZC * ZC) 210 X = XC + XS 212 Y = YC + YS 214 Z = ZC + ZS 220 REM CALCULATE EC = EARTH - COMET DISTANCE 230 REM AC = RIGHT ASCENSION OF COMET; DC = DECLINATION OF COMET 240 EC = SQR(X * X + Y * Y + Z * Z) 245 IF KY = 1 THEN 250 248 DM = .005772 * EC * N M = M - DM KY = 1 GOTO 180 250 SA = Z / EC CA = SQR(1 - SA * SA) GOSUB 8500 260 DC = A3 SA = Y / (EC * CA) CA = X / (EC * CA) GOSUB 8500 270 AC = A3 280 REM TH = ANGULAR SEPARATION SUN - COMET 290 CA = SIN(DC) * SIN(DS) + COS(AS1 - AC) * COS(DC) * COS(DS) 300 SA = SQR(1 - CA * CA) GOSUB 8500 310 TH = A3 320 GOSUB 9000 REM FIND SIDEREAL TIME 330 HA = ST - AC REM HA = HOUR ANGLE 332 IF HA < 0 THEN HA = HA + TPI 340 X = -COS(DC) * SIN(HA) 350 Y = SIN(DC) * COS(LA) - COS(DC) * COS(HA) * SIN(LA) 360 Z = SIN(DC) * SIN(LA) + COS(DC) * COS(HA) * COS(LA) 370 SA = Z CA = SQR(1! - SA * SA) GOSUB 8500 380 EL = A3 SA = X / CA CA = Y / CA GOSUB 8500 385 IF EL > 4.72 THEN EL = EL - TP 390 AZ = A3 400 IF JD > 2446470.5# THEN 430 410 MT = 5.47 + (5 * LOG(EC) + 11.1 * LOG(SC)) / 2.3026 420 GOTO 440 430 MT = 4.94 + (5 * LOG(EC) + 7.68 * LOG(SC)) / 2.3026 440 MN = 14.1 + (5 * LOG(EC) + 5 * LOG(SC)) / 2.3026 2000 PRINT CHR$(12); PRINT " POSITION OF COMET HALLEY:" 2010 PRINT " DATE "; MM; "/"; DD; "/"; YY; " JULIAN DAY "; JD PRINT 2020 RA = 3.81971863# REM CONVERT RADIANS TO HOURS 2022 ST = ST * RA H = INT(ST) M2 = INT(60 * (ST - H)) 2024 PRINT " LOCAL TIME: SOLAR= "; HH; " "; M1; " SIDEREAL= "; H; " "; M2 2026 PRINT 2029 XS = INT(10000 * XS) / 10000 YS = INT(10000 * YS) / 10000 ZS = INT(10000 * ZS) / 10000 2030 PRINT " X= "; XS; " Y= "; YS; " Z= "; ZS 2038 PRINT " COORDINATES OF THE SUN: " 2040 AS1 = AS1 * RA H = INT(AS1) M2 = INT(60 * (AS1 - H)) 2050 DS = DS * RD D = INT(DS) M3 = INT(60 * (DS - D)) IF D <= 90 THEN 2070 2060 D = D - 359 M3 = 60 - M3 2070 PRINT " R.A.= "; H; " "; M2; " DEC= "; D; " "; M3 2075 PRINT " COORDINATES OF THE COMET: " 2078 XC = INT(10000 * XC) / 10000 YC = INT(10000 * YC) / 10000 ZC = INT(10000 * ZC) / 10000 2080 PRINT " X= "; XC; " Y= "; YC; " Z= "; ZC 2090 AC = AC * RA H = INT(AC) M2 = INT(60 * (AC - H)) 2100 DC = DC * RD D = INT(DC) M3 = INT(60 * (DC - D)) IF D <= 90 THEN 2120 2110 D = D - 359 M3 = 60 - M3 2120 PRINT " R.A.= "; H; " "; M2; " DEC= "; D; " "; M3 2130 PRINT " R= "; INT(100 * SC + .5) / 100; 2140 PRINT " DELTA= "; INT(100 * EC + .5) / 100 2145 TH = INT(100 * TH * RD + .5) / 100 2150 PRINT " SUN-COMET ANGLE= "; TH 2160 HA = HA * RA X = 0 IF HA > 12 THEN X = 1 HA = 24 - HA 2170 H = INT(HA) M2 = INT(60 * (HA - H)) 2180 PRINT " HOUR ANGLE= "; H; " "; M2; 2190 IF X = 0 THEN PRINT " WEST" 2200 IF X = 1 THEN PRINT " EAST" 2210 AZ = AZ * RD D = INT(AZ) M3 = INT(60 * (AZ - D)) 2220 PRINT " AZIMUTH= "; D; " "; M3 2230 EL = EL * RD D = INT(EL) M3 = INT(60 * (EL - D)) 2240 PRINT " ELEVATION= "; D; " "; M3 2250 MT = INT(10 * MT) / 10 MN = INT(10 * MN) / 10 2260 PRINT PRINT " NUCLEAR MAGNITUDE= "; MN 2270 PRINT " TOTAL MAGNITUDE= "; MT 4990 INPUT " DO YOU WANT ANOTHER DATE? (Y/N) "; C$ 4992 IF C$ = "Y" THEN PRINT CHR$(12); IF C$ = "Y" THEN GOTO 130 4993 PRINT CHR$(12); PRINT " ************* DONE ***************" 4999 GOTO 11000 5000 REM ** CALCULATES JULIAN DAY ** 5010 IF MM <= 2 THEN X = 1 5020 IF MM > 2 THEN X = 0 5030 IF S = 0 THEN C = 2 5040 IF S = 1 THEN C = INT((YY - X) / 100) 5050 A1 = INT(INT(365.25 * (YY - X)) - .75 * C) 5060 B = INT(367 * ((MM - 2) / 12 + X)) 5070 IF B < 0 THEN B = B + 1 5080 JD = 1721088.5# + DD + A1 + B 5090 RETURN 6000 REM 6010 REM *** SETS UP VECTORS TO CONVERT FROM ORBITAL TO ECLIPTIC COORD.** 6020 PX = COS(O2) * COS(O1) - SIN(O2) * SIN(O1) * COS(IN) 6030 PY = COS(O2) * SIN(O1) + SIN(O2) * COS(O1) * COS(IN) 6040 PZ = SIN(O2) * SIN(IN) 6050 QX = -SIN(O2) * COS(O1) - COS(O2) * SIN(O1) * COS(IN) 6060 QY = -SIN(O2) * SIN(O1) + COS(O2) * COS(O1) * COS(IN) 6070 QZ = COS(O2) * SIN(IN) 6080 RETURN 7000 REM ***************************** 7010 REM * ITERATIVE SOLUTION OF * 7020 REM * KEPLER'S EQUATION * 7030 REM ***************************** 7060 X = M / TPI 7070 M = (X - INT(X)) * TPI 7080 EA = M 7090 EB = EA - (EA - M - E * SIN(EA)) / (1 - E * COS(EA)) 7100 IF ABS(EB - EA) / EB < 1E-08 THEN 7130 7110 EA = EB 7120 GOTO 7090 7130 RETURN 7500 REM POSITION OF THE SUN 7510 REM G1=MEAN ANOMALY; OP=PERIHELION LONG; L1-SUN'S LONGITUDE 7520 D1 = (JD - 2415020) + JH 7530 X1 = (358.4758 + .9856 * D1) / 360 7550 G1 = (X1 - INT(X1)) * 360 / RD 7560 OP = (281.2208 + .000047 * D1) / RD 7570 L1 = G1 + .03344 * SIN(G1) + .000349 * SIN(2 * G1) + OP 7580 X1 = L1 / TPI 7585 REM CALCULATE FINAL L1 AND PRECESS 7586 REM THEN SUN AND COMET REFERRED TO SAME EPOCH 7590 L1 = (X1 - INT(X1)) * TPI + .0002437 * (2433282.5# - JD) / 365.25 7595 REM SE=SUN-EARTH DISTANCE;XS,YS,ZS =RECT. COORD. OF SUN 7600 SE = .99972 / (1 + .01675 * COS(L1 - OP)) 7610 XS = SE * COS(L1) 7620 YS = SE * SIN(L1) * COS(OB) 7630 ZS = SE * SIN(L1) * SIN(OB) 7640 REM --- AS1=RIGHT ASCENSION OF SUN --- 7645 REM --- DS=DECLINATION OF SUN --- 7650 SA = ZS / SE CA = SQR(1 - SA * SA) GOSUB 8500 7660 DS = A3 SA = YS / (SE * CA) CA = XS / (SE * CA) GOSUB 8500 7670 AS1 = A3 7800 RETURN 8000 REM CALCULATES COORDINATES IN 8010 REM ECLIPTIC SYSTEM 8020 X1 = A * (COS(EA) - E) 8030 Y1 = A * SQR(1 - E * E) * SIN(EA) 8040 X = X1 * PX + Y1 * QX 8050 Y = X1 * PY + Y1 * QY 8060 Z = X1 * PZ + Y1 * QZ 8065 REM *** ROTATE TO EQUATORIAL SYSTEM *** 8070 XC = X YC = Y * COS(OB) - Z * SIN(OB) ZC = Y * SIN(OB) + Z * COS(OB) 8080 RETURN 8500 REM TWO ARGUMENT ARCTAN 8510 REM CA=COS(A);SA=SIN(A);RETURNS A IN PROPER QUADRANT 8520 IF ABS(CA) < 1E-08 THEN 8540 8530 A3 = ATN(SA / CA) GOTO 8560 8540 IF SA > 0 THEN A3 = .5 * PI 8550 IF SA < 0 THEN A3 = 1.5 * PI 8560 IF CA < 0 THEN 8590 8570 IF CA > 0 AND SA < 0 THEN 8600 8580 GOTO 8610 8590 A3 = A3 + PI GOTO 8610 8600 A3 = A3 + 2 * PI 8610 RETURN 9000 REM CALCULATE SIDEREAL TIME 9010 T1 = (JD - 2415020) / 36525 REM CENTURIES SINCE 1900 9020 X1 = (18.64606 + 2400.0513# * T1) / 24 + .5 + JH 9030 ST = (X1 - INT(X1)) * TP REM SIDEREAL TIME RADIANS 9040 RETURN 10000 T2 = ((1986 + YY) / 2 - 1900) / 100 10010 X3 = 3.07234 + (.00186 * T2) Y3 = 20.0468 - (8.500001E-03 * T2) Z3 = Y3 / 15 10020 T3 = (JD - PD) / 365.25 10030 X4 = 7.2722E-05 * T3 * (X3 + (Z3 * SIN(AC) * TAN(DC))) 10040 X5 = 4.8481E-06 * T3 * Y3 * COS(AC) 10050 AC = AC + X4 DC = DC + X5 10060 RETURN 11000 END