//
// Program: osc_ec.cxx
//
// Description: Solve the the armonic oscillator equation using the Euler
//              and Euler-Cromer numerical methods
//
// Author: A. Garfagnini
//
// Date: 9 October 2012
//

#include <iostream>
#include <fstream>
#include <unistd.h>
#include <libgen.h>
#include <stdlib.h>

enum algorithm {EULER, EULER_CROMER};

struct oscillator {
    double x;
    double v;
    double t;
    double energy;
};

using namespace std;

void usage(const char * pname)
{
    cerr << "Usage: " << pname << " [-o file_name]" << endl;
    cerr << "\n  Solve the simple armonic oascillator equation using the"
         << "\n  Euler and Euler-Cromer schemes" << endl;
    cerr << "Options:\n";
    cerr << " -o file_name : save output in a user definable file name\n";
    cerr << " -e           : solve using the Euler scheme (default)\n";
    cerr << " -c           : solve using the Euler-Cromer scheme\n";
}

void print_scheme(ostream & out, algorithm scheme)
{
    switch (scheme)
    {
        case EULER :
            out << "# Scheme: EULER" << endl;
            break;

        case EULER_CROMER :
            out << "# Scheme: EULER CROMER" << endl;
            break;
    }
}

void print_step(ostream & out, oscillator & z, bool print_header=false)
{
    if (print_header) {
        out << "# Time  x  v  Energy" << endl;
    }
    out << z.t << " " << z.x << " " << z.v << " " << z.energy << endl;
}

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

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

        switch (opt) {

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

        case 'e':
            scheme = EULER;
            break;

        case 'c':
            scheme = EULER_CROMER;
            break;

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

    cout << "Time step h: ";
    double h;
    cin >> h;
    if (out_file.is_open()) {
        out_file << "# Time step h: " << h << endl;
    }

    cout << "T_final: ";
    double t_f;
    cin >> t_f;
    if (out_file.is_open()) {
        out_file << "# T_final: " << t_f << endl;
    }

    if (out_file.is_open()) {
        print_scheme(out_file, scheme);
    } else {
        print_scheme(cout, scheme);
    }

    // - Init the system: set the starting conditions
    oscillator z;
    z.x = 1.0;
    z.v = 0.0;
    z.t = 0.0;
    z.energy = z.x * z.x + z.v * z.v;


    if (out_file.is_open()) {
        print_step(out_file, z, true);
    } else {
        print_step(cout, z, true);
    }

    double v_temp;
    while (z.t < t_f) {

        z.t += h;

        switch (scheme)
        {
            case EULER :
                v_temp = z.v;
                z.v -=  h * z.x;
                z.x +=  h * v_temp;
                break;

            case EULER_CROMER :
                z.v -=  h * z.x;
                z.x +=  h * z.v;
                break;
        }

        // Evaluate the energy of the oscillator
        z.energy = z.x * z.x  + z.v * z.v;

        if (out_file.is_open()) {
            print_step(out_file, z);
        } else {
            print_step(cout, z);
        }

    }

    if (out_file.is_open()) {
        out_file.close();
    }

    return 0;
}