- #1
Silviu
- 624
- 11
Hi! I am trying to simulate the rotation of a planet around a star, using the 2nd order Runge-Kutta method (I am starting with this and I will try 4th order later, I am new to this topic). I have this code but it doesn't work. When I plot y(x) I should get a circle, but I don't get it. Any idea what am I doing wrong? Thank you!
<<Mentor's note: please use code tags when posting code.>>
Code:
#include<iostream>
#include <vector>
#include<math.h>
#include <fstream>
usingnamespace std;
class Body{
private:
double G= 1;
double rx;
double ry;
double rz;
double vx;
double vy;
double vz;
double mass;
double dt = 0.001;
double dist,dx,dy;
double x1,vx1,ax1,ax2;
double y1,vy1,ay1,ay2;
public:
Body(double rx, double ry, double rz, double vx, double vy, double vz, double mass){
this->rx=rx;
this->ry=ry;
this->rz=rz;
this->vx=vx;
this->vy=vy;
this->vz=vz;
this->mass=mass;
}
void update(Body b){
//Step 1
dx=b.rx-rx;
dy=b.ry-ry;
dist = sqrt(dx*dx+dy*dy);
ax1=G*b.mass*dx/(dist*dist*dist);
ay1=G*b.mass*dy/(dist*dist*dist);
//Step 2
x1=rx+vx*dt*0.5;
vx1=vx+ax1*dt*0.5;
y1=ry+vy*dt*0.5;
vy1=vy+ay1*dt*0.5;
dx=b.rx-x1;
dy=b.ry-y1;
dist = sqrt(dx*dx+dy*dy);
ax2=G*b.mass*dx/(dist*dist*dist);
ay2=G*b.mass*dy/(dist*dist*dist);
rx=rx+0.5*(vx1+vx)*dt;
ry=ry+0.5*(vy1+vy)*dt;
vx=vx+0.5*(ax1+ax2)*dt;
vy=vy+0.5*(ay1+ay2)*dt;
}
double get_x(){
return rx;
}
double get_y(){
return ry;
}
};
int main(){
Body body1(0,0,0,0,0,0,1000), body2(10,0,0,0,10,0,10);
ofstream pos;
pos.open ("Position.txt");
int N=100000;
for(int i; i<N;i++){
body2.update(body1);
pos<<body2.get_x()<<" "<<body2.get_y()<<endl;
}
pos.close();
}
<<Mentor's note: please use code tags when posting code.>>
Attachments
Last edited by a moderator: