#include <fstream>
#include <cmath>
#include "tops.h"

using namespace std;

DOUBLE dt = 0.1; 
DOUBLE finalt = 12.0; 

void testPropagation(Top & rotor, const char* filename)
{
    Vector omega0 ( 1.0, -1.0, 0.5 );
    Matrix A0     ( 1, 0, 0,
                    0, 1, 0,
                    0, 0, 1 );

    cerr << "# Writing " << filename << endl;

    ofstream  f(filename);
    Vector    omega(omega0);
    Matrix    A(A0);

    for (DOUBLE t = 0.0; t < finalt; t += dt) {

        f << t 
          << ' ' << A.row(0)
          << ' ' << A.row(1)
          << ' ' << A.row(2)
          << ' ' << omega 
          << endl;

        rotor.Propagation(dt, omega, A);
    }
}

void testEvolution(Top & rotor, const char* filename)
{
    Vector omega0 ( 1.0, -1.0, 0.5 );
    Matrix A0     ( 1,0,0,
                    0,1,0,
                    0,0,1 );

    cerr << "# Writing " << filename << endl;

    rotor.Initialization( omega0, A0 );
  
    ofstream f( filename );
  
    for ( DOUBLE t = 0.0; t < finalt; t += dt ) {
  
        Vector omega;
        Matrix A;

        rotor.Evolution(t, omega, A);

        f << t 
          << ' ' << A.row(0)
          << ' ' << A.row(1)
          << ' ' << A.row(2)
          << ' ' << omega 
          << endl;
    }
}

int main()
{
    DOUBLE Ix = 1.0, Iy = 1.5, Iz = 2.0;

    TopSpherical  sphere (Ix);          // Ix = Iy = Iz
    TopProlate    prolate(Ix, Iz);      // Ix < Iy = Iz
    TopOblate     oblate (Ix, Iz);      // Ix = Iy < Iz
    TopAsymmetric asymtop(Ix, Iy, Iz);  // Ix < Iy < Iz
    TopRecur      recurse(Ix, Iy, Iz);  // Ix < Iy < Iz, faster

    testPropagation(sphere,  "spherePrp.dat");
    testPropagation(prolate, "prolatePrp.dat");
    testPropagation(oblate,  "oblatePrp.dat");
    testPropagation(asymtop, "asymtopPrp.dat");
    testPropagation(recurse, "recursePrp.dat");

    testEvolution(sphere,  "sphereEvl.dat");
    testEvolution(prolate, "prolateEvl.dat");
    testEvolution(oblate,  "oblateEvl.dat");
    testEvolution(asymtop, "asymtopEvl.dat");
    testEvolution(recurse, "recurseEvl.dat");

    return 0;
}