/* 
   
   Copyright (C) Bernard Piette 
                 University of Durham (UK) 
                 Department of Mathematical Sciences
   
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
   
   You can modify or use any part of the source code but only under the  
   condition that this Copyright notice is kept in the resulting source file.
*/
#include <pde_RK4.h>
#include <iostream>

using namespace std;

// Constructor for RK4 pde 
// Allocated 2 extra arrays     
//******************************
pde_RK4::pde_RK4(double xmin,double xmax,int nx) 
   : pde(xmin,xmax,nx)
{ tmp1_ = new double[nx_];
  tmp2_ = new double[nx_];    
  tmp3_ = new double[nx_];    
}

// Destructor
//***********
pde_RK4::~pde_RK4()
{ delete [] tmp1_;
  delete [] tmp2_;
  delete [] tmp3_;
}

/**********************************************************/
/* Perform a 4th order Runge-Kutta integration step on    */
/* f(n+1) = fn + 1/6(k1 +2 k2 + 2 k3 + k4)                */
/* k1 = dt F(t,fn)                                        */
/* k2 = dt F(t+dt/2,fn+k1/2)                              */
/* k3 = dt F(t+dt/2,fn+k2/2)                              */
/* k4 = dt F(t+dt,fn+k3)                                  */
/* the grid. Update the time in pde.                      */
/* Eqn(Force) : a function which computes the right       */
/*     hand side of the equation.                         */
/*       Force : vector where to store the result.        */
/*       i_    : lattice index                            */   
/* Return : the updated time                              */     
/**********************************************************/
double pde_RK4::integrate_1step()
{ int i;

    // tmp1 is used to compute fn + 1/6(k1 + 2 k2 + 2 k3 + k4) 
    // We initialise it to 0                                
  for(i=0; i < nx_; i++)
  { tmp1_[i] = 0;
  }

    // tmp1 = 0
  sub_RK4_1step(Fcts_,Fcts_,tmp1_,tmp2_,0.5,1.0);

    // tmp1 = k1
    // tmp2 = fn + 1/2 k1
  t_ += 0.5*dt_; 
  sub_RK4_1step(Fcts_,tmp2_,tmp1_,tmp3_,0.5,2.0);

    // tmp1 = fn + 1/6 (k1 + 2 k2 )
    // tmp3 = fn + 1/2 k2 */
  sub_RK4_1step(Fcts_,tmp3_,tmp1_,tmp2_,1.0,2.0);

    // tmp1 = fn + 1/6 (k1 + 2 k2 + 2 k3)
    // tmp3 = fn + k3 */
  t_ += 0.5*dt_; 
  sub_RK4_1step(Fcts_,tmp2_,tmp1_,0,0.0,1.0);

    // tmp1 = (k1 + 2 k2 + 2 k3 + k4) 
    // we update fn => fn + tmp1 / 6  
  for(i=0; i < nx_; i++)
  { Fcts_[i] += tmp1_[i]/6.0;
  }

    // Time has been updated before the last step 
  return(t_);
}

/**************************************************************/
/* Perform one RK intermediate step                           */
/* f : f(n) the fcts at time tn                               */
/* ef : estimated f for this step.                            */
/* mf : modified f : used to add the contribution of each step*/
/* nef : estimated f for the next step, must be computed here */
/* Cnef : coeficiant used to computed nef                     */ 
/* Cmf  : coeficiant used to computed mf                      */
/**************************************************************/
void pde_RK4::sub_RK4_1step(double *f,double *ef,double *mf,double *nef,
                             double Cnef,double Cmf)
{ int i;

  for(i=0; i < nx_; i++)
  { set_neighbour_fcts(ef,i);

        // Computes the force for point i and  
        // stores it in Force                 
    Eqn(Force_); 

        // Update each fct using Force : f(n+1) = f(n) + Force*dt 
    { mf[i] += dt_ * Force_ * Cmf;

      if(nef) // nef will be 0 for the last step 
      { nef[i] = f[i] + dt_ * Force_ * Cnef;
      }
    }
  }
}