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

  woohoo - looks like i can do gravity using fft

  I've had to negate the gravity strength and swap the x and the y lookups for some unknown reason

*/





//---- Render buffer

PATTERN_RGB buf;





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


int offset_x=0, offset_y=0;


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


  //  TODO: Palette

  int mass_to_rgb_r[ 100000 ];
  int mass_to_rgb_g[ 100000 ];
  int mass_to_rgb_b[ 100000 ];

  int mouse_x=0, mouse_y=0;
  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; i+=2 )
    {
		  for ( j=0; j<n; j++ )
      {
        //  Complex multiplication: (a+ib)(x+iy) = ax-by + i( bx+ay )

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

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

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

    //- 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[i+n*j] = f * u0[i+(n+2)*j];  //  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/3; j<2*DIM/3; j++) for(i=DIM/3; i<2*DIM/3; i++)
    for(j=0; j<DIM; j++) for(i=0; i<DIM; i++)
    {
      mbuf[0][j][i] = mbuf[1][j][i]  = AMBIENT_STUFF + 0.1*(rand()%((int)(10.0*(STUFF_MAX-AMBIENT_STUFF))));
      dx0[j][i]     = dx1[j][i]      = 0*((rand()%7)-3);
      dy0[j][i]     = dy1[j][i]      = 0*((rand()%7)-3);
    }
    



    //- Gravity precalc

	  int x, y;
	
	  //  calculate the distance from the centre point (SIZE,SIZE) in (SIZE*2,SIZE*2) array


    //  If the file dirxDIM.dat exists use that, otherwise calculate the lookup table and store in a dat file
    FILE *fp;
    char fname1[257];
    char fname2[257];
    sprintf( fname1, "dirx%d.dat", DIM );
    sprintf( fname2, "diry%d.dat", DIM );

    if( fp = fopen( fname1, "rb" ) )
    {
      fread( (void*)dirx_lookup_small, 1, sizeof( dirx_lookup_small ), fp );
      fclose( fp );

      fp = fopen( fname2, "rb" );
      fread( (void*)diry_lookup_small, 1, sizeof( diry_lookup_small ), fp );
      fclose( fp );
    }
    else
    {
	    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!!
//          if( (x % DIM) || (y % DIM) )  //  was && - that's wrong!!
	        {
            i=j=0; // clear i and j in case the for loops are commented out


            for( j=-160; j<=160; j++ ) for( i=-160; i<=160; 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

          //  NEGATIVE - I don't know why...

          dirx_lookup[y][x] *= -1;//GRAVITY_STRENGTH;
          diry_lookup[y][x] *= -1;//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];

  //        buf.buf[y][DIM+x][2] = dirx_lookup_small [ y ][ x ];
  //        buf.buf[y][DIM+x][1] = diry_lookup_small [ y ][ x ];
        }
      }


      fp = fopen( fname1, "wb" );
      fwrite( (void*)dirx_lookup_small, 1, sizeof( dirx_lookup_small ), fp );
      fclose( fp );

      fp = fopen( fname2, "wb" );
      fwrite( (void*)diry_lookup_small, 1, sizeof( diry_lookup_small ), fp );
      fclose( fp );


    }





      //  scale by gravity strength
      //  done seperately to the rest of the calculations
      //  to allow the things to be loaded from file or calculated

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


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




    //  palette

    double rr, gg, bb, r;

    for( i=0; i<100000; i++ )
    {

      double r = 1.0 - exp(-(i*0.01)*log(2.0));

//      if( j < 4 ) r = (i+0.0)/(DIM+0.0);  //  put the palette

      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*0.5 );

      bb += r*255;

      if( r > 0.1 )
      {
        bb += (0.1) * 254.99;
        rr += (0.1) * 128;
      }
      else if( r > 0.05 )
      {
        bb += (r-0.05) * 254.99;
        rr += (r-0.05) * 128;
      }

      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;

      if( r > 0.9 && r < 0.95 ) bb -= 100;  //  hint of yellow near the end of the spectrum

      if( bb < 0 ) bb = 0;

      mass_to_rgb_r[i] = rr;
      mass_to_rgb_g[i] = gg;
      mass_to_rgb_b[i] = bb;
    }



    //  TEMP: Exit after creating lookup table
//    exit( 1 );

  }












  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;
    


    int x, y;


    //  Decrease gravitational force over time if necessary

