#include<stdio.h>
#include<stdlib.h>
#include<math.h>

	/*  Zadeklarujemy sobie nowy typ danych punkt. Posiada dwie składowe typu double:
	    x1 i x2. Moglibyśmy użyć tablicy dwuelementowej
	*/

struct punkt{
	double x1,x2;
};

	/* Następująca funkcja oblicza pochodną w punkcie (prawą stronę równania). Jako parametry
	   przyjmuje czas t i wartość rozwiązania y(t), zwraca pochodną y'(t).
	*/

struct punkt pochodna(double t, struct punkt y){
	struct punkt tymcz;
	tymcz.x1=y.x2;
	tymcz.x2=-sin(y.x1);
	return tymcz;
}
int main(){
	struct punkt y,tymcz,k1,k2,k3,k4;
	double czas_biezacy, czas_koncowy,krok,polkrok;
	int ilosc_krokow;
	printf("Program służy do numerycznego rozwiązywania równań różniczkowych\n");
	printf("zwyczajnych metodą Rungego-Kutty 4 rzędu.\n\n");
	printf("Czas startowy przyjmujemy jako 0.\n");
	printf("Proszę podać wychylenie początkowe w radianach y(0)= ");
	scanf("%lf",&(y.x1));
	printf("\nTeraz proszę podać początkową prędkość kątową w radianach na sekundę y'(0)= ");
	scanf("%lf",&(y.x2));
	printf("\nProszę podać czas koncowy obliczeń T= ");
	scanf("%lf",&czas_koncowy);
	printf("\nProszę podać liczbę kroków N= ");
	scanf("%d",&ilosc_krokow);
	krok=czas_koncowy/ilosc_krokow;
	polkrok=krok/2;
	czas_biezacy=0;
	for(;ilosc_krokow&>;=0;--ilosc_krokow){
		printf("| %2.6lf | %2.6lf |\n",czas_biezacy,y.x1);

	/* Następujące obliczenia są zapisane w sposób znacznie bardziej skomplikowany, niż
	   mogłyby być. Wynika to z faktu, że nie zdefiniowaliśmy mnożenia naszej struktury typu punkt
	   przez liczby. Musimy więc każdorazowo mnożyć "po współrzędnych".
	*/

		tymcz=pochodna(czas_biezacy,y);
		k1.x1=polkrok*tymcz.x1;
		k1.x2=polkrok*tymcz.x2;
		tymcz.x1=y.x1+k1.x1;
		tymcz.x2=y.x2+k1.x2;
		czas_biezacy=czas_biezacy+polkrok;
		tymcz=pochodna(czas_biezacy, tymcz);
		k2.x1=polkrok*tymcz.x1;
		k2.x2=polkrok*tymcz.x2;
		tymcz.x1=y.x1+k2.x1;
		tymcz.x2=y.x2+k2.x2;
		tymcz=pochodna(czas_biezacy, tymcz);
		k3.x1=krok*tymcz.x1;
		k3.x2=krok*tymcz.x2;
		czas_biezacy=czas_biezacy+polkrok;
		tymcz.x1=y.x1+k3.x1;
		tymcz.x2=y.x2+k3.x2;
		tymcz=pochodna(czas_biezacy, tymcz);
		k4.x1=krok*tymcz.x1;
		k4.x2=krok*tymcz.x2;
		y.x1=y.x1+k1.x1/3+2*k2.x1/3+k3.x1/3+k4.x1/6;
		y.x2=y.x2+k1.x2/3+2*k2.x2/3+k3.x2/3+k4.x2/6;
	};
	return EXIT_SUCCESS;
}