#include <iostream>
#include <cmath>
#include <fstream>



 using namespace std;

 const float g = 10;
 const float l = 50;
 const float gamma = 0.1 ;
const float omegad=10;
const float A=0;

 struct state
 {
 float theta;
 float epsilon;
 };

 struct derivative
 {
 float dtheta;
 float depsilon;

 derivative& operator*(float& num)
 {
	 dtheta *= num;
	 depsilon *= num;
	 return *this;
 }
 };

 derivative evaluate(state actual, float t, float h, derivative d)
 {
 state s;
 derivative out;

 s.theta = actual.theta + d.dtheta*h;
 s.epsilon = actual.epsilon + d.depsilon*h;

 out.dtheta = s.epsilon;
 out.depsilon = A*cos(omegad*t)-gamma * s.epsilon-g/l*sin(s.theta);

 return out;
 }

 state rk4(state s, float t, float h)
 {
 derivative d0;
 float delta_theta, delta_epsilon;

 d0.depsilon = 0;
 d0.dtheta = 0;

 derivative k1 = evaluate(s, t, 0, d0);
 derivative k2 = evaluate(s, t, h*0.5f, k1*0.5);
 derivative k3 = evaluate(s, t, h*0.5f, k2*0.5);
 derivative k4 = evaluate(s, t, h, k3);

 delta_theta = 1.0f/6.0f * h * (k1.dtheta + 2.0f*k2.dtheta + 2.0f*k3.dtheta + k4.dtheta);
 delta_epsilon = 1.0f/6.0f * h * (k1.depsilon + 2.0f*k2.depsilon + 2.0f*k3.depsilon + k4.depsilon);

 s.theta = s.theta + delta_theta;
 s.epsilon = s.epsilon + delta_epsilon;

 return s;
 }

 int main()
 {
     fstream zapis;
     zapis.open("zapis.txt",ios::out);
 float t = 0;
 float h = 0.1f;
 state s;

 s.theta = 15;
 s.epsilon = 0;

 for(t=0; t<50; t+=h)
 {
 s = rk4(s, t, h);
 zapis<<t<<" "<<s.theta<<endl;
 }
zapis.close();
 return 0;
 }