/* 
Water plugin for Blender 
All the complicated maths for the effect.
*/
	
/* Scot Scriven: Numbers for trig tables*/
#define FSINMAX 2047
#define SINFIX 16
#define FSINBITS 16
#define PI 3.14159265358979323846
static int FSinTab[2048];
static int FCosTab[2048];
static int FSin(int angle);
static int FCos(int angle);

/*How wide and large my array is */
static int gBufferWidth;
static int gBufferLen;

/*The arrays and pointers to them */
static int *gArray1=NULL, *gArray2=NULL;
static int *gSwapper, *gDest, *gSource;

/* vars for trails and swirl*/
static int gXang, gYang;  
static int gSwirlAng;

/*************************/
/* Simple Noise function */
/*************************/
#define MaxNoise 33      /* must be 2^x+1 for and-mask */
#define MNMask 31 /* MaxNoise-2 */
#define NK1 33 /* MaxNoise */
#define NK2 1089 /* NK1*NK1 */
#define NCUBE 35937 /* NK1*NK2 */
#define xyz1 1
#define xy1z 33 /* NK1 */
#define xy1z1 34 /* xy1z+1 */
#define x1yz 1089 /* NK2 */
#define x1yz1 1090 /* x1yz+1 */
#define x1y1z 1122 /* NK2+NK1 */
#define x1y1z1 1123 /* x1y1z+1 */

/* Fast Noise Cube Array */
static float NoiseMTX[NCUBE];

void InitNoise()
{
int x, y, z, i, j, k, idx;
	for (x=0;x<MaxNoise;x++)
		for (y=0;y<MaxNoise;y++)
			for (z=0;z<MaxNoise;z++) {
				idx = x*NK2 + y*NK1 + z;
				NoiseMTX[idx] = rand()/(float)RAND_MAX;
				if ((i=x)==(MaxNoise-1)) i=0;
				if ((j=y)==(MaxNoise-1)) j=0;
				if ((k=z)==(MaxNoise-1)) k=0;
				NoiseMTX[idx] = NoiseMTX[i*NK2 + j*NK1 + k];
			}
}

float FNoise(float x, float y, float z)
{
int ix, iy, iz, AIDX;
float ox, oy, oz, p00, p01, p10, p11, p0, p1;
	x = fabs(x);
	y = fabs(y);
	z = fabs(z);
	ix = ((int)x) & MNMask;
	iy = ((int)y) & MNMask;
	iz = ((int)z) & MNMask;
	ox = x - floor(x);
	oy = y - floor(y);
	oz = z - floor(z);
	AIDX = ix*NK2 + iy*NK1 +iz;
	ox = ox * ox * (3.f - 2.f * ox);
	p00 = (NoiseMTX[AIDX + x1yz] - NoiseMTX[AIDX]) * ox + NoiseMTX[AIDX];
	p01 = (NoiseMTX[AIDX + x1yz1] - NoiseMTX[AIDX + xyz1]) * ox + NoiseMTX[AIDX + xyz1];
	p10 = (NoiseMTX[AIDX + x1y1z] - NoiseMTX[AIDX + xy1z]) * ox + NoiseMTX[AIDX + xy1z];
	p11 = (NoiseMTX[AIDX + x1y1z1] - NoiseMTX[AIDX + xy1z1]) * ox + NoiseMTX[AIDX + xy1z1];
	oy = oy * oy * (3.f - 2.f * oy);
	p0 = (p10 - p00) * oy + p00;
	p1 = (p11 - p01) * oy + p01;
	return ((p1-p0) * (oz * oz * (3.f - 2.f * oz)) + p0);
}


/*Scott Scriven: Setup the trig tables*/
int FSin(int angle)
{
	return FSinTab[angle&FSINMAX];
}

int FCos(int angle)
{
	return FCosTab[angle&FSINMAX];
}

