/*   []----------------------------------------------[]
      |              ClassicalStates.cpp             |
     []----------------------------------------------[]
      |                                              |
      |   AUTHOR:    MFSomers 2005.                  |
      |   USE:       Defines a classical dynamics    |
      |              context 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"

/* ----------------------------------------------- ClassicalStates ------- */

ClassicalStates::ClassicalStates( int NrOfStates, unsigned int TheId ) : 
		 TheParticleListList( NrOfStates = ( NrOfStates < 2 ? 2 : NrOfStates ) ), 
		 TheInteractionListList( NrOfStates = ( NrOfStates < 2 ? 2 : NrOfStates ) )
{
 int n;
 
 Id = TheId;
 TheInternalFlag = 0;
 Type = GeneralClassicalStates;
 
 for( int s = 0; s < NrOfStates; ++s )
 {
  TheParticleListList[ s ] = new ParticleList( DEFAULT_START_NR_PARTICLES );
  if( TheParticleListList[ s ] != NULL )
   for( n = 0; n < TheParticleListList[ s ] -> size(); ++n )
    (*(TheParticleListList[ s ]))[ n ] = NULL;
 }
 
 for( int s = 0; s < NrOfStates; ++s )
 {
  TheInteractionListList[ s ] = new InteractionList( DEFAULT_START_NR_INTERACTIONS ); 
  if( TheInteractionListList[ s ] != NULL )
   for( n = 0; n < TheInteractionListList[ s ] -> size(); ++n )
    (*(TheInteractionListList[ s ]))[ n ] = NULL;
 }
  
}

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

ClassicalStates::ClassicalStates( ClassicalStates& D ) : TheParticleListList( D.TheParticleListList.size() ),
                                                         TheInteractionListList( D.TheInteractionListList.size() )
{
 int s, n;
 
 Id = D.Id;
 TheInternalFlag = D.TheInternalFlag;
 Type = D.Type;
  
 // start by copying by particles from D to our own lists... keep in mind they still point
 // to interactions of D rather than to interactions of ourselfs...
 for( s = 0; s < D.TheParticleListList.size(); ++s )
 {
 
  // copy the particle list ...
  if( D.TheParticleListList[ s ] != NULL )
  {
   TheParticleListList[ s ] = new ParticleList( D.TheParticleListList[ s ] -> size() );
   
   if( TheParticleListList[ s ] != NULL )
    for( n = 0; n < D.TheParticleListList[ s ] -> size(); ++n )
     if( (*(D.TheParticleListList[ s ]))[ n ] != NULL )
     {
      (*(TheParticleListList[ s ]))[ n ] = (*(D.TheParticleListList[ s ]))[ n ] -> Clone(); 
      if( (*(TheParticleListList[ s ]))[ n ] != NULL )
       (*(TheParticleListList[ s ]))[ n ] -> DynamicsContext = this;
     } 
     else
      (*(TheParticleListList[ s ]))[ n ] = NULL; 
  }
  else
   TheParticleListList[ s ] = NULL;
   
  // create interaction state lists... but init them to NULL...
  if( D.TheInteractionListList[ s ] != NULL )
  {
   TheInteractionListList[ s ] = new InteractionList( D.TheInteractionListList[ s ] -> size() );
   
   if( TheInteractionListList[ s ] != NULL )
    for( n = 0; n < D.TheInteractionListList[ s ] -> size(); ++n )
     (*(TheInteractionListList[ s ]))[ n ] = NULL;
  } 
  else 
   TheInteractionListList[ s ] = NULL;
}
 
 // we have now copied the all the particles and initialized the interaction lists to NULL.
 // now start including the interactions of D into ourself. we include them so that the
 // particles and interactions will point to eachother within this dynamical context rather
 // than in the ones we are copying...
 
 for( s = 0; s < D.TheInteractionListList.size(); ++s )   // find a non-NULL entry in the D states
  if( D.TheInteractionListList[ s ] != NULL )
   break;
 
 if( s == D.TheInteractionListList.size() ) return;  
   
 for( n = 0; n < D.TheInteractionListList[ s ] -> size(); ++n )
  if( (*(D.TheInteractionListList[ s ]))[ n ] != NULL )
   IncludeInteraction( (*((*(D.TheInteractionListList[ s ]))[ n ])) );
   
}

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

