seedpoints

6 files changed, 1156 insertions, 1123 deletions

diff --git a/Smoke/Week 2.ncb b/Smoke/Week 2.ncb
index a105095..73d9e5f 100644
--- a/Smoke/Week 2.ncb
+++ b/Smoke/Week 2.ncb
diff --git a/Smoke/Week 2.suo b/Smoke/Week 2.suo
index cda6f03..3b8e114 100644
--- a/Smoke/Week 2.suo
+++ b/Smoke/Week 2.suo
diff --git a/Smoke/Week 2.vcproj b/Smoke/Week 2.vcproj
index 1c78012..e0d727c 100644
--- a/Smoke/Week 2.vcproj
+++ b/Smoke/Week 2.vcproj
@@ -739,6 +739,10 @@
 				>

 			</File>

 			<File

+				RelativePath=".\seedpoint.c"

+				>

+			</File>

+			<File

 				RelativePath="fftw-2.1.3\fftw\timer.c"

 				>

 			</File>

@@ -783,6 +787,10 @@
 				RelativePath="fftw-2.1.3\rfftw\rfftw.h"

 				>

 			</File>

+			<File

+				RelativePath=".\seedpoint.h"

+				>

+			</File>

 		</Filter>

 		<Filter

 			Name="Resource Files"

diff --git a/Smoke/fluids.c b/Smoke/fluids.c
index f249425..b40b507 100644
--- a/Smoke/fluids.c
+++ b/Smoke/fluids.c
@@ -1,1122 +1,1125 @@
-// Usage: Drag with the mouse to add smoke to the fluid. This will also move a "rotor" that disturbs 

-//        the velocity field at the mouse location. Press the indicated keys to change options

-//-------------------------------------------------------------------------------------------------- 

-

-#define max(a, b) (a > b ? a : b)

-#define min(a, b) (a < b ? a : b)

-#define round(x) (int)(x < 0 ? x - 0.5 : x + 0.5)

-

-#define MAX_INFINITE (1.#INF000)

-#define MIN_INFINITE (-1.#INF000)

-

-#define vec_len(x,y) (sqrt((x*x) + (y*y)))

-

-#define FALSE 0

-#define TRUE !FALSE

-

-#ifdef G_OS_WIN32

-#define WIN32_LEAN_AND_MEAN 1

-#include <windows.h>

-#endif

-

-#include <math.h>

-#include <rfftw.h>              //the numerical simulation FFTW library

-#include <stdio.h>              //for printing the help text

-

-#include <GL/gl.h>

-#include <GL/glu.h>

-

-#include "fluids.h"

-

-//--- SIMULATION PARAMETERS ------------------------------------------------------------------------

-const int DIM = 50;             //size of simulation grid

-int var_dims = 50;

-double dt = 0.4;				        //simulation time step

-float visc = 0.001f;				    //fluid viscosity

-fftw_real *vx, *vy;             //(vx,vy)   = velocity field at the current moment

-fftw_real *vx0, *vy0;           //(vx0,vy0) = velocity field at the previous moment

-fftw_real *fx, *fy;	            //(fx,fy)   = user-controlled simulation forces, steered with the mouse 

-fftw_real *rho, *rho0;		    	//smoke density at the current (rho) and previous (rho0) moment 

-fftw_real *hight_array;         //used for hight plot

-int *frame_hist;

-int frame_index = 0;

-rfftwnd_plan plan_rc, plan_cr;  //simulation domain discretization

-

-int   winWidth, winHeight;      //size of the graphics window, in pixels

-int   color_dir = 0;            //use direction color-coding or not

-float vec_scale = 1000;			    //scaling of hedgehogs

-int   draw_smoke = 0;           //draw the smoke or not

-int   draw_vecs = 1;            //draw the vector field or not

-int   scalar_col = 0;           //method for scalar coloring

-int   frozen = 0;               //toggles on/off the animation

-int   olivers_color = 256;

-float clamp_min = 0;

-float clamp_max = 0.9999f;

-int autoscale = FALSE;

-float scale_min = 0;

-float scale_max = 1;

-int vis_dataset = DATASET_RHO;

-int mousebutton;

-int mousebuttonstate;

-int active_slider = 0;

-int use_glyphs = GLYPHS_OFF;

-int glyph_scalar = SCALAR_RHO;

-int glyph_vector = VECTOR_VEL;

-int glyph_sort = GLYPH_CYLINDERS;

-int draw_options = FALSE;

-GLuint startList;

-float threshold1 = 1.8f;

-float threshold2 = 2.0f;

-int isolines_nr = 1;

-int prev_mx = 0;

-int prev_my = 0;

-

-int rotate_plane = 0;

-float x_rot = 0.0f;

-float y_rot = 0.0f;

-float z_rot = 0.0f;

-float x_pos = 0.0f;

-float y_pos = 0.0f;

-float z_pos = -1000.0f;

-

-struct color4f {

-  float r;

-  float g;

-  float b;

-  float a;

-};

-

-struct point {

-  float x;

-  float y;

-  float z;

-};

-

-//------ SIMULATION CODE STARTS HERE -----------------------------------------------------------------

-

-//init_simulation: Initialize simulation data structures as a function of the grid size 'n'. 

-//                 Although the simulation takes place on a 2D grid, we allocate all data structures as 1D arrays,

-//                 for compatibility with the FFTW numerical library.

-void init_simulation(int n)				

-{

-  float LightAmbient[]  = { 0.0f, 0.0f, 0.0f, 2.0f }; // Ambient light values

-  float LightPosition[] = { 0.0f, 0.0f, -900.0f, 1.0f }; // Position of the light source

-	int i; size_t dim;

-  GLUquadricObj *qobj;

-	

-	dim     = n * 2*(n/2+1)*sizeof(fftw_real);        //Allocate data structures

-	vx       = (fftw_real*) malloc(dim); 

-	vy       = (fftw_real*) malloc(dim);

-	vx0      = (fftw_real*) malloc(dim); 

-	vy0      = (fftw_real*) malloc(dim);

-	dim     = n * n * sizeof(fftw_real);

-	fx      = (fftw_real*) malloc(dim);

-	fy      = (fftw_real*) malloc(dim);

-	rho     = (fftw_real*) malloc(dim); 

-	rho0    = (fftw_real*) malloc(dim);

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

-  hight_array = (fftw_real*) malloc(dim);

-  frame_hist = (int*) malloc(n * sizeof(int *));

-

-  for (i = 0; i <= n; i++)

-  {

-      frame_hist[i] = (fftw_real*)malloc(n * 2*(n/2+1)*sizeof(fftw_real));

-  }

-	

-	for (i = 0; i < n * n; i++)                      //Initialize data structures to 0

-	{ vx[i] = vy[i] = vx0[i] = vy0[i] = fx[i] = fy[i] = rho[i] = rho0[i] = hight_array[i] = 0.0f; }

-

-  glShadeModel(GL_SMOOTH);                           // Enable smooth shading

-  glClearColor(0.0f, 0.0f, 0.0f, 0.0f);              // Black background

-  glClearDepth(1.0f);                                // Depth buffer setup

-  glDepthFunc(GL_LESS);                              // The type of depth testing to do

-  glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST); // Really nice perspective calculations

-  glHint(GL_POINT_SMOOTH_HINT, GL_NICEST);           // Really nice point smoothing

-  glEnable(GL_BLEND);

-  glBlendFunc(GL_SRC_ALPHA, GL_ONE);

-  glEnable(GL_DEPTH_TEST);

-  //glEnable(GL_LIGHTING);

-  glEnable(GL_LIGHT0);

-  glEnable(GL_COLOR_MATERIAL);

-

-  glLightfv(GL_LIGHT0, GL_AMBIENT, LightAmbient);

-  glLightfv(GL_LIGHT0, GL_POSITION, LightPosition);

-

-  startList = glGenLists(2);

-  qobj = gluNewQuadric();

-

-  gluQuadricDrawStyle(qobj, GLU_FILL);

-  gluQuadricNormals(qobj, GLU_SMOOTH);

-  glNewList(startList, GL_COMPILE);

-    gluCylinder(qobj, 10, 0.2, 8, 8, 8);

-  glEndList();

-}

-

-int rescale_to_winwidth(float value)

-{

-	return round(value *winWidth);

-}

-

-//FFT: Execute the Fast Fourier Transform on the dataset 'vx'.

-//     'dirfection' indicates if we do the direct (1) or inverse (-1) Fourier Transform

-void FFT(int direction,void* vx)

-{

-	if(direction==1) rfftwnd_one_real_to_complex(plan_rc,(fftw_real*)vx,(fftw_complex*)vx);

-	else             rfftwnd_one_complex_to_real(plan_cr,(fftw_complex*)vx,(fftw_real*)vx);

-}

-

-int clamp(float x) 

-{ return ((x)>=0.0?((int)(x)):(-((int)(1-(x))))); }

-

-//solve: Solve (compute) one step of the fluid flow simulation

-void solve(int n, fftw_real* vx, fftw_real* vy, fftw_real* vx0, fftw_real* vy0, fftw_real visc, fftw_real dt) 

-{

-	fftw_real x, y, x0, y0, f, r, U[2], V[2], s, t;

-	int i, j, i0, j0, i1, j1;

-

-	for (i=0;i<n*n;i++) 

-	{ vx[i] += dt*vx0[i]; vx0[i] = vx[i]; vy[i] += dt*vy0[i]; vy0[i] = vy[i]; }    

-

-	for ( x=0.5f/n,i=0 ; i<n ; i++,x+=1.0f/n ) 

-	   for ( y=0.5f/n,j=0 ; j<n ; j++,y+=1.0f/n ) 

-	   {

-	      x0 = n*(x-dt*vx0[i+n*j])-0.5f; 

-	      y0 = n*(y-dt*vy0[i+n*j])-0.5f;

-	      i0 = clamp((float)x0); s = x0-i0;

-	      i0 = (n+(i0%n))%n;

-	      i1 = (i0+1)%n;

-	      j0 = clamp((float)y0); t = y0-j0;

-	      j0 = (n+(j0%n))%n;

-	      j1 = (j0+1)%n;

-	      vx[i+n*j] = (1-s)*((1-t)*vx0[i0+n*j0]+t*vx0[i0+n*j1])+s*((1-t)*vx0[i1+n*j0]+t*vx0[i1+n*j1]);

-	      vy[i+n*j] = (1-s)*((1-t)*vy0[i0+n*j0]+t*vy0[i0+n*j1])+s*((1-t)*vy0[i1+n*j0]+t*vy0[i1+n*j1]);

-	   }     

-	

-	for(i=0; i<n; i++)

-	  for(j=0; j<n; j++) 

-	  {  vx0[i+(n+2)*j] = vx[i+n*j]; vy0[i+(n+2)*j] = vy[i+n*j]; }

-

-	FFT(1,vx0);

-	FFT(1,vy0);

-

-	for (i=0;i<=n;i+=2) 

-	{

-	   x = 0.5f*i;

-	   for (j=0;j<n;j++) 

-	   {

-	      y = j<=n/2 ? (fftw_real)j : (fftw_real)j-n;

-	      r = x*x+y*y;

-	      if ( r==0.0f ) continue;

-	      f = (fftw_real)exp(-r*dt*visc);

-	      U[0] = vx0[i  +(n+2)*j]; V[0] = vy0[i  +(n+2)*j];

-	      U[1] = vx0[i+1+(n+2)*j]; V[1] = vy0[i+1+(n+2)*j];

-

-	      vx0[i  +(n+2)*j] = f*((1-x*x/r)*U[0]     -x*y/r *V[0]);

-	      vx0[i+1+(n+2)*j] = f*((1-x*x/r)*U[1]     -x*y/r *V[1]);

-	      vy0[i+  (n+2)*j] = f*(  -y*x/r *U[0] + (1-y*y/r)*V[0]);

-	      vy0[i+1+(n+2)*j] = f*(  -y*x/r *U[1] + (1-y*y/r)*V[1]);

-	   }

-	}

-

-	FFT(-1,vx0); 

-	FFT(-1,vy0);

-

-	f = 1.0/(n*n);

- 	for (i=0;i<n;i++)

-	   for (j=0;j<n;j++) 

-	   { vx[i+n*j] = f*vx0[i+(n+2)*j]; vy[i+n*j] = f*vy0[i+(n+2)*j]; }

-} 

