// This program draws the current position of planets in the solar system,
// with necessary exaggerated scale.

use planets.frink
use hashColor.frink

if length[ARGS] > 0
   time = parseDate[ARGS@0]
else
   time = now[]

g = new graphics
g.antialiased[false]

planet = "neptune"
planet = input["Clip to planet:", planet]

// Set up bounds of recursion and clipping rectangle.
b = 1.1 neptunedist
punit = unit[lc[planet] + "dist"]
if punit != undef
{
   clip = 1.1 punit/au
   g.clipRectSides[-clip, -clip, clip, clip]
   b = 1.1 punit
} else
{
   println["Don't understand planet $planet, defaulting to Neptune"]
}


g.color[.5,.5,.5]
g.font["SansSerif", b/50/au]

planetpos = new array

for planet = Planet.planets
{   
   [L, B, R] = planet.getCoordinates[time]
   px = R cos[B] cos[L]
   py = R cos[B] sin[L]
   pc = hashColor[planet.name + "16"]
   //pc = new color[randomFloat[0, 1], randomFloat[0,1], randomFloat[0,1]]
   planetpos.push[[planet, px, py, pc, R]]
}

// Perform the actual recursive plotting of gravitational domains.
testRect[-b, b, -b, b, g, planetpos, 12]

for [planet, px, py, pc, R] = planetpos
{
   g.color[.5,.5,.5]
   // Draw approximation of orbit
   g.drawEllipseCenter[0,0,2 R/au,2 R/au]

   // Draw planets (each 20 times bigger than Jupiter, and that's still not
   // nearly enough, so make an even bigger circle around them.)
   g.color[0,0,0]
   g.text[" " + planet.getName[], px/au, -py/au, "left", "bottom"]
   g.fillEllipseCenter[px/au,-py/au,10 jupiterradius/au, 10 jupiterradius/au]
   g.drawEllipseCenter[px/au,-py/au,b/100/au, b/100/au]
}

// Draw the sun (20 times normal size)
g.color[.5,.5,.5]
g.fillEllipseCenter[0,0,20 sunradius/au, 20 sunradius/au]

g.show[]
//g.print[]
g.write["solarSystemGravity.svg",1024,undef,1]
g.write["solarSystemGravity.png",1024,undef,1]



/** Calculate acceleration due to planet at x, y.  The cool thing here is that
    x and y are *intervals* defining a rectangle.  This allows testing an
    infinity of points within the rectangle and finding the minimum and maximum
    acceleration due to that planet's gravity within the rectangle.
*/
accel[planetA, x, y] :=
{
   [pa, pax, pay] = planetA
   ad2 = (pax-x)^2 + (pay-y)^2
   ad2 = max[ad2, pa.radius^2] // Limit to radius of planet.  This is the magic.
   a = G pa.mass / ad2

   return a
}

/** Recursive function to test an interval containing the specified bounds.
    If a strongest solution exists, the recursion halts.  If the strongest
    acceleration is ambiguous, this breaks it down into 4 sub-rectangles and
    tests each of them recursively.  level is the maximum number of levels to
    split, so the total resolution of the final graph will be 2^level.
*/
testRect[x1, x2, y1, y2, g, planetpos, level] :=
{
   nextLevel = level - 1

   x = new interval[x1, x2]
   y = new interval[y1, y2]
   
   accvec = new array
   for p = planetpos
   {
      acc = accel[p, x, y]
      accvec.push[[p@0, acc, p@3]]
   }

   // Sort accelerations by supremum (strongest possible gravity in rectangle)
   sort[accvec, {|a,b| supremum[b@1] <=> supremum[a@1]}]

   // Test if first planet is certainly greater than second planet
   if accvec@0@1 CGT accvec@1@1
   {
      // First planet is certainly strongest, plot rect
      g.color[accvec@0@2]
      g.fillRectSides[x1/au, -y1/au, x2/au, -y2/au]
   } else
   {
      // Undecided which is stronger; subdivide rectangle
      if (nextLevel >= 0)
      {
         cx = (x1 + x2)/2
         cy = (y1 + y2)/2
         testRect[x1, cx, y1, cy, g, planetpos, nextLevel]
         testRect[cx, x2, y1, cy, g, planetpos, nextLevel]
         testRect[x1, cx, cy, y2, g, planetpos, nextLevel]
         testRect[cx, x2, cy, y2, g, planetpos, nextLevel]
      } else
      {
         // At smallest rectangle, plot as planet with largest supremum
         g.color[accvec@0@2]
         g.fillRectSides[x1/au, -y1/au, x2/au, -y2/au]
      }
   }
}