Elliptical Orbit Kepler

ofey

Hi,

I have just finished a sketch which simulates Kepler orbits. It uses a simple iteration to calculate the distance from the planet to a center and then uses this to calculate the horizontal and vertical components of the acceleration.

I am happy with it but the path doesn't always repeat itself exactly. There is a bit of precessing!

I don't have a question, I just wanted to share, in case in was useful for anyone. There were some posts on this in the old forum, but I can't post there now.

//Written by Shane Gibney
PVector pos, vel, v0, acc;
float radius = 2.0;
PVector CM;
PVector r;
float m, Const = 1000;
void setup() {
  size(600, 600);
  background(0);
  smooth();
  frameRate(24);
  pos = new PVector(150.0, 150.0);//initial position
  v0 = new PVector(1.8, 0.2);//initial velocity
  vel = new PVector(v0.x, v0.y);//velocity
  acc = new PVector(0, 0);
  CM = new PVector(height/2, width/2);
  r = new PVector(0, 0);
}
void draw() {
  if (mousePressed == false) {//hold down mouse to see path
    background(0);
  }
  PVector r = PVector.sub(pos, CM);
  m = r.mag();
  acc.x =(Const*cos(r.heading()))/pow(m, 2);// 1 < cos(r.heading()) < -1
  acc.y =(Const*sin(r.heading()))/pow(m, 2);
  vel.x -= acc.x;
  pos.x += vel.x;
  vel.y -= acc.y;
  pos.y += vel.y; 
  stroke(#6BE1EA);
  fill(#0541B4);
  ellipse(pos.x, pos.y, radius, radius);//black the projectile itself!
  stroke(#FEFF00);
  fill(#FEFF00);
  ellipse(CM.x, CM.y, 20, 20);
  /*println("pos.x = ", pos.x);
   println("pos.y = ", pos.y);
   println("vel.x = ", vel.x);
   println("vel.y = ", vel.y);
   println("acc.x = ", acc.x);
   println("acc.y = ", acc.y);
   println("pos.heading() = ", degrees(pos.heading()));
   println("m = ", m);*/
}
// Anykey to pause and unpause animation
boolean bStop;
void keyPressed()
{
  bStop = !bStop;
  if (bStop)
    noLoop();
  else
    loop();
} 

Thanks, Shane

TechWiz777

Cool! Mind if I copy it into my computer and enjoy it as well?

:)

quark

Hi Shane have moved this to the correct category.

I suspect that the precessing is due to rounding errors introduced by the float data type, might be more accurate using doubles.

Anyway very neat implementation - well done.

ofey

Hi,

TechWiz777 - I'd be delighted if you would use it.

quark - thanks, I will try doubles.... as I happens I can't use doubles as it is, as pow() only takes floats.

If different initial values for velocity and position are used it makes the ellipse more or less accurate.

The following combination of variables for inital velocity and position does something a little unusual... I won't say what. Try it!

PVector pos, vel, v0, acc, CM, r;
float m, Const = 1000, radius = 5.0;
void setup() {
  size(600, 400);
  background(255);
  smooth();
  frameRate(24);
  pos = new PVector(170.0, 130.0);//initial position
  v0 = new PVector(1.8, 0.2);//initial velocity
  vel = new PVector(v0.x, v0.y);//velocity
  acc = new PVector(0, 0);
  CM = new PVector(width/2 + 30, height/2);
  r = new PVector(0, 0);
}
void draw() {
  if (mousePressed == false) {//hold down mouse to see path
    background(255);
  }
  PVector r = PVector.sub(pos, CM);
  m = r.mag();
  acc.x =(Const*cos(r.heading()))/pow(m, 2);// 1 < cos(r.heading()) < -1
  acc.y =(Const*sin(r.heading()))/pow(m, 2);
  vel.x -= acc.x;
  pos.x += vel.x;
  vel.y -= acc.y;
  pos.y += vel.y; 
  stroke(0);
  //line(CM.x, CM.y, pos.x, pos.y);
  stroke(0);
  fill(#0006F7);//blue
  ellipse(pos.x, pos.y, radius, radius);//black the projectile itself!
  stroke(#FFFBB2);//yellow
  fill(#FEFF00);//yellow
  ellipse(CM.x, CM.y, 25, 25);
  /*println("pos.x = ", pos.x);
   println("pos.y = ", pos.y);
   println("vel.x = ", vel.x);
   println("vel.y = ", vel.y);
   println("acc.x = ", acc.x);
   println("acc.y = ", acc.y);
   println("pos.heading() = ", degrees(pos.heading()));
   println("m = ", m);*/
}
// Anykey to pause and unpause animation
boolean bStop;
void keyPressed()
{
  bStop = !bStop;
  if (bStop)
    noLoop();
  else
    loop();
} 

Next I want it to sweep out two areas, demonstrating that equal areas are swept out in equal time. Kepler's second law.

Thanks, Shane

_vk

That's the problem with doubles, Processing functions won't like it... But for math stuff you can use java's Math class:

edited with GoToLoop's obs:

double d = 1.11111111d;
double e = 1.00d;
while(e < 10.00d){
println(Math.pow(d, e));
e++;
}

GoToLoop

As mentioned, Processing's math functions are nothing more than wrappers for Java's Math class:
http://docs.oracle.com/javase/8/docs/api/java/lang/Math.html

Beware that Processing's pre-processor slaps in the f suffix in every literal that contains either a dot . or use the scientific notation:

• .5 becomes 0.5f
• 10. becomes 10.0f
• 1e3 becomes 1000.0f

In order to counter that interference, we gotta put our own d suffix on them -> .5d, 10.d, 1e3d! >-)

Another problem is that Processing's PVector uses 3 float fields -> x, y, z.
Java AWT's Point class can be a replacement for a double 2D vector though:
docs.oracle.com/javase/8/docs/api/java/awt/Point.html

P.S.: @_vk, take notice that classes which belong to package "java.lang" got auto import already! ;;)
What's common usage in the Java world for Math is: import static java.lang.Math.*; <):)

_vk

@GoToLoop thanks for the clarifications. :) I edited the code above, is it ok now?

hydrodog

General comment: If you are trying to code numerically, it is worth your while to do the whole thing in double, and only convert to float when you draw at the end. You start with too few digits accuracy in float, and can often lose the entire computation.

For squaring a number, x*x is about 50 times faster than Math.pow(x,2) as well as being more accurate, so avoid pow for integer powers.

In this case, you are using cosine and sine so you are probably ok. Incidentally, Kepler's law is misnamed. It's not really a law at all, more of an empirical formula. The real underlying law is Newton: F=ma. So if you want to make a full gravity simulator (which I am assigning my students right now by coincidence) you would use:

F = Gm1m1/r^2

Since acceleration = F/m, this means that the acceleration on any body due to another body of mass m is:

a = Gm/r^2

where r is the distance between them.

I set up a solver using just high school math. It's not the most numerically stable, but it works. You know a (though you need it in the direction of each body, a vector).

then:

x = x + (0.5 * ax * dt + vx)*dt; 
y = y + (0.5 * ay * dt + vy)*dt; 
vx = vx + ax * dt
vy = vy + ay * dt

you can add a z component if you want. The above works as long as you pick a relatively small dt (3600 second= 1 hour) or smaller.

ofey

Thanks you all for the comments, doubles and Java Maths class and the like is a scary path for me to go down. But i will take the plunge when I get a chance.

ofey

If you are trying to code numerically, it is worth your while to do the whole thing in double, and only convert to float when you draw at the end.

I agree.

For squaring a number, x*x is about 50 times faster than Math.pow(x,2) as well as being more accurate, so avoid pow for integer powers.

Ok i didn't know that. So I suppose then the pow() function is really only for high powers where it would be inconvenient to write them all out.

Kepler's law is misnamed. It's not really a law at all, more of an empirical formula.

True.

The real underlying law is Newton: F=ma. So if you want to make a full gravity simulator (which I am assigning my students right now by coincidence) you would use:

F = Gm1m1/r^2

Since acceleration = F/m, this means that the acceleration on any body due to another body of mass m is:

a = Gm/r^2

where r is the distance between them.

That is what I did actually. I was thinking of using the Kepler equation and trying to find solutions so that I could plot an (x,y) coordinate. But it can not be solved analytically. But then I thought well actually the particle is just accelerating toward a point, and only needs to be given an initial position and velocity. I didn't assign a value to the mass being orbited or a value for the universal gravitational constant. This is all taken care of by the constant called 'const' in my sketch. I had to do a little calculation to find out what this constant should be for this program. Previously I had written a sketch demonstrating projectile motion, with a vertical acceleration of 0.045. So this sketch for an elliptic orbit was similar except that now the acceleration is not uniformly down, but always pointing toward a constant point, the focus of the ellipse. (Well it is not always an ellipse, depends on the eccentricity, which depends on the initial position from the focus and initial velocity.) Of course in this new sketch the acceleration had components horizontal and vertical. The acceleration itself is an inverse-square function of the distance between the particle at any given time and the focus, where I've draw a Sun. Of course the Sun shouldn't really be at the focus. The center-of-mass of the system should be there. I decided that a const = 1000 would best normalise (possible not the correct mathematical word) the acceleration between the most distant point in the window, i.e. the corners and where the acceleration is greatest at the 'Sun', where it is actually infinitely large.

I set up a solver using just high school math. It's not the most numerically stable, but it works. You know a (though you need it in the direction of each body, a vector).

then: x = x + (0.5 * ax * dt + vx)dt; y = y + (0.5 * ay * dt + vy)dt; vx = vx + ax * dt vy = vy + ay * dt

This is interesting... I see that you are using the equations of motion for x and y directions. What does dt actually mean in a processing sketch? I know it is a time step size, but how does it actually relate to real time, or does it? It seems to me that a stepsize in time is just the time it takes the program to step through one loop of draw().

Besides sweeping out equal areas in equal times, as mentioned above, it would be really nice to take this sketch and be able to interact so that the masses could be increased or decreased. If the were equal they should orbit each other in a binary star like system. Anyway just some ideas.

Thanks, Shane

GoToLoop

An online animation using dt time. Although it's called eTime there: /:)
http://studio.processingtogether.com/sp/pad/export/ro.9jAqLpHlv-8K3/latest

ofey

Hi GoToLoop,

Thanks for the link to your sketch. It is great by the way. I see that it calculates the time for one loop of draw(). So as I guessed correctly that the time for one loop is a step size dt. But you've actually measured it. So now all your velocities are in actual units of pixels/second and so accelerations can be in pixels/second squared, when the time is real time. Very useful.

Only one thing I did not understand in you sketch was the line with >> in it.

What does the double greater than symbol mean?

Thanks,

Shane

GoToLoop

Only one thing I did not understand in you sketch was the line with >> in it.

That is the right bitshift operator: http://processing.org/reference/rightshift.html
In that y = height >> 1; case, it's just quickly integer-dividing height by 2. ;)

ofey

Hi,

I have been working an this sketch a bit and now it nicely sweeps out a constant area. I used Daniel Shiffman's example of a snake following mouse positions into an array and then shift them along that array. Works well but there is a strange white line that appears across the area when the particle comes closest to the so called Sun. Does anyone have any idea how to correct this?

PVector pos, vel, v0, acc, CM, r;
float m, Const = 1000, radius = 5.0;
int num = 30;
float[] xpos = new float[num]; 
float[] ypos = new float[num];
int x = 20, y = 10, w = 80, h = 20;
int onColour = #D4D6D8;
int offColour = #EBEEF2;
boolean reset = false;
boolean path = false;
boolean pause = false;
void setup() {
  size(620, 400);
  background(255);
  smooth();
  textSize(15);
  pos = new PVector(140.0, 130.0);//initial position B
  v0 = new PVector(1.8, -0.5);//normal with pos B near perfect ellipse
  vel = new PVector(v0.x, v0.y);//velocity
  acc = new PVector(0, 0);
  CM = new PVector(width/2 + 30, height/2);
  r = new PVector(0, 0);
  for (int i = 0; i < xpos.length; i ++ ) {
    //xpos[i] = 0; 
    //ypos[i] = 0;
    xpos[i] = pos.x; 
    ypos[i] = pos.y;
  }
}
void draw() {
  buttons();
  for (int i = 0; i < xpos.length-1; i ++ ) {
    xpos[i] = xpos[i+1];
    ypos[i] = ypos[i+1];
  }
  stroke(5);
  fill(200);
  PVector r = PVector.sub(pos, CM);
  m = r.mag();
  acc.x =(Const*cos(r.heading()))/(m*m);// 1 < cos(r.heading()) < -1
  acc.y =(Const*sin(r.heading()))/(m*m);
  vel.x -= acc.x;
  pos.x += vel.x;
  vel.y -= acc.y;
  pos.y += vel.y; 
  xpos[xpos.length-1] = pos.x;
  ypos[ypos.length-1] = pos.y;
  fill(#A7D8E5);
  beginShape();
  curveVertex(xpos[0], ypos[0]);
  for (int i = 0; i < xpos.length -1; i++) {
    curveVertex(xpos[i], ypos[i]);
  }
  curveVertex(xpos[num-1], ypos[num-1]);
  endShape();
  noStroke();
  beginShape();
  vertex(xpos[0], ypos[0]);
  vertex(xpos[num-1], ypos[num-1]);
  vertex(CM.x, CM.y);
  endShape();
  stroke(0);
  fill(#0006F7);//blue
  ellipse(pos.x, pos.y, radius, radius);//black the projectile itself!
  stroke(#FFFBB2);//yellow
  fill(#FEFF00);//yellow
  ellipse(CM.x, CM.y, 25, 25);
  /*println("pos.x = ", pos.x);
   println("pos.y = ", pos.y);
   println("vel.x = ", vel.x);
   println("vel.y = ", vel.y);
   println("acc.x = ", acc.x);
   println("acc.y = ", acc.y);
   println("pos.heading() = ", degrees(pos.heading()));*/
  println("m = ", m);
}
void buttons() {
  if (!path) {//hold down mouse to see path
    background(255);
  }
  if (path) {
    fill(onColour);
  } else {
    fill(offColour);
  }
  rect(20, y, w, h);
  fill(0);
  text("Show path", 20 + 3, y + 15);
  if (reset) {
    fill(onColour);
  } else {
    fill(offColour);
  }
  rect(120, y, w, h);
  fill(0);
  text("Reset", 120 + 3, y + 15);
  if (pause) {
    fill(onColour);
  } else {
    fill(offColour);
  }
  rect(220, y, w, h);
  fill(0);
  text("Pause", 220 + 3, y + 15);
  if (pause) {
    noLoop();
  } else {
    loop();
  }
  if (!pause) {
    fill(offColour);
  }
  rect(220, y, w, h);
  fill(0);
  text("Pause", 220 + 3, y + 15);
  fill(255);
}
void mousePressed() {
  if (mouseX > 20 && mouseX < 20+w && mouseY > y && mouseY < y+h) {
    path = !path;
  }
  if (mouseX > 120 && mouseX < 120+w && mouseY > y && mouseY < y+h) {
    reset = !reset;
    setup();
  }
  if (mouseX > 220 && mouseX < 220+w && mouseY > y && mouseY < y+h) {
    fill(offColour);
    rect(220, y, w, h);
    fill(0);
    text("Pause", 220 + 3, y + 15);
    pause = !pause;
    if (pause)
      noLoop();
    else
      loop();
  }
}

Thanks,

Shane

ofey

Hi,

I made a few changes and the area is now swept out correctly,

PVector pos, vel, v0, acc, CM, r;
float m, Const = 1000, radius = 5.0;
int num = 30;
float[] xpos = new float[num]; 
float[] ypos = new float[num];
int x = 20, y = 10, w = 80, h = 20;
int onColour = #D4D6D8;
int offColour = #EBEEF2;
boolean reset = false;
boolean path = false;
boolean pause = false;
void setup() {
  size(620, 400);
  background(255);
  smooth();
  textSize(15);
  pos = new PVector(140.0, 130.0);//initial position B
  v0 = new PVector(1.8, -0.5);//normal with pos B near perfect ellipse
  vel = new PVector(v0.x, v0.y);//velocity
  acc = new PVector(0, 0);
  CM = new PVector(width/2 + 30, height/2);
  r = new PVector(0, 0);
  for (int i = 0; i < xpos.length; i ++ ) {
    xpos[i] = pos.x; 
    ypos[i] = pos.y;
  }
}
void draw() {
  buttons();
  for (int i = 0; i < xpos.length-1; i ++ ) {
    xpos[i] = xpos[i+1];
    ypos[i] = ypos[i+1];
  }
  stroke(5);
  fill(200);
  PVector r = PVector.sub(pos, CM);
  m = r.mag();
  acc.x =(Const*cos(r.heading()))/(m*m);// 1 < cos(r.heading()) < -1
  acc.y =(Const*sin(r.heading()))/(m*m);
  vel.x -= acc.x;
  pos.x += vel.x;
  vel.y -= acc.y;
  pos.y += vel.y; 
  xpos[xpos.length-1] = pos.x;
  ypos[ypos.length-1] = pos.y;
  fill(#A7D8E5);
  noStroke();
  beginShape();
  for (int i = 0; i < xpos.length -1; i++) {
    vertex(xpos[i], ypos[i]);
    vertex(xpos[i+1], ypos[i+1]);
    vertex(CM.x, CM.y);
  }
  endShape();
  stroke(0);
  fill(#0006F7);//blue
  ellipse(xpos[num-1], ypos[num-1], radius, radius);//black the projectile itself!
  stroke(#FFFBB2);//yellow
  fill(#FEFF00);//yellow
  ellipse(CM.x, CM.y, 25, 25);
  /*println("pos.x = ", pos.x);
   println("pos.y = ", pos.y);
   println("vel.x = ", vel.x);
   println("vel.y = ", vel.y);
   println("acc.x = ", acc.x);
   println("acc.y = ", acc.y);
   println("pos.heading() = ", degrees(pos.heading()));*/
  println("m = ", m);
}
void buttons() {
  if (!path) {//hold down mouse to see path
    background(255);
  }
  if (path) {
    fill(onColour);
  } else {
    fill(offColour);
  }
  rect(20, y, w, h);
  fill(0);
  text("Show path", 20 + 3, y + 15);
  if (reset) {
    fill(onColour);
  } else {
    fill(offColour);
  }
  rect(120, y, w, h);
  fill(0);
  text("Reset", 120 + 3, y + 15);
  if (pause) {
    fill(onColour);
  } else {
    fill(offColour);
  }
  rect(220, y, w, h);
  fill(0);
  text("Pause", 220 + 3, y + 15);
  if (pause) {
    noLoop();
  } else {
    loop();
  }
  if (!pause) {
    fill(offColour);
  }
  rect(220, y, w, h);
  fill(0);
  text("Pause", 220 + 3, y + 15);
  fill(255);
}
void mousePressed() {
  if (mouseX > 20 && mouseX < 20+w && mouseY > y && mouseY < y+h) {
    path = !path;
  }
  if (mouseX > 120 && mouseX < 120+w && mouseY > y && mouseY < y+h) {
    reset = !reset;
    setup();
  }
  if (mouseX > 220 && mouseX < 220+w && mouseY > y && mouseY < y+h) {
    fill(offColour);
    rect(220, y, w, h);
    fill(0);
    text("Pause", 220 + 3, y + 15);
    pause = !pause;
    if (pause)
      noLoop();
    else
      loop();
  }
}

Thanks,

Shane

ofey

Hi,

The latest version of this code is a bit easier to use,

//Written by Shane Gibney
PVector pos, vel, v0, acc, CM, r;
float m, Const = 1000, radius = 5.0;
int num = 30;
float[] xpos = new float[num]; 
float[] ypos = new float[num];
int x = 20, y = 10, w = 80, h = 20;
boolean sweep = false;
boolean path = false;
boolean pause = false;
void setup() {
  size(620, 400);
  background(255);
  smooth();
  textSize(13);
  pos = new PVector(140.0, 130.0);//initial position 
  v0 = new PVector(1.8, -0.5);//initial velocity near perfect ellipse
  vel = new PVector(v0.x, v0.y);//velocity
  acc = new PVector(0, 0);
  CM = new PVector(width/2 + 30, height/2);
  r = new PVector(0, 0);
  for (int i = 0; i < xpos.length; i ++ ) {
    xpos[i] = pos.x; 
    ypos[i] = pos.y;
  }
}
void draw() {
  buttons();
  for (int i = 0; i < xpos.length-1; i ++ ) {
    xpos[i] = xpos[i+1];
    ypos[i] = ypos[i+1];
  }
  stroke(5);
  fill(200);
  PVector r = PVector.sub(pos, CM);
  m = r.mag();
  acc.x =(Const*cos(r.heading()))/(m*m);// 1 < cos(r.heading()) < -1
  acc.y =(Const*sin(r.heading()))/(m*m);
  vel.x -= acc.x;
  pos.x += vel.x;
  vel.y -= acc.y;
  pos.y += vel.y; 
  xpos[xpos.length-1] = pos.x;
  ypos[ypos.length-1] = pos.y;
  if (sweep) {
    fill(#A7DAFF);
    noStroke();
    beginShape();
    for (int i = 0; i < xpos.length -1; i++) {
      vertex(xpos[i], ypos[i]);
      vertex(xpos[i+1], ypos[i+1]);
      vertex(CM.x, CM.y);
    }
    endShape();
  }
  stroke(0);
  fill(#0006F7);//blue
  ellipse(xpos[num-1], ypos[num-1], radius, radius);//the particle itself!
  stroke(#FFFBB2);//yellow
  fill(#FEFF00);//yellow
  ellipse(CM.x, CM.y, 25, 25);
  /*println("pos.x = ", pos.x);
   println("pos.y = ", pos.y);
   println("vel.x = ", vel.x);
   println("vel.y = ", vel.y);
   println("acc.x = ", acc.x);
   println("acc.y = ", acc.y);
   println("pos.heading() = ", degrees(pos.heading()));
   println("m = ", m);*/
}
void buttons() {
  stroke(#999999);
  if (!path) {//hold down mouse to see path
    background(255);
  }
  if (path) {
    fill(#BBBBBB);//on
    rect(20, y, w, h);//path
    fill(255);
    text("Show path", 20 + 3, y + 15);
    fill(#DDDDDD);
  } else {
    fill(#DDDDDD);//off
    rect(20, y, w, h);//path
    fill(0, 0, 255);
    text("Show path", 20 + 3, y + 15);
    fill(#DDDDDD);
  }
  if (sweep) {
    fill(#BBBBBB);//on
    rect(120, y, w, h);//path
    fill(255);
    text("Area", 120 + 3, y + 15);
    fill(#DDDDDD);
  } else {
    fill(#DDDDDD);//off
    rect(120, y, w, h);//path
    fill(0, 0, 255);
    text("Area", 120 + 3, y + 15);
    fill(#DDDDDD);
  }
  fill(#DDDDDD);//off
  rect(220, y, w, h);//path
  fill(0, 0, 255);
  text("Pause", 220 + 3, y + 15);
  fill(#DDDDDD);
}
void mousePressed() {
  if (mouseX > 20 && mouseX < 20 + w && mouseY > y && mouseY < y + h) {
    path = !path;
  }
  if (mouseX > 120 && mouseX < 120 + w && mouseY > y && mouseY < y + h) {
    sweep = !sweep;
  }
  if (mouseX > 220 && mouseX < 220 + w && mouseY > y && mouseY < y + h) {
    if (!pause) {
      stroke(#999999);
      fill(#BBBBBB);//on
      rect(220, y, w, h);//path
      fill(255);
      text("Pause", 220 + 3, y + 15);
      fill(#DDDDDD);
    }
    pause = !pause;
    if (pause)
      noLoop();
    else
      loop();
  }
}

Thanks, Shane

ofey

This code is on github.com

https://github.com/shanegibney/KeplerOrbitProcessing