-

-

-// diffuse_matter: This function diffuses matter that has been placed in the velocity field. It's almost identical to the

-// velocity diffusion step in the function above. The input matter densities are in rho0 and the result is written into rho.

-void diffuse_matter(int n, fftw_real *vx, fftw_real *vy, fftw_real *rho, fftw_real *rho0, fftw_real dt) 

-{

-	fftw_real x, y, x0, y0, s, t;

-	int i, j, i0, j0, i1, j1;

-

-	for ( x=0.5f/n,i=0 ; i<n ; i++,x+=1.0f/n )

-		for ( y=0.5f/n,j=0 ; j<n ; j++,y+=1.0f/n ) 

-		{

-			x0 = n*(x-dt*vx[i+n*j])-0.5f; 

-			y0 = n*(y-dt*vy[i+n*j])-0.5f;

-			i0 = clamp((float)x0);

-			s = x0-i0;

-			i0 = (n+(i0%n))%n;

-			i1 = (i0+1)%n;

-			j0 = clamp((float)y0);

-			t = y0-j0;

-			j0 = (n+(j0%n))%n;

-			j1 = (j0+1)%n;

-			rho[i+n*j] = (1-s)*((1-t)*rho0[i0+n*j0]+t*rho0[i0+n*j1])+s*((1-t)*rho0[i1+n*j0]+t*rho0[i1+n*j1]);

-		}    

-}

-

-//set_forces: copy user-controlled forces to the force vectors that are sent to the solver. 

-//            Also dampen forces and matter density to get a stable simulation.

-void set_forces(void) 

-{

-	int i;

-	for (i = 0; i < DIM * DIM; i++) 

-	{

-        rho0[i]  = 0.995 * rho[i];

-        fx[i] *= 0.85; 

-        fy[i] *= 0.85;

-        vx0[i]    = fx[i]; 

-        vy0[i]    = fy[i];

-	}

-}

-

-

-void calculate_hight_plot(void)

-{

-  int i;

-	for (i = 0; i < DIM * DIM; i++) 

-	{

-    switch (vis_dataset)

-    {

-		  case DATASET_FORCE:

-			  hight_array[i] = sqrt(fx[i] * fx[i] + fy[i] * fy[i]) * 100;

-			  break;

-		  case DATASET_VEL:

-			  hight_array[i] = sqrt(vx[i] * vx[i] + vy[i] * vy[i]) * 2000;

-			  break;

-		  case DATASET_RHO:

-			  hight_array[i] = rho[i] * 10;

-			  break;

-      default:

-		  case DATASET_DIVV:

-		  case DATASET_DIVF:

-        hight_array[i] = 0.0f;

-			  break;

-	  }

-  }

-}

-

-void copy_frame(void)

-{

-  switch (vis_dataset)

-  {

-	  case DATASET_FORCE:

-		//  hight_array[i] = sqrt(fx[i] * fx[i] + fy[i] * fy[i]) * 100;

-		  break;

-	  case DATASET_VEL:

-		 // hight_array[i] = sqrt(vx[i] * vx[i] + vy[i] * vy[i]) * 2000;

-		  break;

-	  case DATASET_RHO:

-		  // hight_array[i] = rho[i] * 10;

-      memcpy(frame_hist[frame_index], rho, sizeof(DIM * 2 * (DIM / 2 + 1) * sizeof(fftw_real)));

-      break;

-    default:

-	  case DATASET_DIVV:

-	  case DATASET_DIVF:

-     // hight_array[i] = 0.0f;

-		  break;

-  }

-  frame_index = (frame_index + 1) % DIM;

-}

-

-//do_one_simulation_step: Do one complete cycle of the simulation:

-//      - set_forces:       

-//      - solve:            read forces from the user

-//      - diffuse_matter:   compute a new set of velocities

-//      - gluPostRedisplay: draw a new visualization frame

-

-

-void calculate_one_simulation_step(void)

-{

-    set_forces();

-	  solve(DIM, vx, vy, vx0, vy0, visc, dt);

-	  diffuse_matter(DIM, vx, vy, rho, rho0, dt);

-    calculate_hight_plot();

-    copy_frame();

-}

-

-

-//------ VISUALIZATION CODE STARTS HERE -----------------------------------------------------------------

-

-

-//rainbow: Implements a color palette, mapping the scalar 'value' to a rainbow color RGB

-void rainbow(float value,float* R,float* G,float* B)

-{                          

-   const float dx=0.8f; 

-   if (value<0) value=0; if (value>1) value=1;

-   value = (6-2*dx)*value+dx;

-   *R = (float)max(0.0f,(3-fabs(value-4.0f)-fabs(value-5.0f))/2.0f);

-   *G = (float)max(0.0f,(4-fabs(value-2.0f)-fabs(value-4.0f))/2.0f);

-   *B = (float)max(0.0f,(3-fabs(value-1.0f)-fabs(value-2.0f))/2.0f);

-}

-

-

-void colormap_fire(float value,float* R,float* G,float* B)

-{

-	/* Colormap Fire

-	 * A fire effect deals with two parts, first a drop from red to yellow (halfway)

-	 * during which time the Red component remains full e.g. 1. The Green component is

-	 * slowly added to turn the red into orange, and then yellow (R & G =Y). After this

-	 * point, the Red and Green component (e.g. yellow) have to drop simulataniously

-	 * to go from yellow down to black.

-	 */

-	*B = 0;

-	*G = 0;

-	*R = 0;

-

-	if (value <= (0.01)) {

-		/* whilst value is 0 - 0.5 both red and green equally change to create yellow*/

-		*R = *G = value;

-	} else {

-		/* whilst value is 0.5 - 1 Red is always fully on while the Green component is

-		 * added in steps to go from red to orange to yellow.

-		 */

-		*G = 0.9f - value;

-		*R = 0.8f; // not 1, makes red deeper, more intense

-	}

-}

-

-float remap(float value)

-{

-	value -= scale_min;

-	value /= (scale_max - scale_min);

-

-	//value = (value < 0) ? 0 : (value > 1) ? 1 : value;

-

-	return value;

-}

-

-//set_colormap: Sets three different types of colormaps

-struct color4f set_colormap(float vy, int draw_bar, float alpha)

-{

-   float R, G, B;

-   struct color4f return_value;

-

-   if (autoscale)

-   {

-	   vy = remap(vy);

-   }

-

-   if (!(draw_bar))

-   {

-     if (vy < clamp_min) vy = clamp_min;

-     if (vy > clamp_max) vy = clamp_max;

-   } else {

-     

-   }

-

-   vy *= olivers_color;

-   vy = (float)(int)(vy);

-   vy /= olivers_color;

-

-   if (scalar_col==COLOR_BLUE_GREEN_RED) {

-     if  (vy < -0.1) {

-       R = G = 0;

-       vy -= -0.1;

-       vy /= 0.9;

-       B = -vy;

-     } else if (vy < 0.1) {

-       R = B = 0;

-       vy += 0.1;

-       vy /= 0.2;

-       G = vy;

-     } else {

-       vy -= 0.1;

-       vy /= 0.9;

-       R = vy;

-       G = B = 0;

-     }

-   } else if (scalar_col==COLOR_BLACKWHITE) {

-       R = G = B = vy;

-   }

-   else if (scalar_col==COLOR_WILRIK) {

-	   colormap_fire(vy, &R, &G, &B);

-   }

-   else if (scalar_col==COLOR_OLIVER) {

-	  rainbow(vy,&R,&G,&B);   

-   }

-   else if (scalar_col==COLOR_RAINBOW)

-       rainbow(vy,&R,&G,&B); 

-   else if (scalar_col==COLOR_BANDS)

-       {  

-          const int NLEVELS = 7;

-          vy *= NLEVELS; vy = (float)(int)(vy); vy/= NLEVELS; 

-	      rainbow(vy,&R,&G,&B);   

-	   }

-

-   glColor4f(R,G,B,alpha);

-

-   return_value.r = R;

-   return_value.g = G;

-   return_value.b = B;

-   return_value.a = alpha;

-   return return_value;

-}

-

-

-//direction_to_color: Set the current color by mapping a direction vector (x,y), using

-//                    the color mapping method 'method'. If method==1, map the vector direction

-//                    using a rainbow colormap. If method==0, simply use the white color

-void direction_to_color(float x, float y, int method)

-{

-	float r,g,b,f;

-

-	switch (method) {

-	  case 3:

-		g = 1;

-		break;

-	  case 2:

-		r = 1;

-		g = b = 0;

-		break;

-	  case 1:

-		f = (float)(atan2((double)y, (double)x) / 3.1415927 + 1);

-		r = f;

-		if(r > 1) r = 2 - r;

-	  	g = f + 0.66667f;

-		if(g > 2) g -= 2;

-		if(g > 1) g = 2 - g;

-		b = f + 2.0f * 0.66667f;

-		if(b > 2) b -= 2;

-		if(b > 1) b = 2 - b;

-		break;

-	  case 0:

-	  default: r = g = b = 1;

-	  break;

-	}

-	glColor3f(r,g,b);

-}

-

-float get_dataset(int index)

-{

-  fftw_real cell_x = (fftw_real)winWidth  / (fftw_real)(DIM + 1); // Grid cell width

-	fftw_real cell_y = (fftw_real)winHeight / (fftw_real)(DIM + 1); // Grid cell heigh

-	float return_value, par_der_x, par_der_y;

-  static int klote = 0;

-

-	return_value = 0;

-

-	switch (vis_dataset)

-  {

-		case DATASET_FORCE:

-			return_value = (float)sqrt((fx[index] * fx[index]) + (fy[index] * fy[index]));

-			break;

-		case DATASET_VEL:

-			return_value = (float)sqrt((vx[index] * vx[index]) + (vy[index] * vy[index]));

-			break;

-		case DATASET_RHO:

-		default:

-			return_value = (float)rho[index];

-			break;

-		case DATASET_DIVV:

-      par_der_x = (float)vx[index + 1] - (vx[index] / cell_x);

-      par_der_y = (float)vy[index + 1] - (vy[index] / cell_y);

-      return_value = (float)(par_der_x + par_der_y);

-			break;

-		case DATASET_DIVF:

-      par_der_x = (float)fx[index + 1] - (fx[index] / cell_x);

-      par_der_y = (float)fy[index + 1] - (fy[index] / cell_y);

-      return_value = (float)(par_der_x + par_der_y);

-			break;

-	}

-	

-	return return_value;

-}

