two body problem

waerx

so im trying to draw two body problem in 2D, it does quite good with out the PVectors, i tried to make the code look better with using vectors and the program doesnt run good for small time steps like it did before (same initial values in both case) how could i make this program better to run faster and calculate precisely?

  float G=0.00000667408, m1=10, m2=900000000, t0=0, h;
  float ti=0;
  PVector x10, x20, x1, x2, v10, v20, v1, v2, r1e, r2e, v1t, v2t, xx1, xx2;
  float r1, r2;
void setup() {  // setup() runs once
  size(1150, 900);
  x10=new PVector(450.0,250.0);
  x20=new PVector(600.0,450.0);
  v10=new PVector(0.53,-0.15);
  v20=new PVector(-0.005,0.006);
  x1=x10.copy();
  x2=x20.copy();
  xx1=x10.copy();
  xx2=x20.copy();
  v1=v10.copy();
  v2=v20.copy();
  v1t=v1.copy();
  v2t=v2.copy();
  h=0.99;
  frameRate(500);  
}
void draw() {
  background(204); 
  ti=t0+h; //here if h is 1 its going "good" but for small values like 0.7 and below they crush
  r1e=PVector.sub(x2, x1);
  r2e=PVector.sub(x1, x2);
  r1=x1.mag();
  r2=x2.mag();
  v1=PVector.add(v1t, r1e.mult(G*m2*(ti-t0)/(r1*r1*r1)));
  v2=PVector.add(v2t, r2e.mult(G*m1*(ti-t0)/(r2*r2*r2)));
  x1=PVector.add(xx1, v1.mult(ti-t0));
  x2=PVector.add(xx2, v2.mult(ti-t0));
  xx1=x1.copy();
  xx2=x2.copy();
  v1t=v1.copy();
  v2t=v2.copy();
  t0=ti;
  fill(50);
  ellipse(x1.x, x1.y, 15, 15);
  ellipse(x2.x, x2.y, 45, 45);
}

waerx

here is the original program which not have the PVectors and works for smaller time steps:

  float order=1; //2nd order if its 1 else first order only (Euler method)
  float G=0.0000667408*1.5, m1=1000, m2=9000000, t0=0;
  float ti=0; 
  float h;
  float x10=350; //initial values
  float y10=250;

  float x20=600;
  float y20=450;
  
  float x1=x10; //variables
  float y1=y10;
  
  float x2=x20;
  float y2=y20;

  float vx10=1.3; //initial values
  float vy10=-0.6;

  float vx20=-0.05;
  float vy20=-0.06;
  
  float r1x=x20-x10; //starting directions
  float r1y=y20-y10;

  float r2x=x10-x20;
  float r2y=y10-y20;
  
  float v2tx=vx20; //i-1 component storage for the variable_(i) - variable_(i-1) calculations
  float v2ty=vy20;
  
  float v1tx=vx10;
  float v1ty=vy10;

  float v1ttx=vx10;
  float v1tty=vy10;
 
  float v2ttx=vx20;
  float v2tty=vy20;

  float xx1=x10;
  float yy1=y10;

  float xx2=x20;
  float yy2=y20;
  
  float r1=0; //length of the vectors will be stored in here
  float r2=0;
  //G=0.0000000000667408; //original G value
void setup() {  // setup() runs once
  size(1150, 900);
  frameRate(35000);  
}
 