/* Pre-calculates sin and cos values to speed later use */
void FCreateSines()
{
int i;
double angle;
	for(i=0; i<2048; i++)
	{
		angle = (float)i * (PI/1024.0);
		FSinTab[i] = (int)(sin(angle) * (float)0x10000);
		FCosTab[i] = (int)(cos(angle) * (float)0x10000);
	}
}

/* +:+:+:+:+:+:+:+:+:+:+: BLURRING ROUTINE +:+:+:+:+:+:+:+:+:+:+: */
/* Code to perform the blurring that fakes the water */
void doWave(Cast *cast)
{
	int newh;
	int x, y;   
	int xm, ym, xp, yp, yc;
	int s;
	
	/* E: calculate filter so it wraps around making the texture butSeamless */
	/* D: The trick here is to simply wrap the variables around
		the edges if they go out of bounds - so the effect is continued from the other
		side - very cunning!
	*/
	for (s = 0; s < cast->speed ; s++) {/* "Speed" is really just doing the sim multiple times */
		if (cast->butSeamless)
		{
			for (y = 0; y < gBufferWidth; y++) {
				if ((ym = y-1)==-1) ym = gBufferWidth-1;
				if ((yp = y+1)==gBufferWidth) yp = 0;
				ym *= gBufferWidth;
				yp *= gBufferWidth;
				yc = y*gBufferWidth;
				for (x=0;x<gBufferWidth;x++) {
					if ((xm = x-1)==-1) xm = gBufferWidth-1;
					if ((xp = x+1)==gBufferWidth) xp = 0;
					newh = ( (gSource[ym+xm] + gSource[yc+xm] + gSource[yp+xm]
							+ gSource[ym+x]  + gSource[yp+x]
							+ gSource[ym+xp] + gSource[yc+xp] + gSource[yp+xp])
							>> 2) - gDest[yc+x];
					gDest[yc+x] = newh - (newh >> cast->density);
				}
			}
		} else
		{/* not seamless */
			for (y=1; y < gBufferWidth-2; y++) {
				ym = y-1;
				yp = y+1;
				ym *= gBufferWidth;
				yp *= gBufferWidth;
				yc = y*gBufferWidth;
				for (x=1; x<gBufferWidth-2; x++) {
					xm = x-1;
					xp = x+1;
					newh = ( (gSource[ym+xm] + gSource[yc+xm] + gSource[yp+xm]
							+ gSource[ym+x]  + gSource[yp+x]
							+ gSource[ym+xp] + gSource[yc+xp] + gSource[yp+xp])
							>> 2) - gDest[yc+x];
					gDest[yc+x] = newh - (newh >> cast->density);
				}
			}  
		}
		/*Swap pointers to buffers! Thanks Mike! */
		gSwapper = gDest;
		gDest = gSource;
		gSource = gSwapper;
	} /* loop */
}

/* Draws a hard-edged, heavy drop into the buffer */
void heavyDrop(int x, int y, int radius, int height)
{
/* Scott Scriven's code: */
int rquad;
int cx, cy, cyq;
int left, top, right, bottom;
	
	rquad = radius * radius;
	
	/* Make a randomly-placed blob... */
	if(x<0) x = 1 + radius + (rand() % (gBufferWidth-2*radius-1));
	if(y<0) y = 1 + radius + (rand() % (gBufferWidth-2*radius-1));
	
	left=-radius;  right = radius;
	top=-radius;  bottom = radius;
	
	/* Perform edge clipping... */
	if(x - radius < 1) left -= (x-radius-1);
	if(y - radius < 1) top  -= (y-radius-1);
	if(x + radius > gBufferWidth-1) right -= (x+radius-gBufferWidth+1);
	if(y + radius > gBufferWidth-1) bottom-= (y+radius-gBufferWidth+1);
	
	for(cy = top; cy < bottom; cy++)
	{
		cyq = cy*cy;
		for(cx = left; cx < right; cx++)
		{
		if(cx*cx + cyq < rquad)
			gSource[gBufferWidth*(cy+y) + (cx+x)] += height;
		}
	}

}