-

-void set_autoscaling(void)

-{

-  int k;

-  float value;

-

-  scale_min = scale_max = get_dataset(0);

-

-  for (k = 1; k < DIM * DIM; k++)

-  {

-    value = get_dataset(k);

-

-    if (scale_min > value) { scale_min = value; }

-    if (scale_max < value) { scale_max = value; }

-  }

-

-  //rho_threshold = (scale_min + scale_max) / 2;

-}

-

-void render_glyph(float x_value, float y_value, float i, float j)

-{

-  float x0, y0, z0, x1, y1, z1, x_dev, y_dev, z_dev, size, length;

-  float scale = (float)((float)DIM / (float)var_dims) / 6;

-  double theta, in_prod;

-  fftw_real wn, hn;

-

-  size = sqrt((x_value * x_value * 20) + (y_value * y_value * 20)) * 3 * scale;

-

-  wn = (fftw_real)winWidth / (fftw_real)(var_dims + 1);  // Grid cell width

-	hn = (fftw_real)winHeight / (fftw_real)(var_dims + 1); // Grid cell heigh

-

-  x0 = wn + (fftw_real)i * wn;

-  y0 = hn + (fftw_real)j * hn;

-  z0 = 0.0f;

-  

-  x1 = x0 + (vec_scale * x_value)/4;

-  y1 = y0 + (vec_scale * y_value)/4;

-  z1 = 0.0f;

-

-  // inner product

-  x_dev = x1 - x0;

-  y_dev = y1 - y0;

-  length = sqrt(x_dev * x_dev + y_dev * y_dev);

-  

-  x_dev = (length == 0) ? 0 : x_dev / length;

-  y_dev = (length == 0) ? 1 : y_dev / length;

-  

-  in_prod = x_dev * 0 + y_dev * 1;

-  theta = acos(in_prod) * (180/3.141592654);

-  if (x1 > x0) { theta *= -1; }

-

-  switch(glyph_sort)

-  {

-    default:

-    case GLYPH_LINES:

-      glBegin(GL_LINES);

-      glVertex3f(x0, y0, z0);

-      glVertex3f(x1, y1, z0);

-      glEnd();

-      break;

-

-    case GLYPH_ARROWS:

-      glPushMatrix();

-      glTranslatef(x0, y0, z0);

-      glRotatef(theta, 0.0, 0.0, 1.0);

-      glTranslatef(-x0, -y0, -z0);

-

-      glBegin(GL_TRIANGLES);

-        glVertex2d(-10 * size + x0, -25 * size + y0);

-        glVertex2d( 0  * size + x0,  25 * size + y0);

-		    glVertex2d( 10 * size + x0, -25 * size + y0);

-      glEnd();

-

-      glRotatef(-theta, 0.0, 0.0, 1.0);

-      glPopMatrix();

-      break;

-

-    case GLYPH_CYLINDERS:

-      x_dev = x1 - x0;

-      y_dev = y1 - y0;

-      length = sqrt(x_dev * x_dev + y_dev * y_dev) / 8;

-

-      glPushMatrix();

-      

-        glTranslatef(x0, y0, -1.0);

-        glRotatef(theta, 0.0, 0.0, 1.0);

-        glRotatef(-90, 1.0, 0.0, 0.0);

-        glScalef(length, length, length);

-        glCallList(startList);

-

-      glPopMatrix();

-      break;

-

-    case GLYPH_SPHERES:

-      glPushMatrix();

-      

-        glTranslatef(x0, y0, -1.0);

-        glRotatef(theta, 0.0, 0.0, 1.0);

-        glRotatef(-90, 1.0, 0.0, 0.0);

-        glScalef(length, length, length);

-        glCallList(startList);

-

-      glPopMatrix();

-      break;

-  }

-}

-

-void draw_glyphs(void)

-{

-  int i, j, idx;

-  float value, scale, idxcf, idxrf;

-

-  if (use_glyphs)

-  {

-    scale = (float)((float)DIM / (float)var_dims);

-    for (i = 0; i < var_dims; i++)

-    {

-	    for (j = 0; j < var_dims; j++)

-      {

-        idxcf = round(j * scale) * DIM;

-        idxrf = round(i * scale);

-        idx = idxcf + idxrf;

-

-        switch (glyph_scalar)

-        {

-          default:

-          case SCALAR_RHO:

-            value = rho[idx];

-            break;

-

-          case SCALAR_VEL:

-            value = vec_len(vx[idx], vy[idx]);

-            break;

-

-          case SCALAR_FORCE:

-            value = vec_len(fx[idx], fy[idx]);

-            break;

-        }

-

-        set_colormap(value, FALSE, 0.5f);

-        switch (glyph_vector)

-        {

-          default:

-          case VECTOR_VEL:

-            render_glyph(vx[idx], vy[idx], i, j);

-            break;

-

-          case VECTOR_FORCE:

-            render_glyph(fx[idx], fy[idx], i, j);

-            break;

-        }

-      }

-    }

-  }

-}

-

-#define percentage(a, b) (min(a, b) / max(a, b))

-

-void draw_isolines(float threshold)

-{

-  fftw_real wn = (fftw_real)winWidth  / (fftw_real)(DIM + 1); // Grid cell width

-	fftw_real hn = (fftw_real)winHeight / (fftw_real)(DIM + 1); // Grid cell height

-

-  int idx;

-  int i, j;

-  float v0, v1, v2, v3;

-  float x0, y0, x1, y1, x_offset, y_offset;

-  int state = 0;

-  static int prev_state = 0;

-

-  v0 = v1 = v2 = v3 = 0.0f;

-  x0 = y0 = x1 = y1 = 0.0f;

-

-  if (isolines_nr == 1) glLineWidth(2.0f);

-  else glLineWidth(2.0f);

-  glBegin(GL_LINES);

-

-  for (i = 0; i < DIM - 1; i++)

-  {

-    for (j = 0; j < DIM - 1; j++)

-    {

-      state = 0;

-      idx = (j * DIM) + i;

-      set_colormap(get_dataset(idx), 0, 1.0f);

-

-      v0 = get_dataset(idx + DIM);

-      v1 = get_dataset(idx + 1 + DIM);

-      v2 = get_dataset(idx + 1);

-      v3 = get_dataset(idx);

-

-      if (v0 >= threshold) { state += 1; }

-      if (v1 >= threshold) { state += 2; }

-      if (v2 >= threshold) { state += 4; }

-      if (v3 >= threshold) { state += 8; }

-

-      x_offset = wn + (fftw_real)i * wn;

-      y_offset = wn + (fftw_real)j * hn;

-

-      switch (state)

-      {

-        case 5:

-        case 10:

-          x0 = 0;

-          y0 = percentage(fabs(v3 - v0), threshold) * hn;

-          x1 = percentage(fabs(v1 - v0), threshold) * wn;

-          y1 = hn;

-          if (x0 != y0 != x1 != y1 != 0.0f) {

-            glVertex3i(x_offset + x0, y_offset + y0, 0.0f);

-            glVertex3i(x_offset + x1, y_offset + y1, 0.0f);

-          }

-          // no break !!

-        case 4:

-        case 11:

-          x0 = percentage(fabs(v3 - v2), threshold) * wn;

-          y0 = 0;

-          x1 = wn;

-          y1 = percentage(fabs(v2 - v1), threshold) * hn;

-          break;

-        case 6:

-        case 9:

-          x0 = percentage(fabs(v3 - v2), threshold) * wn;

-          y0 = 0;

-          x1 = percentage(fabs(v1 - v0), threshold) * wn;

-          y1 = hn;

-          break;

-        case 7:

-        case 8:

-          x0 = percentage(fabs(v3 - v2), threshold) * wn;

-          y0 = 0;

-          x1 = 0;

-          y1 = percentage(fabs(v3 - v0), threshold) * hn;

-          break;

-        case 3:

-        case 12:

-          x0 = 0;

-          y0 = percentage(fabs(v3 - v0), threshold) * hn;

-          x1 = wn;

-          y1 = percentage(fabs(v2 - v1), threshold) * hn;

-          break;

-        case 2:

-        case 13:

-          x0 = percentage(fabs(v1 - v0), threshold) * wn;

-          y0 = hn;

-          x1 = wn;

-          y1 = percentage(fabs(v2 - v1), threshold) * hn;

-          break;

-        case 1:

-        case 14:

-          x0 = 0;

-          y0 = percentage(fabs(v3 - v0), threshold) * hn;

-          x1 = percentage(fabs(v1 - v0), threshold) * wn;

-          y1 = hn;

-          break;

-        case 0: case 15: default:

-          x0 = y0 = x1 = y1 = 0.0f;

-          break;

-      }

-

-      // draw line

-      if (x0 != y0 != x1 != y1 != 0.0f)

-      {

-        glVertex3i(x_offset + x0, y_offset + y0, 0.0f);

-        glVertex3i(x_offset + x1, y_offset + y1, 0.0f);

-      }

-

-    }

-  }

-

-  glEnd();

-}

-

-//W

-

-//visualize: This is the main visualization function

-void visualize(void)

