RPK.h

#ifndef __RPK_RPK_H__
#define __RPK_RPK_H__

#include "../Constants/Constants.h"
#include "../InfoStruct/InfoStruct.h"
#include "../Math/Functions.h"
#include "../Particle/Particle.h"
#include "../Units/Units.h"
#include "../Vector/Vector.hpp"
#include "PressureGauge.h"

#include "rapidjson/document.h"
#include "rapidjson/prettywriter.h"

#include <boost/numeric/odeint.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <math.h>
#include <stdio.h>

namespace MDSimulation
{
    const int DIM = 2;

    namespace odeint = boost::numeric::odeint;

    typedef boost::numeric::ublas::vector<double> state_type;
    typedef boost::numeric::ublas::matrix<double> matrix_type;
    typedef odeint::runge_kutta_fehlberg78<state_type> error_stepper_type;
    typedef odeint::controlled_runge_kutta<error_stepper_type> controlled_stepper_type;

    // Numerical solution to Rayleigh-Plesset equation for bubble collapse
    class RPK
    {
    private:

        const double& CurrentTime;

        double kR;       // Interpolation coefficient, aka velocity = R'

        state_type Time;
        state_type PrevRadius;
        state_type NextRadius;

        std::vector<double> Times;
        std::vector<double> Radii;
        std::vector<double> Velocities;

        controlled_stepper_type ControlledSolver;

        // Constants computed using the InfoStruct data. These are kept around for
        // performance reasons.
        double R03;
        double a3;
        double P0;
        double Density;
        double Viscosity;
        double SurfaceTension;
        double SpeedOfSound;
        double Omega;

        // Coefficients for the linear interpolation of the polytropic exponent.
        double A;
        double B;

        double TimeStep;

        // Pressure interpolation coefficients
        double ComputedPressure;    // Pressure used by the RPK function (ideal or pgauge)
        double ComputedDpdt;

        PressureGauge PGauge;

        InfoStruct Info;

        double StepSize;

        FILE* Actual;

        void Init(const std::string& outDir = "");

    public:

        RPK(const double& currentTime, const InfoStruct& info, const std::string& outDir = "");

        ~RPK();

        double NextStep();

        void AddSample(const double time, const double pressure);

        double R(const double time) const;

        double V(const double time) const;

        double VSmart(const double time) const;

        double d2Rdt2() const;

        double NextStepTime();

        void ReadPressureGauge();

        double FindCrossingTime(const Particle& n) const;

        template <typename Writer>
        void Serialize(Writer& writer) const;

        // void Read(const boost::property_tree::ptree& tree);

        state_type operator()(const state_type& x, state_type& dxdt, const double t) const;

        double AcousticPressure(const double time) const;

        double AcousticPressureDerivative(const double time) const;
    };

    template <typename Writer>
    void RPK::Serialize(Writer& writer) const
    {

        writer.String("kR")             , writer.Double(kR);
        writer.String("Time")           , writer.Double(Time);
        writer.String("Time")           , writer.Double(Time);
        writer.String("TimeStep")       , writer.Double(TimeStep);
        writer.String("A")              , writer.Double(A);
        writer.String("B")              , writer.Double(B);
        writer.String("R03")            , writer.Double(R03);
        writer.String("a3")             , writer.Double(a3);
        writer.String("P0")             , writer.Double(P0);
        writer.String("Density")        , writer.Double(Density);
        writer.String("Viscosity")      , writer.Double(Viscosity);
        writer.String("SurfaceTension") , writer.Double(SurfaceTension);
        writer.String("SpeedOfSound")   , writer.Double(SpeedOfSound);
        writer.String("Omega")          , writer.Double(Omega);

        writer.String("PGauge");
        PGauge.Serialize<Writer>(writer);
    }
}

#endif