## 12 Replies - 889 Views - Last Post: 04 August 2018 - 09:47 AMRate Topic:     //<![CDATA[ rating = new ipb.rating( 'topic_rate_', { url: 'https://www.dreamincode.net/forums/index.php?app=forums&module=ajax&section=topics&do=rateTopic&t=412009&amp;s=5272e8b913dd2cbf1c877b8f07e09898&md5check=' + ipb.vars['secure_hash'], cur_rating: 0, rated: 0, allow_rate: 0, multi_rate: 1, show_rate_text: true } ); //]]>

### #1 PGR Reputation: 13
• Posts: 65
• Joined: 10-January 11

# Plotting planet orbits by time

Posted 01 August 2018 - 11:08 PM

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.
It's the OrbitControl.java calculateOrbitElements() and calculateOrbitElements2() I'm stuck on.

Is This A Good Question/Topic? 0

## Replies To: Plotting planet orbits by time

### #2 modi123_1 • • Suitor #2
•    Reputation: 15357
• Posts: 61,571
• Joined: 12-June 08

## 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?

### #3 PGR Reputation: 13
• Posts: 65
• Joined: 10-January 11

## 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.
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 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];
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 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);
}

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.Ω = 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);

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], orbit[i], orbit[i], orbit[i], orbit[i], orbit[i], orbit[i]);

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];
}
for (int j = 0;j < testData.length;j++){
for (int i = 0;i < testData.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 {
InputCapsule in = im.getCapsule(this);
//TODO: load properties of this Control, e.g.
}

@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 ndc85430 Reputation: 981
• Posts: 3,870
• Joined: 13-June 14

## 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?

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 baavgai • • Dreaming Coder
•   Reputation: 7501
• Posts: 15,544
• Joined: 16-October 07

## 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

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)

## 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?

### #7 PGR Reputation: 13
• Posts: 65
• Joined: 10-January 11

## 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.

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

### #8 PGR Reputation: 13
• Posts: 65
• Joined: 10-January 11

## Re: Plotting planet orbits by time

Posted 02 August 2018 - 08:43 PM modi123_1, on 01 August 2018 - 11:19 PM, said:

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

Language is Java.
No, it's not a paid job. It's just to test a (fixed) orbit by time.

### #9 modi123_1 • • Suitor #2
•    Reputation: 15357
• Posts: 61,571
• Joined: 12-June 08

## 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 PGR Reputation: 13
• Posts: 65
• Joined: 10-January 11

## 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.Ω = 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 modi123_1 • • Suitor #2
•    Reputation: 15357
• Posts: 61,571
• Joined: 12-June 08

## Re: Plotting planet orbits by time

Posted 03 August 2018 - 11:19 AM

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

### #12 PGR Reputation: 13
• Posts: 65
• Joined: 10-January 11

## Re: Plotting planet orbits by time

Posted 03 August 2018 - 10:06 PM

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

### #13 modi123_1 • • Suitor #2
•    Reputation: 15357
• Posts: 61,571
• Joined: 12-June 08

## 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

Page 1 of 1

 .related ul { list-style-type: circle; font-size: 12px; font-weight: bold; } .related li { margin-bottom: 5px; background-position: left 7px !important; margin-left: -35px; } .related h2 { font-size: 18px; font-weight: bold; } .related a { color: blue; }