I'm supposed to integrate differential equations for $r$ and $\theta$ in order to simulate orbital motion. The differential equation I used for $r$ is second order and $\theta$ is first order. The end goal is to graph $r(t)$, $\theta(t)$, which I'll probably do with matplotlib in python. I was wondering if anyone with prior knowledge of orbits and Runge-Kutta could let me know where I went wrong with my code? It compiles without error but isn't returning values that I expected. This is my code with parameters describing Earth's orbit:
#include <stdio.h>
#include <math.h>
#define N 3 /* number of first order equations */
#define dist 0.00001 /* stepsize in t*/
#define MAX 1 /* max for t */
#define MAXTHETA 2*M_PI
FILE *output; /* internal filename */
void runge4(double x, double y[], double step);
double f(double x, double y[], int i);
int main()
{
double t, y[N];
int j;
void runge4(double x, double y[], double step); /* Runge-Kutta function */
double f(double x, double y[], int i); /* function for derivatives */
output=fopen("kepler.rtf", "w"); /* external filename */
y[0]= 1.02343552273569; /* initial position */
y[1]= 6.18068911; /* initial velocity */
y[2] = 0; /* initial angle */
fprintf(output, "0\t%f\t%f\n", y[0], y[2]); /* prints time and position */
for (j=1; j*dist<=MAX ;j++) /* time loop */
{
t=j*dist;
runge4(t, y, dist);
fprintf(output, "%f\t%f\t%f\n", t, y[0], y[2]); /* prints time, position, theta */
}
fclose(output);
}
void runge4(double x, double y[], double step)
{
double h=step/2.0, /* the midpoint */
t1[N], t2[N], t3[N], /* temporary storage arrays */
k1[N], k2[N], k3[N],k4[N]; /* for Runge-Kutta */
int i;
for (i=0;i<N;i++) t1[i]=y[i]+0.5*(k1[i]=step*f(x, y, i));
for (i=0;i<N;i++) t2[i]=y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
for (i=0;i<N;i++) t3[i]=y[i]+ (k3[i]=step*f(x+h, t2, i));
for (i=0;i<N;i++) k4[i]= step*f(x+step, t3, i);
for (i=0;i<N;i++) y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
double f(double x, double y[], int i)
{
double k = 39.4784176044, m = 0.000003003, h = 0.00000300299, l = 0.0000188781;
if (i==0) return(y[1]); /* derivative of first equation (dr/dt)*/
if (i==1) return((-k/(m*y[0]*y[0])) + (l*l)/(m*y[0]*y[0]*y[0])); /* derivative of second equation (d^2r/dt^2) */
else return(h/(y[0]*y[0])); /* derivative of theta (d(theta)/dt) */
}
The second order differential equation I used for the position, $r$, is:
$$\ddot{r} = \frac{-k}{mr^2} + \frac{l^2}{mr^3}$$
where the parameter $l$ is the angular momentum and $k$ is the gravitational constant multiplied my the mass of the object that the planet is orbiting, in this case the sun. The first order differential equation I used for $\theta$ is:
$$\dot{\theta} = \frac{h}{r^2}$$
$$h = l\left(\frac{m_1 m_2}{m_1 + m_2}\right)^{-1}$$
where $m_1$ and $m_2$ are the masses of the two bodies in the system, and the mass term in the equation for $h$ is the reduced mass. All of the units of parameters are in AU, years, and solar masses. If anybody could point me in the right direction that would be awesome.