I'm (still) working on a virtual solar system. I barely finished college algebra (so the symbols confuse me. What the heck is α? I don't have a button for a fish on my keyboard! ;-)). I look at so many sites (DreanInCode, (Game Dev)StackOverflow, Wiki) to get help with the math; I find it, but can't understand it. Code... code... half my kingdom for code!

Another way to say it is, their are dozens (hundreds, thousands) of people doing the same thing - solar system/space sims; people keep asking the same questions. I know... we will break the internet if we solved everyone's physics/coding questions.

If we had a tab on the physics pages with solved coded problems, any student/person could go their and find code. Yes, students should work out their own solutions but, us recreational coders could make something without re-inventing the wheel. This:

https://www.math.ubc...-01a/orbits.pdf

is what I'm working on now. I found info hear:

https://en.wikipedia...unction_of_time

2/3rds of the way down is what I'm looking for but, solving it and turning it onto code is hard for me.

Hummmm... I'll make a post in questions with my current problem (and I uploaded my test application to Github). Oh, I got a link to that already. I just feel that i'm asking the same question again.

https://github.com/Mipada/OrbitTest

It's the OrbitControl.java calculateOrbitElements() and calculateOrbitElements2() I'm stuck on.

# Plotting planet orbits by time

Page 1 of 1## 12 Replies - 813 Views - Last Post: 04 August 2018 - 09:47 AM

##
**Replies To:** Plotting planet orbits by time

### #2

## Re: Plotting planet orbits by time

Posted 01 August 2018 - 11:19 PM

What language?

Is this a for pay job, or is there a question?

Is this a for pay job, or is there a question?

### #3

## Re: Plotting planet orbits by time

Posted 01 August 2018 - 11:44 PM

I found this great article about how to plot a planet's orbit in 3D by time. I understand most of it but get lost right after 'Plotting Planets' (Hommage to the old TV show Get Smart. Chief: Did you get all of that. Max: Not all of it Chief. Chief: What part didn't you get? Max: The part after, listen to this.) Only thing is, I barely finished college algebra and can't figure some things out.

https://www.math.ubc...-01a/orbits.pdf

So, I'm stuck converting formulas to code.

Is

v = w * x

Vector3d v = new Vector3d(w.x * x.x, w.y * x.y, w.z * x.z);//

And is

xu + yv + zw

float f = x.x * u.x + y.y * v.y + z.z * w.z;//a float?

for example. javax.vecmath.Vector3d doesn't have a vector multiply function.

I have a link to it and a test program on Github.

https://github.com/Mipada/OrbitTest

It's in JavaMonkeyEngine so the control gets attached to the node and the code gets called every cycle. It's Y is up/down.

I'm stuck at calculateOrbitElements(), calculateOrbitElements2() and getTimeTranslation().

Thanks for the replies.

Can someone help me turn the expressions into code. I can write functions so, how to I get the initial value of T, the mean anomaly and the eccentric anomaly in code?

Ack, so sorry. I shouldn't have posted so late in the evening. I added the link to the source of the equations.

'''

https://www.math.ubc...-01a/orbits.pdf

So, I'm stuck converting formulas to code.

Is

v = w * x

Vector3d v = new Vector3d(w.x * x.x, w.y * x.y, w.z * x.z);//

And is

xu + yv + zw

float f = x.x * u.x + y.y * v.y + z.z * w.z;//a float?

for example. javax.vecmath.Vector3d doesn't have a vector multiply function.

I have a link to it and a test program on Github.

https://github.com/Mipada/OrbitTest

It's in JavaMonkeyEngine so the control gets attached to the node and the code gets called every cycle. It's Y is up/down.

I'm stuck at calculateOrbitElements(), calculateOrbitElements2() and getTimeTranslation().

Thanks for the replies.

Can someone help me turn the expressions into code. I can write functions so, how to I get the initial value of T, the mean anomaly and the eccentric anomaly in code?

Ack, so sorry. I shouldn't have posted so late in the evening. I added the link to the source of the equations.

'''

