/*   []----------------------------------------------[]
      |                   Particle.cpp               |
     []----------------------------------------------[]
      |                                              |
      |   AUTHOR:    MFSomers 2005.                  |
      |   USE:       Defines a particle to be used   |
      |              in a classical dynamics         |
      |              framework...                    |
      |                                              |
     []----------------------------------------------[]
*/

// 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"

/* ---------------------------------------------------- Particle --------- */

Particle::Particle( unsigned int TheID ) : Interactions( 1 )
{
 Position.X = Position.Y = Position.Z = 0.0;
 Momentum.X = Momentum.Y = Momentum.Z = 0.0;
 Mass = Charge = Radius = 0.0;
 Id = TheID;
 Interactions[ 0 ] = NULL;
 ColorIndex = TheInternalFlag = 0;
 Type = AGeneralParticle; 
 DynamicsContext = NULL;
}

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

Particle::Particle( Vector ThePosition, Vector TheMomentum, double TheMass, double TheCharge,
                    double TheRadius, int TheColorIndex, unsigned int TheID ) : Interactions( 1 )
{
 Position = ThePosition;
 Momentum = TheMomentum;
 Mass = TheMass;
 Charge = TheCharge;
 Radius = TheRadius;
 ColorIndex = TheColorIndex;
 Id = TheID;
 Interactions[ 0 ] = NULL;
 TheInternalFlag = 0;
 Type = AGeneralParticle;
 DynamicsContext = NULL;
}

/* ----------------------------------------------------------------------- */ 
		 
Particle::Particle( Particle& P ) : Interactions( P.Interactions.size() )
{
  Position = P.Position;
  Momentum = P.Momentum;
  Mass = P.Mass;
  Charge = P.Charge;
  Radius = P.Radius;
  Id = P.Id;
  TheInternalFlag = P.TheInternalFlag;
  Type = P.Type;  
  ColorIndex = P.ColorIndex;
  DynamicsContext = NULL;
  
  for( int N = 0; N < P.Interactions.size(); ++N )
   Interactions[ N ] = P.Interactions[ N ];
}

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

void Particle::IncludeInteraction( Interaction& I )
{
 int                EmptySpot, Pos;
 InteractionPointer NullPointer = NULL;
 
 /* loop through array and find empty spot or return if interaction already in list */
 for( EmptySpot = -1, Pos = 0; Pos < Interactions.size(); Pos++ )
 {
  if( Interactions[ Pos ] == NULL ) 
  {
   if( EmptySpot <= -1 )
    EmptySpot = Pos;
  }  
  else 
   if( Interactions[ Pos ] -> Id == I.Id )
   {
    if( Interactions[ Pos ] != &I )
    {
     Interactions[ Pos ] = &I;             // refresh links however to this I, see ClassicalDynamics !
     I.IncludeParticle( (*this) );
    } 
    return;
   } 
 }
   
 if( EmptySpot <= -1 )     /* if no empty spots where found */
 {
  Interactions.insert( Interactions.end(), 5, NullPointer );
  EmptySpot = Pos;
 } 
 
 Interactions[ EmptySpot ] = &I;
 I.IncludeParticle( (*this) );     
 
/* now if this particle is part of a dynamical context, and this interaction is not
   yet included into the same context, include it into the context. if it is included
   however, all the copy particles should also have it included */
   
 if( DynamicsContext != NULL )
  DynamicsContext -> IncludeInteractionInParticle( I, (*this) );
   
}

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

Particle& Particle::operator=( Particle& P )
{
 InteractionPointer NullPointer = NULL;
 int N, n;
 
 if( &P == &(*this) ) return( (*this) );
 
 Position = P.Position;
 Momentum = P.Momentum;
 Mass = P.Mass;
 Charge = P.Charge;
 Radius = P.Radius;
 Id = P.Id;
 ColorIndex = P.ColorIndex;
 TheInternalFlag = P.TheInternalFlag;
 Type = P.Type;

 N = P.Interactions.size() - Interactions.size();
 if( N > 0 ) Interactions.insert( Interactions.end(), N, NullPointer );

 for( n = 0; n < P.Interactions.size(); ++n )          // copy table
  Interactions[ n ] = P.Interactions[ n ];
   
 for( ; n < Interactions.size(); ++n )                 // and set rest to NULL 
   Interactions[ n ] = NULL;    

 return( (*this) );
}

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