-{

-  float dataset;

-	static float fadein = 0;

-	int        i, j, idx; double px,py;

-	fftw_real  wn = (fftw_real)winWidth / (fftw_real)(DIM + 1);   // Grid cell width

-	fftw_real  hn = (fftw_real)winHeight / (fftw_real)(DIM + 1);  // Grid cell height

-  static int delay = 0;

-  float hight_value = 0.0f;

-

-  glTranslatef(-298, -295, -725);

-

-  if (delay == 10) {

-    set_autoscaling();

-    delay = 0;

-  } else delay++;

-

-  for (i = 0; i < winWidth; i++)

-  {

-		float value, clamp_scaled_min, clamp_scaled_max, scale_scaled_min, scale_scaled_max;

-

-		value = (float)((float)i/winWidth);

-

-		set_colormap(value, TRUE, 1.0f);

-

-		clamp_scaled_min = clamp_min *winWidth;

-		clamp_scaled_max = clamp_max *winWidth;

-

-		scale_scaled_min = scale_min *winWidth;

-		scale_scaled_max = scale_max *winWidth;

-

-		glBegin(GL_LINES);

-			if (!(i %20)) {

-				glColor3f(0.7, 0.7, 0.7);

-				glVertex2i(i -1, winHeight -6);

-				glVertex2i(i -1, winHeight -22);

-				glVertex2i(i, winHeight -6);

-				glVertex2i(i, winHeight -22);

-				glVertex2i(i +1, winHeight -6);

-				glVertex2i(i +1, winHeight -22);

-			} else {

-				glVertex2i(i, winHeight -6);

-				glVertex2i(i, winHeight -18);

-			}

-			glColor3f(0.2, 0.2, 0.2);

-			glVertex2i(clamp_scaled_min, winHeight -6);

-			glVertex2i(clamp_scaled_min, winHeight -18);

-			glVertex2i(clamp_scaled_max, winHeight -6);

-			glVertex2i(clamp_scaled_max, winHeight -18);

-			if (autoscale) {

-				glColor3f(1, 0, 0);

-				glVertex2i(scale_scaled_min, winHeight -4);

-				glVertex2i(scale_scaled_min, winHeight -18);

-				glVertex2i(scale_scaled_max, winHeight -4);

-				glVertex2i(scale_scaled_max, winHeight -18);

-			}

-		glEnd();

-

-		glBegin(GL_TRIANGLES);

-			glColor3f(0.9, 0.9, 0.9);

-			glVertex2i(clamp_scaled_min, winHeight -18);

-			glVertex2i(clamp_scaled_min -4, winHeight -25);

-			glVertex2i(clamp_scaled_min +4, winHeight -25);

-			glVertex2i(clamp_scaled_max, winHeight -18);

-			glVertex2i(clamp_scaled_max -4, winHeight -25);

-			glVertex2i(clamp_scaled_max +4, winHeight -25);

-			if (autoscale) {

-				glColor3f(1, 0, 0);

-				glVertex2i(scale_scaled_min, winHeight -6);

-				glVertex2i(scale_scaled_min -4, winHeight);

-				glVertex2i(scale_scaled_min +4, winHeight);

-				glVertex2i(scale_scaled_max, winHeight -6);

-				glVertex2i(scale_scaled_max -4, winHeight);

-				glVertex2i(scale_scaled_max +4, winHeight);

-			}

-		glEnd();

-	}

-

-  // Rotate field

-  glLoadIdentity();

-  glTranslatef(x_pos, y_pos, z_pos);

-  glRotatef(x_rot, 1.0f, 0.0f, 0.0f);

-  glRotatef(y_rot, 0.0f, 1.0f, 0.0f);

-  glRotatef(z_rot, 0.0f, 0.0f, 1.0f);

-  

-  glTranslatef(-(wn*DIM)/2, -(hn*DIM)/2, 0.0f);

-

-  // draw isolines

-  if (glyph_scalar == SCALAR_RHO)

-  {

-    int count;

-    float iso_scale;

-

-    if (isolines_nr) iso_scale = fabs(threshold1 - threshold2) / isolines_nr;

-    else iso_scale = 0.0f;

-

-    for (count = 0; count < isolines_nr; count++)

-    {

-      draw_isolines(min(threshold1, threshold2) + count * iso_scale);

-    }

-  }

-

-	// draw the smoke

-  if (draw_smoke)

-  {

-	  glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);

-

-	  for (j = 0; j < DIM - 1; j++)			//draw smoke

-	  {

-		  glBegin(GL_TRIANGLE_STRIP);

-

-		  i = 0;

-		  px = wn + (fftw_real)i * wn;

-		  py = hn + (fftw_real)j * hn;

-		  idx = (j * DIM) + i;

-  		

-		  glColor3f(rho[idx],rho[idx],rho[idx]);

-		  glVertex3f(px,py, hight_array[idx]);

-  						

-		  for (i = 0; i < DIM - 1; i++) 

-		  {

-			  px = wn + (fftw_real)i * wn;

-			  py = hn + (fftw_real)(j + 1) * hn;

-			  idx = ((j + 1) * DIM) + i;

-        set_colormap(get_dataset(idx), FALSE, 1.0f);

-        glVertex3f(px, py, hight_array[idx]);

-			  px = wn + (fftw_real)(i + 1) * wn;

-			  py = hn + (fftw_real)j * hn;

-			  idx = (j * DIM) + (i + 1);

-			  set_colormap(get_dataset(idx), FALSE, 1.0f);

-			  glVertex3f(px, py, hight_array[idx]);

-		  }

-

-		  px = wn + (fftw_real)(DIM - 1) * wn;

-		  py = hn + (fftw_real)(j + 1) * hn;

-		  idx = ((j + 1) * DIM) + (DIM - 1);

-		  set_colormap(get_dataset(idx), FALSE, 1.0f);

-		  glVertex3f(px, py, hight_array[idx]);

-		  glEnd();

-    }

-  }

-

-  draw_glyphs();

-

-	/* Draw the options screen if enabled, hide otherwise */

-	if (draw_options) {

-		glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);

-		glBegin(GL_QUADS);

-			glColor4f(0, 0, 0, 0.2);

-			glVertex2d(40, 40);

-			glVertex2d(40, winHeight -40);

-			glVertex2d(winWidth -40, winHeight -40);

-			glVertex2d(winWidth -40, 40);

-		glEnd();

-		glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);

-		glBegin(GL_QUADS);

-			glColor4f(0.9, 0.9, 0.9, 0.2);

-			glVertex2d(40, 40);

-			glVertex2d(40, winHeight -40);

-			glVertex2d(winWidth -40, winHeight -40);

-			glVertex2d(winWidth -40, 40);

-		glEnd();

-		glBegin(GL_QUADS);

-			glColor3f(1, 1, 1);

-			glVertex2d(50, 50);

-			glVertex2d(50, 70);

-			glVertex2d(70, 70);

-			glVertex2d(70, 50);

-		glEnd();

-		if (use_glyphs) {

-			glBegin(GL_LINES);

-				glVertex2d(50, 50);

-				glVertex2d(70, 70);

-				glVertex2d(50, 70);

-				glVertex2d(70, 50);

-			glEnd();

-		}

-	} else {

-		fadein = 0;

-	}

-}

-

-

-//------ INTERACTION CODE STARTS HERE -----------------------------------------------------------------

-

-

-//keyboard: Handle key presses

-void keyboard(unsigned char key, int x, int y) 

-{

-	switch (key) 

-	{

-	  case 't': dt -= 0.001; break;

-	  case 'T': dt += 0.001; break;

-	  case 'c': color_dir = (color_dir +1) %3; break;

-	  case 'S': vec_scale *= 1.2; break;

-	  case 's': vec_scale *= 0.8; break;

-	  case 'V': visc *= 5; break;

-	  case 'v': visc *= 0.2; break;

-    case 'y': toggle_glyph_usage(); break;

-    case 'x': toggle_smoke_usage(); break;

-	  case 'm': scalar_col++; if (scalar_col>COLOR_OLIVER) scalar_col=COLOR_BLACKWHITE; break;

-	  case 'a': frozen = 1-frozen; break;

-	  case 'q': exit(0); break;

-    case 'r': toggle_mouse_rotate(); break;

-    default: break;

-	}

-}

-

-void click(int button, int state, int mx, int my)

-{

-  /* First check wether we are on the top half, or bottom half of the bar, if neither then

-   * the mouseinput is for the smoke

-   */

-  if (!state)

-  {

-    prev_mx = mx;

-    prev_my = my;

-  }

-

-  if (my <25) {

-    /* Click received on button bar */

-    if (my <13) {

-      /* Upper half */

-      /* check wether mx is min or max slider */

-      if (mx <(rescale_to_winwidth(scale_min +((scale_max - scale_min) /2)))) {

-        active_slider = MOUSE_SCALE_MIN;

-      } else if (mx >(rescale_to_winwidth(scale_max -((scale_max - scale_min) /2)))) {

-        active_slider = MOUSE_SCALE_MAX;

-      }

-    } else {

-      /* Bottom half */

-      /* check wether mx is min or max slider */

-      if (mx <(rescale_to_winwidth(clamp_min +((clamp_max - clamp_min) /2)))) {

-        active_slider = MOUSE_CLAMP_MIN;

-      } else if (mx >(rescale_to_winwidth(clamp_max -((clamp_max - clamp_min) /2)))) {

-        active_slider = MOUSE_CLAMP_MAX;

-      }

-    }

-	}

-}

-

-void mouse_rotate(int mx, int my)

-{

-  x_rot -= ((float)(my - prev_my) / 5.0f);

-  y_rot -= ((float)(mx - prev_mx) / 5.0f);

-

-  prev_mx = mx;

-  prev_my = my;

-}

-

-// drag: When the user drags with the mouse, add a force that corresponds to the direction of the mouse

-//       cursor movement. Also inject some new matter into the field at the mouse location.

-void drag(int mx, int my) 

-{

-		int xi,yi,X,Y; double  dx, dy, len;

-		int clamp_scaled_min, clamp_scaled_max;

-		static int lmx=0,lmy=0;				//remembers last mouse location

-

-		/* Density calculations etc are only computed on area's below the slider */

-    if (rotate_plane)

-    {

-      mouse_rotate(mx, my);

-      return;

-    }

-

-		if (my > 25) {

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

-			xi = (int)clamp((double)(DIM + 1) * ((double)mx / (double)winWidth));

-			yi = (int)clamp((double)(DIM + 1) * ((double)(winHeight - my) / (double)winHeight));

-

-			X = xi; Y = yi;

-

-			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 *= 0.1 / len; dy *= 0.1 / len; }

-			fx[Y * DIM + X] += dx; 

-			fy[Y * DIM + X] += dy;

-			rho[Y * DIM + X] = 10.0f;

-			lmx = mx; lmy = my;

-		} else {

-			switch (active_slider) {

-				case MOUSE_SCALE_MIN:

-					scale_min = (float)mx /winWidth;

-					break;

-				case MOUSE_SCALE_MAX:

-					scale_max = (float)mx /winWidth;

-					break;

-				case MOUSE_CLAMP_MIN:

-					clamp_min = (float)mx /winWidth;

-					break;

-				case MOUSE_CLAMP_MAX:

-					clamp_max = (float)mx /winWidth;

-					break;

-				case MOUSE_SMOKE:

-				default:

-					break;

-			}

-	 }

-}

-

-void zoom_in(void) {

-  z_pos += 20;

-}

-

-void zoom_out(void) {

-  z_pos -= 20;

-}

-

-int get_glyph_usage(void) { return use_glyphs; }

-

-void selectColor(int arg) { scalar_col = arg; }

-void selectOliverscolor(int arg) { olivers_color = arg; }

-void select_dataset(int arg) { vis_dataset = arg; }

-void toggle_autoscale(void) { autoscale = (autoscale) ? FALSE : TRUE ; }

-void toggle_clamping(void) { }

-void show_options(void)     { draw_options = (draw_options) ? FALSE : TRUE; }

-void toggle_glyph_usage(void)  { use_glyphs = (use_glyphs) ? FALSE : TRUE; }

-void toggle_smoke_usage(void)  { draw_smoke = (draw_smoke) ? FALSE : TRUE; }

-void toggle_mouse_rotate(void)  { rotate_plane = (rotate_plane) ? FALSE : TRUE; }

-void set_glyph_scalar(int arg) { glyph_scalar = arg; }

-void set_glyph_vector(int arg) { glyph_vector = arg; }

