code orbit

// planet
// understanding newton's laws

#include <cstdlib>
#include <iostream>
#include <cstdio>
#include <iomanip>
#include <cmath>

using namespace std;

#define D 2 // number of dimensions
struct particle {
double x[D] ; // (x,y) coordinates
double v[D] ; // velocity
double a[D] ; // accelleration
};

void showstate(particle &a, particle &b, double t);
double displacement_sun(particle a);
double displacement_earth(particle a, particle b);
void positionstep(particle &a, double dt);
void velocitystep(particle &a, double dt);
void accelerationa(particle &a);
void accelerationb(particle &a, particle &b);



int main()
{
particle a, b;
double r, rs, re, t=0.0, dt=0.001;

a.x[0]=1;
a.x[1]=0;
a.v[0]=0;
a.v[1]=1;
a.a[0]=0;
a.a[1]=0;

b.x[0]=1.1;
b.x[1]=0;
b.v[0]=0;
b.v[1]=1.1;
b.a[0]=0;
b.a[1]=0;

for (int i=1; i<=10000; i++)
{
	positionstep(a, dt*0.5);
	positionstep(b, dt*0.5);
	accelerationa(a);
	accelerationb(a, b);
	velocitystep(a, dt);
	velocitystep(b, dt);
	positionstep(a, dt*0.5);
	positionstep(b, dt*0.5);
	t += dt;
	showstate(a,b,t);
}

return 0;
}

void showstate(particle &a, particle &b, double t)
{
	cout << t << "\t" << a.x[0] << "\t" << a.x[1] << "\t" << a.v[0] << "\t" << a.v[1]
	<< "\t" << b.x[0] << "\t" << b.x[1] << "\t" << b.v[0] << "\t" << b.v[1] << endl;
}

double displacement_sun(particle a)
{
	double r;
	r = sqrt((a.x[0])*(a.x[0])+(a.x[1])*(a.x[1]));
	return r;
}

double displacement_earth(particle a, particle b)
{
	double r;
	r = sqrt(((b.x[0])-(a.x[0]))*((b.x[0])-(a.x[0]))+((b.x[1])-(a.x[1]))*((b.x[1])-(a.x[1])));
	return r;
}

void accelerationa(particle &a)
{
	double G=1, M=1, r;
	r = displacement_sun(a);
	a.a[0]=-G*M*(a.x[0])/(r*r*r);
	a.a[1]=-G*M*(a.x[1])/(r*r*r);
}

void accelerationb(particle &a, particle &b)
{
	int i;
	double G=1, M=1, m=0.1, rs, re;
	rs = displacement_sun(a);
	re = displacement_earth(a, b);
	for (i=0; i<=D; i++)
	{
		if (sqrt(b.x[i]*b.x[i])>sqrt(a.x[i]*a.x[i]))
		{
			b.a[i]=-G*M*(b.x[i])/(rs*rs*rs)-G*m*(b.x[i]-a.x[i])/(re*re*re);
		}
		if (sqrt(b.x[i]*b.x[i])<sqrt(a.x[i]*a.x[i]))
		{
			b.a[i]=-G*M*(b.x[i])/(rs*rs*rs)+G*m*(b.x[i]-a.x[i])/(re*re*re);
		}
	}

}

void velocitystep(particle &a, double dt)
{
	a.v[0]+=(a.a[0])*dt;
	a.v[1]+=(a.a[1])*dt;
}

void positionstep(particle &a, double dt)
{
	a.x[0]+=+(a.v[0])*dt;
	a.x[1]+=+(a.v[1])*dt;
}