Download or view fitPlutoOrbit.frink in plain text format
/** This is an attempt to find generalized ellipse equation that matches
Pluto's orbit. For now, we're ignoring the latitude and only plotting
the longitude.
We take Pluto's position at 5 points using planets.frink and fit an
ellipse to those points. We then analyze the ellipse to find easy ways
to plot it.
See https://www.wolframalpha.com/input?i=x%5E2+%2B+0.9354427275709442482+y%5E2+%2B+-0.0046031330525519022+x+y+%2B+-13.348467654240820399+x+%2B+-12.492529944210118768+y+%2B+-1309.2827048558484214+%3D+0
https://www.anirdesh.com/math/algebra/rotation-of-conics.php
https://www.anirdesh.com/math/algebra/ellipse-rotation.php
*/
use planets.frink
use fitConicSection.frink
use functionUtils.frink
use allTransforms.frink
start = #1886#
end = #2098#
timestep = (end-start) / 4
args = new array
t = start
while t <= end
{
[L,B,R] = Planet.Pluto.getCoordinates[t]
[x,y,z] = heliocentricToXYZ[L,B,R]
args.push[x]
args.push[y]
t = t + timestep
}
//println[args]
// All dimensions are dividided by au to give dimensionless numbers.
[x1,y1,x2,y2,x3,y3,x4,y4,x5,y5] = divide[args, au]
[B,C,D,E,F] = fitPoints[x1,y1,x2,y2,x3,y3,x4,y4,x5,y5]
A = 1
// Give the equation of the ellipse in general text form
// A x^2 + B x y + C y^2 + D x + E y + F === 0
println["Pluto orbit equation: (units are au)"]
println[fitPointsToText[x1,y1,x2,y2,x3,y3,x4,y4,x5,y5]]
// Find functions for plotting the curves given x.
println["\nFunctions in x:"]
funcs = fitPointsToYFunctions[x1,y1,x2,y2,x3,y3,x4,y4,x5,y5]
println[formatTable[mapList["inputForm", funcs], "left"]]
// Find functions for plotting the curves given y.
println["\nFunctions in y:"]
funcs = fitPointsToXFunctions[x1,y1,x2,y2,x3,y3,x4,y4,x5,y5]
println[formatTable[mapList["inputForm", funcs], "left"]]
d = now[]
// d = # 1992 October 13 TD # // Testing against Meeus example 37.a
// Find current coordinates of Pluto.
[long,lat,radius] = Planet.Pluto.getCoordinates[d]
// Display current heliocentric coordinates (long, lat, radius) also called
// (L,B,R)
// Compare against http://cosinekitty.com/solar_system.html
println["\nCurrent position (heliocentric coordinates):"]
println["longitude= " + (long->"deg") + "\nlatitude = " + (lat->"deg") + "\nradius = " + (radius->"au")]
[x,y,z] = divide[heliocentricToXYZ[long,lat,radius], au]
// Find current heliocentric cartesian (x,y,z) coordinates
println["\nCurrent Cartesian coordinates (au)"]
println["x=$x, y=$y, z=$z"]
// Characterize the ellipse, finding its rotation angle, center, eccentricity,
// and an equation for translated conic section (centered around 0,0)
[theta, cx, cy, eccentricity, A1, B1, C1, D1, E1, F1, a, b] = characterizeEllipse[A,B,C,D,E,F]
println["\ntheta is " + (theta->"deg") + " counterclockwise"]
println["\nCenter is $cx, $cy"]
println["\nEccentricity is $eccentricity"]
println["\nTranslated coefficients: A=$A1, B=$B1, C=$C1, D=$D1, E=$E1, F=$F1"]
println["\nx axis radius: $a"]
println["y axis radius: $b"]
// Automatically translate conic section (to center it around 0,0)
[A2,B2,C2,D2,E2,F2] = centerConic[A,B,C,D,E,F]
println["\nCentered coefficients: A=$A2, B=$B2, C=$C2, D=$D2, E=$E2, F=$F2"]
// Derotate centered conic given the rotation angle found above
[A3,B3,C3,D3,E3,F3] = rotateConic[A2,B2,C2,D2,E2,F2,-theta]
//println["\nDerotated coefficients (manually): A=$A3, B=$B3, C=$C3, D=$D3, E=$E3, F=$F3"]
// Automatically derotate the conic
[A4,B4,C4,D4,E4,F4] = derotateConic[A2,B2,C2,D2,E2,F2]
println["\nDerotated coefficients: A=$A4, B=$B4, C=$C4, D=$D4, E=$E4, F=$F4"]
// Renormalize the coefficients so that A=1
[A5,B5,C5,D5,E5,F5] = normalizeConic[A4,B4,C4,D4,E4,F4]
println["\nNormalized coefficients: A=$A5, B=$B5, C=$C5, D=$D5, E=$E5, F=$F5"]
g = new graphics
// Use interval arithmetic to plot original best equation
g1 = plotPoints[x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,-31,45,-32,46,true]
g.add[g1]
// Use interval arithmetic to Draw Pluto's orbit as if it were centered.
g.color[1,0,0]
g2 = plotConic[A2,B2,C2,D2,E2,F2,-31,45,-32,45]
g.add[g2]
// Use interval arithmetic to raw pluto's orbit as if it were centered,
// derotated, and normalized
g.color[1,1,0]
g5 = plotConic[A5,B5,C5,D5,E5,F5,-31,45,-32,45]
g.add[g5]
// Make a simple algorithm (using only translations, rotations, and an ellipse)
// to draw Pluto's orbit.
g6 = new graphics
g6.saveTransform[]
g6.translate[cx, -cy] // y is down in Frink 2-D graphics so invert y
g6.rotate[0,0,theta] // This is clockwise rotation which negates the ccw
g6.color[1,0,0]
g6.drawEllipseCenter[0,0,2a,2b]
g6.restoreTransform[]
g.add[g6]
// Draw other planet orbits
g.color[0,0,0]
pnames = ["mercury", "venus", "earth", "mars", "jupiter", "saturn", "uranus", "neptune"]
for name = pnames
{
diam = 2 unit["${name}dist"] / au
g.drawEllipseCenter[0,0,diam,diam]
}
g.color[0,1,0,.7]
g.fillEllipseCenter[x, -y, 1, 1]
// Sun
g.color[1,1,0]
g.fillEllipseCenter[0, 0, 1/2, 1/2]
// Draw to files
g.write["fitPlutoOrbit.png", 1000, undef]
g.write["fitPlutoOrbit.svg", 1000, undef]
g.write["fitPlutoOrbit.html", 1000, undef]
g.invertGrays[].show[]
//g.print[]
Download or view fitPlutoOrbit.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen, eliasen@mindspring.com