#include "particle.hpp"

Particle::Particle()
// Constructor sets default values
{
   t=0.0;
   x=0.0;
   v=0.0;
   m=1.0;
}


void Particle::SetPos(double x_)
{
   x=x_;
}

double Particle::GetPos()
{
  return x;
}

void Particle::SetVel(double v_)
{
   v=v_;
}

double Particle::GetVel()
{
   return v;
}

void Particle::SetMass(double m_)
{
   m=m_;
}

double Particle::GetMass()
{
   return m;
}

double Particle::GetTime()
{
   return t;
}

void Particle::IntegrateEL(double dt, double (*f)(double,double,double))
// Euler method  f is a function that evalues the force for a given x,v,t
{
    double xnext=x+v*dt;
    double vnext=v+f(x,v,t)*dt/m;

    x=xnext;
    v=vnext;
    t+=dt;
}


void Particle::IntegrateMP(double dt, double (*f)(double,double,double))
// Midpoint method  f is a function that evalues the force for a given x,v,t
{
    double xm=x+v*dt/2;
    double vm=v+f(x,v,t)*dt/(2*m);

    double xnext = x + v*dt + f(x,v,t)*dt*dt/(2*m);
    double vnext = v + f(xm, vm, t+dt/2)*dt/m;

    x=xnext;
    v=vnext;
    t+=dt;
}

void Particle::IntegrateRK4(double dt, double (*f)(double,double,double))
// Runge-Kutta RK4
{
// derivatives at t
   double f1x=v;
   double f1v=f(x,v,t)/m;

// First estimate of derivatives at t+=dt/2
   double f2x=v+dt/2*f1x;
   double f2v=f(x+dt/2*f1x,v+dt/2*f1v,t+dt/2)/m;

// Corrected estimates at t+=dt/2
   double f3x=v+dt/2*f2x;
   double f3v=f(x+dt/2*f2x,v+dt/2*f2v,t+dt/2)/m;

// Final estimates at t+=dt
   double f4x=v+dt*f3x;
   double f4v=f(x+dt*f3x,v+dt*f3v,t+dt)/m;

   double xnext=x+dt*(f1x+2.0*f2x+2.0*f3x+f4x)/6.0;
   double vnext=v+dt*(f1v+2.0*f2v+2.0*f3v+f4v)/6.0;

   x=xnext;
   v=vnext;
   t+=dt;
}