万有引力を考慮した質点の運動

実現したいこと

万有引力を考慮した物体の運動をプロットしたい

前提

大学の講義でc言語を用いて万有引力を用いた質点の運動を表すコードを書いているのですが、うまくいきません。初学者なので汚いコードですみません。

発生している問題・エラーメッセージ

エラーコードは出ないのですが、実際にプロットしてみると一次関数のような形になってしまいます。

該当のソースコード

c

1#include <stdio.h>2#include <math.h>3 4int main(){5 double v0,v,vx,vy,x0,y0,x,y,dt,a,r,G,M,ax,ay,rad;6 int i;7 8 v0 = 5;//初期速度9 v = v0;//速度vの初期化10 11 dt = 0.0001;//微小時間変化12 13 G = 6.67*pow(10,-11);//万有引力定数14 15 x0 = 1;16 y0 = 3;17 18 rad = atan2(y0,x0);19 vx = v * cos(rad);//x軸方向の速度20 vy = v * sin(rad);//y軸方向の速度21 22 M = pow(10,14);//中心にある物体の質量23 24 x = x0;//xの初期化25 y = y0;//yの初期化26 27 for(i = 0;i <= 100000;i++){28 rad = atan2(y,x);//偏角を計算29 r = sqrt(pow(x,2)+pow(y,2));//原点からの距離30 31 a = G*M/pow(r,2);//中心方向に対する加速度32 ax = -a * cos(rad);//x軸方向に対する加速度33 ay = -a * sin(rad);//y軸方向に対する加速度34 35 vx = vx + ax * dt;//x軸に対する速度をオイラー法で計算36 vy = vy + ay * dt;//y軸に対する速度をオイラー法で計算37 38 x = x + vx * dt;//xをオイラー法で計算39 y = y + vy * dt;//yをオイラー法で計算40 printf("%f %f\n",x,y);//xとyを出力41 }42}

試したこと

物理的な考え方が間違っているのかと思い、sinやcosをいじったり、radにM_PIを足すなどしてみましたがやはりうまくいきませんでした。