int Particle::NrOfInteractionsIncluded( void )
{
 int i, I;
 
 for( i = I = 0; i < Interactions.size(); ++i )
  if( Interactions[ i ] != NULL ) ++I;
  
 return( I ); 
}

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

void Particle::ExcludeInteraction( unsigned int ID )
{
 InteractionPointer I;

 /* loop through array remove interactions with given ID */
 for( int Pos = 0; Pos < Interactions.size(); Pos++ )
 {
  I = Interactions[ Pos ];
  
  if( I != NULL )
   if( I -> Id == ID )
   {
    Interactions[ Pos ] = NULL;
    
    I -> ExcludeParticle( Id );
    
    /* now if we are part of a dynamical context, make sure that the interaction is also
    excluded from all the copy particles within the same dynamical context */
      
    if( DynamicsContext != NULL ) DynamicsContext -> ExcludeInteractionFromParticle( ID, Id ); 
   } 
 }
 
}

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

void Particle::DeleteInteraction( unsigned int ID )
{
 InteractionPointer I;

 /* now if we are part of a dynamical context, make sure that the interaction is
    deleted from all the copy particles within the same dynamical context */
      
 if( DynamicsContext != NULL ) DynamicsContext -> DeleteInteraction( ID );
 
 /* loop through array remove interactions with given ID */
 for( int Pos = 0; Pos < Interactions.size(); Pos++ )
 {
  I = Interactions[ Pos ];
  
  if( I != NULL )
   if( I -> Id == ID )
   {
    Interactions[ Pos ] = NULL;
    delete I;
   } 
 }
 
}

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

void Particle::ExcludeInteraction( Interaction& I )
{
 ExcludeInteraction( I.Id );
} 

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

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

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

int Particle::IsInteractionIncluded( unsigned int ID )
{
 int                Size = Interactions.size();
 InteractionPointer I;
 
 for( int Pos = 0; Pos < Size; Pos++ )
 {
  I = Interactions[ Pos ];
  
  if( I != NULL )
   if( I -> Id == ID )
    return( Pos );
 }
 
 return( -1 );
}

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

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

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

void Particle::ExcludeAllInteractions( void )
{
 InteractionPointer I;

 for( int Pos = 0; Pos < Interactions.size(); Pos++ )
 {
  I = Interactions[ Pos ];
  
  Interactions[ Pos ] = NULL;
  
  if( I != NULL )
  {
   I -> ExcludeParticle( Id );
   
   /* now if we are part of a dynamical context, please exclude all interactions from
      all the copy particles also within the same dynamical context */
   
   if( DynamicsContext != NULL )
    DynamicsContext -> ExcludeInteractionFromParticle( I -> Id, Id );
  }
   
 }
 
}

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

void Particle::DeleteAllInteractions( void )
{
 InteractionPointer I;

 for( int Pos = 0; Pos < Interactions.size(); Pos++ )
 {
  I = Interactions[ Pos ];
  
  Interactions[ Pos ] = NULL;
  
  if( I != NULL )
  {
   /* now if we are part of a dynamical context, please delete these interactions
      also within the same context */
  
   if( DynamicsContext != NULL )
    DynamicsContext -> DeleteInteraction( I -> Id );
   else           
    delete I;
  } 
 }
 
}

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

int Particle::GetInteractionListSize( void )
{
 return( Interactions.size() );
}

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

InteractionPointer Particle::GetInteractionPointer( int Index )
{
 return( ( Index < 0 ? NULL : ( Index >= Interactions.size() ? NULL : Interactions[ Index ] ) ) );
}

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

