//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
//   Program electrons
//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*

#include "cip.h"
#include "vector2.h"
#include "nxgraph.h"

//---- physical constants
const double Mele    = 1.0;                 // mass of an electron
const double Eele    = 1.0;                 // charge of an electron
const double Aint    = 1.0;                 // coupling constant of ele-ele interaction
const double Bext    = 2.0;                 // external magnetic field
const double Omega   = Eele*Bext/Mele;      // Larmor anglur frequency

//---- experimental settings
const double dT      = 1.0/32;             // time slice
const int N          = 4;                   // number of electrons

//---- graphic setting
const int WIN_WIDTH  = 256;
const int WIN_HEIGHT = 256;
const double MAG     =  64.0;

//---- declaration of structure Matter
struct Matter
{
  Vector2 p, q, f;
};

//---- prototypes of functions
void    Init( Matter M[] );
void    Evolve( Matter M[], double dt );
void    EvolveU( Matter M[], double dt );
void    EvolveK( Matter M[], double dt );
void    CalcField( Matter M[] );
double  CalcEnergy( Matter M[] );
void    Draw( Matter M[] );


//---- main function
int main( void )
{
  double T=0.0;
  Matter M[N];

  NXOpenWindow("Dance of Electrons", WIN_WIDTH, WIN_HEIGHT );

  Init( M );

  while( NXNoEvents() ){
    Evolve( M, dT );
    T += dT;
    Draw( M );
  }

  NXCloseWindow();
  return 0;
}

//---- initializes properties of electrons
void Init( Matter M[] )
{
  for( int n=0; n<N; n++ ){
      // locates electrons on circumference
    double ang = 2*M_PI/N*n;
    M[n].p = Vector2( 0.0, 0.0 );
    M[n].q = Vector2( cos(ang), sin(ang) );
  }
}

//---- evolves the system by using symplectic integrator
void Evolve( Matter M[], double dt )
{
  int n;

  for( n=0; n<N; n++ ){
    const double co = cos(Omega*0.5*dt);
    const double si = sin(Omega*0.5*dt);
    
    M[n].q.x += (1.0/(Mele*Omega))*( M[n].p.x*si - M[n].p.y*(1-co) );
    M[n].q.y += (1.0/(Mele*Omega))*( M[n].p.x*(1-co) + M[n].p.y*si );

    const double p_x = M[n].p.x;
    M[n].p.x = p_x*co - M[n].p.y*si;
    M[n].p.y = p_x*si + M[n].p.y*co;
  }

  CalcField( M );
  for( n=0; n<N; n++ ){
    M[n].p += dt * M[n].f;
  }

  for( n=0; n<N; n++ ){
    const double co = cos(Omega*0.5*dt);
    const double si = sin(Omega*0.5*dt);
    
    M[n].q.x += (1.0/(Mele*Omega))*( M[n].p.x*si - M[n].p.y*(1-co) );
    M[n].q.y += (1.0/(Mele*Omega))*( M[n].p.x*(1-co) + M[n].p.y*si );

    const double p_x = M[n].p.x;
    M[n].p.x = p_x*co - M[n].p.y*si;
    M[n].p.y = p_x*si + M[n].p.y*co;
  }
}


//---- calculates electric force amoung two electrons
inline Vector2 Interaction( Matter& M1, Matter& M2 )
{
  const Vector2 vec_r = M1.q - M2.q;
  const double  r     = abs(vec_r);
  return Aint*Eele*Eele/(r*r*r) * vec_r;
}

//---- calculates internal electric field
void CalcField( Matter M[] )
{
  int n1, n2;
    // resets M[].f
  for( n1=0; n1<N; n1++ ){
    M[n1].f = Vector2( 0.0, 0.0 );
  }
    // sums up M[].f
  for( n1=0; n1<N; n1++ ){
    for( n2=n1+1; n2<N; n2++ ){
      Vector2 f = Interaction( M[n1], M[n2] );
      M[n1].f += f;
      M[n2].f -= f;  // Newton's law "The third law of motion"
    }
  }
}


//---- calculates total energy
double CalcEnergy( Matter M[] )
{
  int n1, n2;
  double E;
    // sums up energy
  for( n1=0, E=0.0; n1<N; n1++ ){
    E += 0.5/Mele*norm(M[n1].p);                  // kinetic Energy
    for( n2=n1+1; n2<N; n2++ ){
      E += Aint*Eele*Eele/abs(M[n1].q-M[n2].q);   // potential Energy
    }
  }
  return E;
}

//---- draws traces
void Draw( Matter M[] )
{
  for( int n=0; n<N; n++ ){
    NXDrawPoint( WIN_WIDTH/2  + int(MAG*M[n].q.x),
                 WIN_HEIGHT/2 + int(MAG*M[n].q.y) );
  }
  NXDrawPrintf( 32, 32, "Energy=%20.17lf", CalcEnergy(M) );
}