Articles Hierarchy

Articles Home » Processing » SOLarsystem 2

SOLarsystem 2

and planets elliptic orbit


is now a (TAB) menu page from my sketch_3D_view play tool.

some time ago i play processing.js and try a planet system view but i give up,
as i question my presettings about the real constellation,
meaning a start position for our planets on a specific day, from where my calc would work.
So, even if wrong, still a nice show with some real data and nice relative movements.
a other big question was about the varying speed on a ellipsoid orbit.

Keyboard OPERATION: 3D view SOLAR System
[Tab] for axis; [r] for view reset
for planet info use [m][v][e][M][c][j][s][u][n][p]
[t] for set time = now, [UP][DOWN] for time +- day, [RIGHT][LEFT] for time +- month
to make it all visible, the planets are upsized ( rel. to the orbits), and the sun is downsized ( from the planets), and the moon orbit upsized.

for history, info, limitations ... see original project and original online view here

now all operation also with buttons:

and with the new PAN ( full PTZ operation ) // and recheck on RPI 3B

still use my own buttons, but clean up Main and SOL parameter menu by new button functions.

and planets elliptic orbit

what a poor show:
it is 2018 AD and my planets go on elliptic orbits by a constant angle every day,
( well, better as using circles )

while 1602 AD Kepler solved that problem about elliptic orbits already, more here and here

the Antikythera mechanism is supposed to have a gear for elliptic movements, that was in the year 100..200 BC

and in the aegypt / syria / china / middle america / have been cultures what recorded astronomical data for up to 700 years continuously , and unless that was not only dumb accounting i think these people knew more about astronomy as you and me.

or lunar calender from 32000 BC see here

so i start a test about elliptic movement graph
but will need parameter to play with, like orbital period and major distance Aphelion
and minor distance ( or the eccentricity to calc it from ).
here a play tool doing same, you can step one day and get same angle,
start from Mercury data but can set eccentricity and reduce orbital days.

play this
but this was only the start point for the real math/physics problem i have:

at a different distance the planet has a different speed:
v = sqrt[GM*(2/r - 1/a)]
EarthAt perihelion, Earth's distance from the Sun is r=a(1-e) and at aphelion, it's r=a(1+e).
G=6.673*10^-11 N m2/kg2
M=1.989*10^30 kg
a=1.496*10^11 m
So plugging in the numbers, the speed at perihelion is 30,300 m/s and at aphelion it's 29,300 m/s.
from here

play this
while above formula calculates a speed i have a question to its vector,
when "alpha" is the angle from the center of the ellipsis, "tau" is the angle from the SUN foci point.
the speed vector flight direction (yellow line) is 90 deg off from "alpha" or "tau" ( red or blue ( i used here ) line)?

now i try to calculate ( in mini steps ( 0.5 deg / 720 steps )) position speed,
and integral/sum way and time.
by pack it all in arrays.
the result is impressive,
the way S ( around the ellipsis) ends up as 2.406 AU
the time T ends up at 87.964 d at pa = 0.387098 AU
not bad for Mercury : 87.969 d
+ + with pa , pe from EARTH also get correct info!
now i must think how to show this data.
idea is to select a day ( of T_orbit ) and show the position and the "covered" area ( of last day) in red.

++ usable in a wider range of pa with auto zoom / stepsize
++ add mouse wheel zoom
++ with speed and day graph ( over angle )
+++ also show the day difference to the NOT Kepler movement ( const angle per day )
yes, that looks like a sinus and as i showed in my old spreadsheet can be used as easy correction instead of this real INTEGRAL way of calculation.
( hundred years old mechanical tower clock used that already ) see here and here

but depending on the data / numbers it might not be so accurate,
( the data are accurate, but the show methode expects to show data from selected day next.000
while the 0.5deg / 720 rec array might contain data only for next.000 to next.999 / if at all / so up to one day off )
but how to find the exact day/to/day change if the array contains n floats for one day
or even if the array not contain that day ( because a outer planet need more than one day for 0.5deg )?
a interpolation between 2 records might be needed to find the planets position ... for (new) day.
the actual show version 5 looks ( on "+/-" click ) for a record with the "next" day. as next.xxxx [d]

so i should be "go back in mini-mini" steps between [alpha] and [alpha-1] to find a position what fits better to next.000x [d]
but instead doing interpolation/itteration in every visualisation loop,
( what would make it a impossible sollution for the big SOLarsystem code with all planets ),
i have new idea:
at begin ( and after change of pa "Semi-major axis" or pe "Eccentricity" ) i recalculate the whole array anyway.
i should in a second run repair / overwrite the record containing next.xxxx [d] ( first of an other day )
with correct data ( of x,y,r,alpha,tau,v,S ) for the improved T[but same i/record] next.000x [d] fit.
after this "repair" array same long, but not so aeqidistant ( in angle ) as now
delta alpha[i-1] , alpha[i] , alpha[i+1] not same anymore.
i call this concept:

