Program ODE;

CONST
     h=0.1;

TYPE
    vektor=array[1..4] of real;

PROCEDURE f(tn:real; w:vektor; var fi:vektor);
VAR
	r:real;
	
BEGIN
	r:=w[1]*w[1]+w[3]*w[3];

	fi[1]:=w[2];
	fi[2]:=-w[1]/sqrt(r*r*r);
	fi[3]:=w[4];
	fi[4]:=-w[3]/sqrt(r*r*r);	
END;    

PROCEDURE RKstep(tn:real; var w:vektor);
VAR
	k1,k2,k3,k4,fi,pom:vektor;
	i:integer;
	
BEGIN
	f(tn,w,fi);
	FOR i:=1 TO 4 DO k1[i]:=h*fi[i];
	
	FOR i:=1 TO 4 DO pom[i]:=w[i]+k1[i]/2;
	f(tn+h/2,pom,fi);
	FOR i:=1 TO 4 DO k2[i]:=h*fi[i];

 	FOR i:=1 TO 4 DO pom[i]:=w[i]+k2[i]/2;
	f(tn+h/2,pom,fi);
	FOR i:=1 TO 4 DO k3[i]:=h*fi[i];

 	FOR i:=1 TO 4 DO pom[i]:=w[i]+k3[i];
	f(tn+h,pom,fi);
	FOR i:=1 TO 4 DO k4[i]:=h*fi[i];
	
 	FOR i:=1 TO 4 DO BEGIN
 		w[i]:=w[i]+(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6;
 	END;
	
END;

PROCEDURE Euler(tn:real; var w:vektor);
VAR
	fi:vektor;
	i:integer;
	
BEGIN
	f(tn,w,fi);
	FOR i:=1 TO 4 DO w[i]:=w[i]+h*fi[i];
END;

VAR
	wr,we:vektor;
	i:integer;
	tn:real;
	fileout:text;
	
BEGIN
	ASSIGN(fileout,'rk.dat');
	REWRITE(fileout);
	
	tn:=0;
	wr[1]:=1;
	wr[2]:=0;
	wr[3]:=0;
	wr[4]:=1;

	we[1]:=1;
	we[2]:=0;
	we[3]:=0;
	we[4]:=1;
	
	FOR i:=1 TO 300 DO BEGIN
		tn:=i*h;
		RKstep(tn,wr);
		Euler(tn,we);
		WRITELN(fileout,tn:12:6,' ',wr[1]:12:6,' ',wr[3]:12:6,' ',we[1]:12:6,' ',we[3]:12:6);
	END;

	CLOSE(fileout);
        writeln('Hotovo');
	READLN;
END.