ClassicalStates& ClassicalStates::operator=( ClassicalStates& D )
{
 int s, n;
 ParticleListPointer    ParticleListNullPointer = NULL;
 InteractionListPointer InteractionListNullPointer = NULL;
 
 // check if we need to do something first
 if( &D == &(*this) ) return( (*this) );
 
 Id = D.Id;
 TheInternalFlag = D.TheInternalFlag;
 Type = D.Type;
 
 // then delete all of the lists in our two state lists
 DeleteAllTheParticleLists(); 
 DeleteAllTheInteractionLists(); 
 
 // make sure new state list (list of lists) is equal of size or bigger
 s = D.NrOfStatesIncluded() - TheParticleListList.size();   
 if( s > 0 )                                                               
 {
  TheParticleListList.insert( TheParticleListList.end(), s, ParticleListNullPointer );
  TheInteractionListList.insert( TheInteractionListList.end(), s, InteractionListNullPointer );  
 }
 
  // start by copying by particles from D to our own lists... keep in mind they still point
 // to interactions of D rather than to interactions of ourselfs...
 for( s = 0; s < D.TheParticleListList.size(); ++s )
 {
 
  // copy the particle list ...
  if( D.TheParticleListList[ s ] != NULL )
  {
   TheParticleListList[ s ] = new ParticleList( D.TheParticleListList[ s ] -> size() );
   
   if( TheParticleListList[ s ] != NULL )
    for( n = 0; n < D.TheParticleListList[ s ] -> size(); ++n )
     if( (*(D.TheParticleListList[ s ]))[ n ] != NULL )
     {
      (*(TheParticleListList[ s ]))[ n ] = (*(D.TheParticleListList[ s ]))[ n ] -> Clone(); 
      if( (*(TheParticleListList[ s ]))[ n ] != NULL )
       (*(TheParticleListList[ s ]))[ n ] -> DynamicsContext = this;
     } 
     else
      (*(TheParticleListList[ s ]))[ n ] = NULL; 
  }
  else
   TheParticleListList[ s ] = NULL;
   
  // create interaction state lists... but init them to NULL...
  if( D.TheInteractionListList[ s ] != NULL )
  {
   TheInteractionListList[ s ] = new InteractionList( D.TheInteractionListList[ s ] -> size() );
   
   if( TheInteractionListList[ s ] != NULL )
    for( n = 0; n < D.TheInteractionListList[ s ] -> size(); ++n )
     (*(TheInteractionListList[ s ]))[ n ] = NULL;
  } 
  else 
   TheInteractionListList[ s ] = NULL;
}
 
 // we have now copied the all the particles and initialized the interaction lists to NULL.
 // now start including the interactions of D into ourself. we include them so that the
 // particles and interactions will point to eachother within this dynamical context rather
 // than in the ones we are copying...
 
 for( s = 0; s < D.TheInteractionListList.size(); ++s )   // find a non-NULL entry in the D states
  if( D.TheInteractionListList[ s ] != NULL )
   break;

 if( s == TheInteractionListList.size() ) return( (*this) );  
   
 for( n = 0; n < D.TheInteractionListList[ s ] -> size(); ++n )
  if( (*(D.TheInteractionListList[ s ]))[ n ] != NULL )
   IncludeInteraction( (*((*(D.TheInteractionListList[ s ]))[ n ])) );
   
 // blank rest of states to NULL
 for( s = D.TheParticleListList.size(); s < TheParticleListList.size(); ++s )                     
 {
  TheParticleListList[ s ] = NULL;
  TheInteractionListList[ s ] = NULL;  
 }
 
 return( (*this) );
}

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

