Program ODE;

CONST
     h=0.01;
     pi=3.141592654;
     c=0.1;
     m=1;
     g=1;

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

PROCEDURE f(tn:real; w:vektor; var fi:vektor);
VAR
	d:real;

BEGIN

        d:=c*(w[3]*w[3]);

	fi[1]:=w[3]*cos(w[4]);
	fi[2]:=w[3]*sin(w[4]);
	fi[3]:=-d/m -g*sin(w[4]);
        fi[4]:=-g/w[3]*cos(w[4]);

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;

VAR
	wr:vektor;
	i:integer;
	tn,cil:real;
	u,xf:vektor2;
	s:string[2];
	s2:string[11];
	fileout:text;

BEGIN

	u[1]:=45;
	u[2]:=5;
	u[3]:=(u[1]+u[2])/2;
	cil:=11;

	FOR i:=1 TO 3 DO BEGIN
		Str(i,s);
		s2:='strel'+s+'.dat';
		ASSIGN(fileout,s2);
		REWRITE(fileout);

		tn:=0;
		wr[1]:=0;
		wr[2]:=0;
		wr[3]:=4;
		wr[4]:=u[i]/180*pi;

		WHILE (wr[2]>=0) DO BEGIN
			tn:=tn+h;
			RKstep(tn,wr);
			WRITELN(fileout,tn:5:3,' ',wr[1]:5:3,' ',wr[2]:5:3,' ',wr[3]:5:3,' ',wr[4]:5:3);
		END;
		xf[i]:=wr[1]-cil;
		CLOSE(fileout);

	END;

	FOR i:=4 TO 14 DO BEGIN
		Str(i,s);
		s2:='strel'+s+'.dat';
		ASSIGN(fileout,s2);
		REWRITE(fileout);

		IF(xf[3]<0) THEN BEGIN
			xf[2]:=xf[3];
			u[2]:=u[3];
			u[3]:=(u[1]+u[2])/2;
		END
		ELSE BEGIN
			xf[1]:=xf[3];
			u[1]:=u[3];
			u[3]:=(u[1]+u[2])/2;
		END;

		tn:=0;
		wr[1]:=0;
		wr[2]:=0;
		wr[3]:=4;
		wr[4]:=u[3]/180*pi;

		WHILE (wr[2]>=0) DO BEGIN
			tn:=tn+h;
			RKstep(tn,wr);
			WRITELN(fileout,tn:5:3,' ',wr[1]:5:3,' ',wr[2]:5:3,' ',wr[3]:5:3,' ',wr[4]:5:3);
		END;
		xf[3]:=wr[1]-cil;
		CLOSE(fileout);

	END;

	WRITELN('Hotovo');
	READLN;

END.