|
|
麻烦你多次了，这回不好意思麻烦你了。
我先自己再研究研究.感觉这东西真难作。
总结一下现在状况：
1。核心代码没敢写，几乎完全抄老外的，实际上我还没看懂。
2。fft用的是fftw,参数范围{1,-1},测试计算结果同老外的一样。
3。动画也加进去了，updata步长为1.
4。运行效果确实有海水的感觉了。
现在的问题，海水中出现2条特别光滑的东西？？？？
sevecol老大你不用替我改，我真不好意识麻烦你了。我是”太笨“。
原码；
int printf(const char* pStr, ...){va_list v1;va_start(v1, pStr);CString str;str.FormatV(pStr, v1);con<<str;return 0;}
#include <z9/math/zfft.h>
typedef std::complex<float> ZC;
//////-----------------------------------------------------------------------------------
#define NX 128
#define NY 128
#define NXTIMESNY NX*NY
#define INV_SQRT_TWO (1.0f)/sqrt(2.0f)
#define MAX_WORLD_X 256
#define MAX_WORLD_Y 256
#define BIG_NX (NX+1)
#define BIG_NY (NY+1)
#define STEPX MAX_WORLD_X/NX
#define MINX 10
#define MAXX 100
#define MINY 10
#define MAXY 100
#define GRAV_CONSTANT 9.81f //gravitational constant metric
typedef struct {
float real,imag;
} COMPLEX;
#define PI 3.141592653589793238462643
const double a_global = 0.0008f;
const double wind_global[2] = {30, 10};
COMPLEX mH0[NX][NY];
double hold_horizontal[NX][NY][4];//store k[0],k[1],klen,klen*klen;
COMPLEX mDeltaX[NX][NY];
COMPLEX mDeltaY[NX][NY];
std::vector<ZVector3> gVS;
ZC c[NX][NY];
ZC cc[NX][NY];
ZCFFT2d gfft;
ZVector3 vs[NX*NY];
SXYZN vs1[(NX-1)*(NY*2 - 2)*3];
long gNNN;
void SetTF_(IDirect3DDevice8* pDevice, ULONG uColor, D3DTEXTUREOP op = D3DTOP_SELECTARG2){
pDevice->SetRenderState(D3DRS_TEXTUREFACTOR, uColor);
pDevice->SetTextureStageState(0, D3DTSS_COLOROP, op);
pDevice->SetTextureStageState(0, D3DTSS_COLORARG2, D3DTA_TFACTOR);
pDevice->SetTextureStageState(0, D3DTSS_ALPHAARG1, D3DTA_TFACTOR);
}
void ClearTF_(IDirect3DDevice8* pDevice){
pDevice->SetTextureStageState(0, D3DTSS_COLOROP, D3DTOP_MODULATE);
pDevice->SetTextureStageState(0, D3DTSS_COLORARG2, D3DTA_CURRENT);
pDevice->SetTextureStageState(0, D3DTSS_ALPHAARG1, D3DTA_TEXTURE);
}
float neg1Pow(int k){
return pow(-1,k);
}
double ranf(){
return (double)rand()/(double)RAND_MAX;
}
void gauss(double work[2]){
// Algorithm by Dr. Everett (Skip) Carter, Jr.
double x1, x2, 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; // first gauss random
work[1] = x2 * w; // second gauss random
}
double phillips(double a,double k[2],const double* wind){
double k2 = k[0]*k[0]+k[1]*k[1];
if (k2==0)
return 0;
double v2 = wind[0]*wind[0]+wind[1]*wind[1];
double EL = v2 / GRAV_CONSTANT;
// the factor *exp(-sqrt(k2)*1.0) can get rid of small waves by increasing 1.0
double ret = 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)*0.002);
// this value must be positive since we take the square root later
return ret;
}
void calculate_ho(){
double horizontal[2];
double root_of_phillips;
double gauss_value[2];// this is where gauss generator returns two random values
for (int i=0;i<NX;i++)
{
for (int j=0;j<NY;j++)
{
// hold_horizontal saves these fixed calculations for later use:k[2],klen,klen*klen
horizontal[0]=hold_horizontal[j][0]=2.0*PI*((double)i-.5*NX)/MAX_WORLD_X;// center the origin
horizontal[1]=hold_horizontal[j][1]=2.0*PI*((double)j-.5*NX)/MAX_WORLD_Y;// center the origin
hold_horizontal[j][3]=hold_horizontal[j][0]*hold_horizontal[j][0]+hold_horizontal[j][1]*hold_horizontal[j][1];
hold_horizontal[j][2]=sqrt(hold_horizontal[j][3]);
// horizontal[0]=hold_horizontal[j][0]=2.0*PI*((float)i-.5*NX)/MAX_WORLD_X;
// horizontal[1]=hold_horizontal[j][1]=2.0*PI*((float)j-.5*NX)/MAX_WORLD_Y;
gauss(gauss_value);
root_of_phillips=sqrt(phillips(a_global,horizontal,wind_global));
//mH0[j].real=INV_SQRT_TWO*gauss_value[0]*root_of_phillips;
//mH0[j].imag=INV_SQRT_TWO*gauss_value[1]*root_of_phillips;
mH0[j].real=root_of_phillips;
mH0[j].imag=root_of_phillips;
}
}
}
void ttt_move_water(){
double time_diff, new_time,kvector[2],klength,wkt;;
static double rate=1.0; //controls rate of movement
gNNN++;
time_diff=gNNN; // time passed since program started *rate
int yHalf = NY/2 + 1;
for (int i = 0; i<yHalf; ++i)
{
int yLine = i*NY;
// Mirror the y line index for calculation of -k
// The line below evalutes yLineMirr = y == 0 ? 0 : (mYSize-y)*mXSize;
// by wrapping the heightmap, since it is periodic.
int yLineMirr = ((NY-i)% NY)*NY;
for (int j = 0; j<NX; ++j)
{
//kvector[0]=2.0*PI*((float)i-.5*NX)/NX;//I was using these
//kvector[1]=2.0*PI*((float)j-.5*NY)/NY;
// kvector[0]=2.0*PI*((float)i-.5*NX)/MAX_WORLD_X;//but I think I should use these
// kvector[1]=2.0*PI*((float)j-.5*NY)/MAX_WORLD_Y;
kvector[0]=hold_horizontal[j][0];
kvector[1]=hold_horizontal[j][1];
// klength=sqrt(kvector[0]*kvector[0]+kvector[1]*kvector[1]);
klength=hold_horizontal[j][2];
wkt = sqrt(klength * GRAV_CONSTANT) * time_diff;
// yLineMirr+mXSize-x is the index of -k
int kNegIndex = yLineMirr*NY + ((NY-j)% NY);
// This is h~(K, t) from the Tessendorf paper.
c[j].real( mH0[j].real*cos(wkt)
+ mH0[j].imag*sin(wkt)
+ mH0[NX - i-1][NY - j-1].real*cos(wkt)
- mH0[NX - i-1][NY -j-1].imag*sin(wkt) );
c[j].imag( mH0[j].imag*cos(wkt)
+ mH0[j].real*sin(wkt)
-mH0[NX - i-1][NY - j-1].imag*cos(wkt)
- mH0[NX - i-1][NY -j-1].real*sin(wkt) );
// Store it for later transformation
// h~(-K) = conj(h~(K))
if (i != yHalf-1)
{
c[NX - i-1][NY - j-1].imag( mH0[j].real*cos(wkt)
+ mH0[j].imag*sin(wkt)
+ mH0[NX - i-1][NY - j-1].real*cos(wkt)
- mH0[NX - i-1][NY -j-1].imag*sin(wkt) );
c[NY - i-1][NY - j-1].real( mH0[j].imag*cos(wkt)
+ mH0[j].real*sin(wkt)
-mH0[NX - i-1][NY - j-1].imag*cos(wkt)
- mH0[NX - i-1][NY -j-1].real*sin(wkt) );
}
}
}
gfft.InvFFT_((ZC*)c, (ZC*)c, FALSE);
}
void ttt_render_tri(IDirect3DDevice8* mpDevice){
///_move_h();
float off = 2;
float xS = 0, zS = 0;
float x = xS, z = zS, y=0;
////vs = new ZVector3[gR*gC];
ZVector3* pV = vs;
ZC* pF = (ZC*)(float*)c;
long k = 1;
for(long i=0; i<NX; i++){
x = xS + i*off;
z = zS;
for(long j=0; j<NY; j++){
pV->x = x; pV->z = z;pV->y = (1/30.0f)*pF->real()*(float)neg1Pow(i+j);///(1/50.0f)*sqrt(pF->real()*pF->real() + pF->imag()*pF->imag());
pV++;
z += off;
pF++;
}k=-k;
}
long triCount = (NX-1)*(NY*2 - 2);
///vs1 = new SXYZN[(gR-1)*(gC*2 - 2)*3];
SXYZN* pV1 = vs1;
pV = vs;
D3DMATERIAL8 mat;
memset(&mat, 0, sizeof(mat));
D3DXCOLOR color(0xff0000ff);
mat.Diffuse = color;
mpDevice->SetMaterial(&mat);
mpDevice->SetRenderState(D3DRS_LIGHTING, TRUE);
mpDevice->SetVertexShader(D3DFVF_XYZ|D3DFVF_NORMAL);
////mpDevice->SetRenderState(D3DRS_FILLMODE, D3DFILL_WIREFRAME);
for(i=0; i<NX-1; i++){
for(long j=0; j<NY-1; j++){
ZVector3 v1,v2,v3, v4, v5, v6, vn;
v1 = vs[i*NX + j];
v2 = vs[(i+1)*NX + j+1];
v3 = vs[(i+1)*NX + j];
v4 = vs[i*NX + j+1];
v5 = v3 - v2;
v6 = v1 - v2;
D3DXVec3Cross(&vn, &v5, &v6);
vn.normal();
///if(vn.y<0)vn.y = -vn.y;
memcpy(&pV1->x, &v1, 12);
memcpy(&pV1->nx, &vn, 12);
pV1++;
memcpy(&pV1->x, &v2, 12);
memcpy(&pV1->nx, &vn, 12);
pV1++;
memcpy(&pV1->x, &v3, 12);
memcpy(&pV1->nx, &vn, 12);
pV1++;
ZMatrix4 mm;
mpDevice->DrawPrimitiveUP(D3DPT_TRIANGLELIST, 1, pV1-3, sizeof(SXYZN));
v5 = v2 - v4;
v6 = v1 - v4;
D3DXVec3Cross(&vn, &v5, &v6);
vn.normal();
///if(vn.y<0)vn.y = -vn.y;
///vn.setXYZ(0,1,0);
memcpy(&pV1->x, &v1, 12);
memcpy(&pV1->nx, &vn, 12);
pV1++;
memcpy(&pV1->x, &v4, 12);
memcpy(&pV1->nx, &vn, 12);
pV1++;
memcpy(&pV1->x, &v2, 12);
memcpy(&pV1->nx, &vn, 12);
////printf("x:%0.3f, y:%0.3f, z:%0.3f, nx:%0.3f, ny:%0.3f, nz:%0.3f\r\n", pV1->x, pV1->y, pV1->z, pV1->nx, pV1->ny, pV1->nz);
pV1++;
mpDevice->DrawPrimitiveUP(D3DPT_TRIANGLELIST, 1, pV1-3, sizeof(SXYZN));
}
}
//long rowCount = (gR-1)*2;
//long colTriCount = (gC-1)*2;
//for(i=0; i<rowCount; i++){
//mpDevice->DrawPrimitiveUP(D3DPT_TRIANGLELIST, colTriCount, vs1+i*colTriCount, sizeof(SXYZN));
//}
///mpDevice->DrawPrimitiveUP(D3DPT_TRIANGLELIST, triCount, vs1, sizeof(SXYZN));
mpDevice->SetRenderState(D3DRS_LIGHTING, FALSE);
}
void OutRender_(IDirect3DDevice8* pDevice){
ttt_move_water();
ttt_render_tri(pDevice);
///printf("coutn:%i, %i, %i, %i\r\n", gLineCount, 64*64, pDevice, gVS.begin());
//SetTF_(pDevice, 0xffff0000);
//pDevice->SetVertexShader(D3DFVF_XYZ);
//pDevice->DrawPrimitiveUP(D3DPT_LINELIST, gLineCount, gVS.begin(), sizeof(SXYZ));
//ClearTF_(pDevice);
}
void OutInit_(){
gNNN = ::GetTickCount();
calculate_ho();
gfft.Init_(NX, NY);
}
|
-
|
发表于 2005-7-26 13:45:00