void draw() {  // draw() loops forever, until stopped
  background(204);
  stroke(100);

  ti=t0+0.1; //delta t time step, smaller means better calculation but takes longer
  
  r1x=x2-x1; //directions calculated every step
  r1y=y2-y1;

  r2x=x1-x2;
  r2y=y1-y2;
  
  r1=sqrt(r1x*r1x+r1y*r1y); //length calculated every step
  r2=sqrt(r2x*r2x+r2y*r2y);
  
  
  //euler 2nd order method: y(n+1)=yn+hy'n+h^2/2y''n+O(h^3) where O(h^3) is the error, h=ti-t0
  
  v1tx=G*m2/(r1*r1*r1)*r1x*(ti-t0)+v1ttx-order*(3*G*m2/(r1*r1*r1*r1)*r1x*(ti-t0)*(ti-t0)+(ti-t0)*(ti-t0)*G*m2*(x2-xx2-xx1+x1)/(r1*r1*r1)); //new velocity
  v1ty=G*m2/(r1*r1*r1)*r1y*(ti-t0)+v1tty-order*(3*G*m2/(r1*r1*r1*r1)*r1y*(ti-t0)*(ti-t0)+(ti-t0)*(ti-t0)*G*m2*(y2-yy2-yy1+y1)/(r1*r1*r1));

  v2tx=G*m1/(r2*r2*r2)*r2x*(ti-t0)+v2ttx-order*(3*G*m1/(r2*r2*r2*r2)*r2x*(ti-t0)*(ti-t0)+(ti-t0)*(ti-t0)*G*m1*(x1-xx1-xx2+x2)/(r2*r2*r2));
  v2ty=G*m1/(r2*r2*r2)*r2y*(ti-t0)+v2tty-order*(3*G*m1/(r2*r2*r2*r2)*r2y*(ti-t0)*(ti-t0)+(ti-t0)*(ti-t0)*G*m1*(x1-xx1-xx2+x2)/(r2*r2*r2));
   
  x1=v1tx*(ti-t0)+xx1+order*(ti-t0)*(ti-t0)*G*m2/(r1*r1*r1)*r1x;
  y1=v1ty*(ti-t0)+yy1+order*(ti-t0)*(ti-t0)*G*m2/(r1*r1*r1)*r1y;

  x2=v2tx*(ti-t0)+xx2+order*(ti-t0)*(ti-t0)*G*m1/(r2*r2*r2)*r2x;
  y2=v2ty*(ti-t0)+yy2+order*(ti-t0)*(ti-t0)*G*m1/(r2*r2*r2)*r2y;
  
  v1ttx=v1tx;
  v1tty=v1ty;

  v2ttx=v2tx;
  v2tty=v2ty;
  xx1=x1;
  yy1=y1;

  xx2=x2;
  yy2=y2;
  
  t0=ti;
  ellipse(x1, y1, 20, 20);
  ellipse(x2, y2, 45, 45);
  fill(50);
  text("m1 mass position:",10,15); 
  text("m1 mass:",10,75);
  text(m1,10,85);
  text(x1,10,25); 
  text(y1,10,35); 
  text("m1 mass velocity (x10000):",10,45);
  text(v1tx*10000,10,55); 
  text(v1ty*10000,10,65); 
  
  text("m2 mass position:",500,15); 
  text("m2 mass:",500,75);
  text(m2,500,85);
  text(x2,500,25); 
  text(y2,500,35); 
  text("m2 mass velocity(x10000):",500,45);
  text(v2tx*10000,500,55); 
  text(v2ty*10000,500,65); 
}

koogs

v2tx=G*m1/(r2*r2*r2)*r2x*(ti-t0)+v2ttx-order*(3*G*m1/(r2*r2*r2*r2)*r2x*(ti-t0)*(ti-t0)+(ti-t0)*(ti-t0)*G*m1*(x1-xx1-xx2+x2)/(r2*r2*r2));

does your space bar work?

TfGuy44

Pleasedon'tjokeaboutbrokenspacebars.:-(

koogs

Seriously though, if you want people to debug your code then make it easy for them to read. An unbroken stream of characters and symbols that long, coupled with all your variables being combinations of the same 3 letters, does not help.

koogs

where is the equiv of lines 73-77 in the first bit of code? i don't see it.

waerx

In the upper code i have the problem, where it is first order, so if order=0 in the second they are the same. It means:

v1tx=G m2/(r1 r1 r1) r1x (ti-t0)+v1ttx v1ty=G m2/(r1 r1 r1) r1y (ti-t0)+v1tty

these are the new velocities which we use for the orbital calculation, and my problem is that in the second code i can use small time steps like 0.2 to calculate more precise orbits, and in the first code i need to use times step 1.0 (so h=1 and it works but that means big mistake in the calculation i think)

koogs

