/*
    tops.h
  
    Header file for the top classes for free rotation of general rotors
  
    Ramses van Zon, May 9, 2007
  
    TopSpherical : for spherical rotors,  for which Ix = Iy = Iz. See [1,2].
    TopProlate   : for prolate rotors,    for which Ix < Iy = Iz. See [1,2].
    TopOblate    : for spherical rotors,  for which Ix = Iy < Iz. See [1,2].
    TopAsymmetric: for asymmetric rotors, for which Ix < Iy < Iz. See [2].
    TopRecursive : also for asymmetric rotors, but computed using
                   a recursive scheme. See [3].
  
    [1] Lisandro Hernandez de la Pena, Ramses van Zon, Jeremy Schofield
        and Sheldon B. Opps,
        "Discontinuous molecular dynamics for semiflexible and rigid bodies"
        Journal of Chemical Physics 126, 074105 (2007).
  
    [2] Ramses van Zon and Jeremy Schofield,
        "Numerical implementation of the exact dynamics of free rigid bodies"
        Journal of Computational Physics, doi:10.1016/j.jcp.2006.11.019 (2007).
  
    [3] Ramses van Zon and Jeremy Schofield,
        "Symplectic algorithms for simulations of rigid-body systems using
         the exact solution of free motion", 
        Physical Review E 75, 056701 (2007).
*/
#ifndef _TOPSH_
#define _TOPSH_

#include "vecmat3.h"

// abstract base class Top, derivatives of which compute torque-free rotation

class Top 
{
 public:
  // Initialize
  virtual void Initialization( const Vector& omega, const Matrix& A ) = 0;
  
  // Find the values of the angular momenta and attitude matrix at time t
  virtual void Evolution( double t, Vector& omega, Matrix& A ) = 0;
  
  // Initialization and evolution in one
  void Propagation( double dt, Vector& omega, Matrix& A )
  {
    Initialization( omega, A );
    Evolution( dt, omega, A );
  }  

  virtual ~Top() {}
};


// class TopSpherical computes the torque-free rotation of a spherical rotor

class TopSpherical: public Top
{
 public:
  // Initialize the moments of inertia
  explicit TopSpherical( double _I );

  // To initialize the Top class at time zero, precompute things etc:
  void Initialization( const Vector& omega, const Matrix& A );

  // To find the values of the angular momenta and attitude matrix at time t:
  void Evolution( double t, Vector& omega, Matrix& A );

 private:
  double I;
  Vector omegab;
  Matrix A0;
};


// class TopProlate computes the torque-free rotation of a prolate rotor

class TopProlate: public Top
{
 public:
  // Initialize the moments of inertia
  TopProlate( double _Ix, double _Iz );

  // To initialize the Top class at time zero, precompute things etc:
  void Initialization( const Vector& omega, const Matrix& A );

  // To find the values of the angular momenta and attitude matrix at time t:
  void Evolution( double t, Vector& omega, Matrix& A );

 private:
  double I1;
  double I3;
  Vector omegab0;
  Vector LboverI1;
  Matrix A0;
  double omegap;
};


// class TopOblate computes the torque-free rotation of a oblate rotor

class TopOblate: public Top
{
 public:
  // Initialize the moments of inertia
  TopOblate( double _Ix, double _Iz );

  // To initialize the Top class at time zero, precompute things etc:
  void Initialization( const Vector& omega, const Matrix& A );

  // To find the values of the angular momenta and attitude matrix at time t:
  void Evolution( double t, Vector& omega, Matrix& A );

 private:
  double I1;
  double I3;
  Vector omegab0;
  Vector LboverI1;
  Matrix A0;
  double omegap;
};


// class TopAsymmetric computes the torque-free rotation of a general rotor
//
// Based on:
//    Ramses van Zon and Jeremy Schofield,
//    "Numerical implementation of the exact dynamics of free rigid bodies"
//    J. Comput. Phys. (2007), in press, doi:10.1016/j.jcp.2006.11.019 .

class TopAsymmetric: public Top
{
 public:
  TopAsymmetric( double _Ix, double _Iy, double _Iz );

  // To initialize the Top class at time zero, precompute things etc:
  void Initialization( const Vector& omega, const Matrix& A );

  // To find the values of the angular momenta and attitude matrix at time t:
  void Evolution( double t, Vector& omega, Matrix& A );

  ~TopAsymmetric();
  
 private:
  double Ix;
  double Iy;
  double Iz;
  int orderFlag;
  int numTerms;
  double I1;
  double I2;
  double I3;
  double omega1m;
  double omega2m;
  double omega3m;
  double omegap;
  double relfreq;
  double epsilon;
  double A1;
  double A2;
  double L;
  double m;
  Matrix B;
  int maxNumTerms;
  double* Cqr;
  double* Cqi;
  
};


// class TopRecur also computes the torque-free rotation of a general 
// asymmetric rotor using a recursive formula  for the 'last angle'
//
// Based on:
//    Ramses van Zon and Jeremy Schofield,
//   "Symplectic algorithms for simulations of rigid-body systems using
//    the exact solution of free motion"
//    Phys. Rev. E 75 (2007) 056701

class TopRecur: public Top
{
 public:
  // Initialize the moments of inertia
  TopRecur( double _Ix, double _Iy, double _Iz );
 
  // Initialize the Top class at time zero, precompute things etc:
  void Initialization( const Vector& omega, const Matrix& A );

  // Find the values of the angular momenta and attitude matrix at time t:
  void Evolution( double t, Vector& omega, Matrix& A );

  ~TopRecur();
  
 private:
  int orderFlag;
  double Ix;
  double Iy;
  double Iz;
  double **coeff;
  double *term;
  double initialomega1;
  double initialomega2;
  double initialomega3;
  double I1;             // internal moments of inertia
  double I2;
  double I3;
  double omega1m;
  double omega2m;
  double omega3m;
  double omegap;
  double snepsilon;
  double cnepsilon;
  double dnepsilon;
  double twiceE;
  double L;
  double m; 
  Matrix B;
  double X[2];    // two because of two ordering cases
  double Y[2];
  double*** q[2]; // two three-index arrays for the recursion
  
  void firstOmega(double omega1i, double omega2i, double omega3i,
                  double omega1f, double omega2f, double omega3f,
                  double &Omegai, double &Omegaf);
  void nextOmega(int j, double &Omegai, double &Omegaf);
};

#endif