back swinging door

one of the best ways for this itteration is the methode
dalpha = alpha[i] - alpha[i-1],
alphatest = alpha[i-1]+ dalpha / 2
test on T(alphatest), if dx/2 ok, or decide go up or down and use dx/4 .....
that is not a interpolation as it not uses a simpler ( linear) math between 2 points,
it just try to find a better alpha[i] where T[i] fits next.000x d
the error limit for ( next.xxxx - next.000x ) is adjustable and has great influence in the computational amount.
also a better start as dx/2 might help ( here a little interpolation could help), but is also risky.
more see binary-search ( on sorted lists )

worst case would be a planet with pa = 1.57 AU ( like MARS 1.52 AU )
what has a orbital period of about 720d
where every record needs that optimization ( alpha shift), because possibly could be
T[i-1] 1.5,
T[i] 2.5,... [d]
so actually all that becomes little bit useless.

while i was working on this i lost my windows PC ( harddrive? windows? )
so if the HD is recoverable there is a old ellipsis_v6
and a not saved ellipsis_v6 temp file somewhere ?processing?.
if not i have to start on a other system from ellipsis_v5 again.

in ellipsis 6.3 i try a easy repair, one run with a interpolated new alpha
and do all calculations with that again.
add i make a button to enable / disable this.

also a new button for a autorun PLAY [> true] ( and new day(s) every 0.5 sec )
so when pe >0.7 and autorun true can see the changing speed very nicely.

i think first i stop here and try to put what i have already back into the SOLarsystem code,
worry the array size and the calculation job is too much for all planets.
but sure, as this all is just a processing test, no problem to reduce the array length (720) and live with a bigger integration error.

but i can not copy this code 11 times (planets array) into, i need to make it more smart,
6 times array [11,720] or array [11,6,720] here
or possibly learn about class / structured records first. here and here
start with the class array for the planets data first, as test:

for the big one
[ 11 planets ] [ 720 angle steps ] [ 6 orbital parameter ] // what works as float array
i got stuck at second step after first array from class works,

now i find that info
but actually it was not that difficult,
from class_array_5
orbit_rec[][] O;

