Download or view solarSystemGravity.frink in plain text format
// 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]
}
}
}
Download or view solarSystemGravity.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