Numerical Methods in Pascal: Euler 2D

Program Planetary_Motion;
{https://syeilendrapramuditya.wordpress.com}
{Gerak bulan mengitari bumi karena pengaruh gaya gravitasi}
{Data t[i],x[i],y[i],vx[i],vy[i] ditulis ke dalam file excel (D:\orbit.xls) }
{Data D:\orbit.xls tersebut kemudian diplot menggunakan Excel}
uses crt;
const GM = 3.986344692e14;
Type
 DataPlot = record
            dataT,dataX,dataY,dataVX,dataVY : real;
            end;
 ListData = array[0..1000] of DataPlot;
var t0,x0,y0,vx0,vy0,m,k,b,t,x,y,vx,vy,x12,y12,vx12,vy12,dt : real;
    i,N : integer;
    FileData : text;
    S:ListData;

Begin
clrscr;
writeln;
writeln('      Solusi Numerik PDB Orde 2 Terkopling');
writeln('              Kasus  Gerak  Bulan         ');
writeln;
writeln('Nilai konstanta yang digunakan dalam perhitungan :');
writeln('Massa Bumi  = 5.9742e24 Kg               ');
writeln('Massa Bulan = 7.3500e24 Kg               ');
writeln('d           = 3.844e5 Km                 ');
writeln('T           = 27.322 hari                ');
writeln('G           = 6.6726e-11 N*sqr(m)/sqr(Kg ');
writeln;
writeln('Nilai parameter (default) yang digunakan dalam perhitungan :');
writeln;
writeln('x0  = 3.844e5 ');
writeln('y0  = 0       ');
writeln('vx0 = 0       ');
writeln('vy0 = 32000   ');
writeln('dt  = 1       ');
writeln('N   = 1000    ');
writeln;
write('Masukan nilai x0,y0,vx0,vy0 pada t = 0  : ');readln(x0,y0,vx0,vy0);
write('Masukan parameter iterasi (dt dan N) : ');readln(dt,N);
writeln;
t0 := 0;
Assign(FileData,'D:\orbit.xls');
rewrite(FileData);
 S[0].DataT  := 0;
 S[0].DataX  := x0;
 S[0].DataY  := y0;
 S[0].DataVX := vx0;
 S[0].DataVY := vy0;
 writeln(FileData,'      t               x               y               vx               vy');
 writeln(FileData,s[0].dataT:20:5,s[0].dataX:20:5,s[0].DataY:20:5,s[0].DataVX:20:5,s[0].DataVY:20:5);
for i := 1 to N do
 begin
 x12  := x0 + (dt/2) * vx0;
 vx12 := vx0 - ( (GM * x0) /  ( exp(1.5 * ln( sqr(x0) + sqr(y0) ))) ) * (dt/2);
 y12  := y0 + (dt/2) * vy0;
 vy12 := vy0 - ( (GM * y0) /  ( exp(1.5 * ln( sqr(x0) + sqr(y0) ))) ) * (dt/2);
 x  := x0 + dt * vx12;
 vx := vx0 - ( (GM * x12)  /  (exp(1.5 * ln(sqr(x0) + sqr(y0) )))) * (dt);
 y  := y0 + dt * vy12;
 vy := vy0 - ( (GM * y12)  /  (exp(1.5 * ln(sqr(x0) + sqr(y0) )))) * (dt);
 t := t0 + dt;
 S[i].DataT  :=  t;
 S[i].DataX  :=  x;
 S[i].DataY  :=  y;
 S[i].DataVX := vx;
 S[i].DataVY := vy;
 writeln(FileData,s[i].dataT:20:5,s[i].dataX:20:5,s[i].DataY:20:5,s[i].DataVX:20:5,s[i].DataVY:20:5);
 t0  :=  t;
 x0  :=  x;
 y0  :=  y;
 vx0 := vx;
 vy0 := vy;
 end;
close(FileData);
writeln('Data t,x,y,vx,vy telah disimpan di D:\orbit.xls , silahkan buka file tersebut');
readln;
end.

• programming