#include <fstream.h>
#include <iomanip.h>

// Chapter 8.1, p.375
// Euler
template <class Type>
void euler(Type a, Type b, Type& x, int n)
{
  runtype h, t;
  h = (b - a)/n;
  t = a;
  outfile1 << "k=  0:x(" << setw(5) << t << ") = " << setw(5) << x << endl;

  for (int k = 1; k <= n; k++)
    {
      x += h * f1(t, x);
      t += h;
	  outfile1 << "k=" <<setw(3) << k <<":x("<<setw(5) << t << ") =" << setw(5) << x <<endl;
    }
}

// Chapter 8.1, p.378
// Taylor
template <class Type>
void taylor(Type a, Type b, Type& x, int n)
{
  runtype h, t;
  runtype  x1, x2, x3, x4;
  h = (b - a)/n;
  t = a;
 
  outfile2 << "k=  0:x(" << setw(5) << t << ") = " << setw(5) << x << endl;

  for (int k = 1; k <= n; k++)
  {
    x1 = 1 + x*x + t*t*t;
    x2 = 2*x*x1 + 3*t*t;
    x3 = 2*x*x2 + 2*x1*x1 + 6*t;
    x4 = 2*x*x3 + 6*x1*x2 + 6;
    x += h*(x1 + h/2*(x2 + h/3*(x3 + h/4*x4)));
    t = a + k*h;
	outfile2 << "k=" <<setw(3) << k <<":x("<<setw(5) << t << ") =" << setw(5) << x <<endl;
  }
}
// Chapter 8.2, p.387-38
// Runge-Kutta 4
template <class Type>
void rk4(ofstream& out, Type f(Type, Type), Type t, Type x, Type h, int n)
{
  int j;
  Type K1, K2, K3, K4, ta;

  out << "j=  0: x("<< setw(5) << t <<") =" << setw(5) << x<< endl;
  ta = t;
  for (j = 1; j <= n; j++)
  {
    K1 = h*f(t,x);
    K2 = h*f(t + h/2, x + K1/2);
    K3 = h*f(t + h/2, x + K2/2);
    K4 = h*f(t + h, x + K3);
    x += (K1 + 2*(K2 + K3) + K4)/6;
    t = ta + j*h;
	out<< "j=" << setw(3) << j << ": x("<< setw(5) << t << ") =" << setw(5) <<x<< endl;
  }
}


// Chapter 8.3, p.396-397
// Runge-Kutta-Fehlberg
template <class Type>
void rk45(Type f(Type, Type), Type& t, Type& x, Type h, Type& epsilon)
{

  Type c20 = 0.25, c21 = 0.25;
  Type c30 = 0.375, c31 = 0.09375, c32 = 0.28125;
  Type c40,c41, c42,c43;
  Type c51, c52 = -8.0, c53, c54;
  Type c60 = 0.5, c61, c62 = 2, c63, c64;
  Type c65 = -0.275;
  Type a1, a2 = 0, a3, a4, a5 = -0.2;
  Type b1, b2 = 0, b3, b4, b5= -0.18, b6;
  Type K1, K2, K3, K4, K5, K6, x4;

  c40 = (Type) 12/ (Type) 13;
  c41 = (Type) 1932/(Type) 2197;
  c42 = (Type) -7200/(Type) 2197;
  c43 = (Type) 7296/(Type) 2197;
  c51 = c53 = (Type) 439/ (Type) 216;
  c54 = (Type) -845/(Type) 4104;
  c61 = (Type) -8/(Type) 27;
  c63 = (Type) -3544/(Type) 2565;
  c64 = (Type) 1859/(Type) 4104;
  a1 = (Type) 25/(Type) 216;
  a3 = (Type) 1408/(Type) 2565;
  a4 = (Type) 2197/(Type) 4104;
  b1 = (Type) 16/(Type) 135;
  b3 = (Type) 6656/(Type) 12825;
  b4 = (Type) 28561/(Type) 56430;
  b6 = (Type) 2/(Type) 55;


  K1 = h*f(t, x);
  K2 = h*f(t + c20*h, x + c21*K1);
  K3 = h*f(t + c30*h, x + c31*K1 + c32*K2);
  K4 = h*f(t + c40*h, x + c41*K1 + c42*K2 + c43*K3);
  K5 = h*f(t + h, x + c51*K1 + c52*K2 + c53*K3 + c54*K4 );
  K6 = h*f(t + c60*h, x + c61*K1 + c62*K2 + c63*K3 + c64*K4 + c65*K5);

  x4 = x + a1*K1 + a3*K3 + a4*K4 + a5*K5;
  x += b1*K1 + b3*K3 + b4*K4 + b5*K5 + b6*K6;
  t += h;
  epsilon = fabs(x - x4);

}


// Chapter 8.3, p.398-399
// Adaptive Runge-Kutta


//-----------------------------------------------------------------------------
template <class Type>
void rk45ad (ofstream& out, Type f(Type, Type), Type t, Type x, Type h, Type tb,
           int itmax, Type emin, Type emax, Type hmin, Type hmax)
{
  Type delta, epsilon, d, xsave, tsave;
  int k, iflag;

  out << "n         h          t          x"<< endl;
  out << "0" << setw(10) << h << setw(10) << t << setw(10) << x << endl;

  delta = 0.00005;
  iflag = 1;
  k = 0;
  while (k <= itmax)
  {
    k++;
    if (fabs(h) < hmin)
      h = sign(h)*hmin;
    if (fabs(h) > hmax)
      h = sign(h)*hmax;
    d = fabs(tb - t);
    if (d <= fabs(h) )
    {
      iflag = 0;
      if (d <= delta*max(fabs(tb), fabs(t)))
	break;
      h = sign(h)*d;
    }
    xsave = x;
    tsave = t;
    rk45(f, t, x, h, epsilon);
	out << "k =" << setw(3) << k << "  h =" << setw(10) << h 
		<<"  t =" << setw(10) << t <<"  x =" << setw(10) << x 
		<<"  epsilon =" << setw(10) << epsilon<< endl;
    if (iflag == 0)
      break;
    if (epsilon < emin)
      h = 2*h;
    if (epsilon > emax)
    {
      h = h/2;
      x = xsave;
      t = tsave;
      k--;
    }
  }
  out<< "itmax =" << setw(5) << itmax 
	  <<"  iflag =" << setw(5) << iflag;
}