class orbit_rec {
float pAlpha, pX, pY, pR, pS, pT;
orbit_rec(float talpha, float tx, float ty, float tr, float ts, float tt) {
pAlpha = talpha; pX = tx; pY = ty; pR = tr; pS = ts; pT = tt;
void display() { point(pX * rzoom, pY * rzoom); }

void setup() {
O = new orbit_rec[planets][istep];
for (int p = 0; p < planets; p = p+1) {
for (int i = 0; i < istep; i = i+1) {
O[p][i] = new orbit_rec(0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
} // create and fill with 0.0

for (int p = 1; p < 3; p = p+1) { // here have only planets 1 and 2
for (int i = 0; i < istep; i = i+1) {
O[p][i].pAlpha = alphastep * i;
O[p][i].pX = P[p].pa*(cos(alphastep * i)) + P[p].pf; // AU
O[p][i].pY = -P[p].pb*(sin(alphastep * i)); // AU but y pos for pixel "up"
O[p][i].pR = pow((pow((P[p].pa*P[p].pe + O[p][i].pX),2)+pow(O[p][i].pY,2)),0.5); // AU from SUN to Planet
} // again run over planets and angles to calc the position and more
} // end setup

void draw() {
for (int p = 1; p < 3; p = p+1) { //planets
for (orbit_rec thisPlanet : O[p]) { thisPlanet.display(); } // make a dot for each planet and anglestep
} //planets

the code little bit long and ugly as it contains prior good test code, just to document valid syntax... pls. try this
you might say planet 0 == SUN is a wasted array,
but for me it just means there is spare
for the 2 (or N ) mass movement Barycenter
and more far reaching, the Sunsystem orbit in Milkyway

so, back, but to what?
the tool i want to work on is the
sketch_3D_view esp. the TAB "SOL"
because i now have learned how to use
* * array of class for the
*** P[p] planets data and for the
*** O[p][i] planets orbits over 720 angle steps
* * and i coded the orbits with the position / speed / way / time / integral
acc. to Kepler second law correct.
* * and i have a "one step only" next day correction ( what improves the position indication for the inner planets only )

and that i have to build into sketch_3D_view.

while the download here has the old version V 1.4.3
i have from my corrupted windows hard disk ( and strange behaving Google Drive )
a version what contains same SOL.pde and a SOL(1).pde ( bigger and about 12 hours "newer" )
best i can do is to make a
sketch_3D_view_v150, try to get SOL(1).pde renamed as SOL.pde running and checked, prior to get started.
Now start declaring and filling / initial calculating / the Planets array to replace the old orbital parameter arrays,
and besides that was lots of work, i run into many traps,
- - the index shift ( mercury[0] to sun[0] ) was a bad one.
- - i tried to put the color into the class array but all come out black,
put R,G,B as float in array and found that my calls draw planet / draw moon
now even more bad coding regarding handling the planets data.
+ + but as i build up that array i checked all data again v.s. WIKI
and give more decimales as used in the old arrays.
- - - and see a bad error
i see that the orbits neptun and pluto because of plutos eccentricity nearly touch,
i investigated if that could be correct and yes, so i checked out a date i found when the planets "switched" position ( between 7. Februar 1979 and 11. Februar 1999 Pluto closer to sun )
and checked at my tool, looks like i am 180deg OFF in the start position (J2000) mean anomaly understanding?.
need to check later

in 1.5.0.b i enable planet 1 .. 11 ( spare config as
earth II / 180 deg OFF / but with a ecc 0.7 /
to see what i am up to ( by kepler ellipsis speed )
now earth II passes mercury / venus / earth / mars orbit
i have no place for a add enable button, so i enable it at start and disable it with the reset button!

anyhow that TEST / PLAY planet pay off, when the array calc ( at init ) was running ( v1.5.1 )
there is a problem about the speed calc. a recheck show that with increasing eccentricity the speed change overshoots and even run into a numeric error, like for EARTH II over eccentricity (pe) 0.5
why i not see that in my elliptic test program? a code transfer error?
i not find here, so i check here
Mercury Earthcompare_____________Mercury | Earth
Mean orbital velocity (km/s) 47.36 | 29.78
Max. orbital velocity (km/s) 58.98 | 30.29
Min. orbital velocity (km/s) 38.86 | 29.29
and data for MIN and MAX for Mercury and Earth by my tool correct
and a earth with ecc 0.8 calculates / even i have nothing to compare with.
so dig deep to find the transfer error, i used the focus double for dX [ solved ]

an other way of check is the Kepler 3. law ( semi major axis >> T orbit )
so play with eccentricity should not effect the orbital time
now i list Tsum ( Kepler integral ) and Torg ( from internet )
and even not same number, compare
EARTH _ Tsum 365.283 Torg 365.256 ( 40min error / year ? )
EARTH II Tsum 365.283 Torg 365.256 ( earth with ecc 0.7 and 180deg shift)

now i have the old ( const angle per day ) show
the array with the Kepler / position / speed / way / time / info.
how to show that? can show both at the beginning to compare?

start to show the ellipsis by the pX pY points from array ( in red )
and see how it match
+ after adjust same zoom factor from px [AU] to screen pix
+ change position: translate(+-focus,0),
as old program use elliptic center and a,b BUT
array use SUN center prX, prY.
+ use RED color for the array data

now for using the array data to show the position i first need to find the i ( in istep 720 )
where the day fits to the actual day ( counted from after J2000 )
first i choose the record before ( later i might use the optimization tools from the ellipsis code)
but as i first need to do the J2000 position angle:
what is called Mean Anomaly / a angle ( left wise ) starting at the "pericenter"
what means without it the (all) planet would be at the Perihelion ( the short R1 )
while my array orbital year starts at the Aphelion ( long R2 )
and with this i found my 180deg mistake.

anyhow, the job is to go first to the start position at J2000, and then go to the actual day as days after J2000 ( - x * Torbit ), while the sum of this 2 positions might again be bigger a Torbit.
while working on this i noticed that i could not "go" back J2000 ( the day indication still correct, but the planets not moving / this was never a problem with the old code )

but esp. for earth II with its ecc pe = 0.7 the speed thing is clearly visible,
even the steps ( days (365) or month(12) / by mouse better keyboard [arrow UP or RIGHT]) are not perfect here.

for the earth MOON thing i also need a new code,
in the old program i used a copy of the first draw_planet code (hard coded again for the earth)
and go to its position. from there i call the draw_planet code for the moon.
this time i think it is more easy,
i want while i loop all planet just remember the "thisi" of earth position as (global) "thisi_earth"
so later for the moon i just try to do a

translate(O[earth][thisi_earth].pX,O[earth][thisi_earth].pY); // from sun position jump to actual earth center
translate ( back )
(well it was a little bit more complicated..)

but there is a other trick needed inside the draw_planet_Kepler, because for the moon i need to cheat a bigger orbit, as the real orbit is inside the oversized earth sphere, so would be invisible.

also i need the white line
SUN - earth AND earth - moon
in case earth is selected
end rev 1.5.1

an other old open point is the 3d axis rotation / where prX was 0.0 and rotateX disabled.
want to check on that again, as i forget what the problem was when i wrote that code ?
as i have the EARTH II for play i try how a 90.0deg. angle for prX would look:
that makes sense, but does not fit to "Argument of periapsis" description.
and "Argument of perihelion" like MERCURY 29.124deg also not. well i might even have to switch it with "Longitude of ascending node" ( and rethink the sequence for that three rotations. )
ok, first i collect all that internet data for the planets "Argument of perihelion" into prX array to get a idea / what a disaster

i leave the data in but disabled the rotateX until i know what i am talking about.
The angle picture from WIKI, i see again and again, still not understand, try to find other sources like this
if the "Argument of perihelion" i get from the WIKI ( and i just put into prX for test )
is a rotation of the ellipsis in the direction of the rotation ( so for a circle it would just add up to the position )
i need to change my code and use it in rotateY ( put it into prY )
and use then the "Longitude of ascending node" in prX.
it has nothing to do with astronomical understanding, but now i have the 3 lists of rotation in my array i can have a better look.
the "Longitude of ascending node" now prY are small values ( max 12 deg ), are a much more acceptable disturbance of the usually flat orbits disk.
while the "Argument of perihelion" now prX can have any value of 360deg and might be a added rotation what should happen in my prY ( rotateY ) plane only.
and that is what i test first:
exchange the prX prY columns in the array and start with disable rotate X and Z ( all should be flat in one plane )
after that test i enable the other rotations again and do tests if the sequence of these rotations make any difference.
well, it did not play out, i rotated the data twice and ended up again with a disabled rotateX

thought of the day:HA, i even might have a bigger mistake, actually i do a 180deg correction in the position ( for the J2000 start at the small axis / perhelion. and i added it to the start,
but possibly i should rotate the ellipsis? or you get up and sit on the chair on the other side at the table? or start the same program 6 month later?

it also has to do with the reference?
my ellipsis has the sun as left focal center and the long axis "Aphelion" to the right ( +x ) direction
( what per definition should point where? to the galactic center? NO possibly the "Perihelion" should point there? as the angles seems to be defined for it.)
and the XY plane ( processing screen X(-Y) ) is what? a mean orbital plane for our sun or a absolute plane for the galactic spiral plane?

here i found some text info for the picture from WIKI

as for me also the processing rotation thing is difficult i try this, also want confirm about the sequence for the 3 rotations i asked here

for the OPERATION PTZ i use already a PUSH POP Matrix,
but as i was told at the forum that command can be stacked ( used again inside),
so i would not have to back rotate the coordinate system after drawing each planet orbit.

also in the last play code i packed all the mouse operation PTZ to a extra TAB what cleaned up the main draw loop.
that is a nice revision for the main project "sketch_3D_view"
here 1.5.3 on a RPI:
DEBIAN stretch / Linux version 4.14.62-v7+ / Raspberry Pi 3 Model B Rev 1.2


to show a SOLarsystem you need:

+ SUN and Planets ( with the weight and radius info )

+ + and elliptic orbits for the planets, defined by
pa Semi-major axis
pe Eccentricity
( and the from this calculated focus, Semi-minor axis, R1 R2 .. and even orbital period )
and 3 angles for the xyz rotation of the ellipsis

+ + + and the position of the planet on this elliptic orbit
by a start position ( angle ) at J2000 and the days between NOW and J2000.
the assumption of "constant angle per day" only works for circular orbits,
for a Kepler movement of the planets ( higher speed when close to the SUN )
need to calculate distance ( to SUN ) and speed for each position
and integral dS and dT over the ellipsis, and save all to a ARRAY.

So when you later ask for a certain day what position that planet is?
* find what day of the planets year it is,
* * walk through the array to find a record that has a DAY of the planets year
close to the one you look for,
* * * and use the position info of that record to show it.

(+) if the eccentricity of a orbit is bigger ( for a test set to > 0.7 ), the speedy Kepler movement is good visible.


for the elliptic Kepler orbit tool play this
for the elliptic orientation try this
for the SOLarsystem go back here