#include <windows.h>
#include <math.h>

#include <math.h>
//#include <rfftw.h>

#include "Pattern.h"



/*



Parent version (003):

  use forward projection instead of semi-lagrangian lookups
  combine ood2 lookup with dirx and diry

This version:

  Made gravity travel at the speed of light instead of infitiely fast

  Method:
    Store previous mass frames
    Instead of using current mass frame for each pixel's calculation use a past mass
    frame depending on the distance of the pixel from the one it is affecting.

  I'll call the technique Pass Mass Lookup

*/





//---- Render buffer

PATTERN_RGB buf;


/*

  //  strange amoeba
#define GRAVITY_STRENGTH -0.00333
#define STUFF_MAX      2
#define AMBIENT_STUFF  0
#define SOL 2
#define WEAKER_FORCE
#define DIM 32
#define VELOCITY_SCALE 1
*/

/*
  //  Chaotic galaxy life
#define GRAVITY_STRENGTH 1
#define STUFF_MAX      2
#define AMBIENT_STUFF  0
#define SOL 20000
#define WEAKER_FORCE
#define DIM 32
#define VELOCITY_SCALE 1
*/

/*
//self-reproducing mini-galaxies

#define GRAVITY_STRENGTH 0.0125
#define STUFF_MAX      2
#define AMBIENT_STUFF  0
#define SOL 20000
#define WEAKER_FORCE
#define DIM 32
#define VELOCITY_SCALE 1
*/

/*
//  less galaxies, still self-reproducing
#define GRAVITY_STRENGTH 0.005
#define STUFF_MAX      2
#define AMBIENT_STUFF  0
#define SOL 20000
#define WEAKER_FORCE
#define DIM 32
#define VELOCITY_SCALE 1
*/

/*
//  sub-particle swarm
#define GRAVITY_STRENGTH 0.5
#define STUFF_MAX      0.01
#define AMBIENT_STUFF  0
#define SOL 20000
#define WEAKER_FORCE
#define DIM 32
#define VELOCITY_SCALE 1
*/

/*
//  wispy smoke
#define GRAVITY_STRENGTH 0.1
#define STUFF_MAX      0.125
#define AMBIENT_STUFF  0
#define SOL 20000
#define WEAKER_FORCE
#define SWAP_DX_AND_DY
#define DIM 41
#define VELOCITY_SCALE 1
*/

/*
//  atomic soup
#define GRAVITY_STRENGTH 1
#define STUFF_MAX      0.125
#define AMBIENT_STUFF  0
#define SOL 20000
#define WEAKER_FORCE
#define DIM 41
#define VELOCITY_SCALE 1
*/

/*
//  Galaxy with some nice rotational motion
#define GRAVITY_STRENGTH 0.005
#define STUFF_MAX      15.525
#define AMBIENT_STUFF  0
#define SOL 20000
//#define WEAKER_FORCE
#define DIM 32
#define VELOCITY_SCALE 1
*/

#define GRAVITY_STRENGTH 0.005
#define STUFF_MAX      25
#define AMBIENT_STUFF  0
#define SOL 200000
//#define WEAKER_FORCE
#define DIM 32
#define VELOCITY_SCALE 1
//#define DONT_WRAP_UNVIERSE
//#define SPEED_LIMIT 5.0


//  Speed of Light


//  Ammount of previous mass frames to store (depends on speed of light and universe size)

//  next power of 2 above ((int)( (DIM*1.5) / SOL ))

#define PASTLEN 64


