/** 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[]