/*Draws a fuzzy drop into the buffer */
void fuzzyDrop(int x, int y, int radius, int height)
{
/* Scott Scriven's code: */
int cx, cy;
int left,top,right,bottom;
int square, dist;
int radsquare = radius * radius;
float length = (1024.0/(float)radius)*(1024.0/(float)radius);
	if(x<0) x = 1 + radius + (rand() % (gBufferWidth-2*radius-1));
	if(y<0) y = 1 + radius + (rand() % (gBufferWidth-2*radius-1));
	
	radsquare = (radius*radius);
	
	left=-radius;  right = radius;
	top=-radius;  bottom = radius;
	
	/* Perform edge clipping... */
	if(x - radius < 1) left -= (x-radius-1);
	if(y - radius < 1) top  -= (y-radius-1);
	if(x + radius > gBufferWidth-1) right -= (x+radius-gBufferWidth+1);
	if(y + radius > gBufferWidth-1) bottom-= (y+radius-gBufferWidth+1);
	
	for(cy = top; cy < bottom; cy++)
	{
		for(cx = left; cx < right; cx++)
		{
		square = cy*cy + cx*cx;
		if(square < radsquare)
		{
			dist = sqrt(square*length);
			gSource[gBufferWidth*(cy+y) + cx+x] += (int)((FCos(dist)+0xffff)*(height)) >> 19;
		}
		}
	}
}

/*Draw a large drop in the buffer */
void warpDrop(int x, int y, int radius, int height)
{
/* Scott Scriven's code: */
int cx, cy;
int left,top,right,bottom;
int square;
int radsquare = radius * radius;
	radsquare = (radius*radius);
	
	height /= 64;
	
	left=-radius; right = radius;
	top=-radius; bottom = radius;
	
	/* Perform edge clipping... */
	if(x - radius < 1) left -= (x-radius-1);
	if(y - radius < 1) top  -= (y-radius-1);
	if(x + radius > gBufferWidth-1) right -= (x+radius-gBufferWidth+1);
	if(y + radius > gBufferWidth-1) bottom-= (y+radius-gBufferWidth+1);
	
	for(cy = top; cy < bottom; cy++)
	{
		for(cx = left; cx < right; cx++)
		{
		square = cy*cy + cx*cx;
		if(square < radsquare)
		{
			gSource[gBufferWidth*(cy+y) + cx+x] += (radius-sqrt(square))*(float)(height);
		}
		}
	}
}

/* Draw a box-shaped drop in the buffer */
void boxDrop (int x, int y, int radius, int height)
{
/* Scott Scriven's code: */
int cx, cy;
int left, top, right, bottom;
	if(x<0) x = 1+radius+ rand()%(gBufferWidth-2*radius-1);
	if(y<0) y = 1+radius+ rand()%(gBufferWidth-2*radius-1);
	
	left=-radius; right = radius;
	top=-radius; bottom = radius;
	
	/* Perform edge clipping... */
	if(x - radius < 1) left -= (x-radius-1);
	if(y - radius < 1) top  -= (y-radius-1);
	if(x + radius > gBufferWidth-1) right -= (x+radius-gBufferWidth+1);
	if(y + radius > gBufferWidth-1) bottom-= (y+radius-gBufferWidth+1);
	
	for(cy = top; cy < bottom; cy++)
	{
		for(cx = left; cx < right; cx++)
		{
		gSource[gBufferWidth*(cy+y) + (cx+x)] = height;
		}
	}
}

