#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
/*A program designed to plot the trajectory of a rocket on a mission to land on the moon*/
 
int main(int argc, char *argv[])
{
const float pii=3.1415; /*This defines the value of pi*/
const float u=0.0066; /*This defines the intial velocity of rocket and is equivalent to v0*/
const float angle=((89.9*pii)/180.0); /*This defines the angle of projection*/
const float xe=0, ye=0; /*These are the Earth's coords: (0, 0) - ie. where it's centre lies*/
const float xm=0, ym=225.7; /*These are the Moon's coords: (0, 225.7)*/
const float Mm=1.0, Me=83.3; /*These are the masses of the Moon & Earth Respectively*/
const float Rm=1.0, Re=3.7; /*These are the radii of the Moon & Earth Respectively*/
const float G=(9.81*pow(10, -7)); /*This is the Gravitational constant)*/
 
int dt=10; /*This is the step size delta.t*/
float xr=0, yr=3.7; /*These are the coords of the Rocket)*/
float vx, vy; /*These are the current velocities in both x & y directions*/
float ax=0, ay=0; /*These are the accelerations in both x & y directions*/
float rrx, rry; /*These are the components of the vector from Earth to Rocket*/
float rmx, rmy; /*These are the components of the vector from earth to Moon*/
/*Do NOT confuse with Re & Rm which are radii*/
 
vx = u*cos(angle);
vy = u*sin(angle);
 
int t=0;
/*Using a while loop*/
while(pow((pow(xm-xr,2)) + (pow(ym-yr,2)),1/2)>=Rm)
/*When the distance to the centre of the moon < radius of the moon,*/
/*then the rocket has landed and the calculation can be stopped*/
{
t=t+dt; /*Time increases by 10 each time*/
rmx = xm - xe;
rmy = ym - ye;
rrx = xr - xe;
rry = yr - ye;
ax = (((-G*Me)*xr)/(pow(((rrx*rrx)+(rry*rry)),(3/2)))) - ((G*Mm*(xr-xm))/(pow((pow(((rrx-rmx)),2)+(pow(((rry-rmy)),2))),(3/2)))); /*Acceleration in x-direction*/
ay = (((-G*Me)*yr)/(pow(((rrx*rrx)+(rry*rry)),(3/2)))) - ((G*Mm*(yr-ym))/(pow((pow(((rrx-rmx)),2)+(pow(((rry-rmy)),2))),(3/2)))); /*Acceleration in y-direction*/
vx=vx+(dt*ax); /*Step 4 of E