Microsoft Small Basic

Program Listing: CLJ212-0

List Program

' Orbit
' Version 0.2
' Copyright © 2016 Nonki Takahashi. The MIT License.
' Last update 2016-06-06
' Program ID CLJ212-0
' Reference:
' 1) JPL, Keplerian Elements for Approximate Positions of the Major Planets
' http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf
' 2) 国立天文台, 理科年表 平成28年 (Chronological Scientific Tables 2016), 2015
' http://www.rikanenpyo.jp
'
GraphicsWindow.Title = "Orbit 0.2"
debug = "False"
Init()
DrawSun()
Today = Clock.ElapsedMilliseconds / 1000 / 3600 / 24 + 2415020.5 ' JD[day]
For T_eph = Today - 365 / 2 To Today + 365 / 2
ShowDate()
If date = "2016-05-31" Then
GraphicsWindow.BrushColor = "LightGray"
txtSuperMars = Shapes.AddText("Super Mars 2016")
For xx = gw To -100 Step -1
Shapes.Move(txtSuperMars, xx , gh - 50)
Program.Delay(5)
EndFor
Shapes.Remove(txtSuperMars)
EndIf
For iP = 1 To nP - 1
planet = index[iP]
CalcPlanetPos()
DrawPlanet()
If debug Then
GraphicsWindow.MouseDown = OnMouseDown
mouseDown = "False"
While Not[mouseDown]
Program.Delay(200)
EndWhile
EndIf
EndFor
Program.Delay(30)
EndFor
Sub OnMouseDown
mouseDown = "True"
EndSub
Sub ShowDate
' param T_eph - Julian Day [day]
d2000 = T_eph - 2451545.5
year = 2000 + Math.Floor(d2000 / 365.25)
If Math.Remainder(year, 4) = 0 Then
dom[2] = 29
EndIf
day = Math.Floor(Math.Remainder(d2000, 365.25))
For month = 1 To 12
If day < dom[month] Then
Goto break
Else
day = day - dom[month]
EndIf
EndFor
break:
If month < 10 Then
month = Text.Append(0, month)
EndIf
day = day + 1
If day < 10 Then
day = Text.Append(0, day)
EndIf
date = year + "-" + month + "-" + day
Shapes.SetText(txtDate, date)
EndSub
Sub CalcPlanetPos
' param planet
' param T_eph
' In order to obtain the coordinates of one of the planets as a given Julian Ephemeris Date, T_eph,
buf = "T_eph=" + T_eph + "[year]" + CRLF
' 1. Compute the value of each of that planet's six elements: a = ao + Δa T,, etc., where T,
' the number of centuries past J2000.0, is T = (T_eph - 2451545.0) / 36525.
T = (T_eph - 2451545.0) / 36525
buf = buf + "T=" + T + "[century]" + CRLF
a = tb[planet]["ao"] + tb[planet]["Δa"] * T
buf = buf + "a=" + a + "[au]" + CRLF
e = tb[planet]["eo"] + tb[planet]["Δe"] * T
buf = buf + "e=" + e + CRLF
I = Math.Remainder(tb[planet]["Io"] + tb[planet]["ΔI"] * T, 360)
buf = buf + "I=" + I + "[degree]" + CRLF
L = Math.Remainder(tb[planet]["Lo"] + tb[planet]["ΔL"] * T, 360)
buf = buf + "L=" + L + "[degree]" + CRLF
ϖ = Math.Remainder(tb[planet]["ϖo"] + tb[planet]["Δϖ"] * T, 360)
buf = buf + "ϖ=" + ϖ + "[degree]" + CRLF
Ω = Math.Remainder(tb[planet]["Ωo"] + tb[planet]["ΔΩ"] * T, 360)
' 2. Compute the argument of perihelion, ω, and the mean anomaly, M :
' ω = ϖ - Ω ; M = L - ϖ
ωp = ϖ - Ω
buf = buf + "Ω=" + Ω + "[degree]" + CRLF
buf = buf + "ω=" + ωp + "[degree]" + CRLF
M = Math.Remainder(L - ϖ + 180, 360) - 180
' M = L - ϖ
buf = buf + "M=" + M + "[degree]" + CRLF
' 3. Modulus the mean anomaly so that -180° ≦ M < +180° and then obtain the eccentric
' anomaly, E, from the solution of Kepler's equation (see below):
' M = E - e_deg sin E,
' where e_deg = 180/π e = 57.29578 e.
e_deg = Math.GetDegrees(e)
buf = buf + "e_deg=" + e_deg + "[degree]" + CRLF
_M = Math.GetRadians(M)
EA = M + e_deg * Math.Sin(_M)
tol = "1E-6"
Calc3Equ()
While Math.Abs(ΔEA) > tol
Calc3Equ()
EndWhile
buf = buf + "E" + "=" + EA + "[degree]" + CRLF
' 4. Compute the planet's heliocentric coodinates in its orbital plane, r', with the x'-axis
' aligned from the focus to the perihelion:
' x' = a(cos E - e) ; y' = a √(1 - e^2) sin E ; z' = 0
_EA = Math.GetRadians(EA)
xh = a * (Math.Cos(_EA) - e)
yh = a * Math.SquareRoot(1 - e * e) * Math.Sin(_EA)
zh = 0
' 5. Compute the coordinates, r_ecl, in the J2000 ecliptic plane, with the x-axis aligned
' toward the equinox:
' r_ecl = _M r' ≡ _Rz (-Ω) _Rx (-I) _Rz (-ω) r'
' so that
' x_ecl = (cos ω cos Ω - sin ω sin Ω cos I) x' + (- sin ω cos Ω - cos ω sin Ω cos I) y'
' y_ecl = (cos ω sin Ω + sin ω cos Ω cos I) x' + (- sin ω sin Ω + cos ω cos Ω cos I) y'
' z_ecl = (sin ω sin I) x' + (cos ω sin I) y'
_ωp = Math.GetRadians(ωp)
_Ω = Math.GetRadians(Ω)
_I = Math.GetRadians(I)
x_ecl = (Math.Cos(_ωp) * Math.Cos(_Ω) - Math.Sin(_ωp) * Math.Sin(_Ω) * Math.Cos(_I)) * xh
x_ecl = x_ecl + (-Math.Sin(_ωp) * Math.Cos(_Ω) - Math.Cos(_ωp) * Math.Sin(_Ω) * Math.Cos(_I)) * yh
y_ecl = (Math.Cos(_ωp) * Math.Sin(_Ω) + Math.Sin(_ωp) * Math.Cos(_Ω) * Math.Cos(_I)) * xh
y_ecl = y_ecl + (-Math.Sin(_ωp) * Math.Sin(_Ω) + Math.Cos(_ωp) * Math.Cos(_Ω) * Math.Cos(_I)) * yh
z_ecl = (Math.Sin(_ωp) * Math.Sin(_I)) * xh
z_ecl = z_ecl + (Math.Cos(_ωp) * Math.Sin(_I)) * yh
buf = buf + "x_ecl=" + x_ecl + "[au]" + CRLF
buf = buf + "y_ecl=" + y_ecl + "[au]" + CRLF
buf = buf + "z_ecl=" + z_ecl + "[au]" + CRLF
If debug Then
Shapes.SetText(txtDbg, buf)
EndIf
' 6. If desired, obtain the equatorial coodinates in the "ICRF" or "J2000 frame", r_eq :
' x_eq = x_ecl
' y_eq = + cos ε y_ecl - sin ε z_ecl
' z_eq = + sin ε y_ecl + cos ε z_ecl
' wherre the obliquility at J2000 is ε = 23°.43928.
EndSub
Sub Calc3Equ
' Solution of Kepler's Equation, M = E - e_deg sin E
' Given the mean anomaly, M, and the eccentricity, e_deg, both in degrees, start with
' E0 = M + e_deg sin M
' and iterate the following three equations, with n = 0, 1, 2, ..., until |ΔE| ≦ tol (nothing
' that e_deg is in degrees; e is in radians):
' ΔM = M - (En - e_deg sin En) ; ΔE = ΔM / (1 - e cos En) ; En+1 = En + ΔE.
' For the aproximate formulae in this present context, tol = 10^-6 degrees in sufficient.
_EA = Math.GetRadians(EA)
ΔM = M - (EA - e_deg * Math.Sin(_EA))
ΔEA = ΔM / (1 - e * Math.Cos(_EA))
EA = EA + ΔEA
buf = buf + "ΔM=" + ΔM + CRLF
buf = buf + "ΔE=" + ΔEA + CRLF
buf = buf + "E" + "=" + EA + CRLF
EndSub
Sub DrawSun
xo = gw / 2
yo = gh / 2
pO = "Ma"
rO = yo - 10
r = tb[index[nP]]["r"] / (tb[pO]["ao"] * au) * rO * 10
GraphicsWindow.BrushColor = tb[index[nP]]["color"]
GraphicsWindow.FillEllipse(xo - r, yo - r, 2 * r, 2 * r)
For iP = 1 To nP
If index[iP] = pO Then
nP = iP + 1
EndIf
EndFor
EndSub
Sub DrawPlanet
' param planet - planet index
' param a - semi-mejor axis [au]
' param e - eccenricity
' param ϖ - longitude of perihelion [degree]
_ϖ = Math.GetRadians(ϖ)
a = a / tb[pO]["ao"] * rO ' semi-major axis
b = Math.SquareRoot(1 - e * e) * a ' semi-minor axis
If tb[planet]["orbit"] = "" Then
GraphicsWindow.PenWidth = 1
GraphicsWindow.PenColor = "#333333"
GraphicsWindow.BrushColor = "Transparent"
tb[planet]["orbit"] = Shapes.AddEllipse(2 * a, 2 * b)
xc = xo + e * a * Math.Cos(_ϖ + PI)
yc = yo - e * a * Math.Sin(_ϖ + PI)
x = xc - a
y = yc - b
Shapes.Move(tb[planet]["orbit"], x, y)
Shapes.Rotate(tb[planet]["orbit"], -ϖ)
EndIf
If tb[planet]["planet"] <> "" Then
Shapes.Remove(tb[planet]["planet"])
EndIf
GraphicsWindow.PenColor = "Transparent"
GraphicsWindow.BrushColor = tb[planet]["color"]
r = tb[planet]["r"] / (tb[pO]["ao"] * au) * rO * 500
tb[planet]["planet"] = Shapes.AddEllipse(2 * r, 2 * r)
x = xo + x_ecl / tb[pO]["ao"] * rO
y = yo - y_ecl / tb[pO]["ao"] * rO
Shapes.Move(tb[planet]["planet"], x - r / 2, y - r / 2)
EndSub
Sub Init
CRLF = Text.GetCharacter(13) + Text.GetCharacter(10)
Not = "False=True;True=False;"
PI = Math.Pi
dom = "1=31;2=28;3=31;4=30;5=31;6=30;7=31;8=31;9=30;10=31;11=30;12=31;"
gw = 598
gh = 428
GraphicsWindow.Width = gw
GraphicsWindow.Height = gh
GraphicsWindow.BackgroundColor = "Black"
au = "1.4959787E8" ' [km]
' tb[][] - planet elements table
' Me - Mercury
' V - Venus
' E - Earth
' Ma - Mars
' J - Jupter
' Sa - Saturn
' U - Uranus
' N - Neptune
'
' Table 1 1)
' Keplerian elements and their rates, with respect to the mean ecliptic
' and equinox of J2000, valid for the time-interval 1800 AD - 2050 AD.
' ao, Δa - semi-major axis [au, au/century]
' eo, Δe - eccentricity [ , /century]
' Io, ΔI - inclination [degree, degree/century]
' Lo, ΔL - mean longitude [degree, degree/century]
' ϖo, Δϖ - longitude of perihelion [degree, degree/century]
' Ωo, ΔΩ - longitude of the ascending node [degree, degree/century]
tb["Me"] = "ao= 0.38709927;eo= 0.20563593;Io= 7.00497902;Lo= 252.25032350;ϖo= 77.45779628;Ωo= 48.33076593;"
tb["Me"] = tb["Me"] + "Δa= 0.00000037;Δe= 0.00001906;ΔI=-0.00594749;ΔL=149472.67411175;Δϖ= 0.16047689;ΔΩ= -0.12534081;"
tb["V"] = "ao= 0.72333566;eo= 0.00677672;Io= 3.39467605;Lo= 181.97909950;ϖo=131.60246718;Ωo= 76.67984255;"
tb["V"] = tb["V"] + "Δa= 0.00000390;Δe=-0.00004107;ΔI=-0.00078890;ΔL= 58517.81538729;Δϖ= 0.00268329;ΔΩ= -0.27769418;"
tb["E"] = "ao= 1.00000261;eo= 0.01671123;Io=-0.00001531;Lo= 100.46457166;ϖo=102.93768193;Ωo= 0.0;"
tb["E"] = tb["E"] + "Δa= 0.00000562;Δe=-0.00004392;ΔI=-0.01294668;ΔL= 35999.37244981;Δϖ= 0.32327364;ΔΩ= 0.0;"
tb["Ma"] = "ao= 1.52371034;eo= 0.09339410;Io= 1.84969142;Lo= -4.55343205;ϖo=-23.94362959;Ωo= 49.55953891;"
tb["Ma"] = tb["Ma"] + "Δa= 0.00001847;Δe= 0.00007882;ΔI=-0.00813131;ΔL= 19140.30268499;Δϖ= 0.44441088;ΔΩ= -0.29257343;"
tb["J"] = "ao= 5.20288700;eo= 0.04838624;Io= 1.30439695;Lo= 34.39644051;ϖo= 14.72847983;Ωo=100.47390909;"
tb["J"] = tb["J"] + "Δa=-0.00011607;Δe=-0.00013253;ΔI=-0.00183714;ΔL= 3034.74612775;Δϖ= 0.21252668;ΔΩ= 0.20469106;"
tb["Sa"] = "ao= 9.53667594;eo= 0.05386179;Io= 2.48599187;Lo= 49.95424423;ϖo= 92.59887831;Ωo=113.66242448;"
tb["Sa"] = tb["Sa"] + "Δa=-0.00125060;Δe=-0.00050991;ΔI= 0.00193609;ΔL= 1222.49362201;Δϖ= -0.41897216;ΔΩ= -0.28867794;"
tb["U"] = "ao=19.18916464;eo= 0.04725744;Io= 0.77263783;Lo= 313.23810451;ϖo=170.95427630;Ωo= 74.01692503;"
tb["U"] = tb["U"] + "Δa=-0.00196176;Δe=-0.00004397;ΔI=-0.00242939;ΔL= 428.48202785;Δϖ= 0.40805281;ΔΩ= 0.04240589;"
tb["N"] = "ao=30.06992276;eo= 0.00859048;Io= 1.77004347;Lo= -55.12002969;ϖo= 44.96476227;Ωo=131.78422574;"
tb["N"] = tb["N"] + "Δa= 0.00026291;Δe= 0.00005105;ΔI= 0.00035372;ΔL= 218.45945325;Δϖ= -0.32241464;ΔΩ= -0.00508664;"
' r - equatorial radius of the planet [km] 2)
' color - color for the planet
tb["Su"] = "r=696000;color=#FFFFCC;"
tb["Me"] = tb["Me"] + "r= 2440;color=LightGray;"
tb["V"] = tb["V"] + "r= 6052;color=CornSilk;"
tb["E"] = tb["E"] + "r= 6378;color=#336699;"
tb["Ma"] = tb["Ma"] + "r= 3396;color=#FF6633;"
tb["J"] = tb["J"] + "r=71492;color=Tan;"
tb["Sa"] = tb["Sa"] + "r=60268;color=DarkKhaki;"
tb["U"] = tb["U"] + "r=25559;color=SkyBlue;"
tb["N"] = tb["N"] + "r=24764;color=RoyalBlue;"
nP = Array.GetItemCount(tb)
index = Array.GetAllIndices(tb)
GraphicsWindow.BrushColor = "LightGray"
txtDate = Shapes.AddText("Date")
Shapes.Move(txtDate, gw - 100, 10)
If debug Then
txtDbg = Shapes.AddText("Debug")
Shapes.Move(txtDbg, 10, 10)
EndIf
EndSub