#define FRAMES_BEFORE_PAUSE -1





  //---- Galaxy data

  double m0[ DIM ][ DIM ];  //  Mass double buffer 0
  double m1[ DIM ][ DIM ];  //  Mass double buffer 1

  double mbuf[ PASTLEN ][ DIM ][ DIM ];  //  Mass buffers, storing previous frames

  double dx0[ DIM ][ DIM ], dx1[ DIM ][ DIM ];  //  X-velocity double buffer
  double dy0[ DIM ][ DIM ], dy1[ DIM ][ DIM ];  //  Y-velocity double buffer

  double last_px[ DIM ][ DIM ], last_py[ DIM ][ DIM ];  //  Storage of gravity field (un-necessary for calc, just for display)

  double dirx_lookup[ DIM*2 ][ DIM*2 ];
  double diry_lookup[ DIM*2 ][ DIM*2 ];


  //- Lookup tables to enable past mass lookup

  int    past_distance_i[ DIM*2 ][ DIM*2 ];
  double past_distance_d[ DIM*2 ][ DIM*2 ];



  //- Animation data

  int frame = 0;
  int time = 0;



  int mouse_x, mouse_y;
  int lmx=0, lmy=0;  //  last known mouse pos




    double (*current_m )[DIM], (*new_m )[DIM];  //  Current and new mass
    double (*current_dx)[DIM], (*new_dx)[DIM];  //  Current and new x speed
    double (*current_dy)[DIM], (*new_dy)[DIM];  //  Current and new y speed








  //--- Startup function

  void Pattern_Startup()
  {
    int i, j;

    //  Initial empty mass and speed buffers

    for(j=0; j<DIM; j++) for(i=0; i<DIM; i++)
    {
      m0[j][i]  = m1[j][i]  = 0;
      dx0[j][i] = dx1[j][i] = 0;
      dy0[j][i] = dy1[j][i] = 0;
    }

    for(j=0; j<DIM-0; j++) for(i=0; i<DIM-0; i++)
    {
      mbuf[0][j][i] = mbuf[1][j][i]  = 10;//(rand()%10);
      dx0[j][i]     = dx1[j][i]      = 0;//(rand()%7)-3;
      dy0[j][i]     = dy1[j][i]      = 0;//(rand()%7)-3;
    }

    for(j=0; j<DIM-0; j++) for(i=0; i<DIM-0; i++)
    {
      mbuf[0][j][i] = mbuf[1][j][i]  = AMBIENT_STUFF + 0.001*(rand()%((int)(1000.0*(STUFF_MAX-AMBIENT_STUFF))));
      dx0[j][i]     = dx1[j][i]      = 1*((rand()%7)-3);
      dy0[j][i]     = dy1[j][i]      = 1*((rand()%7)-3);
    }
    
    //- Gravity precalc

	  int x, y;
	
	  //  calculate the distance from the centre point (SIZE,SIZE) in (SIZE*2,SIZE*2) array
	
	  double dx, dy, d2, d;

	  for( y=0; y<DIM*2; y++ )
	  {
	    for( x=0; x<DIM*2; x++ )
	    {
        dirx_lookup[ y ][ x ] = 0;

        if( x != DIM || y != DIM )  //  was && - that's wrong!!
	      {
          for( j=-20; j<=20; j++ ) for( i=-20; i<=20; i++ )
          {

#ifdef SWAP_DX_AND_DY

	          dy =   (x-DIM+i*DIM);
	          dx = - (y-DIM+j*DIM);
#else

	          dx = (x-DIM+i*DIM);
	          dy = (y-DIM+j*DIM);
#endif

	          
	          d2 = dx*dx + dy*dy;   //  Distance squared
	  
	          d = sqrt( d2 );  //  Distance

            #ifdef WEAKER_FORCE

              dirx_lookup[ y ][ x ] += dx / d2 / d2;  //  X Influence of a unit of mass
	            diry_lookup[ y ][ x ] += dy / d2 / d2;  //  Y Influence of a unit of mass

            #else

              dirx_lookup[ y ][ x ] += dx / d / d2;  //  X Influence of a unit of mass
	            diry_lookup[ y ][ x ] += dy / d / d2;  //  Y Influence of a unit of mass

            #endif
          }

          //- Do speed of light lookup
//
	        dx = (x-DIM);
	        dy = (y-DIM);
	        
	        d2 = dx*dx + dy*dy;   //  Distance squared

	        d = sqrt( d2 );  //  Distance
//
          d /= SOL;

          past_distance_i[ y ][ x ] = (int)d;      //  Integer part of how far into the past to look
          past_distance_d[ y ][ x ] = d - (int)d;  //  Fractional part of how far into the past to look
	      }
	      else  //  Special case for centre point
	      {
	        dirx_lookup[ y ][ x ] = 0;
	        diry_lookup[ y ][ x ] = 0;

          past_distance_i[ y ][ x ] = 0;
          past_distance_d[ y ][ x ] = 0;
	      }
	    }
	  }

  }












  void recalc()
  {
    if( FRAMES_BEFORE_PAUSE >= 0 && frame >= FRAMES_BEFORE_PAUSE ) return;

    int i, j;

    double newi, newj;
    double mm;
    double px, py;
    double  dx, dy;
    

    //- Add mouse generated matter
/*
    if( mouse_drag )
    {
      int x, y;
      x = mouse_x;
      y = mouse_y;
      for(j=y-2; j<=y+2; j++) for(i=x-2; i<=x+2; i++)
        if(i >= 0 && i < DIM && j >= 0 && j < DIM)
        {
          current_dx[j][i] = 10*mouse_dx * 1.345 / (10+current_m[j][i]);// + (j-y)*0.0121;
          current_dy[j][i] = 10*mouse_dy * 1.345 / (10+current_m[j][i]);// + (x-i)*0.0121;
          current_m[j][i] += 10;
        }
    }
*/  
    
    int ii, jj;
    int ii1, jj1;

    //- Calculate new speeds
    
    int x, y;

    int    past_i;
    double past_d;
    double m, past_mbuf0, past_mbuf1;

    for(j=0; j<DIM; j++) for(i=0; i<DIM; i++)
    {
	    px=0;
	    py=0;
	    
	    //  Calculate force for (i,j) by calculating the influence of every pixel (x,y)
	
      for( y=0; y<DIM; y++ )
      {
        for( x=0; x<DIM; x++ )
        {
          //  The current_m here needs to be a linear interpolation of previous masses
          //  depending on the distance

          past_i = past_distance_i[ DIM+y-j ][ DIM+x-i ];  //  Integer how far to look back
          past_d = past_distance_d[ DIM+y-j ][ DIM+x-i ];  //  Decimal how far to look back

          past_mbuf0 = mbuf[ (frame-past_i)&31 ][ y ][ x ];    //  Get mass from past_i
          past_mbuf1 = mbuf[ (frame-past_i-1)&31 ][ y ][ x ];  //  Get mass from frame before past_i

          //  now this may just work, but how do I tell?

          m = past_mbuf0 * (1.0-past_d) + past_mbuf1 * past_d;  //  perceived mass of the pixel

          //  TODO: Past mass lookup

          px += dirx_lookup[DIM+y-j][DIM+x-i] * m;
          py += diry_lookup[DIM+y-j][DIM+x-i] * m;
        }
      }

      px *= GRAVITY_STRENGTH;
      py *= GRAVITY_STRENGTH;
      
      last_px[j][i] = px;
      last_py[j][i] = py;
      
      //  Update V - acceleration of pixel is force (px,py) / mass

      if( current_m[j][i] > 0 )
      {
        current_dx[j][i] += px / current_m[j][i];
        current_dy[j][i] += py / current_m[j][i];

        #ifdef SPEED_LIMIT
          px = current_dx[j][i]*current_dx[j][i] + current_dy[j][i]*current_dy[j][i];
          if( px > SPEED_LIMIT*SPEED_LIMIT )
          {
            px = sqrt( px / SPEED_LIMIT );
            current_dx[j][i] /= px;
            current_dy[j][i] /= px;
          }
        #endif
      }
      else
      {
        current_dx[j][i] = 0;
        current_dy[j][i] = 0;
      }
    }



    //- Clear destination arrays

    for(j=0; j<DIM; j++) for(i=0; i<DIM; i++)
    {
      new_m[j][i] = 0;
      new_dx[j][i] = 0;
      new_dy[j][i] = 0;
    }



    //- Move matter and velocities by sub-pixel correct forward projection

    double m00, m01, m10, m11;

    for(j=0; j<DIM; j++) for(i=0; i<DIM; i++)
    {
      mm = current_m [j][i];   //  Mass of this pixel last time
      if( mm == 0 ) continue;  //  No point doing anything if there is no stuff

      //  New position

      newi = (double)(i) + current_dx[j][i] * VELOCITY_SCALE;
      newj = (double)(j) + current_dy[j][i] * VELOCITY_SCALE;

      //  Wrap around
/*
      if( newi <  0   ) newi += DIM;
      if( newi >= DIM ) newi -= DIM;

      if( newj <  0   ) newj += DIM;
      if( newj >= DIM ) newj -= DIM;
*/

      //  Simply discard if it falls out of the array
#ifdef DONT_WRAP_UNVIERSE
      if( newi < 0    ) continue;
      if( newj < 0    ) continue;
      if( newi > DIM-1 ) continue;
      if( newj > DIM-1 ) continue;
#endif


      //  Integer positions with (un-necessary) wrap-around
      
      px = newi - floor(newi);  //  Fractional parts of position
      py = newj - floor(newj);

      if( px < 0 ) px += 1;
      if( py < 0 ) py += 1;

      ii = newi - floor( newi / DIM )*DIM;
      jj = newj - floor( newj / DIM )*DIM;
      
//      if( px < 0 ) px += 1;
//      if( py < 0 ) py += 1;

      //  Add loads of DIM's on so that the mod works properly

      ii = ii % DIM;
      jj = jj % DIM;

      if( ii < 0 ) ii += DIM;
      if( jj < 0 ) jj += DIM;

      ii1 = (ii+1) % DIM;
      jj1 = (jj+1) % DIM;

      if( jj < 0 || jj >= DIM || ii < 0 || ii >= DIM || px < 0 || py < 0 )
      {
        // random breakpoint - should never happen but does if you don't have (A) and (B)
        int a = DIM;
        a+=3;
        //  TEMP: Just move onto next pixel; if it's going too fast mass just gets destroyed for now
        continue;
      }

      //  Move mass
            
      new_m [ jj  ][ ii  ] += m11 = mm * (1-px) * (1-py);
      new_m [ jj  ][ ii1 ] += m01 = mm * (  px) * (1-py);
      new_m [ jj1 ][ ii  ] += m10 = mm * (1-px) * (  py);
      new_m [ jj1 ][ ii1 ] += m00 = mm * (  px) * (  py);
      
      //  Move speed

      dx = current_dx[j][i];
      dy = current_dy[j][i];

      new_dx[ jj  ][ ii  ] += dx * m11;  //  Speed is weiDIMted by the portion of the
      new_dx[ jj  ][ ii1 ] += dx * m01;  //  mass it is bringing into each square
      new_dx[ jj1 ][ ii  ] += dx * m10;
      new_dx[ jj1 ][ ii1 ] += dx * m00;

      new_dy[ jj  ][ ii  ] += dy * m11;
      new_dy[ jj  ][ ii1 ] += dy * m01;
      new_dy[ jj1 ][ ii  ] += dy * m10;
      new_dy[ jj1 ][ ii1 ] += dy * m00;
    }
    
    //  Scale down the velocities (the new velocity is a weighted average of the velocities passed in)
    //  array currently holds total of the velocities, needs to be divided by the mass now
    
    for(j=0; j<DIM; j++) for(i=0; i<DIM; i++)
    {
      if( new_m[j][i] > 0 )
      {
        new_dx[j][i] /= new_m[j][i];
        new_dy[j][i] /= new_m[j][i];
      }
    }    


  }
