void ClassicalStates::IncludeParticle( Particle& P )
{
 int                EmptySpot, Pos, I, s;
 ParticlePointer    NullPointer = NULL;
 InteractionPointer AnInteraction;

 for( s = 0; s < TheParticleListList.size(); ++s )   // find a non-NULL entry in the states
  if( TheParticleListList[ s ] != NULL )
   break;
 
 if( s == TheParticleListList.size() ) return;  
   
 /* loop through first array and find empty spot or return if particle already in first list */
 for( EmptySpot = -1, Pos = 0; Pos < TheParticleListList[ s ] -> size(); Pos++ )
 {
  if( (*(TheParticleListList[ s ]))[ Pos ] == NULL )
  {
   if( EmptySpot <= -1 )
    EmptySpot = Pos;
  } 
  else 
   if( (*(TheParticleListList[ s ]))[ Pos ] -> Id == P.Id )
    return;
 }

 /* if no empty spots where found, increase all lists with NULL pointers */ 
 if( EmptySpot <= -1 )     
 {
  for( s = 0; s < TheParticleListList.size(); ++s )
   if( TheParticleListList[ s ] != NULL )
    TheParticleListList[ s ] -> insert( TheParticleListList[ s ] -> end(), DEFAULT_DELTA_NR_PARTICLES, NullPointer );
  EmptySpot = Pos;
 } 
 
 /* make a copy of P and add it into each list */
 for( s = 0; s < TheParticleListList.size(); ++s )
  if( TheParticleListList[ s ] != NULL )
  {
   (*(TheParticleListList[ s ]))[ EmptySpot ] = P.Clone();   
   if( (*(TheParticleListList[ s ]))[ EmptySpot ] != NULL )
    (*(TheParticleListList[ s ]))[ EmptySpot ] -> DynamicsContext = this;
  }  
 
 /* we have now a copy of particle P in our lists, the copy however points 
    to interactions, which in turn points to the original particle etc. we
    need to make sure that the copy particles point to copy interactions, 
    which in turn point correctly to the copy particles. 
       
    if one of the interactions of the particle is not yet included, automatically
    include that interaction and thus make copy interactions. 
   
    tricky buisness but we need it ;-).
 */
   
 for( I = P.Interactions.size() - 1; I >= 0; --I )      // loop over each interaction of P
 {
  AnInteraction = P.Interactions[ I ];                    // get pointer to the interaction to get ID
  if( AnInteraction == NULL ) continue;
    
  Pos = IsInteractionIncluded( AnInteraction -> Id );     // see if interaction is already present
  
  if( Pos <= -1 ) 
  {
   // the interaction of particle P, of which we made copies, is not yet
   // included into our own internal lists. we therefore need to include
   // it and automatically make copies of it...
   
   IncludeInteraction( (*AnInteraction) );
  }
  else
  {
   // the interaction of particle P was found in our own list. we now need
   // to make sure that the copy particles point to the copy interactions and
   // the copy interactions point to the copy particles... it is at this point
   // where the copy particles IncludeInteraction shoud update the interaction
   // pointer to our own copied interaction from the list...
   
   for( s = 0; s < TheParticleListList.size(); ++s )
    if( ( TheParticleListList[ s ] != NULL ) && ( TheInteractionListList[ s ] != NULL ) )
     (*(TheParticleListList[ s ]))[ EmptySpot ] -> IncludeInteraction( (*((*(TheInteractionListList[ s ]))[ Pos ])) );
   
  } 
 } 
   
  
}

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

