我应该积分微分方程$r$和$\theta$为了模拟轨道运动。我用的微分方程$r$是二阶并且$\theta$是第一顺序。最终目标是绘制图表$r(t)$,$\theta(t)$，我可能会在 python 中使用 matplotlib。我想知道是否有事先了解轨道和 Runge-Kutta 的人可以让我知道我的代码哪里出错了？它编译没有错误，但没有返回我期望的值。这是我的代码，其中包含描述地球轨道的参数：

#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) */
}

我用于位置的二阶微分方程，$r$， 是：

$$\ddot{r} = \frac{-k}{mr^2} + \frac{l^2}{mr^3}$$

参数在哪里$l$是角动量和$k$是引力常数乘以行星运行的物体的质量，在这种情况下是太阳。我用于的一阶微分方程$\theta$是：

$$\dot{\theta} = \frac{h}{r^2}$$

$$h = l\left(\frac{m_1 m_2}{m_1 + m_2}\right)^{-1}$$

在哪里$m_1$和$m_2$是系统中两个物体的质量，方程中的质量项为$h$是减少的质量。所有参数的单位都是 AU、年和太阳质量。如果有人能指出我正确的方向，那就太棒了。