//---- Pattern_Generate ------------------------------------------------------//



PATTERN_RGB *Pattern_Generate()
{

  //- Decide which speed buffer is primary and which is secondary

  if( (frame&1) == 0 )
  {
    current_dx = dx0; new_dx = dx1;
    current_dy = dy0; new_dy = dy1;
  }
  else
  {
    current_dx = dx1; new_dx = dx0;
    current_dy = dy1; new_dy = dy0;
  }



  //- Next mass buffer

  current_m = mbuf[ (frame  )&31 ];
  new_m     = mbuf[ (frame+1)&31 ];



  //- Animate

  recalc();



  //- Render

  int i, j, k, idx, c;

      int x, y;

    int    past_i;
    double past_d;
    double m, past_mbuf0, past_mbuf1;

    int ii, jj;

    //  plot lookup table

  for (j = 0; j < DIM*2; j++)
  {
		for (i = 0; i < DIM*2; i++, idx++)
    {
      c = 128+dirx_lookup[j][i]*13256;
      if( c < 0 ) c = 0; if( c > 255 ) c = 255;
      buf.buf[j][i+DIM*2][0] = c;
      
      c = 128+diry_lookup[j][i]*13256;
      if( c < 0 ) c = 0; if( c > 255 ) c = 255;
      buf.buf[j][i+DIM*2][1] = c;

      buf.buf[j][i+DIM*2][2] = 0;
    }
  }


  //  plot the universe

//  int x, y;

	for (j = 0; j < DIM; j++)
  {
		idx = j * DIM;
						
		for (i = 0; i < DIM; i++, idx++)
    {
      //  mass

      c = (int)(new_m[j][i] * 25 );

      if( c < 0 ) c = 0;
      if( c > 255 ) c = 255;

			buf.buf[j][i][0] = c;
			buf.buf[j][i][1] = c;
			buf.buf[j][i][2] = c;



      //- relative view

      ii = FRAMES_BEFORE_PAUSE >= 0 && frame < FRAMES_BEFORE_PAUSE ? DIM/2 : mouse_x-DIM;
      jj = FRAMES_BEFORE_PAUSE >= 0 && frame < FRAMES_BEFORE_PAUSE ? DIM/2 : mouse_y;
      x = i;
      y = j;

      //  The current_m here needs to be a linear interpolation of previous masses
      //  depending on the distance

      past_i = past_distance_i[ DIM+y-jj ][ DIM+x-ii ];  //  Integer how far to look back
      past_d = past_distance_d[ DIM+y-jj ][ DIM+x-ii ];  //  Decimal how far to look back

      past_mbuf0 = mbuf[ (frame-past_i)&31 ][ y ][ x ];    //  Get mass from past_i
      past_mbuf1 = mbuf[ (frame-past_i-1)&31 ][ y ][ x ];  //  Get mass from frame before past_i

      //  now this may just work, but how do I tell?

      m = past_mbuf0 * (1.0-past_d) + past_mbuf1 * past_d;  //  perceived mass of the pixel

      //  Put pixel

      c = (int)(m * 25 );

      if( c < 0 ) c = 0;
      if( c > 255 ) c = 255;

      buf.buf[j][i+DIM][0] = c;
			buf.buf[j][i+DIM][1] = c;
			buf.buf[j][i+DIM][2] = c;



      // x velocity

      c = (int)(last_px[j][i] * 128.0 + 128.0);

      if( c < 0 ) c = 0;
      if( c > 255 ) c = 255;

			buf.buf[j+DIM][i][0] = c;
			buf.buf[j+DIM][i][1] = c;
			buf.buf[j+DIM][i][2] = c;

      //  y velocity

      c = (int)(last_py[j][i] * 128.0 + 128.0);

      if( c < 0 ) c = 0;
      if( c > 255 ) c = 255;

			buf.buf[j+DIM][i+DIM][0] = c;
			buf.buf[j+DIM][i+DIM][1] = c;
			buf.buf[j+DIM][i+DIM][2] = c;



      //  tiled universe

      //  r has to be from 0..1 as m goes from 0..infinity

      //  r = 1 - 2^(-m)
      //  y^z = exp( z * ln( y ) )

      if( DIM <= 42 )
      {
  //      double r = 1.0 - 1.0/(new_m[j][i]+1.0);
        double r = 1.0 - exp(-new_m[j][i]*log(2.0));

        int rr =  255.0 * r + r * 255.0 * sin( 3.14159*0.66 + r*3.14159*1.5 );
        int gg =  255.0 * r + r * 255.0 * sin( 3.14159*1.33 + r*3.14159*1.5 );
        int bb =  255.0 * r + r * 255.0 * sin( 3.14159* r*3.14159*1.5 );

        if( rr < 0 ) rr = 0;
        if( gg < 0 ) gg = 0;
        if( bb < 0 ) bb = 0;

        if( rr > 255 ) rr = 255;
        if( gg > 255 ) gg = 255;
        if( bb > 255 ) bb = 255;

        for( y=0; y<2; y++ ) for( x=0; x<2; x++ )
        {
			    buf.buf[ j*2 + DIM*2 + DIM*y*2     ][ i*2 + DIM*x*2     ][ 2 ] = rr;
			    buf.buf[ j*2 + DIM*2 + DIM*y*2     ][ i*2 + DIM*x*2     ][ 1 ] = gg;
			    buf.buf[ j*2 + DIM*2 + DIM*y*2     ][ i*2 + DIM*x*2     ][ 0 ] = bb;

          buf.buf[ j*2 + DIM*2 + DIM*y*2 + 1 ][ i*2 + DIM*x*2     ][ 2 ] = rr;
			    buf.buf[ j*2 + DIM*2 + DIM*y*2 + 1 ][ i*2 + DIM*x*2     ][ 1 ] = gg;
			    buf.buf[ j*2 + DIM*2 + DIM*y*2 + 1 ][ i*2 + DIM*x*2     ][ 0 ] = bb;

          buf.buf[ j*2 + DIM*2 + DIM*y*2     ][ i*2 + DIM*x*2 + 1 ][ 2 ] = rr;
			    buf.buf[ j*2 + DIM*2 + DIM*y*2     ][ i*2 + DIM*x*2 + 1 ][ 1 ] = gg;
			    buf.buf[ j*2 + DIM*2 + DIM*y*2     ][ i*2 + DIM*x*2 + 1 ][ 0 ] = bb;

          buf.buf[ j*2 + DIM*2 + DIM*y*2 + 1 ][ i*2 + DIM*x*2 + 1 ][ 2 ] = rr;
			    buf.buf[ j*2 + DIM*2 + DIM*y*2 + 1 ][ i*2 + DIM*x*2 + 1 ][ 1 ] = gg;
			    buf.buf[ j*2 + DIM*2 + DIM*y*2 + 1 ][ i*2 + DIM*x*2 + 1 ][ 0 ] = bb;
        }
      }

    }

	}

  //- Draw mini past frames

  for( k=0; k<PASTLEN; k++ )
  {
    for( j=0; j<8; j++ )
    {
      for( i=0; i<8; i++ )
      {
        c = (int)( mbuf[ k ][ j*DIM/8 ][ i*DIM/8 ] * 25 );

        if( c < 32 ) c = 32;
        if( c > 255 ) c = 255;

        buf.buf[ j+DIM*2 ][ k*8+i ][ 0 ] = c;
        buf.buf[ j+DIM*2 ][ k*8+i ][ 1 ] = c;
        buf.buf[ j+DIM*2 ][ k*8+i ][ 2 ] = c;
      }
    }
  }



  //- Move on

  if( FRAMES_BEFORE_PAUSE < 0 || frame < FRAMES_BEFORE_PAUSE ) frame ++;
  time ++;



  //- Return buffer for displaying

  return &buf;
}