void ClassicalStates::IncludeInteraction( Interaction& I )
{
 int                EmptySpot, Pos, P, s;
 InteractionPointer NullPointer = NULL;
 ParticlePointer    AParticle;

 for( s = 0; s < TheInteractionListList.size(); ++s )   // find a non-NULL entry in the states
  if( TheInteractionListList[ s ] != NULL )
   break;

 if( s == TheInteractionListList.size() ) return;  
 
 /* loop through first array and find empty spot or return if interaction already in first list */
 for( EmptySpot = -1, Pos = 0; Pos < TheInteractionListList[ s ] -> size(); Pos++ )
 {
  if( (*(TheInteractionListList[ s ]))[ Pos ] == NULL )
  {
   if( EmptySpot <= -1 )
    EmptySpot = Pos;
  }
  else 
   if( (*(TheInteractionListList[ s ]))[ Pos ] -> Id == I.Id )
    return;
 }
 
 /* if no empty spots where found, increase all lists with NULL pointers */ 
 if( EmptySpot <= -1 )     
 {
  for( s = 0; s < TheInteractionListList.size(); ++s )
   if( TheInteractionListList[ s ] != NULL )
    TheInteractionListList[ s ] -> insert( TheInteractionListList[ s ] -> end(), DEFAULT_DELTA_NR_INTERACTIONS, NullPointer );
  EmptySpot = Pos;
 } 

 /* make a copy of I and add it into each list */
 for( s = 0; s < TheInteractionListList.size(); ++s )
  if( TheInteractionListList[ s ] != NULL )
  {
   (*(TheInteractionListList[ s ]))[ EmptySpot ] = I.Clone();  
   if( (*(TheInteractionListList[ s ]))[ EmptySpot ] != NULL )
    (*(TheInteractionListList[ s ]))[ EmptySpot ] -> DynamicsContext = this;
  }  

 /* we have now a copy of interaction I in our lists, the copy however points 
    to particles, which in turn points to the original interaction etc. we
    need to make sure that the copy interaction point to copy particles, 
    which in turn point correctly to the copy interaction. 
       
    if one of the particles of the interaction is not yet included, automatically
    include that particle and thus make copy particles.
    
    tricky buisness but we need it ;-).
 */

 for( P = I.Particles.size() - 1; P >= 0; --P )      // loop over each particle of I
 {
  AParticle = I.Particles[ P ];                      // get pointer to the particle to get ID
  if( AParticle == NULL ) continue;
    
  Pos = IsParticleIncluded( AParticle -> Id );         // see if particle is already present
  
  if( Pos <= -1 ) 
  {
   // the particle of interaction I, of which we made copies, is not yet
   // included into our own internal lists. we therefore need to include
   // it and automatically make copies of it...
   
   IncludeParticle( (*AParticle) );
  }
  else
  {
   // the particle of interaction I was found in our own list. we now need
   // to make sure that the copy interactions point to the copy particles and
   // the copy particles point to the copy interactions... it is at this point
   // where the copy interactions IncludeParticle shoud update the particle
   // pointer to our own copied particle from the list...
   
   for( s = 0; s < TheInteractionListList.size(); ++s )
    if( ( TheInteractionListList[ s ] != NULL ) && ( TheParticleListList[ s ] != NULL ) )
     (*(TheInteractionListList[ s ]))[ EmptySpot ] -> IncludeParticle( (*((*(TheParticleListList[ s ]))[ Pos ])) );
   
  } 
 } 

}

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

void ClassicalStates::DeleteParticle( Particle& P )
{
 DeleteParticle( P.Id );
} 

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

void ClassicalStates::DeleteInteraction( Interaction& I )
{
 DeleteInteraction( I.Id );
}

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

void ClassicalStates::DeleteParticle( unsigned int ID )
{
 int             Pos, s, ss;
 ParticlePointer P;

 for( ss = 0; ss < TheParticleListList.size(); ++ss )   // find a non-NULL entry in the states
  if( TheParticleListList[ ss ] != NULL )
   break;
 
 if( ss == TheParticleListList.size() ) return;  
 
 /* loop through the first list and find particles with given ID */
 for( Pos = 0; Pos < TheParticleListList[ ss ] -> size(); Pos++ )
 {
  P = (*(TheParticleListList[ ss ]))[ Pos ];
  
  if( P != NULL )                       
   if( P -> Id == ID )
   {
   
    for( s = 0; s < TheParticleListList.size(); ++s )    // if particle found, remove it from all lists
    {
     if( TheParticleListList[ s ] == NULL ) continue;
     
     P = (*(TheParticleListList[ s ]))[ Pos ];   
     
     (*(TheParticleListList[ s ]))[ Pos ] = NULL;
     
     if( P != NULL )
      delete P;
    } 
    
   } 
 }
 
}

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

