//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
//   Program fluid
//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*

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

//---- definition of fundermental constants
//---- experimental settings
const double RE          = 100.0;           // Reynolds number
const double IRE         = 1.0/RE;
const double dT          = 0.0625;            // temporal step
const int    XI_NUM      = 24;              // radial cells
const int    TH_NUM      = 48;              // angular cells
const double dXI         = 0.125;           // radial step
const double dTH         = 2*M_PI/(TH_NUM-2); // angular step
const double IdXI        = 1.0/dXI;
const double IdTH        = 1.0/dTH;

//---- index indicators
const int    XI_SUR      = 0;               // index of surface
const int    XI_FAR      = XI_NUM-1;        // index of farthest
const int    TH_STA      = 1;               // index of start angle
const int    TH_END      = TH_NUM-2;        // index of end angle

//---- graphics settings
const int    WIN_WIDTH   = 256;             // width of window
const int    WIN_HEIGHT  = 256;             // height of window
const int    CYLINDER_R  = 24;              // radius of cylinder
const double V_MAG       = 16.0;            // magnify factor of velocity

//---- prototype of functions
void InitFlow( void );
void SetBorder( void );
void ModifyPsi( void );
void UpdateOmg( void );
void DrawFlow( void );


//---- global vairables
//---- field variables
double Psi[XI_NUM][TH_NUM];  // Flow Function field
double Omg[XI_NUM][TH_NUM];  // Vorticity field 

//---- tables

class ExpTable
{
  double expo[4*XI_NUM];
public:
  ExpTable( void ){
    for( int i=-2*XI_NUM ; i<2*XI_NUM ; i++ ){
      expo[2*XI_NUM+i] = exp( i*dXI );
    }
  }
  inline double operator () ( int i ){
    return expo[2*XI_NUM+i];
  }
} Exp;

class SinTable
{
  double sine[TH_NUM];
public:
  SinTable( void ){
    for( int j=TH_STA ; j<=TH_END ; j++ ){
      sine[j] = sin( j*dTH );
    }
  }
  inline double operator () ( int j ){
    return sine[j];
  }
} Sin;

class CosTable
{
  double cosi[TH_NUM];
public:
  CosTable( void ){
    for( int j=TH_STA ; j<=TH_END ; j++ ){
      cosi[j] = cos( j*dTH );
    }
  }
  inline double operator () ( int j ){
    return cosi[j];
  }
} Cos;

class XTable
{
  int x[XI_NUM][TH_NUM];
public:
  XTable( void ){
    for( int i=XI_SUR ; i<=XI_FAR ; i++ ){
      for( int j=TH_STA ; j<=TH_END ; j++ ){
	x[i][j] = WIN_WIDTH/2  + int(CYLINDER_R*Exp(i)*Cos(j));
      }
    }
  }
  inline int operator () ( int i, int j ){
    return x[i][j];
  }
} X;

class YTable
{
  int y[XI_NUM][TH_NUM];
public:
  YTable( void ){
    for( int i=XI_SUR ; i<=XI_FAR ; i++ ){
      for( int j=TH_STA ; j<=TH_END ; j++ ){
	y[i][j] = WIN_HEIGHT/2 + int(CYLINDER_R*Exp(i)*Sin(j));
      }
    }
  }
  inline int operator () ( int i, int j ){
    return y[i][j];
  }
} Y;








//---- main function
int main( void )
{
  NXOpenWindow( "2D flow around a cylinder", WIN_WIDTH, WIN_HEIGHT );
  NXEnable( NX_DOUBLEBUFFER );

  InitFlow();           // inttializes flow as uniformal flow

  while( NXNoEvents() ){
    SetBorder();        // sets border conditions
    ModifyPsi();        // solves phi by eq13
    UpdateOmg();        // calculates next omg by eq10
    DrawFlow();         // draws velocity filed
  }

  NXCloseWindow();
  return 0;
}





//---- intialize flow as uniformal flow
void InitFlow( void )
{
  int i, j;
  for( i=XI_SUR ; i<=XI_FAR ; i++ ){
    for( j=TH_STA ; j<=TH_END ; j++ ){
      Psi[i][j] = -Exp(i)*Sin(j);
      Omg[i][j] = 0.0;
    }
  }
}


//---- set border conditions on surface and far away
void SetBorder( void )
{
  int j;
  for( j=TH_STA ; j<=TH_END ; j++ ){
      // surface
    Omg[XI_SUR][j] = (-2.0*IdXI*IdXI) * Psi[XI_SUR+1][j];
    Psi[XI_SUR][j] = 0.0;
      // far away
    Psi[XI_FAR][j] = -Exp(XI_FAR)*Sin(j);
    Omg[XI_FAR][j] = 0.0;
  }
}


//---- derive flow function field from Laplace equation
void ModifyPsi( void )
{
  int i, j, loop=8;

  while( loop-- ){
    for( i=XI_SUR+1 ; i<=XI_FAR-1 ; i++ ){
      Psi[i][TH_STA-1] = Psi[i][TH_END];   // cyclic condition
      Psi[i][TH_END+1] = Psi[i][TH_STA];
      for( j=TH_STA ; j<=TH_END ; j++ ){
        Psi[i][j] =
          (0.5/(IdXI*IdXI+IdTH*IdTH))*(
            +(IdXI*IdXI)*(Psi[i+1][j]+Psi[i-1][j])
            +(IdTH*IdTH)*(Psi[i][j-1]+Psi[i][j+1])
            +Exp(2*i)*Omg[i][j]
          );
      }
    }
  }
}


//---- update vorticity field
void UpdateOmg( void )
{
  int i, j;
  for( i=XI_SUR+1 ; i<=XI_FAR-1 ; i++ ){
    Omg[i][TH_STA-1] = Omg[i][TH_END];   // periodic condition
    Omg[i][TH_END+1] = Omg[i][TH_STA];
    for( j=TH_STA ; j<=TH_END ; j++ ){
      Omg[i][j] +=
        dT*Exp(-2*i)*(
          + (IdXI*0.5) * (Psi[i+1][j]-Psi[i-1][j])
          * (IdTH*0.5) * (Omg[i][j+1]-Omg[i][j-1])
          - (IdTH*0.5) * (Psi[i][j+1]-Psi[i][j-1])
          * (IdXI*0.5) * (Omg[i+1][j]-Omg[i-1][j])
        )
        +
        dT*Exp(-2*i)*IRE*(
          + (IdXI*IdXI) * (Omg[i+1][j]-2*Omg[i][j]+Omg[i-1][j])
          + (IdTH*IdTH) * (Omg[i][j+1]-2*Omg[i][j]+Omg[i][j-1])
        );
    }
  }
}


//---- derive velocity field and draw them 
void DrawFlow( void )
{
  int i, j;
  double Vxi, Vth, Vx, Vy;
  
  NXClearWindow();
  for( i=XI_SUR+1 ; i<=XI_FAR-1 ; i++ ){
    for( j=TH_STA ; j<=TH_END ; j++ ){
      Vxi = (+0.5*IdTH) * Exp(-i)*(Psi[i][j+1]-Psi[i][j-1]);
      Vth = (-0.5*IdXI) * Exp(-i)*(Psi[i+1][j]-Psi[i-1][j]);
      Vx  = Vxi*Cos(j) - Vth*Sin(j);
      Vy  = Vxi*Sin(j) + Vth*Cos(j);

      NXDrawMoveto( X(i,j), Y(i,j) );
      NXDrawLinerel( int(V_MAG*Vx), int(V_MAG*Vy) );
    }
  }
  NXFlush();
}