void MouseMove( int x, int y )
{
  mouse_x = x;
  mouse_y = y;
}

/*
 * When the user drags with the mouse, add a force
 * that corresponds to the direction of the mouse
 * cursor movement. Also inject som new matter into
 * the field at the mouse location.
 *
 */
void drag(int mx, int my)
{
	double  dx, dy, len;
	int     X, Y;


	/* Compute the array index that corresponds to the
	 * cursor location */

	X = mx;
	Y = my;

	if (X > (DIM - 1)) {
		X = DIM - 1;
	}
	if (Y > (DIM - 1)) {
		Y = DIM - 1;
	}
	if (X < 0) {
		X = 0;
	}
	if (Y < 0) {
		Y = 0;
	}

	/* Add force at the cursor location */

//	my = winHeight - my;
/*
	dx = mx - lmx;
	dy = my - lmy;
	len = sqrt(dx * dx + dy * dy);
  if (len != 0.0) { 
    dx *= 1.f / len;
    dy *= 1.f / len;
  }
	u_u0[Y * DIM + X] += dx;
	u_v0[Y * DIM + X] += dy;
*/
	/* Increase matter density at the cursor location  */

	current_m[Y][X] += 10.0f;
	new_m[Y][X] += 10.0f;

	current_dx[Y][X] = new_dx[Y][X] = current_dy[Y][X] = new_dy[Y][X] = 0;

	lmx = mx;
	lmy = my;
}