12 Replies - 734 Views - Last Post: 04 August 2018 - 09:47 AM Rate Topic: -----

#1 PGR   User is offline

  • D.I.C Head

Reputation: 13
  • View blog
  • 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.
https://github.com/Mipada/OrbitTest
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   User is offline

  • Suitor #2
  • member icon



Reputation: 14929
  • View blog
  • Posts: 59,607
  • 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?
Was This Post Helpful? 0
  • +
  • -

#3 PGR   User is offline

  • D.I.C Head

Reputation: 13
  • View blog
  • 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.
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

Was This Post Helpful? 0
  • +
  • -

#4 ndc85430   User is offline

  • I think you'll find it's "Dr"
  • member icon

Reputation: 976
  • View blog
  • Posts: 3,849
  • Joined: 13-June 14

Re: Plotting planet orbits by time

Posted 02 August 2018 - 12:27 AM

View PostPGR, 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".
Was This Post Helpful? 0
  • +
  • -

#5 baavgai   User is offline

  • Dreaming Coder
  • member icon


Reputation: 7420
  • View blog
  • Posts: 15,377
  • Joined: 16-October 07

Re: Plotting planet orbits by time

Posted 02 August 2018 - 09:07 AM

View PostPGR, 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)

Was This Post Helpful? 0
  • +
  • -

#6 victorman   User is offline

  • New D.I.C Head

Reputation: 2
  • View blog
  • Posts: 19
  • Joined: 03-July 18

Re: Plotting planet orbits by time

Posted 02 August 2018 - 02:53 PM

View Postndc85430, 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/
Was This Post Helpful? 0
  • +
  • -

#7 PGR   User is offline

  • D.I.C Head

Reputation: 13
  • View blog
  • 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

Was This Post Helpful? 0
  • +
  • -

#8 PGR   User is offline

  • D.I.C Head

Reputation: 13
  • View blog
  • Posts: 65
  • Joined: 10-January 11

Re: Plotting planet orbits by time

Posted 02 August 2018 - 08:43 PM

View Postmodi123_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.
Was This Post Helpful? 0
  • +
  • -

#9 modi123_1   User is offline

  • Suitor #2
  • member icon



Reputation: 14929
  • View blog
  • Posts: 59,607
  • 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.
Was This Post Helpful? 0
  • +
  • -

#10 PGR   User is offline

  • D.I.C Head

Reputation: 13
  • View blog
  • 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.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);
        //
        
    }


Was This Post Helpful? 0
  • +
  • -

#11 modi123_1   User is offline

  • Suitor #2
  • member icon



Reputation: 14929
  • View blog
  • Posts: 59,607
  • 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?
Was This Post Helpful? 0
  • +
  • -

#12 PGR   User is offline

  • D.I.C Head

Reputation: 13
  • View blog
  • 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);
Was This Post Helpful? 0
  • +
  • -

#13 modi123_1   User is offline

  • Suitor #2
  • member icon



Reputation: 14929
  • View blog
  • Posts: 59,607
  • 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
Was This Post Helpful? 0
  • +
  • -

Page 1 of 1