double Particle::TheTotalInteractionPotential( void )
{
 int                Size = Interactions.size();
 double             TotalInteraction = 0.0;
 InteractionPointer I;

 for( int Index = 0; Index < Size; Index++ )
 {
  I = Interactions[ Index ];
  
  if( I != NULL )
   TotalInteraction += I -> InteractionPotential();
 }  
 
 return( TotalInteraction );  
}

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

Vector Particle::TheTotalForce( void )
{
 int                Size = Interactions.size();
 Vector             TotalForce;
 InteractionPointer I;
 
 TotalForce.X = TotalForce.Y = TotalForce.Z = 0.0;
 
 for( int Index = 0; Index < Size; Index++ )
 {
  I = Interactions[ Index ];
  
  if( I != NULL )
   TotalForce += I -> InteractionForce( Id );
 }  
 
 return( TotalForce );  
}

/* --------------------------------------------------- Functions --------- */

Vector CentreOfMass( const ParticleList& TheParticles )
{
 int    Size = TheParticles.size();
 double TotalMass;
 Vector TheCentreOfMass;
 
 TotalMass = TheCentreOfMass.X = TheCentreOfMass.Y = TheCentreOfMass.Z = 0.0;

 for( int P = 0; P < Size; ++P )
  if( TheParticles[ P ] != NULL )
  {
   if( TheParticles[ P ] -> Mass != 0.0 )
   {
    TheCentreOfMass += TheParticles[ P ] -> Position / TheParticles[ P ] -> Mass;
    TotalMass += TheParticles[ P ] -> Mass;
   } 
  } 
   
 return( TheCentreOfMass / TotalMass ); 
}

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

Vector DipoleMoment( const ParticleList& TheParticles )
{
 int    Size = TheParticles.size();
 Vector TheMoment;
 
 TheMoment.X = TheMoment.Y = TheMoment.Z = 0.0;
 
 for( int P = 0; P < Size; ++P )
  if( TheParticles[ P ] != NULL )
   TheMoment += TheParticles[ P ] -> Charge * TheParticles[ P ] -> Position;
   
 return( TheMoment ); 
}

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

double TotalKineticEnergy( const ParticleList& TheParticles )
{
 int    Size = TheParticles.size();
 double E = 0.0;

 for( int P = 0; P < Size; ++P )
  if( TheParticles[ P ] != NULL )
   E += TheParticles[ P ] -> TheKineticEnergy();
   
 return( E );  
}

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

double TotalMass( const ParticleList& TheParticles )
{
 int    Size = TheParticles.size();
 double Mass = 0.0;

 for( int P = 0; P < Size; ++P )
  if( TheParticles[ P ] != NULL )
   Mass += TheParticles[ P ] -> Mass;
   
 return( Mass );  
}

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

double TotalCharge( const ParticleList& TheParticles )
{
 int    Size = TheParticles.size();
 double Charge = 0.0;
 
 for( int P = 0; P < Size; ++P )
  if( TheParticles[ P ] != NULL )
   Charge += TheParticles[ P ] -> Charge;
   
 return( Charge );  
}

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

double Virial( const ParticleList& TheParticles, const Vector& TheBox )
{
 int    Size = TheParticles.size();
 double TheVirial = 0.0;
 
 for( int P = 0; P < Size; ++P )
  if( TheParticles[ P ] != NULL )
   TheVirial += Dot( TheParticles[ P ] -> TheTotalForce(), MimimumSystemImageConventionVector( TheParticles[ P ] -> Position, TheBox ) );
     
 return( TheVirial );   
}

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

double MeanDistanceSquared( const ParticleList& TheParticles )
{
 int    Size = TheParticles.size();
 double TheMeanDistanceSquared = 0.0;
 int    NrOfParticles = 0;
 
 for( int P = 0; P < Size; ++P )
  if( TheParticles[ P ] != NULL )
  {
   TheMeanDistanceSquared += AbsSqr( TheParticles[ P ] -> Position );
   ++NrOfParticles;
  } 
     
 return( NrOfParticles ? TheMeanDistanceSquared / NrOfParticles : 0.0 );   
}

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