We are about to switch to a new forum software. Until then we have removed the registration on this forum.

- All Categories 25.7K
- Announcements & Guidelines 13
- Common Questions 30
- Using Processing 22.1K
- Programming Questions 12.2K
- Questions about Code 6.4K
- How To... 4.2K
- Hello Processing 72
- GLSL / Shaders 292
- Library Questions 4K
- Hardware, Integration & Other Languages 2.7K
- Kinect 668
- Arduino 1K
- Raspberry PI 188
- Questions about Modes 2K
- Android Mode 1.3K
- JavaScript Mode 413
- Python Mode 205
- Questions about Tools 100
- Espanol 5
- Developing Processing 548
- Create & Announce Libraries 211
- Create & Announce Modes 19
- Create & Announce Tools 29
- Summer of Code 2018 93
- Rails Girls Summer of Code 2017 3
- Summer of Code 2017 49
- Summer of Code 2016 4
- Summer of Code 2015 40
- Summer of Code 2014 22
- p5.js 1.6K
- p5.js Programming Questions 947
- p5.js Library Questions 315
- p5.js Development Questions 31
- General 1.4K
- Events & Opportunities 288
- General Discussion 365

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

## Comments

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

:)

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.

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!

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

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:

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.*;`

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

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:

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.

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.

I agree.

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.

True.

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.

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

An online animation using

dttime. Although it's calledeTimethere: /:)http://studio.processingtogether.com/sp/pad/export/ro.9jAqLpHlv-8K3/latest

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

That is the right bitshift operator: http://processing.org/reference/rightshift.html

In that

`y = height >> 1;`

case, it's just quickly integer-dividingheightby 2. ;)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?

Thanks,

Shane

Hi,

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

Thanks,

Shane

Hi,

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

Thanks, Shane

This code is on github.com

https://github.com/shanegibney/KeplerOrbitProcessing