//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
//   Program phonon
//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*

#include <stdio.h>
#include <math.h>
#define  ANIMATION
#include "nxgraph.h"

//---- physical setting
#define N          (128)
#define K          (1.0)
#define M          (1.0)
#define dT         (1.0/16)

//---- graphic setting
#define WIN_WIDTH  (256)
#define WIN_HEIGHT (128)


void Init( double x[N], double p[N] )
{
  for( int i=0 ; i<N ; i++ ){
    x[i] = exp(-0.01*(i-N/2)*(i-N/2));
    p[i] = 0.0;
  }
}

void Edge( double x[N], double p[N] )
{
  x[0] = x[N-1] = 0.0;
  p[0] = p[N-1] = 0.0;
}

void Evolve( double x[N], double p[N] )
{
  int i;

  for( i=0 ; i<=N-1 ; i++ ){
    x[i] += p[i]/M*dT;
  }
  p[0] += K*( x[1] - x[0] )*dT;
  for( i=1 ; i<=N-2 ; i++ ){
    p[i] += K*( x[i+1] + x[i-1] - 2*x[i] )*dT;
  }
  p[N-1] += K*( x[N-2] - x[N-1] )*dT;
}

void Draw( double x[N] )
{
  NXClearWindow();
  NXDrawMoveto( 2*0, WIN_HEIGHT/2 - (int)(64*x[0]) );
  for( int i=1 ; i<N ; i++ ){
    NXDrawLineto( 2*i, WIN_HEIGHT/2 - (int)(64*x[i]) );
  }
  NXFlush();
}


//---- main function
int main(void)
{
  XEvent ev;
  NXOpenWindow("Phonon Wave", WIN_WIDTH, WIN_HEIGHT );

  double x[N], p[N];

  Init( x, p );
  do{
    Edge( x, p );

    for( int interval=0 ; interval<16 ; interval++ ){
      Evolve( x, p );
    }
    Draw( x );
    NXCheckEvent( NX_NOWAIT, ev );
  }while( ev.type != KeyPress );

  NXCloseWindow();

  return(0);
}