/*
 Solver for second order Ordinary Differential Equations

 x"(t) = f(t,x,x')    equation
 x(ti) = xi           initial value 1
 x'(ti)= vi           initial value 2

 the second order ODE is solved as a set of two first order ODEs
 x'(t) = f1(t,x)
 x"(t) = f2(t,x,x')

 Methods (select one by a key)
 key = 0; simple Euler
 key = 1; modified Euler (predictor-corrector)
 key = 2; 4-th order Runge-Kutta

 Alex Godunov: Last revision March 2007
*/
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <string>

using namespace std;

/* function prototypes */
double f1(double, double, double);
double f2(double, double, double);
double euler2d(double(*)(double, double, double),
               double(*)(double, double, double),
               double, double, double, double,
               double&, double&);
double euler2m(double(*)(double, double, double),
               double(*)(double, double, double),
               double, double, double, double,
               double&, double&);
double rk4_2nd(double(*)(double, double, double),
               double(*)(double, double, double),
               double, double, double, double,
               double&, double&);

int main()
{
    double ti, xi, vi, tf, xf, vf, dt, tmax;
    double energy;
    int key;
    const string method[3] = {"simple Euler","modified Euler","4th order Runge-Kutta"};

/* output: file and formats */
    ofstream file;
    file.open ("ode02.dat");
    file.precision(6);
    file.setf(ios::fixed | ios::showpoint);
    cout.precision(6);
    cout.setf(ios::fixed | ios::showpoint);

/* initial information */
    key =  2;             // select a method (key = 0, 1, 2)
    ti = 0.0;             // initial value for variable
    xi = 0.0;             // initial value for function x(t)
    vi = 1.0;             // initial
    dt = 0.1;             // step size for integration
    tmax = 12.0;          // integrate from ti till tmax
/* end of initial information */

    file << setw(30) << method[key] << endl;
    file << setw(12) << "t" << setw(12) <<"x"<< setw(12) << "x'" << endl;
    file << setw(12) << ti << setw(12) << xi << setw(12) << vi   << endl;

/* integration of ODE */
    while (ti <= tmax)
    {
        tf = ti + dt;
        if(key == 0) euler2d(f1,f2,ti,xi,vi,tf,xf,vf);
        if(key == 1) euler2m(f1,f2,ti,xi,vi,tf,xf,vf);
        if(key == 2) rk4_2nd(f1,f2,ti,xi,vi,tf,xf,vf);

        file << setw(12) << tf << setw(12) << xf  << setw(12) << vf << endl;
        ti = tf;
        xi = xf;
        vi = vf;
    }
    system ("pause");
    return 0;
}

/*
  Definition of the x'(t) = f1(t,x,x') = x' by the definition
*/
    double f1(double t, double x, double v)
{
    double d1x;
    d1x = v;
    return d1x;
}
/*
 *  Definition of the x"(t) = f2(t,x,x')
*/
    double f2(double t, double x, double v)
{
    double d2x;
    d2x = (-1.0)*x;  //simple harmonic oscillator
    return d2x;
}

/*-----------------------------------------------------------------------
 Program to solve 1st order Ordinary Differential Equations
 x'(t) = d1x(t,x)
 x"(t) = d2x(t,x)
 method:    simple Euler method
 input ...
 d1x(t,x) - first order derivative: supplied by a user
 d2x(t,x) - second order derivative: supplied by a user
 ti  - initial value for an independent variable (t)
 xi  - initial value for x(t)
 vi  - initial value for x'(t)
 tf  - find solution for this point t
 output ...
 xf and vf - solutions at point tf, i.e. x(tf) and x'(tf)
-----------------------------------------------------------------------*/
double euler2d(double(*d1x)(double, double, double),
               double(*d2x)(double, double, double),
               double ti, double xi, double vi, double tf,
               double& xf, double& vf)
{
    xf = xi + d1x(ti,xi,vi)*(tf-ti);
    vf = vi + d2x(ti,xi,vi)*(tf-ti);
   return 0.0;
}