void ClassicalStates::DeleteInteraction( unsigned int ID )
{
 int                Pos, s, ss;
 InteractionPointer I;
 
 for( ss = 0; ss < TheInteractionListList.size(); ++ss )   // find a non-NULL entry in the states
  if( TheInteractionListList[ ss ] != NULL )
   break;
 
 if( ss == TheInteractionListList.size() ) return;  

 /* loop through the first list and find interactions with given ID */
 for( Pos = 0; Pos < TheInteractionListList[ ss ] -> size(); Pos++ )
 {
  I = (*(TheInteractionListList[ ss ]))[ Pos ];
  
  if( I != NULL )                       
   if( I -> Id == ID )
   {
   
    for( s = 0; s < TheInteractionListList.size(); ++s )    // if interaction found, remove it from all lists
    {
     if( TheInteractionListList[ s ] == NULL ) continue;
     
     I = (*(TheInteractionListList[ s ]))[ Pos ];   
     
     (*(TheInteractionListList[ s ]))[ Pos ] = NULL;
     
     if( I != NULL )
      delete I;
    } 
    
   } 
 }
 
}

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

void ClassicalStates::DeleteAllParticles( void )
{
 int             Pos, s;
 ParticlePointer P;
 
 for( s = 0; s < TheParticleListList.size(); ++s )   
 {
  if( TheParticleListList[ s ] == NULL ) continue; 
  
  for( Pos = 0; Pos < TheParticleListList[ s ] -> size(); ++Pos )
  {
   P = (*(TheParticleListList[ s ]))[ Pos ];
   
   (*(TheParticleListList[ s ]))[ Pos ] = NULL;
   
   if( P != NULL )
    delete P;
  } 
 }
}

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

void ClassicalStates::DeleteAllInteractions( void )
{
 int                Pos, s;
 InteractionPointer I;
 
 for( s = 0; s < TheInteractionListList.size(); ++s )   
 {
  if( TheInteractionListList[ s ] == NULL ) continue;
  
  for( Pos = 0; Pos < TheInteractionListList[ s ] -> size(); ++Pos )
  {
   I = (*(TheInteractionListList[ s ]))[ Pos ];
   
   (*(TheInteractionListList[ s ]))[ Pos ] = NULL;
   
   if( I != NULL )
    delete I;
  } 
 }
}

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

void ClassicalStates::DeleteAllTheParticleLists( void )
{
 DeleteAllParticles();
 for( int s = 0; s < TheParticleListList.size(); ++s )
 {
  if( TheParticleListList[ s ] != NULL )
   delete TheParticleListList[ s ];
  TheParticleListList[ s ] = NULL;
 }
}


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

void ClassicalStates::DeleteAllTheInteractionLists( void )
{
 DeleteAllInteractions();
 for( int s = 0; s < TheInteractionListList.size(); ++s )
 {
  if( TheInteractionListList[ s ] != NULL )
   delete TheInteractionListList[ s ];
  TheInteractionListList[ s ] = NULL;
 }
}

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

int  ClassicalStates::IsInteractionIncluded( unsigned int ID )
{
 int Pos, s;

 for( s = 0; s < TheInteractionListList.size(); ++s )   // find a non-NULL entry in the states
  if( TheInteractionListList[ s ] != NULL )
   break;
 
 if( s == TheInteractionListList.size() ) return( -1 );   
 
 for( Pos = 0; Pos < TheInteractionListList[ s ] -> size(); Pos++ )
 {
  if( (*(TheInteractionListList[ s ]))[ Pos ] != NULL )
   if( (*(TheInteractionListList[ s ]))[ Pos ] -> Id == ID )
    return( Pos );
 }
 
 return( -1 );
}

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

int  ClassicalStates::IsInteractionIncluded( Interaction& I )
{ 
 return( IsInteractionIncluded( I.Id ) );
}

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

