/*
   tops.h

   Header file for computation torque-free rotation of general rigid bodies
  
   Ramses van Zon, May 9, 2007

   References:

   [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 _CTOPSH_
#define _CTOPSH_

struct Vector 
{
  double x;
  double y;
  double z;
};
typedef struct Vector Vector; 

struct Matrix
{
  double xx;  double xy;  double xz;
  double yx;  double yy;  double yz;
  double zx;  double zy;  double zz;
};
typedef struct Matrix Matrix; 

enum TopType
{
  SphericalTop,   /* Ix = Iy = Iz */
  ProlateTop,     /* Ix < Iy = Iz */
  OblateTop,      /* Ix = Iy < Iz */
  AsymmetricTop,  /* Ix < Iy < Iz */
  RecurTop        /* Ix < Iy < Iz */
};

struct Top
{
  enum TopType type;    /* Spherical  Oblate  Prolate  Asymmetric  Recur */
  double Ix;            /*    X                            X         X   */
  double Iy;            /*                                 X         X   */
  double Iz;            /*                                 X         X   */
  int orderFlag;        /*                                 X         X   */
  int numTerms;         /*                                 X             */
  double I1;            /*               X       X         X         X   */
  double I2;            /*                                 X         X   */
  double I3;            /*               X       X         X         X   */
  double omega1m;       /*                                 X         X   */
  double omega2m;       /*                                 X         X   */
  double omega3m;       /*                                 X         X   */
  double omegap;        /*               X       X         X         X   */
  double relfreq;       /*                                 X             */
  double epsilon;       /*                                 X             */
  double A1;            /*                                 X             */
  double A2;            /*                                 X             */
  double L;             /*                                 X         X   */
  double m;             /*                                 X         X   */
  struct Matrix B;      /*     X         A0      A0        X         X   */
  int maxNumTerms;      /*                                 X             */
  double* Cqr;          /*                                 X             */
  double* Cqi;          /*                                 X             */
  struct Vector omegab0;/*     X         X       X                       */
  struct Vector LboverI1;/*              X       X                       */
  double **coeff ;      /*                                           X   */
  double *term;         /*                                           X   */
  double initialomega1; /*                                           X   */
  double initialomega2; /*                                           X   */
  double initialomega3; /*                                           X   */
  double snepsilon;     /*                                           X   */
  double cnepsilon;     /*                                           X   */
  double dnepsilon;     /*                                           X   */
  double twiceE;        /*                                           X   */
  double*** q[2];       /*                                           X   */
};
typedef struct Top Top;

/* fill appropriate Top structure */
void DefineTop( struct Top * top,
                enum TopType type,
                double Ix, double Iy, double Iz );

/* Initialize */
void Initialization( struct Top * top,
                     struct Vector * omega,
                     struct Matrix * A );

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

/* To find the values of the angular momenta and attitude after at time d */
void Propagation( struct Top * top,
                  double dt,
                  struct Vector * omega,
                  struct Matrix * A );

/* release memory */
void FreeTop( struct Top * top );


#endif