/*-----------------------------------------------------------------------
 Program to solve 1st order Ordinary Differential Equations
 x'(t) = d1x(t,x)
 x"(t) = d2x(t,x)
 method: modified Euler method (predictor-corrector)
 input ...
 d1x(t,x) - first order derivative: supplied by a user
 d2x(t,x) - second order derivative: supplied by a user
 ti  - initial value for an independent variable (t)
 xi  - initial value for x(t)
 vi  - initial value for x'(t)
 tf  - find solution for this point t
 output ...
 xf and vf - solutions at point tf, i.e. x(tf) and x'(tf)
-----------------------------------------------------------------------*/
double euler2m(double(*d1x)(double, double, double),
               double(*d2x)(double, double, double),
               double ti, double xi, double vi, double tf,
               double& xf, double& vf)
{
    xf = xi + d1x(ti,xi,vi)*(tf-ti);
    vf = vi + d2x(ti,xi,vi)*(tf-ti);
/* correction */
    xf = xi + (d1x(ti,xi,vi)+d1x(ti,xf,vf))*0.5*(tf-ti);
    vf = vi + (d2x(ti,xi,vi)+d2x(ti,xf,vf))*0.5*(tf-ti);
   return 0.0;
}

/*-----------------------------------------------------------------------
 Program to solve 1st order Ordinary Differential Equations
 x'(t) = d1x(t,x)
 x"(t) = d2x(t,x)
 method: 4th-order Runge-Kutta method
 input ...
 d1x(t,x) - first order derivative: supplied by a user
 d2x(t,x) - second order derivative: supplied by a user
 ti  - initial value for an independent variable (t)
 xi  - initial value for x(t)
 vi  - initial value for x'(t)
 tf  - find solution for this point t
 output ...
 xf and vf - solutions at point tf, i.e. x(tf) and x'(tf)
-----------------------------------------------------------------------*/
double rk4_2nd(double(*d1x)(double, double, double),
               double(*d2x)(double, double, double),
               double ti, double xi, double vi, double tf,
               double& xf, double& vf)
{
      double h,t,k1x,k2x,k3x,k4x,k1v,k2v,k3v,k4v;

      h = tf-ti;
      t = ti;

      k1x = h*d1x(t,xi,vi);
      k1v = h*d2x(t,xi,vi);

      k2x = h*d1x(t+h/2.0,xi+k1x/2.0,vi+k1v/2.0);
      k2v = h*d2x(t+h/2.0,xi+k1x/2.0,vi+k1v/2.0);

      k3x = h*d1x(t+h/2.0,xi+k2x/2.0,vi+k2v/2.0);
      k3v = h*d2x(t+h/2.0,xi+k2x/2.0,vi+k2v/2.0);

      k4x = h*d1x(t+h,xi+k3x,vi+k3v);
      k4v = h*d2x(t+h,xi+k3x,vi+k3v);

      xf = xi + (k1x + 2.0*(k2x+k3x) + k4x)/6.0;
      vf = vi + (k1v + 2.0*(k2v+k3v) + k4v)/6.0;

      return 0.0;
}