/* * To change this license header, choose License Headers in Project Properties. * To change this template file, choose Tools | Templates * and open the template in the editor. */ //package lib.my.app.controls; package my.app.scene.controls; import java.text.DecimalFormat; import javax.swing.text.NumberFormatter; import javax.vecmath.Vector3d; //import javax.vecmath.Vector3f; import my.app.nodes.MassNode; //import lib.my.app.utilities.Utils; import my.app.nodes.Node; /** * * @author Me */ //a semi-major axis //b semi-minor-axis //e eccentricity //Om + w = II (argument of Perihelion) //L(t) = II + M(t) //i inclination //t time public class OrbitControl extends AbstractControl { //Any local variables should be encapsulated by getters/setters so they //appear in the SDK properties window and can be edited. //Right-click a local variable to encapsulate it with getters and setters. static public int X = 0; static public int Y = 1; static public int Z = 2; static int count = 0; int index = 0; static int a_ = 0; static int e_ = 1; static int Ω_ = 2; static int II_ = 3; static int L_ = 4; static int i_ = 5; String name = ""; static int num = 10;//number of bodies static String[] bodyName = new String[num]; static double[][] orbit = new double[num][6]; static int UP = 1;//to adjusst for different 3D libs //arguements private double a = 0d;//semiMajorAxis private double e = 0d;//eccentricity private double T = 1L;//time at perihelion //calculate private double b = 0d;//semiMinorAxis private double f = 0d;//the distance from the centre of the orbit to the Sun private double P = 365.25d;//period (days?) private double i = 0d;//inclination private double Ω = 0d;//longitude of ascending node private double ω = 0d;//angle between ascending node and perihelion private double II = 0d;//II = Ω + ω (0 for Earth) private double L = 0d;//L(t) = II + M(t) //for 3D protected Vector3d localTranslation = new Vector3d();//world trans protected Vector3d worldTranslation = new Vector3d();//world trans protected boolean changed = false; protected Vector3d u = new Vector3d(1d, 0d, 0d);//perihelion (right) protected Vector3d w = new Vector3d(0d, 1d, 0d);//up protected Vector3d v = new Vector3d(0d, 0d, 0d);//positive y/z axis protected Vector3d α = new Vector3d(1d, 0d, 0d); // Vector3d VStar = new Vector3d(); Vector3d V0 = new Vector3d(); Vector3d V1 = new Vector3d(); double M0 = 0d; double M = 0d;//??????????? double E = 0d;//??????????????? // protected double angleOfAscendingNode = 0d;//Autumnal equinox protected Vector3d ascendingNode = new Vector3d(0d, 0d, 0d); protected double angleOfDescendingNode = 180d;//Vernal equinox protected Vector3d descendingNode = new Vector3d(0d, 0d, 0d); protected double angleOfPariapsis = 180d; protected Vector3d pariapsis = new Vector3d(0d, 0d, 0d); protected double angleOfApoapsis = 180d; protected Vector3d apoapsis = new Vector3d(0d, 0d, 0d); DecimalFormat degreeFormat = new DecimalFormat("000.00◦"); NumberFormatter degreeFormatter = new NumberFormatter(degreeFormat); DecimalFormat distanceFormat = new DecimalFormat("00.00o"); NumberFormatter distanceFormatter = new NumberFormatter(distanceFormat); //private double rad;//angle from parent //private double time = 0d;//current time //plane is x, z (up is sun polor up/down) public OrbitControl(String name){ this(name, 0d, 0d, 0d, 0d, 0d); } //plane is x, z (up is sun polor up/down) public OrbitControl(String name, double semiMajorAxis, double eccentricity, double inclination){ this(name, semiMajorAxis, eccentricity, inclination, 0d, 0d); } public OrbitControl(String[] data){ this(data[0]); } public OrbitControl(String name, double semiMajorAxis, double eccentricity/*degrees*/, double inclination/*degrees*/, double longitudeOfAscendingNode, double argumentOfPerihelion){ this.name = name; bodyName[index] = this.name; this.a = semiMajorAxis; this.e = Math.toRadians(eccentricity);//EDIT: convert degs to rads this.i = Math.toRadians(inclination);//EDIT: convert degs to rads this.Ω = longitudeOfAscendingNode; this.ω = argumentOfPerihelion; // index = count++; calculateOrbitElements(); } //Calculate; //mean anomaly (M) //eccentric anomaly (E) public void calculateOrbitElements(){ b = a * Math.sqrt(1d - (e * e)); f = (a * e);//average distance (kilomters?), AUs for now P = Math.sqrt(Math.pow(a, 3));//Kepler's 3rd law//EDIT: Kepler II = Ω + ω; // T = 0.5d;//EDIT: test value for T doulbe T0 = 0.5d;//EDIT: Set value for T0 M0 = 1d; System.out.println("OrbitControl.getTimeTranslation()" + "(" + name + ")" + ",T=" + T + ",P=" + P); M = ((Math.PI * 2d) * ((T0/T)/P));//mean anomaly ???????????? E = M + 0d;//eccentric anomaly ??????????? calculateOrbitElements2(); } public void updateOrbitElements(){ //set new parameters due to thrust //i.e., set new a, b, e //calculateOrbitElements2(); } //calculate the 3D elements public void calculateOrbitElements2(){ // //vectors if (UP == Y){ u.set(1d, 0d, 0d); w.set(1d, 0d, 0d); α.set(Math.cos(Ω), 0d, Math.sin(Ω)); α.normalize(); v.set(u.x * w.x, u.y * w.y, u.z * w.z); α.cross(V1, VStar); v.cross(α, V1); } else if (UP == Z) {//up = 2 u.set(1d, 0d, 0d); w.set(0d, 1d, 0d); α.set(Math.cos(Ω), Math.sin(Ω), 0d); α.normalize(); v.set(u.x * w.x, u.y * w.y, u.z * w.z); α.cross(V1, VStar); v.cross(α, V1); } else{ System.out.println("bad UP code (only 1 (Y axis) or 2 (Z axis))"); System.exit(1); } System.out.println("OribtControl.calculateElements2() " + "α=" + α + ",v=" + v); // //orbit table //a, e, Ω, II, L, i orbit[index][a_] = a; orbit[index][e_] = e; orbit[index][Ω_] = Ω; orbit[index][II_] = II; orbit[index][L_] = L; orbit[index][i_] = i; } //Plotting in 3D //you can change the return signiture to return Vector3ds or Vector3fs //Solor system plane is x,z (with +y going out the sun's north pole) //Get the translation at a certain time //you can change the return signiture to return Vector3ds or Vector3fs //Solor system plane is x,z (with +y going out the sun's north pole) //@Param time //@Param array for return values (can be changed to Vector3f, Vector3d //@Param pTrans coordinate of parent public Vector3d getTimeTranslation(double T0, Vector3d pTrans){ double x, y, z; //if (name.equals("Sun")){ // return worldTranslation; //} // M0 = 360d * (T0 - T)/P; T = T0 - (M0/360d) * P;//seconds // if (UP == 1){ //α = new Vector3d(Math.cos(Ω), 0d, Math.sin(Ω)); } else{ //α = new Vector3d(Math.cos(Ω), Math.sin(Ω), 0d); } v.x = w.x * u.x; v.y = w.y * u.y; v.z = w.z * u.z; localTranslation.set(v.x, v.y, v.z); worldTranslation.add(pTrans, localTranslation); changed = true; if (name.equals("Earth")) System.out.println("OrbitControl.getTimeTranslation() " + "(" + name + ") " + worldTranslation); return worldTranslation; } static public void printTable(){ System.out.println("Planet a e Ω Π L i"); for (int i = 0;i < testData.length;i++){ String s = getPlanetData(i); System.out.println(s); } } //Class Formatter //◦ static private String getPlanetData(int i){ //Class Formatter //System.out.println("name=" + bodyName[i]); String s = String.format("%-8s %6.3f %5.3f %6.1f◦ %6.2f◦ %5.1f◦ %4.2f◦", testData[i][0], orbit[i][0], orbit[i][1], orbit[i][2], orbit[i][3], orbit[i][4], orbit[i][5]); return s; } public void setIndex(int i){ this.index = i; } static public void setTable(){ for (int j = 0;j < testData.length;j++){ bodyName[j] = testData[j][0]; } for (int j = 0;j < testData.length;j++){ for (int i = 0;i < testData[0].length - 1;i++){ orbit[j][i] = Double.valueOf(testData[j][i + 1]); } } } @Override protected void controlUpdate(float tpf) { //TODO: add code that controls Spatial, //e.g. spatial.rotate(tpf,tpf,tpf); //System.out.println("OribtControl.controlUpdate() " + "(" + name + ")"); if (changed){ ((MassNode)spatial).setTranslation(worldTranslation); changed = false; } } static String[][] testData = { // {"Planet", "a", "e", "Ω", "Π", "L", "i"}, {"Mercury", "0.387", "0.206", "48.3", "77.46", "252.3", "7.00"}, //◦ {"Venus", "0.723", "0.007", "76.7", "131.6", "182.0", "3.39"}, {"Earth", "1.000", "0.017", "00.00", "102.9", "100.5", "0.00"}, {"Mars", "1.524", "0.093", "49.6", "336.1", "355.4", "1.85"}, {"Jupiter", "5.203", "0.048", "100.4", "14.3", "34.4", "1.30"}, {"Saturn", "9.555", "0.056", "113.7", "93.1", "50.1", "2.49"}, {"Uranus", "19.22", "0.046", "74.0", "173.0", "314.1", "0.77"}, {"Neptune", "30.11", "0.009", "131.8", "48.1", "304.3", "1.77"} }; /* Planet a e Ω Π L i Mercury 0.387 0.206 48.3◦ 77.46◦ 252.3◦ 7.00◦ Venus 0.723 0.007 76.7◦ 131.6◦ 182.0◦ 3.39◦ Earth 1.000 0.017 − 102.9◦ 100.5◦ 0.00◦ Mars 1.524 0.093 49.6◦ 336.1◦ 355.4◦ 1.85◦ Jupiter 5.203 0.048 100.4◦ 14.3◦ 34.4◦ 1.30◦ Saturn 9.555 0.056 113.7◦ 93.1◦ 50.1◦ 2.49◦ Uranus 19.22 0.046 74.0◦ 173.0◦ 314.1◦ 0.77◦ Neptune 30.11 0.009 131.8◦ 48.1◦ 304.3◦ 1.77◦ */ /* @Override protected void controlRender(RenderManager rm, ViewPort vp) { //Only needed for rendering-related operations, //not called when spatial is culled. } @Override public Control cloneForSpatial(Spatial spatial) { OrbitControl control = new OrbitControl(name, a, e, i, Ω, ω); //TODO: copy parameters to new Control return control; } @Override public void read(JmeImporter im) throws IOException { super.read(im); InputCapsule in = im.getCapsule(this); //TODO: load properties of this Control, e.g. //this.value = in.readFloat("name", defaultValue); } @Override public void write(JmeExporter ex) throws IOException { super.write(ex); OutputCapsule out = ex.getCapsule(this); //TODO: save properties of this Control, e.g. //out.write(this.value, "name", defaultValue); } */ }

