//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
//   Program stars
//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*

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

//---- physical constants
const double Mstar      = 1.0;                   // mass of a star
const double G          = 1.0;                   // Gravity constant

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

//---- 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;        // momentum, position, force
};

//---- prototypes of functions
void    Init( Matter M[] );
void    Evolve( 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 Stars", WIN_WIDTH, WIN_HEIGHT );

  Init( M );

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

  NXCloseWindow();
  return 0;
}

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

//---- evolves the system
void Evolve( Matter M[], double dt )
{
  int n;
  for( n=0; n<N; n++ ){
    M[n].q += 0.5*dt * (1.0/Mstar) * M[n].p;
  }
  CalcField( M );
  for( n=0; n<N; n++ ){
    M[n].p += dt * M[n].f;
  }
  for( n=0; n<N; n++ ){
    M[n].q += 0.5*dt * (1.0/Mstar) * M[n].p;
  }
}


//---- calculates gravity force amoung two stars
inline Vector2 Interaction( Matter& M1, Matter& M2 )
{
  const Vector2 vec_r = M1.q - M2.q;
  const double  r     = abs(vec_r);
  return -G*Mstar*Mstar/(r*r*r) * vec_r;
}

//---- calculates gravity 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;
    }
  }
}

//---- calculates Total Energy
double CalcEnergy( Matter M[] )
{
  int n1, n2;
  double E;
    // sums up the energy
  for( n1=0, E=0.0; n1<N; n1++ ){
    E += (1.0/(2*Mstar))*norm(M[n1].p);         // kinetic energy
    for( n2=n1+1; n2<N; n2++ ){
      E += -G*Mstar*Mstar/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) );
}