#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;







//  Dimensions of universe

#define DIM 32

//  Speed of Light

#define SOL 3.0

//  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 32






  //---- 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;





    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
















  void recalc()
  {
    int i, j;

    int ii, jj;
    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;
        }
    }
*/  
    

    //- 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 ];
          past_mbuf1 = mbuf[ (frame-past_i-1)&31 ][ y ][ x ];

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

          m = past_mbuf0 * (1.0-past_d) + past_mbuf1 * past_d;

          //  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 *= 0.00333;
      py *= 0.00333;
      
      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];
      }
    }



    //- 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];
      newj = (double)(j) + current_dy[j][i];

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

      //  Integer positions with (un-necessary) wrap-around

      ii = (int)newi;
      jj = (int)newj;
      
      px = newi - (int)newi;  //  Fractional parts of position
      py = newj - (int)newj;
//      if( px < 0 ) px += 1;
//      if( py < 0 ) py += 1;

      //  Move mass
            
      new_m [ jj   ][ ii   ] += m11 = mm * (1-px) * (1-py);
      new_m [ jj   ][ ii+1 ] += m01 = mm * (  px) * (1-py);
      new_m [ jj+1 ][ ii   ] += m10 = mm * (1-px) * (  py);
      new_m [ jj+1 ][ ii+1 ] += 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   ][ ii+1 ] += dx * m01;  //  mass it is bringing into each square
      new_dx[ jj+1 ][ ii   ] += dx * m10;
      new_dx[ jj+1 ][ ii+1 ] += dx * m00;

      new_dy[ jj   ][ ii   ] += dy * m11;
      new_dy[ jj   ][ ii+1 ] += dy * m01;
      new_dy[ jj+1 ][ ii   ] += dy * m10;
      new_dy[ jj+1 ][ ii+1 ] += dy * m00;
    }
    
    //  Scale down the velocities (the new velocity is a weiDIMted average of the velocities passed in)
    
    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];
      }
    }    


  }

















  //--- 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]  = (rand()%10);
      dx0[j][i] = dx1[j][i] = (rand()%7)-3;
      dy0[j][i] = dy1[j][i] = (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++ )
	    {
	      if( x != DIM || y != DIM )  //  was && - that's wrong!!
	      {
	        dx = (x-DIM);
	        dy = (y-DIM);
	        
	        d2 = dx*dx + dy*dy;   //  Distance squared
	
	        d = sqrt( d2 );  //  Distance

          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

          //- Do speed of light lookup

          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;
	      }
	    }
	  }

  }











//---- 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;

	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;

      // 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;

    }

	}

  //- 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

  frame ++;
  time ++;



  //- Return buffer for displaying

  return &buf;
}

















/*
 * 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;

  static int lmx=0, lmy=0;

	/* 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;
}