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

template<typename T>
void exercise1() {
  constexpr int lambda =  -1;
  constexpr T y_0 = 1.0, t_min = 0.0, t_max = 5.0;
  std::vector<std::ofstream> files;
  files.emplace_back(std::ofstream {"zad1skok1.dat"});
  files.emplace_back(std::ofstream {"zad1skok2.dat"});
  files.emplace_back(std::ofstream {"zad1skok3.dat"});
  const std::vector<T> deltas_t = {0.01, 0.1, 1.0};

  for (int i {}; i < 3; ++i) {
    T y = y_0, y1;
    for (T t = t_min; t <= t_max; t += deltas_t[i]){
      y1 = y + deltas_t[i] * lambda * y;
      files[i] << t << " " << y1 << " " << exp(lambda * t) << " " << (y1 - exp(lambda * t)) << std::endl;
      y = y1;
    }
  }
}

template<typename T>
void exercise2() {
  constexpr int lambda =  -1;
  constexpr T y_0 = 1.0, t_min = 0.0, t_max = 5.0;
  std::vector<std::ofstream> files;
  files.emplace_back(std::ofstream {"zad2skok1.dat"});
  files.emplace_back(std::ofstream {"zad2skok2.dat"});
  files.emplace_back(std::ofstream {"zad2skok3.dat"});
  const std::vector<T> deltas_t = {0.01, 0.1, 1.0};

  for (int i {}; i < 3; ++i) {
    T y = y_0, y1;
    for (T t = t_min; t <= t_max; t += deltas_t[i]){
      T k1 = lambda * y;
      T k2 = lambda * (y + deltas_t[i] * k1);
      y1 = y + (k1 + k2) * (deltas_t[i] / 2.0);
      files[i] << t << " " << y1 << " " << exp(lambda * t) << " " << (y1 - exp(lambda * t)) << std::endl;
      y = y1;
    }
  }
}

template<typename T>
void exercise3() {
  constexpr int lambda =  -1;
  constexpr T y_0 = 1.0, t_min = 0.0, t_max = 5.0;
  std::vector<std::ofstream> files;
  files.emplace_back(std::ofstream {"zad3skok1.dat"});
  files.emplace_back(std::ofstream {"zad3skok2.dat"});
  files.emplace_back(std::ofstream {"zad3skok3.dat"});
  const std::vector<T> deltas_t = {0.01, 0.1, 1.0};

  for (int i {}; i < 3; ++i) {
    T y = y_0, y1;
    for (T t = t_min; t <= t_max; t += deltas_t[i]){
      T k1 = lambda * y;
      T k2 = lambda * (y + k1 * deltas_t[i] / 2.0);
      T k3 = lambda * (y + k2 * deltas_t[i] / 2.0);
      T k4 = lambda * (y + deltas_t[i] * k3);
      y1 = y + (k1 + 2 * k2 + 2 * k3 + k4) * (deltas_t[i] / 6.0);
      files[i] << t << " " << y1 << " " << exp(lambda * t) << " " << (y1 - exp(lambda * t)) << std::endl;
      y = y1;
    }
  }
}

template<typename T>
void exercise4() {
  constexpr T  delta_t = 0.0001, t_min = 0;

}

int main()
{
  exercise1<double>();
  exercise2<double>();
  exercise3<double>();

  return 0;
}


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

// float fun (float x){
//  return exp(-x);
// }

// float V (float t, float omega){
//  return 10.*sin(omega*t);
// }

// int main (void){




//  //drugiego rzedu

//  fp1 = fopen("dane41.dat", "w");
//  fp2 = fopen("dane42.dat", "w");
//  fp3 = fopen("dane43.dat", "w");
//  FILE *fp4 = fopen("dane44.dat", "w");

//  float dt = pow(10, -4), R=100., L=0.1, C = 0.001, omega0 =1./sqrt(L*C), T0 = 2*M_PI/omega0;
//  //float t, tmin, tmax;
//  tmin = 0;
//  tmax = 4*T0;

//  float Q0 = 0, Q = Q0, I0 = 0, I = 0, omegaV;
//  float k1Q, k1I, k2Q, k2I, k3Q, k3I, k4Q, k4I;

//  omegaV = 0.5*omega0;