int  ClassicalStates::IsParticleIncluded( unsigned int ID )
{
 int Pos, s;

 for( s = 0; s < TheParticleListList.size(); ++s )   // find a non-NULL entry in the states
  if( TheParticleListList[ s ] != NULL )
   break;
 
 if( s == TheParticleListList.size() ) return( -1 );   
 
 for( Pos = 0; Pos < TheParticleListList[ s ] -> size(); Pos++ )
 {
  if( (*(TheParticleListList[ s ]))[ Pos ] != NULL )
   if( (*(TheParticleListList[ s ]))[ Pos ] -> Id == ID )
    return( Pos );
 }
 
 return( -1 );
}

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

int  ClassicalStates::IsParticleIncluded( Particle& P )
{
 return( IsParticleIncluded( P.Id ) );
}

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

int ClassicalStates::NrOfInteractionsIncluded( void )
{
 int i, I, s;
 
 for( s = 0; s < TheInteractionListList.size(); ++s )   // find a non-NULL entry in the states
  if( TheInteractionListList[ s ] != NULL )
   break;
 
 if( s == TheInteractionListList.size() ) return( 0 );   
 
 for( i = I = 0; i < TheInteractionListList[ s ] -> size(); ++i )
  if( (*(TheInteractionListList[ s ]))[ i ] != NULL ) ++I;
  
 return( I ); 
}

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

int ClassicalStates::NrOfParticlesIncluded( void )
{
 int p, P, s;
 
 for( s = 0; s < TheParticleListList.size(); ++s )   // find a non-NULL entry in the states
  if( TheParticleListList[ s ] != NULL )
   break;
 
 if( s == TheParticleListList.size() ) return( 0 );    
 
 for( p = P = 0; p < TheParticleListList[ 0 ] -> size(); ++p )
  if( (*(TheParticleListList[ 0 ]))[ p ] != NULL ) ++P;
  
 return( P ); 
}

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

int ClassicalStates::NrOfStatesIncluded( void )
{
 int s, S;
 
 for( s = S = 0; s < TheParticleListList.size(); ++s )
  if( TheParticleListList[ s ] != NULL ) S++;
  
 return( S ); 
}

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

ParticleListPointer ClassicalStates::GetStateParticlesPointer( int s )
{
 return( s < 0 ? NULL : ( s >= TheParticleListList.size() ? NULL : TheParticleListList[ s ] ) );
}

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

InteractionListPointer ClassicalStates::GetStateInteractionsPointer( int s )
{
 return( s < 0 ? NULL : ( s >= TheInteractionListList.size() ? NULL : TheInteractionListList[ s ] ) );
}

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

ParticlePointer ClassicalStates::GetParticlePointer( int Index, int s )
{
 ParticleListPointer State = ( s < 0 ? NULL : ( s >= TheParticleListList.size() ? NULL : TheParticleListList[ s ] ) );
 return( State == NULL ? NULL : ( Index < 0 ? NULL : ( Index >= State -> size() ? NULL : (*State)[ Index ] ) ) );
}

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

InteractionPointer ClassicalStates::GetInteractionPointer( int Index, int s )
{
 InteractionListPointer State = ( s < 0 ? NULL : ( s >= TheInteractionListList.size() ? NULL : TheInteractionListList[ s ] ) );
 return( State == NULL ? NULL : ( Index < 0 ? NULL : ( Index >= State -> size() ? NULL : (*State)[ Index ] ) ) );
}

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

int ClassicalStates::GetParticleListSize( void )
{
 int s;
 
 for( s = 0; s < TheParticleListList.size(); ++s )   // find a non-NULL entry in the states
  if( TheParticleListList[ s ] != NULL )
   break;
 
 if( s == TheParticleListList.size() ) return( 0 );   
 
 return( TheParticleListList[ s ] -> size() );
}

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

