Post by danmaz on Apr 4, 2018 19:08:11 GMT
Hi everybody
This is a not-too-long and boring astronomical program.
I tried to add proper comments to the main routines, but in order to
grasp all the stuff one should read a introductory textbook about astronomy calculus.
It calculates interplanetary distances (including Sun and Pluto) using different units of lenght,
including time of light (you must specify).
User can choose among 9 planet + sun ; all distances will be computed from that celestial body.
Comments and questions are welcome
Daniele
'PLANET_DISTANCES 2018 by danmaz(mazzad50@gmail.com)
nomainwin
Rad = 57.29577951 ' (180/pi)
'All angles in decimal degrees. They are converted in radians ( /Rad) inside the funcion sin(), cos(), tan()
'Keplerian elements and their rates, with respect to the mean ecliptic and equinox of J2000.
'Their are assumed valid for the time interval 1800AD - 3000AD
'Courtesy E M Standish, Solar System Dynamics Group, JPL, Caltech
data "Earth" , 1.00000018, 0.01673163, -0.00054346, 100.46691572, 102.93005885, -5.11260389
data -0.00000003,-0.00003661, -0.01337178,35999.37306329, 0.31795260, -0.24123856
data "Mercury", 0.38709843, 0.20563661, 7.00559432, 252.25166724, 77.45771895, 48.33961819
data 0.00000000, 0.00002123,-0.00590158,149472.67486623, 0.15940013, -0.12214182
data "Venus" , 0.72332102, 0.00676399, 3.39777545, 181.97970850, 131.76755713, 76.67261496
data -0.00000026,-0.00005107,-0.00043494, 58517.81560260, 0.05679648, -0.27274174
data "Mars" , 1.52371243, 0.09336511, 1.85181869, -4.56813164, -23.91744784, 49.71320984
data 0.00000097, 0.00009149,-0.00724757,19140.29934243, 0.45223625, -0.26852431
data "Jupiter", 5.20248019, 0.04853590, 1.29861416, 34.33479152, 14.27495244,100.29282654
data -0.00002864,-0.00018026,-0.00322699, 3034.90371757, 0.18199196, 0.13024619
data "Saturn" , 9.54149883, 0.05550825, 2.49424102, 50.07571329, 92.86136063,113.63998702
data -0.00003065,-0.00032044, 0.00451969, 1222.11494724, 0.54179478, -0.25015002
data "Uranus" ,19.18797948, 0.04685740, 0.77298127, 314.20276625, 172.43404441, 73.96250215
data -0.00196176,-0.00001550,-0.00180155, 428.49512595, 0.09266985, 0.05739699
data "Neptune",30.06952752, 0.00895439, 1.77005520, 304.22289287, 46.68158724,131.78635853
data 0.00006447, 0.00000818, 0.00022400, 218.46515314, 0.01009938, -0.00606302
data "Pluto" ,39.48686035, 0.24885238,17.14104260, 238.96535011, 224.09702598,110.30167986
data 0.00449751, 0.00006016, 0.00000501, 145.18042903, -0.00968827, -0.00809981
'Arrays dimensioned
dim planetName$ (10)
dim semiMajorAxis (9,2) 'for each of the 9 planet semiMajorAxis (n,1) and its rate (n,2)
dim eccentricity (9,2) 'for each of the 9 planet eccentricity (n,1) and its rate (n,2)
dim inclination (9,2) 'for each of the 9 planet inclination (n,1) and its rate (n,2)
dim meanLongitude (9,2) 'for each of the 9 planet meanLongitude (n,1) and its rate (n,2)
dim longitudeOfPerihelion (9,2) 'for each of the 9 planet longitudeOfPerihelion (n,1) and its rate (n,2)
dim longitudeOfAscendingNode (9,2)'for each of the 9 planet longitudeOfAscendingNode (n,1) and its rate (n,2)
dim x(10),y(10),z(10),dist(10),unit$(4)
'Arrays filled
for i=1 to 9
read a$:planetName$(i)=a$
for j=1 to 2
read a:semiMajorAxis(i,j)=a
read a:eccentricity(i,j)=a
read a:inclination(i,j)=a
read a:meanLongitude(i,j)=a
read a:longitudeOfPerihelion(i,j)=a
read a:longitudeOfAscendingNode(i,j)=a
next j
next i
planetName$(10) = "SUN":unit$(1) = "Million km":unit$(2) = "Million miles":unit$(3) = "Time of light":unit$(4) = "Astronomical Unit"
'Window is prepared
UpperLeftX = 150: UpperLeftY = 150:WindowWidth = 600: WindowHeight = 640
statictext #m.t0,"Date (mm/dd/yyyy)", 20, 40, 100, 30
textbox #m.tb0, 120, 35, 100, 25
statictext #m.t1,"Greenwich Mean Time"+chr$(13)+"(decimal hours)", 250, 35, 110, 60
textbox #m.tb1, 370, 35, 80, 25
combobox #m.combo, planetName$(), [dumb], 30, 120, 160, 40
statictext #m.t1,"Distance Unit ", 240, 85, 100, 25
combobox #m.combo1, unit$(), [dumb], 360, 80, 120, 40
button #m.calc, " Calc. ", [startCalc], UL, 270, 120, 50, 25
button #m.exit, " Exit ", [quit], UL, 350, 120, 50, 25
statictext #m.t2, "Earth", 200, 180, 60, 30
textbox #m.tb2, 300, 175, 220, 25
statictext #m.t3, "Mercury", 200, 220, 60, 30
textbox #m.tb3, 300, 215, 220, 25
statictext #m.t4, "Venus", 200, 260, 60, 30
textbox #m.tb4, 300, 255, 220, 25
statictext #m.t5, "Mars", 200, 300, 60, 30
textbox #m.tb5, 300, 295, 220, 25
statictext #m.t6, "Jupiter", 200, 340, 60, 30
textbox #m.tb6, 300, 335, 220, 25
statictext #m.t7, "Saturn", 200, 380, 60, 30
textbox #m.tb7, 300, 375, 220, 25
statictext #m.t8, "Uranus", 200, 420, 60, 30
textbox #m.tb8, 300, 415, 220, 25
statictext #m.t9, "Neptune", 200, 460, 60, 30
textbox #m.tb9, 300, 455, 220, 25
statictext #m.t10, "Pluto", 200, 500, 60, 30
textbox #m.tb10, 300, 495, 220, 25
statictext #m.t11, "SUN", 200, 540, 60, 30
textbox #m.tb11, 300, 535, 220, 25
open "SUN & PLANETS RELATIVE DISTANCES valid 1800 - 3000 AD" for window as #m
#m.combo, "font courier_new 10":#m.tb0,"!font courier_new 10":#m.tb1,"!font courier_new 10"
#m.tb3,"!font courier_new 10":#m.tb4,"!font courier_new 10":#m.tb5,"!font courier_new 10"
#m.tb6,"!font courier_new 10":#m.tb7,"!font courier_new 10":#m.tb8,"!font courier_new 10"
#m.tb9,"!font courier_new 10":#m.tb10,"!font courier_new 10":#m.tb2,"!font courier_new 10"
#m.tb11,"!font courier_new 10":#m.combo, "font courier_new 10"
#m.combo,"!--select one--"
#m.combo1,"!--select --"
#m.tb0,date$("mm/dd/yyyy")
#m.tb1,"12.00"
[dumb]
wait
[startCalc]
#m.tb0,"!contents? date1$"
#m.tb1,"!contents? time1$":Hour = val(time1$)
'this is the fundamental epoch J2000 (January,01,2000 12:00 Greenwich Mean Time, GMT)
datJ2000 = date$("01/01/2000")
'number of days (and fractions of day) from fundamental epoch J2000 (Jan,01,2000) to the day of observation
days1 = date$(date1$) - datJ2000 + (Hour-12)/24
'transforms days in Julian Centuries (1 Julian Century = 36525 days)
T = days1/36525
#m.combo, "selectionindex? index"
#m.combo1, "selectionindex? index1"
select case index1
case 1:AUcon = 149.5978707 : d1$ = " *1000000 km" ' 1 AU = millions of km
case 2:AUcon = 92.9558073 : d1$ = " *1000000 miles" ' 1 AU = millions of miles
case 3:AUcon = 8.3167464 : d1$ = " min. light time" ' 1 AU = seconds of light
case 4:AUcon = 1: d1$ = " A.U."
end select
'Sets this fundamental epoch J2000 (January,01,2000 12:00 Greenwich Mean Time, GMT)
datJ2000 = date$("01/01/2000")
'Calculates the number of days (and fractions of day) from fundamental epoch J2000 (Jan,01,2000) to the day of observation
days1 = date$(date1$) - datJ2000 + (Hour-12)/24
'Transforms days in Julian Centuries (1 Julian Century = 36525 days)
T = days1/36525
'Now we compute the value of each of the planet's six elements at the time of observation
for i=1 to 9 ' ------------------------ start cycle for each of the nine Planets ----------------------
' here we compute semi major axis, a1 , of the Planet'orbit at the time (T) of observation
a1 = semiMajorAxis(i,1) + T * semiMajorAxis(i,2)
' here we compute eccentricity, ecc , of the Planet'orbit at the time (T) of observation
ecc = eccentricity(i,1) + T * eccentricity(i,2)
' here we compute inclination, In ,of the orbit of the planet to the plane of Ecliptic
In = inclination(i,1) + T * inclination(i,2)
' here we compute La, the longitude of the ascending node (Omega Capital in most textbooks)
La = longitudeOfAscendingNode(i,1) + T * longitudeOfAscendingNode(i,2)
' here we compute the argument of perihelion, omega (omega lowercase in most texbooks)
' from the longitudeOfPerihelion, defined as the sum of omega and La
' longitudeOfPerihelion = omega +La
' omega = longitudeOfPerihelion - La
' omega and La are not on the same plane, omega being on the orbital plane while La lies on ecliptic plane
omega = longitudeOfPerihelion(i,1) + T * longitudeOfPerihelion(i,2) - La
' here we compute the mean anomaly, MA (M in most textbooks)
' meanLongitude of the planet is indeed given by the sum of longitudeOfPerihelion + mean anomaly (MA)
' meanLongitude = longitudeOfPerihelion + MA = omega + La + MA
' therefore MA = meanLongitude - omega - MA
MA = meanLongitude(i,1) + T * meanLongitude(i,2) - omega - La
' and we check that mean anomaly 'MA' is between -180 and +180
MA = ((MA+180) mod 360 ) - 180
if MA < -180 then MA = MA + 360
' here we solve KEPLER's equation
' MA = Eca - ecc*sin(Eca)
' where ecc = eccentricity MA = mean anomaly EcA = eccentric anomaly
' using decimal degrees thus obtaining eccentric anomaly
' this equation cannot be solved directly
' and so here it is solved by an iterative procedure known as Newton method
exitOK = 0
tolerance = 1e-6
dim Ei(31)
Ei(0) = MA + Rad*ecc*sin(M/Rad)
for j=0 to 30
deltaE = (MA + Rad*ecc*sin(Ei(j)/Rad) - Ei(j))/(1 - ecc*cos(Ei(j)/Rad))
Ei(j+1) = Ei(j) + deltaE
if abs(deltaE) < tolerance then exitOK = 1 :exit for
next j
EcA = Ei(j+1)
if exitOK = 0 then print "--> convergence failed !"
'now compute true anomaly v1 and planet distance from sun, r1 , in AU = 149 598 500 km
x1 = a1*(cos(EcA/Rad) - ecc)
y1 = a1*sqr(1 - ecc*ecc)*sin(EcA/Rad)
z1 = 0
v1 = atn(y1/x1)*Rad
if x1<0 and y1<0 then v1 = v1 + 180
if x1<0 and y1>0 then v1 = v1 + 180
if v1 > 360 then v1 = v1 - 360
if v1 < 0 then v1 = v1 + 360
r1 = sqr(x1*x1 + y1*y1)
'now we compute the ecliptic coordinates (x2 y2 z2) in the J2000 ecliptic plane
'the x-axis aligned toward the vernal equinox
'and origin is on the Sun (heliocentric coordinates) Ref. Montenbruck page 74
n1 = La/Rad 'n1, n2, n3 are auxiliary variables
n2 = (v1+omega)/Rad
n3 = In/Rad
x(i) = r1*(cos(n1)*cos(n2) - sin(n1)*sin(n2)*cos(n3))
y(i) = r1*(sin(n1)*cos(n2) + cos(n1)*sin(n2)*cos(n3))
z(i) = r1*(sin(n2)*sin(n3))
next i '------------------------- end of cycle for each planet ----------------------------------
x(10)=0:y(10)=0:z(10)=0 'Sun lies in the coordinate center
' Final printout
for i = 1 to 10
dist(i) = sqr((x(i)-x(index))^2 + (y(i)-y(index))^2 + (z(i)-z(index))^2)
next i
#m.tb2," ";using("#####.##" , dist(1)*AUcon);d1$
#m.tb3," ";using("#####.##" , dist(2)*AUcon);d1$
#m.tb4," ";using("#####.##" , dist(3)*AUcon);d1$
#m.tb5," ";using("#####.##" , dist(4)*AUcon);d1$
#m.tb6," ";using("#####.##" , dist(5)*AUcon);d1$
#m.tb7," ";using("#####.##" , dist(6)*AUcon);d1$
#m.tb8," ";using("#####.##" , dist(7)*AUcon);d1$
#m.tb9," ";using("#####.##" , dist(8)*AUcon);d1$
#m.tb10," ";using("#####.##", dist(9)*AUcon);d1$
#m.tb11," ";using("#####.##", dist(10)*AUcon);d1$
wait
[quit]
close #m
end