oceanutil.cpp

基于FFT的海面模拟。能够实现海面的复制

#include "OceanUtil.h"


float ranf()
{
	return (double)rand()/(double)RAND_MAX;
}

void ouGauss(double work[2])
{
	double x1;
	double x2;
	double w;

	do 
	{
		x1 = 2.0 * ranf() - 1.0;
		x2 = 2.0 * ranf() - 1.0;
		w = x1 * x1 + x2 * x2; 
	} while ( w >= 1.0 );

	w = sqrt( (-2.0 * log( w ) ) / w );
	work[0] = x1 * w;
	work[1] = x2 * w;
}

double ouPhillips(double a,double k[2],double wind[2])
{
	double k2;
	double v2;
	double EL;
	double Phk;
	k2 = k[0]*k[0]+k[1]*k[1];			 
	v2 = wind[0]*wind[0]+wind[1]*wind[1];
	EL = v2 / worldGravity;		    	 

	if (k2==0)
	{
		return 0;
	}

	else
	{
		Phk = a*(exp(-1/(k2*(EL)*(EL)))/(k2*k2))*((k[0]*wind[0]+k[1]*wind[1])*(k[0]*wind[0]+k[1]*wind[1])/(k2*v2))*exp(-sqrt(k2)*1.0);
		return Phk;
	}
}

void   ouNormalize(double vec[])
{
	float Length;
	Length=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
	if (Length==0.0)
	{
		Length=1.0;
		vec[0] /=Length;
		vec[1] /=Length;
		vec[2] /=Length;
	}
}

void   ouCrossProd(double in1[],double in2[],double resulet[3])
{
	resulet[0]=(in1[1]*in2[2])-(in2[1]*in1[2]);
	resulet[1]=(in1[2]*in2[0])-(in2[2]*in1[0]);
	resulet[2]=(in1[0]*in2[1])-(in2[0]*in1[1]);
	ouNormalize(resulet);
}

bool    ouPowerof2(int n,int *m,int *twopm)
{
	if (n <= 1) 
	{
		*m = 0;
		*twopm = 1;
		return(false);
	}
	*m = 1;
	*twopm = 2;

	do 
	{
		(*m)++;
		(*twopm) *= 2;
	} while (2*(*twopm) <= n);

	if (*twopm != n) 
	{
		return(false);
	}
	else
	{
		return(true);
	}
}

int    ouFFT(int dir,int m,double *x,double *y)
{
		long nn,i,i1,j,k,i2,l,l1,l2;
		float c1,c2,tx,ty,t1,t2,u1,u2,z;

		nn = 1 << m;
		i2 = nn >> 1;
		j = 0;
		for (i=0;i<nn-1;i++)
		{
			if (i < j) 
			{     
				tx = x[i];
				ty = y[i];
				x[i] = x[j];
				y[i] = y[j];
				x[j] = tx;
				y[j] = ty;
			}             
			k = i2;

			while (k <= j)
			{
				j -= k;
				k >>= 1;
			}
			j += k;
		}    

		c1 = -1.0;
		c2 = 0.0;
		l2 = 1;
		for (l=0;l<m;l++)
		{
			l1 = l2;
			l2 <<= 1;
			u1 = 1.0;
			u2 = 0.0;
			for (j=0;j<l1;j++)
			{
				for (i=j;i<nn;i+=l2)
				{
					i1 = i + l1;
					t1 = u1 * x[i1] - u2 * y[i1];
					t2 = u1 * y[i1] + u2 * x[i1];
					x[i1] = x[i] - t1;
					y[i1] = y[i] - t2;
					x[i] += t1;
					y[i] += t2;
				}
				z =  u1 * c1 - u2 * c2;
				u2 = u1 * c2 + u2 * c1;
				u1 = z;
			}
			c2 = sqrt((1.0 - c1) / 2.0);
			if (dir == 1)
				c2 = -c2;
			c1 = sqrt((1.0 + c1) / 2.0);
		}

		if (dir == 1) 
		{
			for (i=0;i<nn;i++)
			{
				x[i] /= (float)nn;
				y[i] /= (float)nn;
			}
		}
		return(true);
}

int ouFFT2D(oceanComplex c[][512],int nx,int ny,int dir)
{
		int m;
		int twopm;
		double *real;
		double *imag;

		real = (double *)malloc(nx * sizeof(double));
		imag = (double *)malloc(nx * sizeof(double));

		if (real == NULL || imag == NULL)
		{
			return(false);
		}

		if (!ouPowerof2(nx,&m,&twopm) || twopm != nx)
		{
			return(false);
		}

		for (int j=0;j<ny;j++) 
		{
			for (int i=0;i<nx;i++) 
			{
				real[i] = c[i][j].real;
				imag[i] = c[i][j].imag;
			}

			ouFFT(dir,m,real,imag);

			for (int i=0;i<nx;i++) 
			{
				c[i][j].real = real[i];
				c[i][j].imag = imag[i];
			}
		}
		free(real);
		free(imag);

		real = (double *)malloc(ny * sizeof(double));
		imag = (double *)malloc(ny * sizeof(double));
		if (real == NULL || imag == NULL)
		{
			return(false);
		}
		if (!ouPowerof2(ny,&m,&twopm) || twopm != ny)
		{
			return(false);
		}

		for (int i=0;i<nx;i++)
		{
			for (int j=0;j<ny;j++) 
			{
				real[j] = c[i][j].real;
				imag[j] = c[i][j].imag;
			}

			ouFFT(dir,m,real,imag);

			for (int j=0;j<ny;j++)
			{
				c[i][j].real = real[j];
				c[i][j].imag = imag[j];
			}
		}
		free(real);
		free(imag);

		return(true);
}

float	ouNeg1Pow(int k)
{
	return  pow(-1.0f,k);
}