int ClassicalStates::GetInteractionListSize( void )
{
 int s;
 
 for( s = 0; s < TheInteractionListList.size(); ++s )   // find a non-NULL entry in the states
  if( TheInteractionListList[ s ] != NULL )
   break;
 
 if( s == TheInteractionListList.size() ) return( 0 ); 
 
 return( TheInteractionListList[ s ] -> size() );
}

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

int ClassicalStates::GetStateListSize( void )
{
 return( TheInteractionListList.size() );
}

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

void ClassicalStates::IncludeParticleInInteraction( Particle& P, Interaction& I )
{
 int PPos, IPos;
 
 // first check if the particle is included into the interaction
 I.IncludeParticle( P );
 
 // then check if interaction is included into context
 IPos = IsInteractionIncluded( I.Id );   
 if( IPos <= -1 ) return;
  
 // now make sure the particle is also included into the context
 IncludeParticle( P );
 
 // then get particles position
 PPos = IsParticleIncluded( P.Id );
 if( PPos <= -1 ) return;
 
 // now we make sure that all copy interactions have the copy particle included, if
 // the particle was present in the context before we included it ourselves
 for( int s = 0; s < TheInteractionListList.size(); ++s )
  if( TheInteractionListList[ s ] != NULL )
   (*(TheInteractionListList[ s ]))[ IPos ] -> IncludeParticle( (*((*(TheParticleListList[ s ]))[ PPos ])) );
  
}

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

void ClassicalStates::IncludeInteractionInParticle( Interaction& I, Particle& P )
{
 int PPos, IPos;
 
 // first check if the interaction is included into the particle
 P.IncludeInteraction( I );
 
 // then check if particle is included into context
 PPos = IsParticleIncluded( P.Id );   
 if( PPos <= -1 ) return;
  
 // now make sure the interaction is also included into the context
 IncludeInteraction( I );
 
 // then get interactions position
 IPos = IsInteractionIncluded( I.Id );
 if( IPos <= -1 ) return;
 
 // now we make sure that all copy particles have the copy interaction included, if
 // the interaction was present in the context before we included it ourselves
 for( int s = 0; s < TheParticleListList.size(); ++s )
  if( TheParticleListList[ s ] != NULL )
   (*(TheParticleListList[ s ]))[ PPos ] -> IncludeInteraction( (*((*(TheInteractionListList[ s ]))[ IPos ])) );

}

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

void ClassicalStates::ExcludeInteractionFromParticle( unsigned int IID, unsigned int PID )
{
 int PPos;

 // check if particle is included into context
 PPos = IsParticleIncluded( PID );   
 if( PPos <= -1 ) return;

 // now exclude the interaction from all copy particles
 for( int s = 0; s < TheParticleListList.size(); ++s )
  if( TheParticleListList[ s ] != NULL )
   (*(TheParticleListList[ s ]))[ PPos ] -> ExcludeInteraction( IID );
   
}

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

void ClassicalStates::ExcludeParticleFromInteraction( unsigned int PID, unsigned int IID )
{
 int IPos;
 
 // check if interaction is included into context
 IPos = IsInteractionIncluded( IID );   
 if( IPos <= -1 ) return;

 // now exclude the particle from all the copy interactions
 for( int s = 0; s < TheInteractionListList.size(); ++s )
  if( TheInteractionListList[ s ] != NULL )
   (*(TheInteractionListList[ s ]))[ IPos ] -> ExcludeParticle( PID );
   
}

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

void ClassicalStates::ExcludeInteractionFromParticle( Interaction& I, Particle& P )
{
 P.ExcludeInteraction( I );
 ExcludeInteractionFromParticle( I.Id, P.Id );
}

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

void ClassicalStates::ExcludeParticleFromInteraction( Particle& P, Interaction& I )
{
 I.ExcludeParticle( P );
 ExcludeParticleFromInteraction( P.Id, I.Id );
}

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

