//
// Program: v_doubler.cxx
//
// Description: Compute the voltage output for an unsymetrical voltage
//              doubler circuit.
//              The differential equations are solved using a
//              fourth order Runge Kutta method.
//
//              The program compiles with the vector.o library
//              g++ v_doubler.cxx vector.o
//
// Author: A. Garfagnini
//
// Date: 19 October 2013
//

#include <iostream>
#include <fstream>
#include <cmath>
#include <libgen.h>

using namespace std;

#include "vector.hpp"
using namespace cpl;

// Global variables
double h = 0.01;

void usage(const char * pname)
{
    cerr << "usage: " << pname << " [-o file_name]" << endl;
}

// - 
// The circuit source Voltage, a alternating voltage function
// -
double v_source(double t)
{
    const double nu = 50.0; // 50 Hz
    return ( 5.0 * sin(2.0 * M_PI * nu * t) );
}

// -
// Evaluate the derivative of the input voltage function
// using a centered first derivative
// -
double dvdt_source(double t, double h)
{
    double t_minus = t - h;
    double t_plus  = t + h;
    return ( (v_source(t_plus) - v_source(t_minus))/(2.0 * h) );
}

// -
// Here are coded the Circuit differential equations
// -
Vector f(Vector u)
{
    const double Is = 1.E-8;   // 10 pA;
    const double C1 = 1.E-3;   // 1 mF
    const double C2 = 1.E-3;   // 1 mF
    const double R  = 1.E4;    // 10 kOhm
    const double VT = 0.0256;  // 25 mV

    double t = u[0];
    double x = u[1];
    double y = u[2];

    double ID1 = Is*(exp(-x/VT) - 1.0);
    double ID2 = Is*(exp((x-y)/VT) - 1.0);

    Vector fv(3);
    fv[0] = 1;
    fv[1] = dvdt_source(t, h) + ID1/C1 - ID2/C1;
    fv[2] = ID2/C2 - y/(R*C2);

    return fv;
}

// -
// Perform a single evolution step using a IV ordr Runge-Kutta method
// -
void RK4Step(Vector & u, double h)
{
    Vector k1 = h * f(u);
    Vector k2 = h * f(u + 0.5*k1);
    Vector k3 = h * f(u + 0.5*k2);
    Vector k4 = h * f(u + k3);

    u += (k1 + 2*k2 + 2*k3 + k4)/6.0;
}

void print_step(ostream & out, Vector & u, bool header=false)
{
    if (header) {
        out << "# Time [ms]   V1 [V]   V2 [V]\n";
    }
    double t = u[0];
    double v_s = v_source(t);

    // - time in ms
    out << t*1000 << " " << v_s << " " << u[1] << " " << u[2] << endl;
}

int main(int argc, char * argv[])
{
   int opt;
   ofstream out_file;

    // Handle command line options
    while((opt = getopt(argc, argv, "o:h?")) != -1) {

        switch (opt) {

        case 'o':
            out_file.open(optarg);
            break;

        case 'h':
        case '?':
        default:
            usage(basename(argv[0]));
            return 1;
        }
    }

    ostream & myout = (out_file.is_open()) ? out_file : cout;

    //
    // We decide to evolve the equation for an integer number
    // of voltage source periods, t = 20 ms
    cout << "Time integration, number of periods [1T = 20 ms]: ";
    double Np;
    cin >> Np;

    cout << "Single time step :";
    cin >> h;

    double t_f = Np * 0.02;
    myout << "# Final integration time (starting from 0s): "
         << t_f << " s\n";

    Vector u(3);

    // Set the initial values: t=0, x = 0, y=0;
    u[0] = u[1] = u[2] = 0;

    // - Write the initial values
    print_step(myout, u, true);

    while (u[0] < t_f) {

        RK4Step(u, h);
        print_step(myout, u);

    }

}