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

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

#include "Pattern.h"



/*

  pprizm (prizm ep) - i miss u already (benbecular label, scotland, techno)


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



  woohoo - looks like i can do gravity using fft


*/





//---- 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      10
#define AMBIENT_STUFF  0
#define SOL 200000
//#define WEAKER_FORCE
#define DIM 64
#define VELOCITY_SCALE 1
//#define DONT_WRAP_UNVIERSE
//#define SPEED_LIMIT 5.0

//#define SWAP_DX_AND_DY

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


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

  double dirx_lookup_small[ DIM ][ DIM ];
  double diry_lookup_small[ DIM ][ DIM ];


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






















  //---- do_convolution

  //  Convolutes the image with the convolution and stores the result in dest
  //  dest may point to the same data as image

  fftw_real u[DIM * 2*(DIM/2+1)];
  fftw_real v[DIM * 2*(DIM/2+1)];

  fftw_real u0[DIM * 2*(DIM/2+1)];
  fftw_real v0[DIM * 2*(DIM/2+1)];  //  dirx convolution fft
  fftw_real v1[DIM * 2*(DIM/2+1)];  //  diry convolution fft

  fftw_real rho[DIM * DIM];
  fftw_real rho0[DIM * DIM];  /* Smoke density */

  fftw_real u_u0[DIM * DIM];
  fftw_real u_v0[DIM * DIM];  /* User-induced forces */

  static rfftwnd_plan plan_rc, plan_cr;







  //---- init_FFT

  //  Sets up the FFT engine

  void init_FFT(int n)
  {
	  plan_rc = rfftw2d_create_plan(n, n, FFTW_REAL_TO_COMPLEX, FFTW_IN_PLACE);
	  plan_cr = rfftw2d_create_plan(n, n, FFTW_COMPLEX_TO_REAL, FFTW_IN_PLACE);
  }



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




  void do_convolution( double* image, double* dest, int n, fftw_real* v0 )
  {
    int i, j;

    //- FFT of original: image -> u0

	  for ( i=0 ; i<n; i++ )
    {
		  for ( j=0 ; j<n; j++ )
      {
			  u0[i+(n+2)*j] = image[ i+n*j ]; 
		  }
    }

    FFT(1,u0);

    //- Convolution of kernel is passed in

    //  in v0 (from the outside world it is v0 or v1)

    //- Multiply the bastards u0 * v0 -> u0

    fftw_real x,y, a,b;

	  for ( i=0 ; i<=(n/2); i++ )
    {
		  for ( j=0 ; j<n; j++ )
      {
        //  Complex multiplication: (a+ib)(x+iy) = ax-by + i( bx+ay )

        a = u0[i*2  +(n+2)*j];
        b = u0[i*2+1+(n+2)*j];

        x = v0[i*2  +(n+2)*j];
        y = v0[i*2+1+(n+2)*j];


        a = a*x - b*y;
        b = b*x + a*y;

//        a *= x;
//        b *= y;


        u0[i*2  +(n+2)*j] = a;
        u0[i*2+1+(n+2)*j] = b;
		  }    
	  }

    //- Inverse FFT of result -> dest

    FFT(-1,u0);

	  double f = 1.0 / (n*n+0.0);

 	  for ( i=0 ; i<n; i++ )
    {
		  for ( j=0 ; j<n; j++ )
      {
			  dest[j*n+i] = f * u0[j+(n+2)*i];  //  u0 : don't ? swap i and j
		  }
    }
    
  }






  void do_convolution2( double* image, double* convolution, double* dest, int n )
  {
    int i, j, x, y;

    for( i=0; i<n; i++ )
      for( j=0; j<n; j++ )
        dest[j*n+i] = 0;

    for( i=0; i<n; i++ )
      for( j=0; j<n; j++ )
        for( x=0; x<n; x++ )
          for( y=0; y<n; y++ )
            dest[j*n+i] += image[((j+y)&(n-1))*n+((i+x)&(n-1))] * convolution[y*n+x];

  }



















  //--- Startup function

  void Pattern_Startup()
  {
    int i, j;

    init_FFT( DIM );

    //  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]  = 0;//(rand()%10);
      dx0[j][i]     = dx1[j][i]      = 0;//(rand()%7)-3;
      dy0[j][i]     = dy1[j][i]      = 0;//(rand()%7)-3;
    }

    for(j=DIM/2-10; j<DIM/2+10; j++) for(i=DIM/2-10; i<DIM/2+10; 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]      = 0;//1*((rand()%7)-3);
      dy0[j][i]     = dy1[j][i]      = 0;//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;
        diry_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;
	      }

        //  scale by gravity strength

        dirx_lookup[y][x] *= GRAVITY_STRENGTH;
        diry_lookup[y][x] *= GRAVITY_STRENGTH;
	    }
	  }

	  for( y=0; y<DIM; y++ )
	  {
	    for( x=0; x<DIM; x++ )
	    {
        dirx_lookup_small [ y ][ x ] = dirx_lookup[y+DIM][x+DIM];
        diry_lookup_small [ y ][ x ] = diry_lookup[y+DIM][x+DIM];
      }
    }





    //- FFT of convolution: convolution x -> v0

    int n = DIM;

	  for ( i=0 ; i<n; i++ )
    {
		  for ( j=0 ; j<n; j++ )
      {
			  v0[i+(n+2)*j] = dirx_lookup_small[i][j];  //  swap i and j
		  }
    }

    FFT(1,v0);

    //- FFT of convolution: convolution y -> v1

	  for ( i=0 ; i<n; i++ )
    {
		  for ( j=0 ; j<n; j++ )
      {
			  v1[i+(n+2)*j] = diry_lookup_small[i][j];  //  swap i and j
		  }
    }

    FFT(1,v1);


  }












  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)&(PASTLEN-1) ][ y ][ x ];    //  Get mass from past_i
          past_mbuf1 = mbuf[ (frame-past_i-1)&(PASTLEN-1) ][ 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 += dirx_lookup_small [ ((DIM+y-j)%DIM) ][ ((DIM+x-i)%DIM) ] * m;
          py += diry_lookup_small [ ((DIM+y-j)%DIM) ][ ((DIM+x-i)%DIM) ] * m;
        }
      }

//      px *= GRAVITY_STRENGTH;
//      py *= GRAVITY_STRENGTH;
      
      last_px[j][i] = px;
      last_py[j][i] = py;
    }

*/
    

    do_convolution( (double*)mbuf[ (frame-past_i)&(PASTLEN-1) ], (double*)last_px, DIM, v0 );
    do_convolution( (double*)mbuf[ (frame-past_i)&(PASTLEN-1) ], (double*)last_py, DIM, v1 );


    //  Update V - acceleration of pixel is force (px,py) / mass

    for(j=0; j<DIM; j++) for(i=0; i<DIM; i++)
    {
      if( current_m[j][i] > 0 )
      {
        px = last_px[j][i];
        py = last_py[j][i];

        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  )&(PASTLEN-1) ];
  new_m     = mbuf[ (frame+1)&(PASTLEN-1) ];



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

      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;

			buf.buf[j][i][0] = bb;
			buf.buf[j][i][1] = gg;
			buf.buf[j][i][2] = rr;



/*
        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] * 1280.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] * 1280.0 + 128.0);

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

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


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

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