1 files changed, 28 insertions, 173 deletions

diff --git a/Smoke/fluids.c b/Smoke/fluids.c
index e9756df..8832a3a 100644
--- a/Smoke/fluids.c
+++ b/Smoke/fluids.c
@@ -11,8 +11,10 @@
 #include <stdio.h>
 
 #include "fluids.h"
+#include "colormap.h"
 #include "glyphs.h"
 #include "seedpoint.h"
+#include "streamlines.h"
 #include "funcs.h"
 
 //--- SIMULATION PARAMETERS ------------------------------------------------------------------------
@@ -36,7 +38,6 @@ 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;
@@ -52,7 +53,7 @@ int glyph_vector = VECTOR_VEL;
 int glyph_sort = GLYPH_CYLINDERS;
 int draw_options = FALSE;
 GLuint startList;
-float threshold1 = 1.8f;
+float threshold1 = 0.2f;
 float threshold2 = 2.0f;
 int isolines_nr = 1;
 int prev_mx = 0;
@@ -75,20 +76,21 @@ 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
+  float LightPosition[] = { 0.0f, 0.0f, -700.0f, 1.0f }; // Position of the light source
 	int i; size_t dim;
 	
 	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);
+	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);
@@ -96,8 +98,9 @@ void init_simulation(int n)
 	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 *));
+  frame_hist  = (int*) malloc(n * sizeof(int *));
 
   for (i = 0; i <= n; i++)
   {
@@ -116,7 +119,7 @@ void init_simulation(int n)
   glEnable(GL_BLEND);
   glBlendFunc(GL_SRC_ALPHA, GL_ONE);
   glEnable(GL_DEPTH_TEST);
-  glEnable(GL_LIGHTING);
+  //glEnable(GL_LIGHTING);
   glEnable(GL_LIGHT0);
   glEnable(GL_COLOR_MATERIAL);
 
@@ -250,10 +253,10 @@ void calculate_hight_plot(void)
     switch (vis_dataset)
     {
 		  case DATASET_FORCE:
-			  hight_array[i] = sqrt(fx[i] * fx[i] + fy[i] * fy[i]) * 100;
+			  hight_array[i] = quake_root(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;
+			  hight_array[i] = quake_root(vx[i] * vx[i] + vy[i] * vy[i]) * 2000;
 			  break;
 		  case DATASET_RHO:
 			  hight_array[i] = rho[i] * 10;
@@ -305,44 +308,6 @@ void calculate_one_simulation_step(void)
 
 //------ 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;
@@ -351,73 +316,6 @@ float remap(float 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
@@ -463,10 +361,10 @@ float get_dataset(int index)
 	switch (vis_dataset)
   {
 		case DATASET_FORCE:
-			return_value = (float)sqrt((fx[index] * fx[index]) + (fy[index] * fy[index]));
+			return_value = (float)quake_root((fx[index] * fx[index]) + (fy[index] * fy[index]));
 			break;
 		case DATASET_VEL:
-			return_value = (float)sqrt((vx[index] * vx[index]) + (vy[index] * vy[index]));
+			return_value = (float)quake_root((vx[index] * vx[index]) + (vy[index] * vy[index]));
 			break;
 		case DATASET_RHO:
 		default:
@@ -502,7 +400,7 @@ void set_autoscaling(void)
     if (scale_max < value) { scale_max = value; }
   }
 
-  //rho_threshold = (scale_min + scale_max) / 2;
+  //threshold = (scale_min + scale_max) / 2;
 }
 
 #define percentage(a, b) (min(a, b) / max(a, b))
@@ -532,7 +430,7 @@ void draw_isolines(float threshold)
     {
       state = 0;
       idx = (j * DIM) + i;
-      set_colormap(get_dataset(idx), 0, 1.0f);
+      set_colormap(scalar_col, get_dataset(idx), 0, 1.0f);
 
       v0 = get_dataset(idx + DIM);
       v1 = get_dataset(idx + 1 + DIM);
@@ -623,7 +521,6 @@ void draw_isolines(float threshold)
 //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
@@ -633,18 +530,13 @@ void visualize(void)
 
   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);
+		set_colormap(scalar_col, value, TRUE, 1.0f);
 
 		clamp_scaled_min = clamp_min *winWidth;
 		clamp_scaled_max = clamp_max *winWidth;
@@ -745,19 +637,19 @@ void visualize(void)
 			  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);
+        set_colormap(scalar_col, 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);
+			  set_colormap(scalar_col, 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);
+		  set_colormap(scalar_col, get_dataset(idx), FALSE, 1.0f);
 		  glVertex3f(px, py, hight_array[idx]);
 		  glEnd();
     }
@@ -765,43 +657,7 @@ void visualize(void)
 
   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;
-	}
+  render_streamlines();
 }
 
 
@@ -816,10 +672,10 @@ void keyboard(unsigned char key, int x, int y)
 	  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 'S': vec_scale *= 1.2f; break;
+	  case 's': vec_scale *= 0.8f; break;
+	  case 'V': visc *= 5.0f; break;
+	  case 'v': visc *= 0.2f; 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;
@@ -906,7 +762,7 @@ void drag(int mx, int my)
 			// Add force at the cursor location 
 			my = winHeight - my;
 			dx = mx - lmx; dy = my - lmy;
-			len = sqrt(dx * dx + dy * dy);
+			len = vec_len(dx, dy);
 			if (len != 0.0) {  dx *= 0.1 / len; dy *= 0.1 / len; }
 			fx[Y * DIM + X] += dx; 
 			fy[Y * DIM + X] += dy;
@@ -946,7 +802,6 @@ void zoom_out(void)
 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) { }