Microsoft Small Basic

Program Listing: CLJ212-3
' Orbit
' Version 0.3
' Copyright © 2016 Nonki Takahashi. The MIT License.
' Last update 2016-06-10
' Program ID CLJ212-3
' Reference:
' 1) NASA/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.3"
debug = "False"
pauseOnSuperMars = "False"
years = 100 ' duration
delay = 30 ' slow
delay = 0 ' fast
pO = "Sa" ' outer orbit to show
Init()
Form()
If debug Then
mouseDown = "False"
GraphicsWindow.MouseDown = OnMouseDown
EndIf
For T_eph = T_start To T_end
param = "slider=" + sliderID + ";pos=" + T_eph + ";"
param["show"] = "False"
Controls_MoveSlider()
ShowDate()
For iP = 1 To nP
planet = index[iP]
CalcPlanetPos()
DrawPlanet()
If debug Then
GraphicsWindow.MouseDown = OnMouseDown
mouseDown = "False"
While Not[mouseDown]
Program.Delay(200)
EndWhile
EndIf
EndFor
param["playButton"] = playButtonID
Controls_GetPlayStatus()
While status = "PAUSE"
If mouseDown Then
Controls_TogglePlayStatus()
mouseDown = "False"
EndIf
Program.Delay(200)
Controls_GetPlayStatus()
EndWhile
If mouseDown Then
Controls_TogglePlayStatus()
mouseDown = "False"
EndIf
Program.Delay(delay)
If pauseOnSuperMars And (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
EndFor
param = "playButton=" + playButtonID + ";"
Controls_PlayPause()
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)
For iP = 1 To nP
If index[iP] = pO Then
nP = iP
EndIf
EndFor
EndSub
Sub Form
xo = gw / 2
yo = gh / 2 - 30
rO = yo - 10
GraphicsWindow.BrushColor = "LightGray"
txtDate = Shapes.AddText("Date")
Shapes.Move(txtDate, 50, gh - 30)
If debug Then
txtDbg = Shapes.AddText("Debug")
Shapes.Move(txtDbg, 10, 10)
EndIf
DrawSun()
Today = Clock.ElapsedMilliseconds / 1000 / 3600 / 24 + 2415020.5 ' JD[day]
T_start = Today - years * 365 / 2
T_end = Today + years * 365 / 2
param = "x=" + 10 + ";y=" + (gh - 45) + ";width=" + (gw - 20) + ";"
param["start"] = T_start
param["end"] = T_end
Controls_AddSlider()
sliderID = id
param = "x=" + 10 + ";y=" + (gh - 30) + ";"
Controls_AddPlayButton()
playButtonID = id
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 DrawPlanet
' param planet - planet index
' param a - semi-mejor axis [au]
' param e - eccenricity
' param ϖ - longitude of perihelion [degree]
_ϖ = Math.GetRadians(ϖ)
If debug Then
buf = "a=" + a + CRLF
buf = buf + "tb["+pO+"][ao]=" + tb[pO]["ao"] + CRLF
buf = buf + "rO=" + rO + CRLF
Shapes.SetText(txtDbg, buf)
EndIf
a = a / tb[pO]["ao"] * rO ' semi-major axis
b = Math.SquareRoot(1 - e * e) * a ' semi-minor axis
If "False" And debug Then
buf = "a=" + a + CRLF
buf = buf + "b=" + b + CRLF
Shapes.SetText(txtDbg, buf)
EndIf
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
GraphicsWindow.PenColor = "Transparent"
GraphicsWindow.BrushColor = tb[planet]["color"]
If nP < 5 Then
r = tb[planet]["r"] / (tb[pO]["ao"] * au) * rO * 500
Else
r = 2
EndIf
If tb[planet]["planet"] = "" Then
tb[planet]["planet"] = Shapes.AddEllipse(2 * r, 2 * r)
EndIf
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 DrawSun
If nP < 5 Then
r = tb["Su"]["r"] / (tb[pO]["ao"] * au) * rO * 10
Else
r = 2.5
EndIf
GraphicsWindow.BrushColor = tb["Su"]["color"]
GraphicsWindow.FillEllipse(xo - r, yo - r, 2 * r, 2 * r)
EndSub
Sub JDToDate
' param jd - Julian Day [day] ≧ -32142.5 (-4801 BC/01/01)
' return date - date formatted as YYYY/MM/DD
If jd < 2299160.5 Then ' Julian calendar
d4801BC = Math.Floor(jd + 0.5) - 0.5 + 32142.5
year = (-4800) + Math.Floor(d4801BC / 365.25)
If Math.Remainder(year, 4) = 0 Then
dom[2] = 29
Else
dom[2] = 28
EndIf
day = Math.Floor(Math.Remainder(d4801BC, 365.25))
If year < 1 Then
year = (Math.Abs(year) + 1) + " BC"
EndIF
Else ' Gregorian calendar
d1200 = Math.Floor(jd + 0.5) - 0.5 - 2159350.5
y1200 = Math.Floor(d1200 / 365.2425)
year = 1200 + y1200
If Math.Remainder(year, 4) = 0 And Math.Remainder(year, 100) <> 0 Or Math.Remainder(year, 400) = 0 Then
dom[2] = 29
Else
dom[2] = 28
EndIf
nol = Math.Ceiling(y1200 / 4) - Math.Ceiling(y1200 / 100) + Math.Ceiling(y1200 / 400) ' number of leap year
day = d1200 - nol
day = Math.Floor(Math.Remainder(d1200 - nol, 365))
EndIf
For mm = 1 To 12
If day < dom[mm] Then
month = mm
mm = 12 ' exit For
Else
day = day - dom[mm]
EndIf
EndFor
If month < 10 Then
month = Text.Append(0, month * 1)
EndIf
day = day + 1
If day < 10 Then
day = Text.Append(0, day * 1)
EndIf
date = year + "/" + month + "/" + day
EndSub
Sub ShowDate
' param T_eph - Julian Day [day] ≧ -32142.5 (-4801 BC/01/01)
jd = T_eph
JDToDate()
Shapes.SetText(txtDate, date)
EndSub
Sub Controls_AddSlider
' param["x"] - left x coordinate
' param["y"] - top y coordinate
' param["width"] - width
' param["start"] - start value
' param["end"] - end value
' return id - slider ID
Stack.PushValue("local", GraphicsWindow.BrushColor)
Stack.PushValue("local", GraphicsWindow.PenWidth)
sl = param
GraphicsWindow.BrushColor = "Gray"
GraphicsWindow.FillRectangle(sl["x"], sl["y"], sl["width"], 4)
GraphicsWindow.PenWidth = 0
GraphicsWindow.BrushColor = "Red"
sl["slider"] = Shapes.AddEllipse(10, 10)
Shapes.HideShape(sl["slider"])
sl["pos"] = sl["start"]
sl["play"] = "False"
numSlider = numSlider + 1
slider[numSlider] = sl
id = "Slider" + numSlider
GraphicsWindow.PenWidth = Stack.PopValue("local")
GraphicsWindow.BrushColor = Stack.PopValue("local")
EndSub
Sub Controls_MoveSlider
' param["slider"] - slider ID
' param["pos"] - slider position
' param["show"] - "True" to show slider
iSlider = Text.GetSubTextToEnd(param["slider"], 7)
sl = slider[iSlider]
If sl["start"] <= param["pos"] And param["pos"] <= sl["end"] Then
sl["pos"] = param["pos"]
If sl["start"] < sl["pos"] Then
Stack.PushValue("local", GraphicsWindow.BrushColor)
width = sl["width"] * (sl["pos"] - sl["start"]) / (sl["end"] - sl["start"])
GraphicsWindow.BrushColor = "Red"
GraphicsWindow.FillRectangle(sl["x"], sl["y"], width, 4)
GraphicsWindow.BrushColor = "Gray"
GraphicsWindow.FillRectangle(sl["x"] + width, sl["y"], sl["width"] - width, 4)
GraphicsWindow.BrushColor = Stack.PopValue("local")
EndIf
slider[iSlider] = sl
If param["show"] Then
width = sl["width"] * (sl["pos"] - sl["start"]) / (sl["end"] - sl["start"])
Shapes.Move(sl["slider"], sl["x"] + width - 5, sl["y"] + 2 - 5)
Shapes.ShowShape(sl["slider"])
Else
Shapes.HideShape(sl["slider"])
EndIf
EndIF
EndSub
Sub Controls_AddPlayButton
' param["x"] - left x coodinate
' param["y"] - top y coodinate
' return id - play button ID
Stack.PushValue("local", GraphicsWindow.BrushColor)
Stack.PushValue("local", GraphicsWindow.PenWidth)
pb = param
GraphicsWindow.PenWidth = 0
GraphicsWindow.BrushColor = "LightGray"
pb["play"] = Shapes.AddTriangle(0, 0, 0, 16, 16, 8)
Shapes.Move(pb["play"], pb["x"], pb["y"])
pb["pause1"] = Shapes.AddRectangle(5, 16)
Shapes.Move(pb["pause1"], pb["x"], pb["y"])
Shapes.HideShape(pb["pause1"])
pb["pause2"] = Shapes.AddRectangle(5, 16)
Shapes.Move(pb["pause2"], pb["x"] + 11, pb["y"])
Shapes.HideShape(pb["pause2"])
numPlayButton = numPlayButton + 1
playButton[numPlayButton] = pb
id = "PlayButton" + numPlayButton
GraphicsWindow.PenWidth = Stack.PopValue("local")
GraphicsWindow.BrushColor = Stack.PopValue("local")
GraphicsWindow.MouseDown = OnMouseDown
EndSub
Sub OnMouseDown
mx = GraphicsWindow.MouseX
my = GraphicsWindow.MouseY
mouseDown = "True"
EndSub
Sub Controls_PlayStart
' param["playButton"] - play button ID
iPlayButton = Text.GetSubTextToEnd(param["playButton"], 11)
pb = playButton[iPlayButton]
Shapes.HideShape(pb["play"])
Shapes.ShowShape(pb["pause1"])
Shapes.ShowShape(pb["pause2"])
pb["playing"] = "True"
playButton[iPlayButton] = pb
EndSub
Sub Controls_PlayPause
' param["playButton"] - play button ID
iPlayButton = Text.GetSubTextToEnd(param["playButton"], 11)
pb = playButton[iPlayButton]
Shapes.HideShape(pb["pause1"])
Shapes.HideShape(pb["pause2"])
Shapes.ShowShape(pb["play"])
pb["playing"] = "False"
playButton[iPlayButton] = pb
EndSub
Sub Controls_GetPlayStatus
iPlayButton = Text.GetSubTextToEnd(param["playButton"], 11)
pb = playButton[iPlayButton]
If pb["playing"] Then
status = "PLAYING"
Else
status = "PAUSE"
EndIf
EndSub
Sub Controls_TogglePlayStatus
For iPlayButton = 1 To numPlayButton
pb = playButton[iPlayButton]
x1 = pb["x"]
y1 = pb["y"]
x2 = pb["x"] + 15
y2 = pb["y"] + 15
If x1 <= mx And mx <= x2 And y1 <= my And my <= y2 Then
param["playButton"] = "PlayButton" + iPlayButton
If pb["playing"] Then
Controls_PlayPause()
Else
Controls_PlayStart()
EndIf
EndIf
EndFor
For iSlider = 1 To numSlider
sl = slider[iSlider]
x1 = sl["x"]
y1 = sl["y"] + 2 - 5
x2 = sl["x"] + sl["width"]
y2 = sl["y"] + 2 + 5
If x1 <= mx And mx <= x2 And y1 <= my And my <= y2 Then
param["slider"] = "Slider" + iSlider
param["pos"] = sl["start"] + (sl["end"] - sl["start"]) * (mx - x1) / (x2 - x1)
param["show"] = "True"
Controls_MoveSlider()
T_eph = param["pos"]
ShowDate()
EndIf
EndFor
EndSub