+// Usage: Drag with the mouse to add smoke to the fluid. This will also move a "rotor" that disturbs 
+//        the velocity field at the mouse location. Press the indicated keys to change options
+//-------------------------------------------------------------------------------------------------- 
+
+#define max(a, b) (a > b ? a : b)
+#define min(a, b) (a < b ? a : b)
+#define round(x) (int)(x < 0 ? x - 0.5 : x + 0.5)
+
+#define MAX_INFINITE (1.#INF000)
+#define MIN_INFINITE (-1.#INF000)
+
+#define vec_len(x,y) (sqrt((x*x) + (y*y)))
+
+#define FALSE 0
+#define TRUE !FALSE
+
+#ifdef G_OS_WIN32
+#define WIN32_LEAN_AND_MEAN 1
+#include <windows.h>
+#endif
+
+#include <math.h>
+#include <rfftw.h>              //the numerical simulation FFTW library
+#include <stdio.h>              //for printing the help text
+
+#include <GL/gl.h>
+#include <GL/glu.h>
+
+#include "fluids.h"
+#include "seedpoint.h"
+
+//--- SIMULATION PARAMETERS ------------------------------------------------------------------------
+const int DIM = 50;             //size of simulation grid
+int var_dims = 50;
+double dt = 0.4;				        //simulation time step
+float visc = 0.001f;				    //fluid viscosity
+fftw_real *vx, *vy;             //(vx,vy)   = velocity field at the current moment
+fftw_real *vx0, *vy0;           //(vx0,vy0) = velocity field at the previous moment
+fftw_real *fx, *fy;	            //(fx,fy)   = user-controlled simulation forces, steered with the mouse 
+fftw_real *rho, *rho0;		    	//smoke density at the current (rho) and previous (rho0) moment 
+fftw_real *hight_array;         //used for hight plot
+int *frame_hist;
+int frame_index = 0;
+rfftwnd_plan plan_rc, plan_cr;  //simulation domain discretization
+
+int   winWidth, winHeight;      //size of the graphics window, in pixels
+int   color_dir = 0;            //use direction color-coding or not
+float vec_scale = 1000;			    //scaling of hedgehogs
+int   draw_smoke = 0;           //draw the smoke or not
+int   draw_vecs = 1;            //draw the vector field or not
+int   scalar_col = 0;           //method for scalar coloring
+int   frozen = 0;               //toggles on/off the animation
+int   olivers_color = 256;
+float clamp_min = 0;
+float clamp_max = 0.9999f;
+int autoscale = FALSE;
+float scale_min = 0;
+float scale_max = 1;
+int vis_dataset = DATASET_RHO;
+int mousebutton;
+int mousebuttonstate;
+int active_slider = 0;
+int use_glyphs = GLYPHS_OFF;
+int glyph_scalar = SCALAR_RHO;
+int glyph_vector = VECTOR_VEL;
+int glyph_sort = GLYPH_CYLINDERS;
+int draw_options = FALSE;
+GLuint startList;
+float threshold1 = 1.8f;
+float threshold2 = 2.0f;
+int isolines_nr = 1;
+int prev_mx = 0;
+int prev_my = 0;
+
+// rotation and translation vars
+int rotate_plane = 0;
+int seed_insert = 0;
+float x_rot = 0.0f;
+float y_rot = 0.0f;
+float z_rot = 0.0f;
+float x_pos = 0.0f;
+float y_pos = 0.0f;
+float z_pos = -1000.0f;
+
+// init pos of sim
+float xPos  = -298.0f;
+float yPos  = -295.0f;
+float zPos  = -725.0f;
+
+//------ SIMULATION CODE STARTS HERE -----------------------------------------------------------------
+
+//init_simulation: Initialize simulation data structures as a function of the grid size 'n'. 
+//                 Although the simulation takes place on a 2D grid, we allocate all data structures as 1D arrays,
+//                 for compatibility with the FFTW numerical library.
+void init_simulation(int n)				
+{
+  float LightAmbient[]  = { 0.25f, 0.25f, 0.25f, 1.0f }; // Ambient light values
+  float LightPosition[] = { 0.0f, 0.0f, -900.0f, 1.0f }; // Position of the light source
+	int i; size_t dim;
+  GLUquadricObj *qobj;
+	
+	dim     = n * 2*(n/2+1)*sizeof(fftw_real);        //Allocate data structures
+	vx       = (fftw_real*) malloc(dim); 
+	vy       = (fftw_real*) malloc(dim);
+	vx0      = (fftw_real*) malloc(dim); 
+	vy0      = (fftw_real*) malloc(dim);
+	dim     = n * n * sizeof(fftw_real);
+	fx      = (fftw_real*) malloc(dim);
+	fy      = (fftw_real*) malloc(dim);
+	rho     = (fftw_real*) malloc(dim); 
+	rho0    = (fftw_real*) malloc(dim);
+	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);
+  hight_array = (fftw_real*) malloc(dim);
+  frame_hist = (int*) malloc(n * sizeof(int *));
+
+  for (i = 0; i <= n; i++)
+  {
+      frame_hist[i] = (fftw_real*)malloc(n * 2*(n/2+1)*sizeof(fftw_real));
+  }
+	
+	for (i = 0; i < n * n; i++)                      //Initialize data structures to 0
+	{ vx[i] = vy[i] = vx0[i] = vy0[i] = fx[i] = fy[i] = rho[i] = rho0[i] = hight_array[i] = 0.0f; }
+
+  glShadeModel(GL_SMOOTH);                           // Enable smooth shading
+  glClearColor(0.0f, 0.0f, 0.0f, 0.0f);              // Black background
+  glClearDepth(1.0f);                                // Depth buffer setup
+  glDepthFunc(GL_LESS);                              // The type of depth testing to do
+  glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST); // Really nice perspective calculations
+  glHint(GL_POINT_SMOOTH_HINT, GL_NICEST);           // Really nice point smoothing
+  glEnable(GL_BLEND);
+  glBlendFunc(GL_SRC_ALPHA, GL_ONE);
+  glEnable(GL_DEPTH_TEST);
+  glEnable(GL_LIGHTING);
+  glEnable(GL_LIGHT0);
+  glEnable(GL_COLOR_MATERIAL);
+
+  glLightfv(GL_LIGHT0, GL_AMBIENT, LightAmbient);
+  glLightfv(GL_LIGHT0, GL_POSITION, LightPosition);
+
+  startList = glGenLists(2);
+  qobj = gluNewQuadric();
+
+  gluQuadricDrawStyle(qobj, GLU_FILL);
+  gluQuadricNormals(qobj, GLU_SMOOTH);
+  glNewList(startList, GL_COMPILE);
+    gluCylinder(qobj, 10, 0.2, 8, 8, 8);
+  glEndList();
+}
+
+int rescale_to_winwidth(float value)
+{
+	return round(value *winWidth);
+}
+
+//FFT: Execute the Fast Fourier Transform on the dataset 'vx'.
+//     'dirfection' indicates if we do the direct (1) or inverse (-1) Fourier Transform
+void FFT(int direction,void* vx)
+{
+	if(direction==1) rfftwnd_one_real_to_complex(plan_rc,(fftw_real*)vx,(fftw_complex*)vx);
+	else             rfftwnd_one_complex_to_real(plan_cr,(fftw_complex*)vx,(fftw_real*)vx);
+}
+
+int clamp(float x) 
+{ return ((x)>=0.0?((int)(x)):(-((int)(1-(x))))); }
+
+//solve: Solve (compute) one step of the fluid flow simulation
+void solve(int n, fftw_real* vx, fftw_real* vy, fftw_real* vx0, fftw_real* vy0, fftw_real visc, fftw_real dt) 
+{
+	fftw_real x, y, x0, y0, f, r, U[2], V[2], s, t;
+	int i, j, i0, j0, i1, j1;
+
+	for (i=0;i<n*n;i++) 
+	{ vx[i] += dt*vx0[i]; vx0[i] = vx[i]; vy[i] += dt*vy0[i]; vy0[i] = vy[i]; }    
+
+	for ( x=0.5f/n,i=0 ; i<n ; i++,x+=1.0f/n ) 
+	   for ( y=0.5f/n,j=0 ; j<n ; j++,y+=1.0f/n ) 
+	   {
+	      x0 = n*(x-dt*vx0[i+n*j])-0.5f; 
+	      y0 = n*(y-dt*vy0[i+n*j])-0.5f;
+	      i0 = clamp((float)x0); s = x0-i0;
+	      i0 = (n+(i0%n))%n;
+	      i1 = (i0+1)%n;
+	      j0 = clamp((float)y0); t = y0-j0;
+	      j0 = (n+(j0%n))%n;
+	      j1 = (j0+1)%n;
+	      vx[i+n*j] = (1-s)*((1-t)*vx0[i0+n*j0]+t*vx0[i0+n*j1])+s*((1-t)*vx0[i1+n*j0]+t*vx0[i1+n*j1]);
+	      vy[i+n*j] = (1-s)*((1-t)*vy0[i0+n*j0]+t*vy0[i0+n*j1])+s*((1-t)*vy0[i1+n*j0]+t*vy0[i1+n*j1]);
+	   }     
+	
+	for(i=0; i<n; i++)
+	  for(j=0; j<n; j++) 
+	  {  vx0[i+(n+2)*j] = vx[i+n*j]; vy0[i+(n+2)*j] = vy[i+n*j]; }
+
+	FFT(1,vx0);
+	FFT(1,vy0);
+
+	for (i=0;i<=n;i+=2) 
+	{
+	   x = 0.5f*i;
+	   for (j=0;j<n;j++) 
+	   {
+	      y = j<=n/2 ? (fftw_real)j : (fftw_real)j-n;
+	      r = x*x+y*y;
+	      if ( r==0.0f ) continue;
+	      f = (fftw_real)exp(-r*dt*visc);
+	      U[0] = vx0[i  +(n+2)*j]; V[0] = vy0[i  +(n+2)*j];
+	      U[1] = vx0[i+1+(n+2)*j]; V[1] = vy0[i+1+(n+2)*j];
+
+	      vx0[i  +(n+2)*j] = f*((1-x*x/r)*U[0]     -x*y/r *V[0]);
+	      vx0[i+1+(n+2)*j] = f*((1-x*x/r)*U[1]     -x*y/r *V[1]);
+	      vy0[i+  (n+2)*j] = f*(  -y*x/r *U[0] + (1-y*y/r)*V[0]);
+	      vy0[i+1+(n+2)*j] = f*(  -y*x/r *U[1] + (1-y*y/r)*V[1]);
+	   }
+	}
+
+	FFT(-1,vx0); 
+	FFT(-1,vy0);
+
+	f = 1.0/(n*n);
+ 	for (i=0;i<n;i++)
+	   for (j=0;j<n;j++) 
+	   { vx[i+n*j] = f*vx0[i+(n+2)*j]; vy[i+n*j] = f*vy0[i+(n+2)*j]; }
+} 
+
+
+// diffuse_matter: This function diffuses matter that has been placed in the velocity field. It's almost identical to the
+// velocity diffusion step in the function above. The input matter densities are in rho0 and the result is written into rho.
+void diffuse_matter(int n, fftw_real *vx, fftw_real *vy, fftw_real *rho, fftw_real *rho0, fftw_real dt) 
+{
+	fftw_real x, y, x0, y0, s, t;
+	int i, j, i0, j0, i1, j1;
+
+	for ( x=0.5f/n,i=0 ; i<n ; i++,x+=1.0f/n )
+		for ( y=0.5f/n,j=0 ; j<n ; j++,y+=1.0f/n ) 
+		{
+			x0 = n*(x-dt*vx[i+n*j])-0.5f; 
+			y0 = n*(y-dt*vy[i+n*j])-0.5f;
+			i0 = clamp((float)x0);
+			s = x0-i0;
+			i0 = (n+(i0%n))%n;
+			i1 = (i0+1)%n;
+			j0 = clamp((float)y0);
+			t = y0-j0;
+			j0 = (n+(j0%n))%n;
+			j1 = (j0+1)%n;
+			rho[i+n*j] = (1-s)*((1-t)*rho0[i0+n*j0]+t*rho0[i0+n*j1])+s*((1-t)*rho0[i1+n*j0]+t*rho0[i1+n*j1]);
+		}    
+}
+
+//set_forces: copy user-controlled forces to the force vectors that are sent to the solver. 
+//            Also dampen forces and matter density to get a stable simulation.
+void set_forces(void) 
+{
+	int i;
+	for (i = 0; i < DIM * DIM; i++) 
+	{
+        rho0[i]  = 0.995 * rho[i];
+        fx[i] *= 0.85; 
+        fy[i] *= 0.85;
+        vx0[i]    = fx[i]; 
+        vy0[i]    = fy[i];
+	}
+}
+
+
+void calculate_hight_plot(void)
+{
+  int i;
+	for (i = 0; i < DIM * DIM; i++) 
+	{
+    switch (vis_dataset)
+    {
+		  case DATASET_FORCE:
+			  hight_array[i] = sqrt(fx[i] * fx[i] + fy[i] * fy[i]) * 100;
+			  break;
+		  case DATASET_VEL:
+			  hight_array[i] = sqrt(vx[i] * vx[i] + vy[i] * vy[i]) * 2000;
+			  break;
+		  case DATASET_RHO:
+			  hight_array[i] = rho[i] * 10;
+			  break;
+      default:
+		  case DATASET_DIVV:
+		  case DATASET_DIVF:
+        hight_array[i] = 0.0f;
+			  break;
+	  }
+  }
+}
+
+void copy_frame(void)
+{
+  switch (vis_dataset)
+  {
+	  case DATASET_FORCE:
+		//  hight_array[i] = sqrt(fx[i] * fx[i] + fy[i] * fy[i]) * 100;
+		  break;
+	  case DATASET_VEL:
+		 // hight_array[i] = sqrt(vx[i] * vx[i] + vy[i] * vy[i]) * 2000;
+		  break;
+	  case DATASET_RHO:
+		  // hight_array[i] = rho[i] * 10;
+      memcpy(frame_hist[frame_index], rho, sizeof(DIM * 2 * (DIM / 2 + 1) * sizeof(fftw_real)));
+      break;
+    default:
+	  case DATASET_DIVV:
+	  case DATASET_DIVF:
+     // hight_array[i] = 0.0f;
+		  break;
+  }
+  frame_index = (frame_index + 1) % DIM;
+}
+
+//do_one_simulation_step: Do one complete cycle of the simulation:
+//      - set_forces:       
+//      - solve:            read forces from the user
+//      - diffuse_matter:   compute a new set of velocities
+//      - gluPostRedisplay: draw a new visualization frame
+
+
+void calculate_one_simulation_step(void)
+{
+    set_forces();
+	  solve(DIM, vx, vy, vx0, vy0, visc, dt);
+	  diffuse_matter(DIM, vx, vy, rho, rho0, dt);
+    calculate_hight_plot();
+    copy_frame();
+}
+
+
+//------ VISUALIZATION CODE STARTS HERE -----------------------------------------------------------------
+
+
+//rainbow: Implements a color palette, mapping the scalar 'value' to a rainbow color RGB
+void rainbow(float value,float* R,float* G,float* B)
+{                          
+   const float dx=0.8f; 
+   if (value<0) value=0; if (value>1) value=1;
+   value = (6-2*dx)*value+dx;
+   *R = (float)max(0.0f,(3-fabs(value-4.0f)-fabs(value-5.0f))/2.0f);
+   *G = (float)max(0.0f,(4-fabs(value-2.0f)-fabs(value-4.0f))/2.0f);
+   *B = (float)max(0.0f,(3-fabs(value-1.0f)-fabs(value-2.0f))/2.0f);
+}
+
+
+void colormap_fire(float value,float* R,float* G,float* B)
+{
+	/* Colormap Fire
+	 * A fire effect deals with two parts, first a drop from red to yellow (halfway)
+	 * during which time the Red component remains full e.g. 1. The Green component is
+	 * slowly added to turn the red into orange, and then yellow (R & G =Y). After this
+	 * point, the Red and Green component (e.g. yellow) have to drop simulataniously
+	 * to go from yellow down to black.
+	 */
+	*B = 0;
+	*G = 0;
+	*R = 0;
+
+	if (value <= (0.01)) {
+		/* whilst value is 0 - 0.5 both red and green equally change to create yellow*/
+		*R = *G = value;
+	} else {
+		/* whilst value is 0.5 - 1 Red is always fully on while the Green component is
+		 * added in steps to go from red to orange to yellow.
+		 */
+		*G = 0.9f - value;
+		*R = 0.8f; // not 1, makes red deeper, more intense
+	}
+}
+
+float remap(float value)
+{
+	value -= scale_min;
+	value /= (scale_max - scale_min);
+
+	//value = (value < 0) ? 0 : (value > 1) ? 1 : value;
+
+	return value;
+}
+
+//set_colormap: Sets three different types of colormaps
+struct color4f set_colormap(float vy, int draw_bar, float alpha)
+{
+   float R, G, B;
+   struct color4f return_value;
+
+   if (autoscale)
+   {
+	   vy = remap(vy);
+   }
+
+   if (!(draw_bar))
+   {
+     if (vy < clamp_min) vy = clamp_min;
+     if (vy > clamp_max) vy = clamp_max;
+   } else {
+     
+   }
+
+   vy *= olivers_color;
+   vy = (float)(int)(vy);
+   vy /= olivers_color;
+
+   if (scalar_col==COLOR_BLUE_GREEN_RED) {
+     if  (vy < -0.1) {
+       R = G = 0;
+       vy -= -0.1;
+       vy /= 0.9;
+       B = -vy;
+     } else if (vy < 0.1) {
+       R = B = 0;
+       vy += 0.1;
+       vy /= 0.2;
+       G = vy;
+     } else {
+       vy -= 0.1;
+       vy /= 0.9;
+       R = vy;
+       G = B = 0;
+     }
+   } else if (scalar_col==COLOR_BLACKWHITE) {
+       R = G = B = vy;
+   }
+   else if (scalar_col==COLOR_WILRIK) {
+	   colormap_fire(vy, &R, &G, &B);
+   }
+   else if (scalar_col==COLOR_OLIVER) {
+	  rainbow(vy,&R,&G,&B);   
+   }
+   else if (scalar_col==COLOR_RAINBOW)
+       rainbow(vy,&R,&G,&B); 
+   else if (scalar_col==COLOR_BANDS)
+       {  
+          const int NLEVELS = 7;
+          vy *= NLEVELS; vy = (float)(int)(vy); vy/= NLEVELS; 
+	      rainbow(vy,&R,&G,&B);   
+	   }
+
+   glColor4f(R,G,B,alpha);
+
+   return_value.r = R;
+   return_value.g = G;
+   return_value.b = B;
+   return_value.a = alpha;
+   return return_value;
+}
+
+
+//direction_to_color: Set the current color by mapping a direction vector (x,y), using
+//                    the color mapping method 'method'. If method==1, map the vector direction
+//                    using a rainbow colormap. If method==0, simply use the white color
+void direction_to_color(float x, float y, int method)
+{
+	float r,g,b,f;
+
+	switch (method) {
+	  case 3:
+		g = 1;
+		break;
+	  case 2:
+		r = 1;
+		g = b = 0;
+		break;
+	  case 1:
+		f = (float)(atan2((double)y, (double)x) / 3.1415927 + 1);
+		r = f;
+		if(r > 1) r = 2 - r;
+	  	g = f + 0.66667f;
+		if(g > 2) g -= 2;
+		if(g > 1) g = 2 - g;
+		b = f + 2.0f * 0.66667f;
+		if(b > 2) b -= 2;
+		if(b > 1) b = 2 - b;
+		break;
+	  case 0:
+	  default: r = g = b = 1;
+	  break;
+	}
+	glColor3f(r,g,b);
+}
+
+float get_dataset(int index)
+{
+  fftw_real cell_x = (fftw_real)winWidth  / (fftw_real)(DIM + 1); // Grid cell width
+	fftw_real cell_y = (fftw_real)winHeight / (fftw_real)(DIM + 1); // Grid cell heigh
+	float return_value, par_der_x, par_der_y;
+  static int klote = 0;
+
+	return_value = 0;
+
+	switch (vis_dataset)
+  {
+		case DATASET_FORCE:
+			return_value = (float)sqrt((fx[index] * fx[index]) + (fy[index] * fy[index]));
+			break;
+		case DATASET_VEL:
+			return_value = (float)sqrt((vx[index] * vx[index]) + (vy[index] * vy[index]));
+			break;
+		case DATASET_RHO:
+		default:
+			return_value = (float)rho[index];
+			break;
+		case DATASET_DIVV:
+      par_der_x = (float)vx[index + 1] - (vx[index] / cell_x);
+      par_der_y = (float)vy[index + 1] - (vy[index] / cell_y);
+      return_value = (float)(par_der_x + par_der_y);
+			break;
+		case DATASET_DIVF:
+      par_der_x = (float)fx[index + 1] - (fx[index] / cell_x);
+      par_der_y = (float)fy[index + 1] - (fy[index] / cell_y);
+      return_value = (float)(par_der_x + par_der_y);
+			break;
+	}
+	
+	return return_value;
+}
+
+void set_autoscaling(void)
+{
+  int k;
+  float value;
+
+  scale_min = scale_max = get_dataset(0);
+
+  for (k = 1; k < DIM * DIM; k++)
+  {
+    value = get_dataset(k);
+
+    if (scale_min > value) { scale_min = value; }
+    if (scale_max < value) { scale_max = value; }
+  }
+
+  //rho_threshold = (scale_min + scale_max) / 2;
+}
+
+void render_glyph(float x_value, float y_value, float i, float j)
+{
+  float x0, y0, z0, x1, y1, z1, x_dev, y_dev, z_dev, size, length;
+  float scale = (float)((float)DIM / (float)var_dims) / 6;
+  double theta, in_prod;
+  fftw_real wn, hn;
+
+  size = sqrt((x_value * x_value * 20) + (y_value * y_value * 20)) * 3 * scale;
+
+  wn = (fftw_real)winWidth / (fftw_real)(var_dims + 1);  // Grid cell width
+	hn = (fftw_real)winHeight / (fftw_real)(var_dims + 1); // Grid cell heigh
+
+  x0 = wn + (fftw_real)i * wn;
+  y0 = hn + (fftw_real)j * hn;
+  z0 = 0.0f;
+  
+  x1 = x0 + (vec_scale * x_value)/4;
+  y1 = y0 + (vec_scale * y_value)/4;
+  z1 = 0.0f;
+
+  // inner product
+  x_dev = x1 - x0;
+  y_dev = y1 - y0;
+  length = sqrt(x_dev * x_dev + y_dev * y_dev);
+  
+  x_dev = (length == 0) ? 0 : x_dev / length;
+  y_dev = (length == 0) ? 1 : y_dev / length;
+  
+  in_prod = x_dev * 0 + y_dev * 1;
+  theta = acos(in_prod) * (180/3.141592654);
+  if (x1 > x0) { theta *= -1; }
+
+  switch(glyph_sort)
+  {
+    default:
+    case GLYPH_LINES:
+      glBegin(GL_LINES);
+      glVertex3f(x0, y0, z0);
+      glVertex3f(x1, y1, z0);
+      glEnd();
+      break;
+
+    case GLYPH_ARROWS:
+      glPushMatrix();
+      glTranslatef(x0, y0, z0);
+      glRotatef(theta, 0.0, 0.0, 1.0);
+      glTranslatef(-x0, -y0, -z0);
+
+      glBegin(GL_TRIANGLES);
+        glVertex2d(-10 * size + x0, -25 * size + y0);
+        glVertex2d( 0  * size + x0,  25 * size + y0);
+		    glVertex2d( 10 * size + x0, -25 * size + y0);
+      glEnd();
+
+      glRotatef(-theta, 0.0, 0.0, 1.0);
+      glPopMatrix();
+      break;
+
+    case GLYPH_CYLINDERS:
+      x_dev = x1 - x0;
+      y_dev = y1 - y0;
+      length = sqrt(x_dev * x_dev + y_dev * y_dev) / 8;
+
+      glPushMatrix();
+      
+        glTranslatef(x0, y0, -1.0);
+        glRotatef(theta, 0.0, 0.0, 1.0);
+        glRotatef(-90, 1.0, 0.0, 0.0);
+        glScalef(length, length, length);
+        glCallList(startList);
+
+      glPopMatrix();
+      break;
+
+    case GLYPH_SPHERES:
+      glPushMatrix();
+      
+        glTranslatef(x0, y0, -1.0);
+        glRotatef(theta, 0.0, 0.0, 1.0);
+        glRotatef(-90, 1.0, 0.0, 0.0);
+        glScalef(length, length, length);
+        glCallList(startList);
+
+      glPopMatrix();
+      break;
+  }
+}
+
+void draw_glyphs(void)
+{
+  int i, j, idx;
+  float value, scale, idxcf, idxrf;
+
+  if (use_glyphs)
+  {
+    scale = (float)((float)DIM / (float)var_dims);
+    for (i = 0; i < var_dims; i++)
+    {
+	    for (j = 0; j < var_dims; j++)
+      {
+        idxcf = round(j * scale) * DIM;
+        idxrf = round(i * scale);
+        idx = idxcf + idxrf;
+
+        switch (glyph_scalar)
+        {
+          default:
+          case SCALAR_RHO:
+            value = rho[idx];
+            break;
+
+          case SCALAR_VEL:
+            value = vec_len(vx[idx], vy[idx]);
+            break;
+
+          case SCALAR_FORCE:
+            value = vec_len(fx[idx], fy[idx]);
+            break;
+        }
+
+        set_colormap(value, FALSE, 0.5f);
+        switch (glyph_vector)
+        {
+          default:
+          case VECTOR_VEL:
+            render_glyph(vx[idx], vy[idx], i, j);
+            break;
+
+          case VECTOR_FORCE:
+            render_glyph(fx[idx], fy[idx], i, j);
+            break;
+        }
+      }
+    }
+  }
+}
+
+#define percentage(a, b) (min(a, b) / max(a, b))
+
+void draw_isolines(float threshold)
+{
+  fftw_real wn = (fftw_real)winWidth  / (fftw_real)(DIM + 1); // Grid cell width
+	fftw_real hn = (fftw_real)winHeight / (fftw_real)(DIM + 1); // Grid cell height
+
+  int idx;
+  int i, j;
+  float v0, v1, v2, v3;
+  float x0, y0, x1, y1, x_offset, y_offset;
+  int state = 0;
+  static int prev_state = 0;
+
+  v0 = v1 = v2 = v3 = 0.0f;
+  x0 = y0 = x1 = y1 = 0.0f;
+
+  if (isolines_nr == 1) glLineWidth(2.0f);
+  else glLineWidth(2.0f);
+  glBegin(GL_LINES);
+
+  for (i = 0; i < DIM - 1; i++)
+  {
+    for (j = 0; j < DIM - 1; j++)
+    {
+      state = 0;
+      idx = (j * DIM) + i;
+      set_colormap(get_dataset(idx), 0, 1.0f);
+
+      v0 = get_dataset(idx + DIM);
+      v1 = get_dataset(idx + 1 + DIM);
+      v2 = get_dataset(idx + 1);
+      v3 = get_dataset(idx);
+
+      if (v0 >= threshold) { state += 1; }
+      if (v1 >= threshold) { state += 2; }
+      if (v2 >= threshold) { state += 4; }
+      if (v3 >= threshold) { state += 8; }
+
+      x_offset = wn + (fftw_real)i * wn;
+      y_offset = wn + (fftw_real)j * hn;
+
+      switch (state)
+      {
+        case 5:
+        case 10:
+          x0 = 0;
+          y0 = percentage(fabs(v3 - v0), threshold) * hn;
+          x1 = percentage(fabs(v1 - v0), threshold) * wn;
+          y1 = hn;
+          if (x0 != y0 != x1 != y1 != 0.0f) {
+            glVertex3i(x_offset + x0, y_offset + y0, 0.0f);
+            glVertex3i(x_offset + x1, y_offset + y1, 0.0f);
+          }
+          // no break !!
+        case 4:
+        case 11:
+          x0 = percentage(fabs(v3 - v2), threshold) * wn;
+          y0 = 0;
+          x1 = wn;
+          y1 = percentage(fabs(v2 - v1), threshold) * hn;
+          break;
+        case 6:
+        case 9:
+          x0 = percentage(fabs(v3 - v2), threshold) * wn;
+          y0 = 0;
+          x1 = percentage(fabs(v1 - v0), threshold) * wn;
+          y1 = hn;
+          break;
+        case 7:
+        case 8:
+          x0 = percentage(fabs(v3 - v2), threshold) * wn;
+          y0 = 0;
+          x1 = 0;
+          y1 = percentage(fabs(v3 - v0), threshold) * hn;
+          break;
+        case 3:
+        case 12:
+          x0 = 0;
+          y0 = percentage(fabs(v3 - v0), threshold) * hn;
+          x1 = wn;
+          y1 = percentage(fabs(v2 - v1), threshold) * hn;
+          break;
+        case 2:
+        case 13:
+          x0 = percentage(fabs(v1 - v0), threshold) * wn;
+          y0 = hn;
+          x1 = wn;
+          y1 = percentage(fabs(v2 - v1), threshold) * hn;
+          break;
+        case 1:
+        case 14:
+          x0 = 0;
+          y0 = percentage(fabs(v3 - v0), threshold) * hn;
+          x1 = percentage(fabs(v1 - v0), threshold) * wn;
+          y1 = hn;
+          break;
+        case 0: case 15: default:
+          x0 = y0 = x1 = y1 = 0.0f;
+          break;
+      }
+
+      // draw line
+      if (x0 != y0 != x1 != y1 != 0.0f)
+      {
+        glVertex3i(x_offset + x0, y_offset + y0, 0.0f);
+        glVertex3i(x_offset + x1, y_offset + y1, 0.0f);
+      }
+
+    }
+  }
+
+  glEnd();
+}
+
+//W
+
+//visualize: This is the main visualization function
+void visualize(void)
+{
+  float dataset;
+	static float fadein = 0;
+	int        i, j, idx; double px,py;
+	fftw_real  wn = (fftw_real)winWidth / (fftw_real)(DIM + 1);   // Grid cell width
+	fftw_real  hn = (fftw_real)winHeight / (fftw_real)(DIM + 1);  // Grid cell height
+  static int delay = 0;
+  float hight_value = 0.0f;
+
+  glTranslatef(xPos, yPos, zPos);
+
+  if (delay == 10) {
+    set_autoscaling();
+    delay = 0;
+  } else delay++;
+
+  for (i = 0; i < winWidth; i++)
+  {
+		float value, clamp_scaled_min, clamp_scaled_max, scale_scaled_min, scale_scaled_max;
+
+		value = (float)((float)i/winWidth);
+
+		set_colormap(value, TRUE, 1.0f);
+
+		clamp_scaled_min = clamp_min *winWidth;
+		clamp_scaled_max = clamp_max *winWidth;
+
+		scale_scaled_min = scale_min *winWidth;
+		scale_scaled_max = scale_max *winWidth;
+
+		glBegin(GL_LINES);
+			if (!(i %20)) {
+				glColor3f(0.7, 0.7, 0.7);
+				glVertex2i(i -1, winHeight -6);
+				glVertex2i(i -1, winHeight -22);
+				glVertex2i(i, winHeight -6);
+				glVertex2i(i, winHeight -22);
+				glVertex2i(i +1, winHeight -6);
+				glVertex2i(i +1, winHeight -22);
+			} else {
+				glVertex2i(i, winHeight -6);
+				glVertex2i(i, winHeight -18);
+			}
+			glColor3f(0.2, 0.2, 0.2);
+			glVertex2i(clamp_scaled_min, winHeight -6);
+			glVertex2i(clamp_scaled_min, winHeight -18);
+			glVertex2i(clamp_scaled_max, winHeight -6);
+			glVertex2i(clamp_scaled_max, winHeight -18);
+			if (autoscale) {
+				glColor3f(1, 0, 0);
+				glVertex2i(scale_scaled_min, winHeight -4);
+				glVertex2i(scale_scaled_min, winHeight -18);
+				glVertex2i(scale_scaled_max, winHeight -4);
+				glVertex2i(scale_scaled_max, winHeight -18);
+			}
+		glEnd();
+
+		glBegin(GL_TRIANGLES);
+			glColor3f(0.9, 0.9, 0.9);
+			glVertex2i(clamp_scaled_min, winHeight -18);
+			glVertex2i(clamp_scaled_min -4, winHeight -25);
+			glVertex2i(clamp_scaled_min +4, winHeight -25);
+			glVertex2i(clamp_scaled_max, winHeight -18);
+			glVertex2i(clamp_scaled_max -4, winHeight -25);
+			glVertex2i(clamp_scaled_max +4, winHeight -25);
+			if (autoscale) {
+				glColor3f(1, 0, 0);
+				glVertex2i(scale_scaled_min, winHeight -6);
+				glVertex2i(scale_scaled_min -4, winHeight);
+				glVertex2i(scale_scaled_min +4, winHeight);
+				glVertex2i(scale_scaled_max, winHeight -6);
+				glVertex2i(scale_scaled_max -4, winHeight);
+				glVertex2i(scale_scaled_max +4, winHeight);
+			}
+		glEnd();
+	}
+
+  // Rotate field
+  glLoadIdentity();
+  glTranslatef(x_pos, y_pos, z_pos);
+  glRotatef(x_rot, 1.0f, 0.0f, 0.0f);
+  glRotatef(y_rot, 0.0f, 1.0f, 0.0f);
+  glRotatef(z_rot, 0.0f, 0.0f, 1.0f);
+  
+  glTranslatef(-(wn*DIM)/2, -(hn*DIM)/2, 0.0f);
+
+  // draw isolines
+  if (glyph_scalar == SCALAR_RHO)
+  {
+    int count;
+    float iso_scale;
+
+    if (isolines_nr) iso_scale = fabs(threshold1 - threshold2) / isolines_nr;
+    else iso_scale = 0.0f;
+
+    for (count = 0; count < isolines_nr; count++)
+    {
+      draw_isolines(min(threshold1, threshold2) + count * iso_scale);
+    }
+  }
+
+	// draw the smoke
+  if (draw_smoke)
+  {
+	  glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
+
+	  for (j = 0; j < DIM - 1; j++)			//draw smoke
+	  {
+		  glBegin(GL_TRIANGLE_STRIP);
+
+		  i = 0;
+		  px = wn + (fftw_real)i * wn;
+		  py = hn + (fftw_real)j * hn;
+		  idx = (j * DIM) + i;
+  		
+		  glColor3f(rho[idx],rho[idx],rho[idx]);
+		  glVertex3f(px,py, hight_array[idx]);
+  						
+		  for (i = 0; i < DIM - 1; i++) 
+		  {
+			  px = wn + (fftw_real)i * wn;
+			  py = hn + (fftw_real)(j + 1) * hn;
+			  idx = ((j + 1) * DIM) + i;
+        set_colormap(get_dataset(idx), FALSE, 1.0f);
+        glVertex3f(px, py, hight_array[idx]);
+			  px = wn + (fftw_real)(i + 1) * wn;
+			  py = hn + (fftw_real)j * hn;
+			  idx = (j * DIM) + (i + 1);
+			  set_colormap(get_dataset(idx), FALSE, 1.0f);
+			  glVertex3f(px, py, hight_array[idx]);
+		  }
+
+		  px = wn + (fftw_real)(DIM - 1) * wn;
+		  py = hn + (fftw_real)(j + 1) * hn;
+		  idx = ((j + 1) * DIM) + (DIM - 1);
+		  set_colormap(get_dataset(idx), FALSE, 1.0f);
+		  glVertex3f(px, py, hight_array[idx]);
+		  glEnd();
+    }
+  }
+
+  draw_glyphs();
+  render_seedpoints();
+
+	/* Draw the options screen if enabled, hide otherwise */
+	if (draw_options) {
+		glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
+		glBegin(GL_QUADS);
+			glColor4f(0, 0, 0, 0.2);
+			glVertex2d(40, 40);
+			glVertex2d(40, winHeight -40);
+			glVertex2d(winWidth -40, winHeight -40);
+			glVertex2d(winWidth -40, 40);
+		glEnd();
+		glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
+		glBegin(GL_QUADS);
+			glColor4f(0.9, 0.9, 0.9, 0.2);
+			glVertex2d(40, 40);
+			glVertex2d(40, winHeight -40);
+			glVertex2d(winWidth -40, winHeight -40);
+			glVertex2d(winWidth -40, 40);
+		glEnd();
+		glBegin(GL_QUADS);
+			glColor3f(1, 1, 1);
+			glVertex2d(50, 50);
+			glVertex2d(50, 70);
+			glVertex2d(70, 70);
+			glVertex2d(70, 50);
+		glEnd();
+		if (use_glyphs) {
+			glBegin(GL_LINES);
+				glVertex2d(50, 50);
+				glVertex2d(70, 70);
+				glVertex2d(50, 70);
+				glVertex2d(70, 50);
+			glEnd();
+		}
+	} else {
+		fadein = 0;
+	}
+}
+
+
+//------ INTERACTION CODE STARTS HERE -----------------------------------------------------------------
+
+
+//keyboard: Handle key presses
+void keyboard(unsigned char key, int x, int y) 
+{
+	switch (key) 
+	{
+	  case 't': dt -= 0.001; break;
+	  case 'T': dt += 0.001; break;
+	  case 'c': color_dir = (color_dir +1) %3; break;
+	  case 'S': vec_scale *= 1.2; break;
+	  case 's': vec_scale *= 0.8; break;
+	  case 'V': visc *= 5; break;
+	  case 'v': visc *= 0.2; break;
+    case 'y': toggle_glyph_usage(); break;
+    case 'x': toggle_smoke_usage(); break;
+	  case 'm': scalar_col++; if (scalar_col>COLOR_OLIVER) scalar_col=COLOR_BLACKWHITE; break;
+	  case 'a': frozen = 1-frozen; break;
+	  case 'q': exit(0); break;
+    case 'r': toggle_mouse_rotate(); break;
+    case 'e': toggle_seed_insert(); break;
+    default: break;
+	}
+}
+
+void click(int button, int state, int mx, int my)
+{
+  /* First check wether we are on the top half, or bottom half of the bar, if neither then
+   * the mouseinput is for the smoke
+   */
+  if (!state)
+  {
+    prev_mx = mx;
+    prev_my = my;
+  }
+
+  if (seed_insert)
+  {
+    create_seedpoint(mx, my);
+  }
+
+  if (my <25) {
+    /* Click received on button bar */
+    if (my <13) {
+      /* Upper half */
+      /* check wether mx is min or max slider */
+      if (mx <(rescale_to_winwidth(scale_min +((scale_max - scale_min) /2)))) {
+        active_slider = MOUSE_SCALE_MIN;
+      } else if (mx >(rescale_to_winwidth(scale_max -((scale_max - scale_min) /2)))) {
+        active_slider = MOUSE_SCALE_MAX;
+      }
+    } else {
+      /* Bottom half */
+      /* check wether mx is min or max slider */
+      if (mx <(rescale_to_winwidth(clamp_min +((clamp_max - clamp_min) /2)))) {
+        active_slider = MOUSE_CLAMP_MIN;
+      } else if (mx >(rescale_to_winwidth(clamp_max -((clamp_max - clamp_min) /2)))) {
+        active_slider = MOUSE_CLAMP_MAX;
+      }
+    }
+	}
+}
+
+void mouse_rotate(int mx, int my)
+{
+  x_rot -= ((float)(my - prev_my) / 5.0f);
+  y_rot -= ((float)(mx - prev_mx) / 5.0f);
+
+  prev_mx = mx;
+  prev_my = my;
+}
+
+// drag: When the user drags with the mouse, add a force that corresponds to the direction of the mouse
+//       cursor movement. Also inject some new matter into the field at the mouse location.
+void drag(int mx, int my) 
+{
+		int xi,yi,X,Y; double  dx, dy, len;
+		int clamp_scaled_min, clamp_scaled_max;
+		static int lmx=0,lmy=0;				//remembers last mouse location
+
+		/* Density calculations etc are only computed on area's below the slider */
+    if (rotate_plane)
+    {
+      mouse_rotate(mx, my);
+      return;
+    }
+
+		if (my > 25) {
+			// Compute the array index that corresponds to the cursor location 
+			xi = (int)clamp((double)(DIM + 1) * ((double)mx / (double)winWidth));
+			yi = (int)clamp((double)(DIM + 1) * ((double)(winHeight - my) / (double)winHeight));
+
+			X = xi; Y = yi;
+
+			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 *= 0.1 / len; dy *= 0.1 / len; }
+			fx[Y * DIM + X] += dx; 
+			fy[Y * DIM + X] += dy;
+			rho[Y * DIM + X] = 10.0f;
+			lmx = mx; lmy = my;
+		} else {
+			switch (active_slider) {
+				case MOUSE_SCALE_MIN:
+					scale_min = (float)mx /winWidth;
+					break;
+				case MOUSE_SCALE_MAX:
+					scale_max = (float)mx /winWidth;
+					break;
+				case MOUSE_CLAMP_MIN:
+					clamp_min = (float)mx /winWidth;
+					break;
+				case MOUSE_CLAMP_MAX:
+					clamp_max = (float)mx /winWidth;
+					break;
+				case MOUSE_SMOKE:
+				default:
+					break;
+			}
+	 }
+}
+
+void zoom_in(void) {
+  z_pos += 20;
+}
+
+void zoom_out(void) {
+  z_pos -= 20;
+}
+
+int get_glyph_usage(void) { return use_glyphs; }
+
+void selectColor(int arg) { scalar_col = arg; }
+void selectOliverscolor(int arg) { olivers_color = arg; }
+void select_dataset(int arg) { vis_dataset = arg; }
+void toggle_autoscale(void) { autoscale = (autoscale) ? FALSE : TRUE ; }
+void toggle_clamping(void) { }
+void show_options(void)     { draw_options = (draw_options) ? FALSE : TRUE; }
+void toggle_glyph_usage(void)  { use_glyphs = (use_glyphs) ? FALSE : TRUE; }
+void toggle_smoke_usage(void)  { draw_smoke = (draw_smoke) ? FALSE : TRUE; }
+void toggle_mouse_rotate(void)  { rotate_plane = (rotate_plane) ? FALSE : TRUE; }
+void toggle_seed_insert(void)  { seed_insert = (seed_insert) ? FALSE : TRUE; }
+void set_glyph_scalar(int arg) { glyph_scalar = arg; }
+void set_glyph_vector(int arg) { glyph_vector = arg; }
diff --git a/Smoke/fluids.h b/Smoke/fluids.h
index 7464e73..7e306f7 100644
--- a/Smoke/fluids.h
+++ b/Smoke/fluids.h
@@ -66,6 +66,10 @@ extern int glyph_scalar;
 extern int glyph_vector;
 extern int draw_options;
 