//  for (t = tmin; t <= tmax; t+=dt){
//      k1Q = I;
//      k1I = V(t, omegaV)/L - Q/(L*C) - R*I/L;
//      k2Q = I + dt*k1I/2.;
//      k2I = V(t + 0.5*dt, omegaV)/L + - (Q + dt*k1Q/2.)/(L*C) - R*(I + dt*k1I/2.)/L;
//      k3Q = I + dt*k2I/2.;
//      k3I = V(t + 0.5*dt, omegaV)/L + - (Q + dt*k2Q/2.)/(L*C) - R*(I + dt*k2I/2.)/L;
//      k4Q = I + dt*k3I;
//      k4I = V(t + dt, omegaV)/L + - (Q + dt*k3Q)/(L*C) - R*(I + dt*k3I)/L;
//      Q = Q + dt/6.*(k1Q + 2*k2Q + 2*k3Q + k4Q);
//      I = I + dt/6.*(k1I + 2*k2I + 2*k3I + k4I);
//      fprintf(fp1, "%f %f %f\n", t, Q, I);
//  }

//  omegaV = 0.8*omega0;

//  for (t = tmin; t <= tmax; t+=dt){
//      k1Q = I;
//      k1I = V(t, omegaV)/L - Q/(L*C) - R*I/L;
//      k2Q = I + dt*k1I/2.;
//      k2I = V(t + 0.5*dt, omegaV)/L + - (Q + dt*k1Q/2.)/(L*C) - R*(I + dt*k1I/2.)/L;
//      k3Q = I + dt*k2I/2.;
//      k3I = V(t + 0.5*dt, omegaV)/L + - (Q + dt*k2Q/2.)/(L*C) - R*(I + dt*k2I/2.)/L;
//      k4Q = I + dt*k3I;
//      k4I = V(t + dt, omegaV)/L + - (Q + dt*k3Q)/(L*C) - R*(I + dt*k3I)/L;
//      Q = Q + dt/6.*(k1Q + 2*k2Q + 2*k3Q + k4Q);
//      I = I + dt/6.*(k1I + 2*k2I + 2*k3I + k4I);
//      fprintf(fp2, "%f %f %f\n", t, Q, I);
//  }

//  omegaV = 1.0*omega0;

//  for (t = tmin; t <= tmax; t+=dt){
//      k1Q = I;
//      k1I = V(t, omegaV)/L - Q/(L*C) - R*I/L;
//      k2Q = I + dt*k1I/2.;
//      k2I = V(t + 0.5*dt, omegaV)/L + - (Q + dt*k1Q/2.)/(L*C) - R*(I + dt*k1I/2.)/L;
//      k3Q = I + dt*k2I/2.;
//      k3I = V(t + 0.5*dt, omegaV)/L + - (Q + dt*k2Q/2.)/(L*C) - R*(I + dt*k2I/2.)/L;
//      k4Q = I + dt*k3I;
//      k4I = V(t + dt, omegaV)/L + - (Q + dt*k3Q)/(L*C) - R*(I + dt*k3I)/L;
//      Q = Q + dt/6.*(k1Q + 2*k2Q + 2*k3Q + k4Q);
//      I = I + dt/6.*(k1I + 2*k2I + 2*k3I + k4I);
//      fprintf(fp3, "%f %f %f\n", t, Q, I);
//  }

//  omegaV = 1.2*omega0;

//  for (t = tmin; t <= tmax; t+=dt){
//      k1Q = I;
//      k1I = V(t, omegaV)/L - Q/(L*C) - R*I/L;
//      k2Q = I + dt*k1I/2.;
//      k2I = V(t + 0.5*dt, omegaV)/L + - (Q + dt*k1Q/2.)/(L*C) - R*(I + dt*k1I/2.)/L;
//      k3Q = I + dt*k2I/2.;
//      k3I = V(t + 0.5*dt, omegaV)/L + - (Q + dt*k2Q/2.)/(L*C) - R*(I + dt*k2I/2.)/L;
//      k4Q = I + dt*k3I;
//      k4I = V(t + dt, omegaV)/L + - (Q + dt*k3Q)/(L*C) - R*(I + dt*k3I)/L;
//      Q = Q + dt/6.*(k1Q + 2*k2Q + 2*k3Q + k4Q);
//      I = I + dt/6.*(k1I + 2*k2I + 2*k3I + k4I);
//      fprintf(fp4, "%f %f %f\n", t, Q, I);
//  }
//  fclose(fp1);
//  fclose(fp2);
//  fclose(fp3);
//  fclose(fp4);

//  return 0;
// }