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

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

#include "Pattern.h"



/*


003:

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

004:

  use fft image convolution to calculate gravity (avoids n^4 calc and reduces it to n^2)



*/



static rfftwnd_plan plan_rc, plan_cr;

#define FFT(s,u)\
	if(s==1) rfftwnd_one_real_to_complex(plan_rc,(fftw_real *)u,(fftw_complex*)u);\
	else rfftwnd_one_complex_to_real(plan_cr,(fftw_complex *)u,(fftw_real *)u)





//---- Render buffer

PATTERN_RGB buf;








#define DIM 64

#define DIM2 (DIM*2)







  //---- Galaxy data

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

  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[ DIM2 ][ DIM2 ];
  double diry_lookup[ DIM2 ][ DIM2 ];

  int frame = 0;
  int time = 0;



//---- double buffer pointers

    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






//---- FFT stuff



fftw_real fft_dx[ 4 * DIM2 * DIM2 ];
fftw_real fft_dy[ 4 * DIM2 * DIM2 ];

fftw_real fft_m [ 4 * DIM2 * DIM2 ];




/*

setup:

  fft_dx <- dirx
  fft_dy <- diry

  FFT( fft_dx, 1 )
  FFT( fft_dy, 1 )



convolution:

  fft_m <- m

  FFT( fft_m, 1 )

  px = fft_mm * fft_dx
  py = fft_mm * fft_dy

  FFT( px, -1 )
  FFT( py, -1 )




*/








  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;

    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++ )
        {
          px += dirx_lookup[DIM+y-j][DIM+x-i] * current_m[y][x];
          py += diry_lookup[DIM+y-j][DIM+x-i] * current_m[y][x];
        }
      }

      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;

    //- Set up FFT plans

  	plan_rc = rfftw2d_create_plan( DIM2, DIM2, FFTW_REAL_TO_COMPLEX, FFTW_IN_PLACE );
	  plan_cr = rfftw2d_create_plan( DIM2, DIM2, FFTW_COMPLEX_TO_REAL, FFTW_IN_PLACE );

    //- 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++)
    {
      m0[j][i]  = m1[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<DIM2; y++ )
	  {
	    for( x=0; x<DIM2; 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 ] = 5.0 * dx / d / d2;
	        diry_lookup[ y ][ x ] = 5.0 * dy / d / d2;
	      }
	      else  //  Special case for centre point
	      {
	        dirx_lookup[ y ][ x ] = 0;
	        diry_lookup[ y ][ x ] = 0;
	      }
	    }
	  }



    //- Set up FFT for convolution

   	for( i=0 ; i<DIM2 ; i++ )
    {
	  	for( j=0 ; j<DIM2 ; j++ )
      {
		  	fft_dx[ i + (DIM2+2) * j ] = dirx_lookup[j][i]; 
		    fft_dy[ i + (DIM2+2) * j ] = dirx_lookup[j][i];
      }
    }

    rfftwnd_one_real_to_complex( plan_rc, (fftw_real *)fft_dx, (fftw_complex*)fft_dx );
    rfftwnd_one_real_to_complex( plan_rc, (fftw_real *)fft_dy, (fftw_complex*)fft_dy );

  }






















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



PATTERN_RGB *Pattern_Generate()
{
    //- Decide which buffer is primary and which is secondary

    if(frame == 0)
    {
      current_m  =  m0; new_m  =  m1;
      current_dx = dx0; new_dx = dx1;
      current_dy = dy0; new_dy = dy1;
    }
    else
    {
      current_m  =  m1; new_m  =  m0;
      current_dx = dx1; new_dx = dx0;
      current_dy = dy1; new_dy = dy0;
    }


  //- Animate

  recalc();



  //- Render

  int i, j, idx, c;

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

      c = (int)(current_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;

    }

	}


  int cc;

	for (j = 0; j < DIM2; j++)
  {
		for (i = 0; i < DIM2; i+=2 )
    {
      //  mass

      c = (int)( 128.0 + fft_dx[ i*2 + (DIM2+2)*j ] * 128.0 );

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

      cc = (int)( 128.0 + fft_dy[ i*2 + 1 + (DIM2+2)*j ] * 128.0 );

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

			buf.buf[ (j+DIM) % DIM2 ][ (i+DIM) % DIM2 ][0] = c;
			buf.buf[ (j+DIM) % DIM2 ][ (i+DIM) % DIM2 ][1] = cc;
			buf.buf[ (j+DIM) % DIM2 ][ (i+DIM) % DIM2 ][2] = 0;
    }
	}


/*

	for (j = 0; j < DIM2; j++)
  {
		for (i = 0; i < DIM2; i++ )
    {
      //  mass

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

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

  */

      //- Move on

    frame ^= 1;
    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;
}