This post has been edited by **PGR**: 02 August 2018 - 07:18 PM

### #4

## Re: Plotting planet orbits by time

Posted 02 August 2018 - 12:27 AM

PGR, on 02 August 2018 - 07:44 AM, said:

v = w * x

for example. javax.vecmath.Vector3d doesn't have a vector multiply function. And

xu + yv + zw

What?

for example. javax.vecmath.Vector3d doesn't have a vector multiply function. And

xu + yv + zw

What?

It's not really clear what your question is. Where do these expressions come from? What are the terms here? Apologies, for I've forgotten pretty much all the orbital mechanics I did as an undergrad (which granted wasn't a lot).

Remember also that there are two kinds of multiplication of vectors: the scalar (or dot) product and the vector (or cross) product. It's not clear from your notation which you're referring to.

Also, line 138: that comment should really say "Kepler's 3rd law".

### #5

## Re: Plotting planet orbits by time

Posted 02 August 2018 - 09:07 AM

PGR, on 02 August 2018 - 01:44 AM, said:

v = w * x

for example. javax.vecmath.Vector3d doesn't have a vector multiply function. And

for example. javax.vecmath.Vector3d doesn't have a vector multiply function. And

Indeed, because Java doesn't override operators...

So, I'm afraid you'll have to learn some math. For your 3D vectors, start here: https://en.wikipedia...tion_of_vectors

Given you want a vector as a result, you probably want a cross product.

So,

Quote

public final void cross(Vector3d v1, Vector3d v2)

Sets this vector to the vector cross product of vectors v1 and v2.

-- https://docs.oracle....ecmath.Vector3d,%20javax.vecmath.Vector3d)

Sets this vector to the vector cross product of vectors v1 and v2.

-- https://docs.oracle....ecmath.Vector3d,%20javax.vecmath.Vector3d)

### #6

## Re: Plotting planet orbits by time

Posted 02 August 2018 - 02:53 PM

ndc85430, on 02 August 2018 - 07:27 AM, said:

Where do these expressions come from? What are the terms here?

OP made another related post in caffeine lounge about this, there OP links to the source of the equations.

https://www.dreaminc...i-pages-please/

### #7

## Re: Plotting planet orbits by time

Posted 02 August 2018 - 05:36 PM

I added the source of the equations. To clarify, I'm having trouble converting the equations to code. For example, how do I solve for the original value for T? If I put in 0 (zero) I get a divide by zero error. I guess 360/2Pi is the same direction (degree/rad) so, I'll try using that to get my value for T.

I'm not sure how to get the mean anomaly or the eccentric anomaly.

I'm not sure how to get the mean anomaly or the eccentric anomaly.

This post has been edited by **PGR**: 02 August 2018 - 05:43 PM

### #8

## Re: Plotting planet orbits by time

Posted 02 August 2018 - 08:43 PM

### #9

## Re: Plotting planet orbits by time

Posted 03 August 2018 - 07:10 AM

Which equations are you trying to solve for? It's not clear.

### #10

## Re: Plotting planet orbits by time

Posted 03 August 2018 - 10:23 AM

I cleaned up all the formulas that I could. I'm still stuck on getting the mean anomaly.

'''

'''

public OrbitControl(String name, double semiMajorAxis, double eccentricity/*degrees*/, double inclination/*degrees*/, double longitudeOfAscendingNode, double argumentOfPerihelion, double T0){ index = count++; this.name = name; bodyName[index] = this.name; this.a = semiMajorAxis; this.e = Math.toRadians(eccentricity); this.i = Math.toRadians(inclination); this.Ω = longitudeOfAscendingNode; this.ω = argumentOfPerihelion; // b = a * Math.sqrt(1d - (e * e)); f = (a * e);//average distance (kilomters?), AUs for now P = Math.sqrt(Math.pow(a, 3));//Newton's 3rd law II = Ω + ω; // calculateOrbitElements(T0); } //Calculate; //mean anomaly (M) //eccentric anomaly (E) public void calculateOrbitElements(double T0){ // double M0 = PI2 * (T0 - T)/P; T = T0 - (M0/PI2) * P;//??????????? //M = E - e * Math.sin(E);//solve for E //E = M + e * Math.sin(E); //E = origin.set(a * Math.cos(E) - f, origin.y, b * Math.sin(E)); // L = II + M0; // //orbit table //a, e, Ω, II, L, i orbit[index][a_] = a; orbit[index][e_] = e; orbit[index][Ω_] = Ω; orbit[index][II_] = II; orbit[index][L_] = L; orbit[index][i_] = i; System.out.println("OrbitControl.getTimeTranslation()" + "(" + name + ") " + "T=" + T + ",P=" + P + ",M=" + M + ",E=" + E); calculateOrbitElements2(); } public void updateOrbitElements(){ //set new parameters due to thrust //i.e., set new a, b, e //calculateOrbitElements2(); } //calculate the 3D elements public void calculateOrbitElements2(){ // //vectors if (UP == Y){ α.set(Math.cos(Ω), 0d, Math.sin(Ω)); α.normalize(); u.set(1d, 0d, 0d);//left w.set(0d, 1d, 0d);//up System.out.println("OrbitControl.calculateOrbitElements2() " + "angle=" + Math.toDegrees(u.angle(w))); v.set(u.x * w.x, u.y * w.y, u.z * w.z); α.cross(V1, VStar); v.cross(α, V1); } else if (UP == Z) {//up = 2 u.set(1d, 0d, 0d); w.set(0d, 0d, 1d); α.set(Math.cos(Ω), Math.sin(Ω), 0d); α.normalize(); v.set(u.x * w.x, u.y * w.y, u.z * w.z); α.cross(V1, VStar); v.cross(α, V1); } else{ System.out.println("bad UP code (only 1 (Y axis) or 2 (Z axis))"); System.exit(1); } System.out.println("OribtControl.calculateElements2() " + "(" + name + ") " + "α=" + α + ",v=" + v); // }

### #11

## Re: Plotting planet orbits by time

Posted 03 August 2018 - 11:19 AM

Which equation is 'mean anomaly' in all of that?

### #12

## Re: Plotting planet orbits by time

Posted 03 August 2018 - 10:06 PM

//solve for E

M = E - e * Math.sin(E);

M = E - e * Math.sin(E);

### #13

## Re: Plotting planet orbits by time

Posted 04 August 2018 - 09:47 AM

Assuming E is eccentric anomaly, perhaps use an equation built for that with information you may have.

https://en.wikipedia...centric_anomaly

http://www.jgiesen.d...ler/kepler.html

https://en.wikipedia...centric_anomaly

http://www.jgiesen.d...ler/kepler.html

Page 1 of 1