/*   []----------------------------------------------[]
      |              ClassicalDynamics.cpp           |
     []----------------------------------------------[]
      |                                              |
      |   AUTHOR:    MFSomers 2005.                  |
      |   USE:       Defines a classical dynamics    |
      |              framework class to be used for  |
      |              different equations of motions  |
      |              integrators...                  |
      |                                              |
     []----------------------------------------------[]
*/

// Copyright (C) 2005 M.F. Somers, Theoretical Chemistry Department, Leiden University
//
// This is free software; you can redistribute it and/or modify it under the terms of
// the GNU Lesser General Public License as published by the Free Software Foundation.
//
// http://www.gnu.org/licenses/lgpl.txt
 
#include "Global.h"
#include "Vector.h"
#include "Constants.h"
#include "ClassicalDynamics.h"

/* ------------------------------------------------ ClassicalDynamics ---- */

int ClassicalDynamics::RunTrajectory( double Start, double Step, int nSteps, ConstraintListPointer Constraints )
{
 double At, IntegratorAt, IntegratorStep;
 int    n, IntegratorNSteps, IntegratorContinueFlag;
 
 if( nSteps <= 0 ) return( 0 );
 
 if( TheParticleListList[ 0 ] == NULL ) return( 0 );         // we need at least two states allways
 if( TheParticleListList[ 1 ] == NULL ) return( 0 );
   
 At = Start;
 n = 0;
 
 while( TrajectoryContinuationCheck( At, Step, n, nSteps ) && ( n < nSteps ) )   // as long as trajectory need to run
 {
  
  IntegratorAt = At;                                     // start the integration for the complete step in a single step
  IntegratorStep = Step;                    
  IntegratorNSteps = 1;
  IntegratorContinueFlag = 0;
  
  do
  {
   
   // do integration step but also check if conservation laws are okay and step size is not too small
   if( ( Integrator( IntegratorStep, IntegratorAt, Constraints ) & IntegratorConservationCheck( IntegratorStep ) ) | IntegratorContinueFlag )   
   {
    CycleTheStates( -1 );                               // cycle the states so we can continue
    IntegratorAt += IntegratorStep;                     // with this step size to finish the
    --IntegratorNSteps;                                 // total step size of the trajectory
   }
   else
   {
    IntegratorStep *= 0.5;                             // we did not acieve integration converegence
    IntegratorNSteps *= 2;                             // so we'll halve the step size and do twice as much steps
    if( ( fabs( IntegratorStep ) < CONST_EPSILON ) ) IntegratorContinueFlag = 1;   // if step becomes to small, flag this
   }
   
  } while( IntegratorNSteps > 0 );    // finish all integration steps of trajectory step 
  
  At += Step;                         // we continue with the next step in the trajectory
  ++n;
 }
  
 return( n );
}

/* ----------------------------------------------------------------------- */