C言語でケプラーの公式
１．惑星は太陽を一つの焦点とする楕円軌道を描いて太陽のまわりをまわるというプログラムは
#include <stdio.h>
#include <math.h>
#define GM 4*M_PI*M_PI
#define M 1.0
#define dT 1.0/256
#define QX0 1.0
#define QY0 0.0
#define PX0 0.0
#define PY0 sqrt(GM)*M
int main(void)
{
FILE*fout;
double qx, qy;
double px, py;
double T;
double qxk1, qyk1, pxk1, pyk1;
double qxk2, qyk2, pxk2, pyk2;
double qxk3, qyk3, pxk3, pyk3;
double qxk4, qyk4, pxk4, pyk4;
double tqx, tqy, tpx, tpy;
double q, qi3;
fout = fopen("data16.dat","w");
qx=QX0;qy=QY0;px=PX0;py=PY0;
for(T=0.0;T<10.00;T+=dT)
{
tqx = qx; tqy = qy; tpx = px; tpy = py;
q = hypot( tqx, tqy ); qi3 = 1.0/(q*q*q);
qxk1 = tpx / M * dT;
qyk1 = tpy / M * dT;
pxk1 = -GM * M * tqx * qi3 * dT;
pyk1 = -GM * M * tqy * qi3 * dT;
tqx = qx+0.5*qxk1;
tqy = qy+0.5*qyk1;
tpx = px+0.5*pxk1;
tpy = py+0.5*pyk1;
q = hypot( tqx, tqy ); qi3 = 1.0/(q*q*q);
qxk2 = tpx / M * dT;
qyk2 = tpy / M * dT;
pxk2 = -GM * M * tqx * qi3 * dT;
pyk2 = -GM * M * tqy * qi3 * dT;
tqx = qx+0.5*qxk2;
tqy = qy+0.5*qyk2;
tpx = px+0.5*pxk2;
tpy = py+0.5*pyk2;
q = hypot( tqx, tqy );
qi3 = 1.0/(q*q*q);
qxk3 = tpx / M * dT;
qyk3 = tpy / M * dT;
pxk3 = -GM * M * tqx * qi3 * dT;
pyk3 = -GM * M * tqy * qi3 * dT;
tqx = qx+qxk3;
tqy = qy+qyk3;
tpx = px+pxk3;
tpy = py+pyk3;
q = hypot( tqx, tqy );
qi3 = 1.0/(q*q*q);
qxk4 = tpx / M * dT;
qyk4 = tpy / M * dT;
pxk4 = -GM * M * tqx * qi3 * dT;
pyk4 = -GM * M * tqy * qi3 * dT;
qx += (qxk1 + 2*qxk2 + 2*qxk3 + qxk4)*(1.0/6);
qy += (qyk1 + 2*qyk2 + 2*qyk3 + qyk4)*(1.0/6);
px += (pxk1 + 2*pxk2 + 2*pxk3 + pxk4)*(1.0/6);
py += (pyk1 + 2*pyk2 + 2*pyk3 + pyk4)*(1.0/6);
printf("qx=%f qy=%f px=%f py=%f \n",qx,qy,px,py);
fprintf(fout,"%f, %f\n",qx,qy);
}
fclose(fout);
return 0;
}
と作れたのですが、
2.面積速度は一定である。
3.惑星の公転周期Tの２乗と軌道長半径aの３乗の比T^2/a^3は太陽系のすべての惑星に共通で一定のあたいになる。
を数値的に示せ。
というプログラムが作れません。
教えてください。
お礼
foomufoomuさん、ありがとうございます。 τxの作用面も長方形ですが、この面は1.5倍しなくてもよいのでしょうか。