#ifdef GRAVITY_CHANGE

    for( y=0; y<DIM; y++ )
	  {
	    for( x=0; x<DIM; x++ )
	    {
        dirx_lookup_small [ y ][ x ] *= GRAVITY_CHANGE;
        diry_lookup_small [ y ][ x ] *= GRAVITY_CHANGE;

//        buf.buf[y][DIM+x][2] = dirx_lookup_small [ y ][ x ];
//        buf.buf[y][DIM+x][1] = diry_lookup_small [ y ][ x ];
      }
    }

    //- Have to redo the ffts because of this

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

#endif







    //- 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    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, v1 );
    do_convolution( (double*)mbuf[ (frame-past_i)&(PASTLEN-1) ], (double*)last_py, DIM, v0 );


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

      //- New in 018 - universal expansion

#ifdef EXPANSION
      newi = (newi - DIM/2) * EXPANSION + DIM/2;
      newj = (newj - DIM/2) * EXPANSION + DIM/2;
#endif

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

/*
      if( newi < 0    ||
       newj < 0    ||
       newi > DIM-1 ||
         newj > DIM-1 ) { newi = rand()%DIM; newj = rand()%DIM; }
*/

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

#ifdef VELOCITY_DAMPING
      current_dx[j][i] *= VELOCITY_DAMPING;
      current_dy[j][i] *= VELOCITY_DAMPING;
#endif

      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];
      }
      else
      {
        new_dx[j][i] = 0;
        new_dy[j][i] = 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  )&(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));

//      if( j < 4 ) r = (i+0.0)/(DIM+0.0);  //  put the palette

      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*0.5 );

      bb += r*255;

      if( r > 0.1 )
      {
        bb += (0.1) * 256;
        rr += (0.1) * 128;
      }
      else if( r > 0.05 )
      {
        bb += (r-0.05) * 256;
        rr += (r-0.05) * 128;
      }

        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;

      if( r > 0.9 && r < 0.95 ) bb -= 100;  //  hint of yellow near the end of the spectrum

			buf.buf[j][i][0] = buf.buf[j][i][0]*0.25 + 0.75*bb;
			buf.buf[j][i][1] = buf.buf[j][i][1]*0.25 + 0.75*gg;
			buf.buf[j][i][2] = buf.buf[j][i][2]*0.25 + 0.75*rr;
*/

/*
      double r = new_m[j][i] * 15.0;
      if( r > 255 ) r = 255;
*/
      double mm = new_m[(j-offset_y)&(DIM-1)][(i-offset_x)&(DIM-1)] * 100.0;

      if( mm > 99999 ) mm = 99999;
      int r = mm;

			buf.buf[j][i][0] = mass_to_rgb_b[r];
			buf.buf[j][i][1] = mass_to_rgb_g[r];
			buf.buf[j][i][2] = mass_to_rgb_r[r];

/*
//  tile
      buf.buf[j][i+DIM][0] = buf.buf[j][i][0];
			buf.buf[j][i+DIM][1] = buf.buf[j][i][1];
			buf.buf[j][i+DIM][2] = buf.buf[j][i][2];

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



/*
        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)((new_dx[j][i] * 1280.0 * r + 128.0)*r);

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

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

      
//			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)((new_dy[j][i] * 1280.0 * r + 128.0)*r);

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

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

            buf.buf[j+DIM][i][0] = 255*r;


//			buf.buf[j+DIM][i+DIM][0] = c;
//			buf.buf[j+DIM][i+DIM][1] = c;
//			buf.buf[j+DIM][i+DIM][2] = c;
*/
      // x force
/*
      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][2] = c;

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

      //  y force

      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)
{
  offset_x += mx-mouse_x;
  offset_y += my-mouse_y;

  mouse_x = mx;
  mouse_y = 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;
}