+extern float xPos;
+extern float yPos;
+extern float zPos;
+
 void init_simulation(int n);
 
 void visualize(void);
@@ -93,5 +97,19 @@ void show_options(void);
 void toggle_glyph_usage(void);
 void toggle_smoke_usage(void);
 void toggle_mouse_rotate(void);
+void toggle_seed_insert(void);
 void set_glyph_scalar(int arg);
 void set_glyph_vector(int arg);
+
+struct point {
+  float x;
+  float y;
+  float z;
+};
+
+struct color4f {
+  float r;
+  float g;
+  float b;
+  float a;
+};
diff --git a/Smoke/gtk.c b/Smoke/gtk.c
index 681b6c8..e9a2462 100644
--- a/Smoke/gtk.c
+++ b/Smoke/gtk.c
@@ -288,7 +288,7 @@ scroll_event (GtkWidget      *widget,
 		    gpointer        data)
 {
   //g_print ("%s: \"scroll_event\": ", gtk_widget_get_name (widget));
-
+  //if (event->state == GDK_ALT) {
   if (event->direction == GDK_SCROLL_UP) {
     zoom_in();
   } else if (event->direction == GDK_SCROLL_DOWN) {
@@ -426,6 +426,10 @@ key_press_event (GtkWidget   *widget,
     case GDK_r:
         toggle_mouse_rotate();
         break;
+    
+    case GDK_e:
+        toggle_seed_insert();
+        break;
 
     default:
 //      g_print("\n");