/*Draw a light-weight drop in the buffer */
void lightDrop(int x, int y, int pheight)
{
int xm,xp,ym,yp;
int v;
	/* eeshlo's code (with tweaks): */
	xm = x-1;
	if (xm < 0) xm=0;
	xp = x+1; 
	if (xp >= gBufferWidth) xp=gBufferWidth-1;
	ym = y-1; 
	if (ym < 0) ym=0;
	yp = y+1; 
	if (yp >= gBufferWidth) yp=gBufferWidth-1;
	y *= gBufferWidth; 
	ym *= gBufferWidth; 
	yp *= gBufferWidth; 
	
	v = rand()%pheight+1; 
	v < 96 ? v=96 : v; /* clamp to be at least visible! */
	gSource[x+y] += v;
	v -= 10; /* Reduce intensity only a little */
	gSource[xm+y] += v;
	gSource[xp+y] += v; 
	gSource[x+ym] += v; 
	gSource[x+yp] += v; 
	v -= 10;
	gSource[xm+ym] += v; 
	gSource[xp+ym] += v; 
	gSource[xm+yp] += v; 
	gSource[xp+yp] += v;         
}            

/* Trail effect */
void calcTrail(int *x, int *y)
{
	/* Works out an x,y based on magical trig stuff that
	I will never understand - I kid you not!
	Scott Scriven again */
	*x = (gBufferWidth/2)
		+ ((
			(
			(FSin( (gXang* 65) >>8) >>8) *
			(FSin( (gXang*349) >>8) >>8)
			) * ((gBufferWidth-8)/2)
		) >> 16);
	*y = (gBufferWidth/2)
		+ ((
			(
			(FSin( (gYang*377) >>8) >>8) *
			(FSin( (gYang* 84) >>8) >>8)
			) * ((gBufferWidth-8)/2)
		) >> 16);
	gXang += 13;
	gYang += 12;
}

/* Swirling effect */
void calcSwirl(int *x, int *y)
{
	/* Scott Scriven's code: */
	*x = (gBufferWidth/2)
	+ ((
		(FCos(gSwirlAng)) * (25)
		) >> 16);
	*y = (gBufferWidth/2)
	+ ((
		(FSin(gSwirlAng)) * (25)
		) >> 16);
	gSwirlAng += 50; 
}


/* E: function that returns the interpolated (averaged) value at (u, v) */
/*  (This returns a value from the relative position within the gSource[ ] array: the buffer.) */
float getIntensity(float u, float v)
{
float fx, fy; /* the fractional index values */
int ix, iy, ix2, iy2; /* the integer index values */
float val1, val2, val3, val4; /* the intensities from four neighbouring cells */
	u = 0.5f*u + 0.5f;
	fx = u - floor(u);
	v = 0.5f*v + 0.5f;
	fy = v - floor(v);
	
	/* The range is now (0, 1), multiply with gBufferWidth to get final index values.
	The integer part of these are the index values ix, iy */
	fx *= (float)gBufferWidth;
	fy *= (float)gBufferWidth;
	ix = (int)fx;
	iy = (int)fy;
	/* Though shouldn't really happen, still make sure that ix, iy is in range */
	if (ix<0) ix=0; else if (ix>=gBufferWidth) ix=gBufferWidth-1;
	if (iy<0) iy=0; else if (iy>=gBufferWidth) iy=gBufferWidth-1;
	/* ix2, iy2 = ix+1, iy+1 */
	if ((ix2 = ix+1) >= gBufferWidth) ix2 = gBufferWidth-1;
	if ((iy2 = iy+1) >= gBufferWidth) iy2 = gBufferWidth-1;
	/* Now perform the bilinear interpolation.
		This just calculates an average value from four neighbouring cells. */
	val1 = gSource[(ix*gBufferWidth)  + iy ]; /* value of cell (ix,   iy)   */
	val2 = gSource[(ix2*gBufferWidth) + iy ]; /* value of cell (ix+1, iy)   */
	val3 = gSource[(ix*gBufferWidth)  + iy2]; /* value of cell (ix,   iy+1) */
	val4 = gSource[(ix2*gBufferWidth) + iy2]; /* value of cell (ix+1, iy+1) */
	/* For the interpolation we only need the fractional part of fx, fy */
	fx -= floor(fx);
	fy -= floor(fy);
	/* Now calculate and return the final intensity value */
	return ((1.f-fx)*(1.f-fy)*val1 + (1.f-fx)*fy*val3 +  fx*(1.f-fy)*val2 + fx*fy*val4);
}