(use blank lines between paragraphs or it'll run them all together. i have fixed the above formatting)

no idea what yo are saying. pretend i know nothing about the two body problem.

waerx

In the original code the 73-77 line is (if we use first order, that means y(n+1)=yn+hy'n+O(h^2) ):

v1tx=G m2/(r1 r1 r1) r1x (ti-t0)+v1tt;
v1ty=G m2/(r1 r1 r1) r1y (ti-t0)+v1tty;

v2tx=G*m1/(r2*r2*r2)*r2x*(ti-t0)+v2ttx;
v2ty=G*m1/(r2*r2*r2)*r2y*(ti-t0)+v2tty;

x1=v1tx*(ti-t0)+xx1;
y1=v1ty*(ti-t0)+yy1;
 
x2=v2tx*(ti-t0)+xx2;
y2=v2ty*(ti-t0)+yy2;
   

In the PVector code its

v1=PVector.add(v1t, r1e.mult(G*m2*(ti-t0)/(r1*r1*r1)));
v2=PVector.add(v2t, r2e.mult(G*m1*(ti-t0)/(r2*r2*r2)));
x1=PVector.add(xx1, v1.mult(ti-t0));
x2=PVector.add(xx2, v2.mult(ti-t0));

here i calculate the velocities and orbits in a 2D vectors, r1e, r2e means the directions. In xx1 and xx2 I store the last (i-1) value to the calculate the i-th value. Id like to use this second compact vector form because i think its easier to read it, but i have problems with the code, its only running if the ti-t0 time step is equal 1.

GoToLoop

https://OpenProcessing.org/sketch/437853

/** 
 * Euler 2 Body Masses (v1.1)
 * Mod: GoToLoop (2017/Jul/04)
 *
 * Forum.Processing.org/two/discussion/23308/
 * two-body-problem#Item_9
 *
 * OpenProcessing.org/sketch/437853
 */

static final boolean JAVA = 1/2 != 1/2.;
static final int FONT_SIZE = 18;

static final String FONT_JAVA = "DejaVu Sans Mono";
static final String FONT = JAVA? FONT_JAVA : "monospace";
static final String LF = "\n", COMMA = ", ", M = "M";

static final String FC = "Frame: ", FR = "  -  FPS: ";
static final String MASS = " Mass: ";
static final String POS  = "Position: ";
static final String VEL  = "Velocity (x10000): ";

static final color BG = 0350, FG = 0100;
static final color C1 = #0000FF, C2 = #FF0000;
static final float BOLD = 3.5;

static final int W = 900, H = 650, FPS = 60, SMOOTH = 3;
static final int DIAM1 = 20, DIAM2 = 45, OFFSET = 50;

static final float G = 6.67408e-5; // 6.67408e-11
static final float T = .5, TT = T * T, SPD = 1e4;

static final float M1 = 1e3, GM1 = G*M1, TGM1 = T*GM1;
static final float TTGM1 = TT*GM1, TTG3M1 = 3*TTGM1;
static final float X1 = W>>2, Y1 = H>>2;
static final float VX1 = 1.3, VY1 = -.6;

static final float M2 = 9e6, GM2 = G*M2, TGM2 = T*GM2;
static final float TTGM2 = TT*GM2, TTG3M2 = 3*TTGM2;
static final float X2 = W>>1, Y2 = H>>1;
static final float VX2 = 0, VY2 = 0;

float x1, y1, x2, y2;
float vx1, vy1, vx2, vy2;

void settings() {
  size(W, H, JAVA? FX2D : JAVA2D);
  smooth(SMOOTH);
}

void setup() {
  if (width != W)  settings();
  frameRate(FPS);

  colorMode(RGB);
  ellipseMode(CENTER);
  strokeWeight(BOLD);

  textMode(MODEL);
  textAlign(LEFT, BASELINE);
  textFont(createFont(FONT, FONT_SIZE, true));

  reset();
}

void draw() {
  final float rx1 = x2 - x1, ry1 = y2 - y1;
  final float rx2 = x1 - x2, ry2 = y1 - y2;

  final float r1 = mag(rx1, ry1), r2 = mag(rx2, ry2);

  final float r13 = r1*r1*r1, r14 = r13*r1;
  final float r23 = r2*r2*r2, r24 = r23*r2;

  final float rx1r13 = rx1/r13, ry1r13 = ry1/r13;
  final float rx2r23 = rx2/r23, ry2r23 = ry2/r23;

  vx1 += TGM2*rx1r13 - TTG3M2*rx1/r14;
  vy1 += TGM2*ry1r13 - TTG3M2*ry1/r14;

  vx2 += TGM1*rx2r23 - TTG3M1*rx2/r24;
  vy2 += TGM1*ry2r23 - TTG3M1*ry2/r24;

  x1 += T*vx1 + TTGM2*rx1r13;
  y1 += T*vy1 + TTGM2*ry1r13;

  x2 += T*vx2 + TTGM1*rx2r23;
  y2 += T*vy2 + TTGM1*ry2r23;

  if (JAVA)  title();

  background(BG);
  fill(FG);

  stroke(C1);
  ellipse(x1, y1, DIAM1, DIAM1);

  stroke(C2);
  ellipse(x2, y2, DIAM2, DIAM2);

  fill(C1);
  hud(1, OFFSET, OFFSET, M1, x1, y1, vx1, vy1);

  fill(C2);
  hud(2, width/2 + OFFSET, OFFSET, M2, x2, y2, vx2, vy2);
}

void mousePressed() {
  if (mouseButton != LEFT)  reset();
}

void keyPressed() {
  if (key == ' ' | key == ENTER)  reset();
}

void reset() {
  x1 = X1;
  y1 = Y1;
  x2 = X2;
  y2 = Y2;

  vx1 = VX1;
  vy1 = VY1;
  vx2 = VX2;
  vy2 = VY2;
}

void title() {
  surface.setTitle(FC + frameCount + FR + round(frameRate));
}

void hud(
  int idx, int x, int y, float m, 
  float px, float py, float vx, float vy) {

  final String info =
    M + idx + MASS + round(m) + LF +
    POS + round(px) + COMMA + round(py) + LF +
    VEL + round(vx*SPD) + COMMA + round(vy*SPD);

  text(info, x, y);
}

kfrajer

I think you have to do this replacement in the PVector code at lines 27 and 28:

  r1=r1e.mag();  //x1
  r2=r2e.mag();  //x2

However, this does not solve your problem. Few suggestions. When you present two codes, ensure they have the initial set of parameters. This will debugging easier. Then, I'd test limits in your functions to see if both behave the same. For example, making your static planet massive. Do both codes behave the same? In your case, probably not. It could be a simple bug. It is not easy to find since you are working directly on the equations.

The solution to your initial problem relies on this difference:

x1=PVector.add(xx1, v1.mult(ti-t0));

vs.

x1=PVector.add(xx1, PVector.mult(v1,(ti-t0)));

Kf

waerx

Yes it solved it, thank you. Do you have any advise how to make it calculate faster? If anyone intrested about the working code here it is:

float G=0.00000667408, m1=10, m2=900000000, t0=0.0, h;
  float ti=0.0;
  PVector x10, x20, x1, x2, v10, v20, v1, v2, r1e, r2e, v1t, v2t, xx1, xx2;
  float r1, r2;
void setup() {
  size(1150, 900);
  x10=new PVector(475.0,200.0);
  x20=new PVector(800.0,700.0);
  v10=new PVector(1.6,-0.3);
  v20=new PVector(-0.000,0.000);
  x1=x10.copy();
  x2=x20.copy();
  xx1=x10.copy();
  xx2=x20.copy();
  v1=v10.copy();
  v2=v20.copy();
  v1t=v1.copy();
  v2t=v2.copy();
  h=0.1;
  frameRate(500);  
}
void draw() {
    background(250);
  ti=t0+h;
  r1e=PVector.sub(x2, x1);
  r2e=PVector.sub(x1, x2);
  r1=r1e.mag();
  r2=r2e.mag();
  v1=PVector.add(v1t, r1e.mult(G*m2*(ti-t0)/(r1*r1*r1)));
  v2=PVector.add(v2t, r2e.mult(G*m1*(ti-t0)/(r2*r2*r2)));
  x1=PVector.add(xx1, PVector.mult(v1,(ti-t0)));
  x2=PVector.add(xx2, PVector.mult(v2,(ti-t0)));
  xx1=x1.copy();
  xx2=x2.copy();
  v1t=v1.copy();
  v2t=v2.copy();
  t0=ti;
  fill(55);
  ellipse(x1.x/2, x1.y/2, 10, 10);
  ellipse(x2.x/2, x2.y/2, 45, 45);
  text("G grav. value: (x10000):",500,45);
  text(G*10000,500,55); 
}

koogs

how to make it calculate faster

it's doing 170 fps here.

prince_polka

class body {
  PVector position;
  PVector speed;
  float mass;
  body(float XPOS,float YPOS,float XSPEED,float YSPEED,float MASS) {
    position = new PVector(XPOS,YPOS);
    speed = new PVector(XSPEED,YSPEED);
    mass=MASS;
  }
  void fall(float pointMass) {
    speed.add(PVector.mult(delta,pointMass*force));
    position.add(speed);
    ellipse(position.x, position.y, 40, 40);
  }
}
body A;
body B;
PVector delta;
float force;
void setup() {
  size(500,500); 
  A = new body(250,250,-.01,0,50);
  B = new body(240,400,0.5,0,-5);
  delta = new PVector(0, 0);
  delta = PVector.sub(A.position, B.position);
}
void draw() {
  background(0);
  fill(180);
  delta = PVector.sub(A.position, B.position) ; 
  force = 1./delta.magSq()/delta.mag();
  A.fall(B.mass);
  B.fall(A.mass);
}