/* test output for simple harmonic motion
   x"(t) = -x
   x(0)  = 0.0
   x'(0) = 1.0

                  simple Euler                       modified Euler                4th order Runge-Kutta
           t           x          x'            t           x          x'            t           x          x'
    0.000000    0.000000    1.000000     0.000000    0.000000    1.000000     0.000000    0.000000    1.000000
    0.100000    0.100000    1.000000     0.100000    0.100000    0.995000     0.100000    0.099833    0.995004
    0.200000    0.200000    0.990000     0.200000    0.199000    0.980050     0.200000    0.198669    0.980067
    0.300000    0.299000    0.970000     0.300000    0.296010    0.955299     0.300000    0.295520    0.955337
    0.400000    0.396000    0.940100     0.400000    0.390060    0.920996     0.400000    0.389418    0.921061
    0.500000    0.490010    0.900500     0.500000    0.480209    0.877483     0.500000    0.479425    0.877583
    0.600000    0.580060    0.851499     0.600000    0.565556    0.825194     0.600000    0.564642    0.825336
    0.700000    0.665210    0.793493     0.700000    0.645248    0.764654     0.700000    0.644217    0.764843
    0.800000    0.744559    0.726972     0.800000    0.718487    0.696467     0.800000    0.717356    0.696707
    0.900000    0.817256    0.652516     0.900000    0.784542    0.621316     0.900000    0.783326    0.621611
    1.000000    0.882508    0.570790     1.000000    0.842750    0.539951     1.000000    0.841470    0.540303
    1.100000    0.939587    0.482540     1.100000    0.892532    0.453187     1.100000    0.891207    0.453597
    1.200000    0.987841    0.388581     1.200000    0.933388    0.361891     1.200000    0.932039    0.362359
    1.300000    1.026699    0.289797     1.300000    0.964910    0.266976     1.300000    0.963558    0.267500
    1.400000    1.055679    0.187127     1.400000    0.986783    0.169392     1.400000    0.985449    0.169968
    1.500000    1.074391    0.081559     1.500000    0.998788    0.070113     1.500000    0.997495    0.070738
    1.600000    1.082547   -0.025880     1.600000    1.000806   -0.029867     1.600000    0.999574   -0.029198
    1.700000    1.079959   -0.134135     1.700000    0.992815   -0.129548     1.700000    0.991665   -0.128843
    1.800000    1.066546   -0.242131     1.800000    0.974896   -0.227933     1.800000    0.973848   -0.227201
    1.900000    1.042333   -0.348785     1.900000    0.947228   -0.324039     1.900000    0.946300   -0.323288
    2.000000    1.007454   -0.453019     2.000000    0.910088   -0.416905     2.000000    0.909298   -0.416145
    2.100000    0.962152   -0.553764     2.100000    0.863847   -0.505602     2.100000    0.863210   -0.504845
    2.200000    0.906776   -0.649979     2.200000    0.808968   -0.589243     2.200000    0.808497   -0.588500
    2.300000    0.841778   -0.740657     2.300000    0.745999   -0.666991     2.300000    0.745706   -0.666274
    2.400000    0.767712   -0.824835     2.400000    0.675570   -0.738070     2.400000    0.675465   -0.737392
    2.500000    0.685229   -0.901606     2.500000    0.598385   -0.801767     2.500000    0.598474   -0.801142
    2.600000    0.595068   -0.970129     2.600000    0.515216   -0.857447     2.600000    0.515503   -0.856887
    2.700000    0.498055   -1.029636     2.700000    0.426895   -0.904553     2.700000    0.427382   -0.904071
    2.800000    0.395092   -1.079441     2.800000    0.334306   -0.942613     2.800000    0.334990   -0.942221
    2.900000    0.287148   -1.118950     2.900000    0.238373   -0.971247     2.900000    0.239252   -0.970957
    3.000000    0.175253   -1.147665     3.000000    0.140056   -0.990168     3.000000    0.141122   -0.989992
    3.100000    0.060486   -1.165190     3.100000    0.040339   -0.999188     3.100000    0.041583   -0.999135
    3.200000   -0.056033   -1.171239     3.200000   -0.059781   -0.998216     3.200000   -0.058371   -0.998295
    3.300000   -0.173157   -1.165636     3.300000   -0.159304   -0.987262     3.300000   -0.157743   -0.987480
    3.400000   -0.289720   -1.148320     3.400000   -0.257234   -0.966435     3.400000   -0.255538   -0.966799
    3.500000   -0.404552   -1.119348     3.500000   -0.352591   -0.935944     3.500000   -0.350780   -0.936457
    3.600000   -0.516487   -1.078893     3.600000   -0.444422   -0.896093     3.600000   -0.442518   -0.896760
    3.700000   -0.624376   -1.027244     3.700000   -0.531810   -0.847281     3.700000   -0.529833   -0.848101
    3.800000   -0.727101   -0.964806     3.800000   -0.613879   -0.789997     3.800000   -0.611855   -0.790969
    3.900000   -0.823582   -0.892096     3.900000   -0.689809   -0.724812     3.900000   -0.687764   -0.725934
    4.000000   -0.912791   -0.809738     4.000000   -0.758841   -0.652380     4.000000   -0.756800   -0.653646
    4.100000   -0.993765   -0.718459     4.100000   -0.820285   -0.573424     4.100000   -0.818275   -0.574827
    4.200000   -1.065611   -0.619083     4.200000   -0.873526   -0.488733     4.200000   -0.871574   -0.490264
    4.300000   -1.127519   -0.512522     4.300000   -0.918032   -0.399155     4.300000   -0.916164   -0.400802
    4.400000   -1.178771   -0.399770     4.400000   -0.953357   -0.305586     4.400000   -0.951601   -0.307336
    4.500000   -1.218748   -0.281892     4.500000   -0.979149   -0.208960     4.500000   -0.977529   -0.210799
    4.600000   -1.246938   -0.160018     4.600000   -0.995149   -0.110246     4.600000   -0.993690   -0.112156
    4.700000   -1.262939   -0.035324     4.700000   -1.001198   -0.010428     4.700000   -0.999923   -0.012393
    4.800000   -1.266472    0.090970     4.800000   -0.997235    0.089493     4.800000   -0.996165    0.087495
    4.900000   -1.257375    0.217617     4.900000   -0.983299    0.188520     4.900000   -0.982453    0.186508
    5.000000   -1.235613    0.343355     5.000000   -0.959531    0.285662     5.000000   -0.958925    0.283658
    5.100000   -1.201278    0.466916     5.100000   -0.926167    0.379946     5.100000   -0.925816    0.377974
    5.200000   -1.154586    0.587044     5.200000   -0.883541    0.470432     5.200000   -0.883456    0.468513
    5.300000   -1.095882    0.702502     5.300000   -0.832081    0.556213     5.300000   -0.832270    0.554370
    5.400000   -1.025631    0.812090     5.400000   -0.772299    0.636432     5.400000   -0.772767    0.634689
    5.500000   -0.944422    0.914654     5.500000   -0.704794    0.710287     5.500000   -0.705543    0.708666
    5.600000   -0.852957    1.009096     5.600000   -0.630242    0.777038     5.600000   -0.631270    0.775563
    5.700000   -0.752047    1.094391     5.700000   -0.549386    0.836020     5.700000   -0.550689    0.834710
    5.800000   -0.642608    1.169596     5.800000   -0.463038    0.886641     5.800000   -0.464606    0.885517
    5.900000   -0.525649    1.233857     5.900000   -0.372058    0.928396     5.900000   -0.373881    0.927476
    6.000000   -0.402263    1.286422     6.000000   -0.277358    0.960867     6.000000   -0.279420    0.960168
    6.100000   -0.273621    1.326648     6.100000   -0.179885    0.983729     6.100000   -0.182167    0.983267
    6.200000   -0.140956    1.354010     6.200000   -0.080613    0.996754     6.200000   -0.083094    0.996541
    6.300000   -0.005555    1.368106     6.300000    0.019466    0.999811     6.300000    0.016809    0.999858
    6.400000    0.131256    1.368661     6.400000    0.119350    0.992870     6.400000    0.116544    0.993185
    6.500000    0.268122    1.355536     6.500000    0.218040    0.976001     6.500000    0.215115    0.976588
    6.600000    0.403675    1.328724     6.600000    0.314550    0.949371     6.600000    0.311536    0.950234
    6.700000    0.536548    1.288356     6.700000    0.407914    0.913248     6.700000    0.404845    0.914385
    6.800000    0.665383    1.234701     6.800000    0.497199    0.867992     6.800000    0.494108    0.869400
    6.900000    0.788854    1.168163     6.900000    0.581513    0.814057     6.900000    0.578435    0.815728
    7.000000    0.905670    1.089278     7.000000    0.660011    0.751981     7.000000    0.656982    0.753906
    7.100000    1.014598    0.998711     7.100000    0.731909    0.682385     7.100000    0.728965    0.684551
    7.200000    1.114469    0.897251     7.200000    0.796488    0.605965     7.200000    0.793664    0.608356
    7.300000    1.204194    0.785804     7.300000    0.853102    0.523485     7.300000    0.850433    0.526082
    7.400000    1.282774    0.665385     7.400000    0.901185    0.435771     7.400000    0.898705    0.438553
    7.500000    1.349313    0.537107     7.500000    0.940256    0.343699     7.500000    0.937997    0.346641
    7.600000    1.403023    0.402176     7.600000    0.969925    0.248190     7.600000    0.967918    0.251266
    7.700000    1.443241    0.261874     7.700000    0.989894    0.150199     7.700000    0.988167    0.153380
    7.800000    1.469428    0.117549     7.800000    0.999964    0.050706     7.800000    0.998542    0.053962
    7.900000    1.481183   -0.029393     7.900000    1.000035   -0.049294     7.900000    0.998941   -0.045996
    8.000000    1.478244   -0.177512     8.000000    0.990106   -0.148801     8.000000    0.989359   -0.145493
    8.100000    1.460493   -0.325336     8.100000    0.970275   -0.246820     8.100000    0.969891   -0.243537
    8.200000    1.427959   -0.471385     8.200000    0.940742   -0.342371     8.200000    0.940732   -0.339148
    8.300000    1.380821   -0.614181     8.300000    0.901801   -0.434498     8.300000    0.902174   -0.431370
    8.400000    1.319402   -0.752263     8.400000    0.853842   -0.522280     8.400000    0.854602   -0.519282
    8.500000    1.244176   -0.884204     8.500000    0.797345   -0.604839     8.500000    0.798491   -0.602006
    8.600000    1.155756   -1.008621     8.600000    0.732874   -0.681350     8.600000    0.734402   -0.678714
    8.700000    1.054894   -1.124197     8.700000    0.661075   -0.751048     8.700000    0.662974   -0.748641
    8.800000    0.942474   -1.229686     8.800000    0.582665   -0.813235     8.800000    0.584923   -0.811088
    8.900000    0.819505   -1.323933     8.900000    0.498428   -0.867289     8.900000    0.501027   -0.865431
    9.000000    0.687112   -1.405884     9.000000    0.409207   -0.912671     9.000000    0.412125   -0.911127
    9.100000    0.546524   -1.474595     9.100000    0.315894   -0.948926     9.100000    0.319105   -0.947719
    9.200000    0.399064   -1.529248     9.200000    0.219422   -0.975692     9.200000    0.222897   -0.974841
    9.300000    0.246139   -1.569154     9.300000    0.120755   -0.992701     9.300000    0.124462   -0.992224
    9.400000    0.089224   -1.593768     9.400000    0.020881   -0.999783     9.400000    0.024783   -0.999692
    9.500000   -0.070153   -1.602690     9.500000   -0.079201   -0.996867     9.500000   -0.075143   -0.997172
    9.600000   -0.230422   -1.595675     9.600000   -0.178492   -0.983982     9.600000   -0.174319   -0.984689
    9.700000   -0.389989   -1.572633     9.700000   -0.275998   -0.961257     9.700000   -0.271753   -0.962366
    9.800000   -0.547253   -1.533634     9.800000   -0.370743   -0.928920     9.800000   -0.366471   -0.930429
    9.900000   -0.700616   -1.478909     9.900000   -0.461782   -0.887294     9.900000   -0.457528   -0.889194
   10.000000   -0.848507   -1.408847    10.000000   -0.548202   -0.836795    10.000000   -0.544014   -0.839075
   10.100000   -0.989392   -1.323996    10.100000   -0.629141   -0.777928    10.100000   -0.625064   -0.780573
   10.200000   -1.121791   -1.225057    10.200000   -0.703788   -0.711281    10.200000   -0.699868   -0.714271
   10.300000   -1.244297   -1.112878    10.300000   -0.771397   -0.637522    10.300000   -0.767680   -0.640833
   10.400000   -1.355585   -0.988448    10.400000   -0.831292   -0.557388    10.400000   -0.827821   -0.560991
   10.500000   -1.454430   -0.852890    10.500000   -0.882874   -0.471679    10.500000   -0.879691   -0.475544
   10.600000   -1.539719   -0.707447    10.600000   -0.925628   -0.381254    10.600000   -0.922771   -0.385346
   10.700000   -1.610463   -0.553475    10.700000   -0.959125   -0.287017    10.700000   -0.956632   -0.291298
   10.800000   -1.665811   -0.392429    10.800000   -0.983031   -0.189909    10.800000   -0.980934   -0.194339
   10.900000   -1.705054   -0.225848    10.900000   -0.997107   -0.090902    10.900000   -0.995435   -0.095438
   11.000000   -1.727638   -0.055342    11.000000   -1.001212    0.009014    11.000000   -0.999989    0.004417
   11.100000   -1.733173    0.117422    11.100000   -0.995304    0.108840    11.100000   -0.994553    0.104227
   11.200000   -1.721430    0.290739    11.200000   -0.979444    0.207577    11.200000   -0.979179    0.202996
   11.300000   -1.692357    0.462882    11.300000   -0.953789    0.304239    11.300000   -0.954021    0.299736
   11.400000   -1.646068    0.632118    11.400000   -0.918596    0.397858    11.400000   -0.919332    0.393482
   11.500000   -1.582857    0.796724    11.500000   -0.874217    0.487499    11.500000   -0.875456    0.483296
   11.600000   -1.503184    0.955010    11.600000   -0.821096    0.572264    11.600000   -0.822833    0.568281
   11.700000   -1.407683    1.105328    11.700000   -0.759764    0.651307    11.700000   -0.761989    0.647588
   11.800000   -1.297150    1.246097    11.800000   -0.690835    0.723837    11.800000   -0.693532    0.720425
   11.900000   -1.172541    1.375812    11.900000   -0.614997    0.789129    11.900000   -0.618144    0.786064
   12.000000   -1.034959    1.493066    12.000000   -0.533009    0.846529    12.000000   -0.536581    0.843848
   12.100000   -0.885653    1.596562    12.100000   -0.445691    0.895464    12.100000   -0.449656    0.893201

*/