void ClassicalStates::ExcludeAllInteractionsFromParticle( unsigned int PID )
{
 int PPos;

 // check if particle is included into context
 PPos = IsParticleIncluded( PID );   
 if( PPos <= -1 ) return;

 // now exclude the interaction from all copy particles
 for( int s = 0; s < TheParticleListList.size(); ++s )
  if( TheParticleListList[ s ] != NULL )
   (*(TheParticleListList[ s ]))[ PPos ] -> ExcludeAllInteractions();

}

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

void ClassicalStates::ExcludeAllParticlesFromInteraction( unsigned int IID )
{
 int IPos;
 
 // check if interaction is included into context
 IPos = IsInteractionIncluded( IID );   
 if( IPos <= -1 ) return;

 // now exclude the particle from all the copy interactions
 for( int s = 0; s < TheInteractionListList.size(); ++s )
  if( TheInteractionListList[ s ] != NULL )
   (*(TheInteractionListList[ s ]))[ IPos ] -> ExcludeAllParticles();
   
}

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

void ClassicalStates::ExcludeAllInteractionsFromParticle( Particle& P )
{
 P.ExcludeAllInteractions();
 ExcludeAllInteractionsFromParticle( P.Id );
}

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

void ClassicalStates::ExcludeAllParticlesFromInteraction( Interaction& I )
{
 I.ExcludeAllParticles();
 ExcludeAllParticlesFromInteraction( I.Id );
}

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

void ClassicalStates::DeleteAllInteractionsOfParticle( unsigned int PID )
{
 int PPos;

 // check if particle is included into context
 PPos = IsParticleIncluded( PID );   
 if( PPos <= -1 ) return;

 // now delete the interaction from all copy particles
 for( int s = 0; s < TheParticleListList.size(); ++s )
  if( TheParticleListList[ s ] != NULL )
   (*(TheParticleListList[ s ]))[ PPos ] -> DeleteAllInteractions();
 
}

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

void ClassicalStates::DeleteAllParticlesOfInteraction( unsigned int IID )
{
 int IPos;
 
 // check if interaction is included into context
 IPos = IsInteractionIncluded( IID );   
 if( IPos <= -1 ) return;

 // now delete the particle from all the copy interactions
 for( int s = 0; s < TheInteractionListList.size(); ++s )
  if( TheInteractionListList[ s ] != NULL )
   (*(TheInteractionListList[ s ]))[ IPos ] -> DeleteAllParticles();

}

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

void ClassicalStates::DeleteAllInteractionsOfParticle( Particle& P )
{
 P.DeleteAllInteractions();
}

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

void ClassicalStates::DeleteAllParticlesOfInteraction( Interaction& I )
{
 I.DeleteAllParticles();
}

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

void ClassicalStates::CycleTheStates( int Steps )
{
 int i, s, Size;
 ParticleListPointer    TempP;
 InteractionListPointer TempL;
 
 if( Steps == 0 ) return;

 Size = NrOfStatesIncluded(); 
 Steps = Steps % Size;
 
 if( Steps < 0 )         /* cycle back in state list */
 {
 
  for( i = 0; i > Steps; --i )
  {
   TempL = TheInteractionListList[ Size - 1 ];
   TempP = TheParticleListList[ Size - 1 ];
   for( s = Size - 1; s > 0; --s )
   {
    TheInteractionListList[ s ] = TheInteractionListList[ s - 1 ];
    TheParticleListList[ s ] = TheParticleListList[ s - 1 ];
   }
   TheInteractionListList[ 0 ] = TempL;
   TheParticleListList[ 0 ] = TempP; 
  }
  
 }
 else                    /* cycle forwards in state list */
 {
 
  for( i = 0; i < Steps; ++i )
  {
   TempL = TheInteractionListList[ 0 ];
   TempP = TheParticleListList[ 0 ];
   for( s = 0; s < Size - 1; ++s )
   {
    TheInteractionListList[ s ] = TheInteractionListList[ s + 1 ];
    TheParticleListList[ s ] = TheParticleListList[ s + 1 ];
   }
   TheInteractionListList[ Size - 1 ] = TempL;
   TheParticleListList[ Size - 1 ] = TempP; 
  }
  
 }
 
}

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