随机表面生成

#include <iostream> #include <cmath> #include <cstdlib> #include <iomanip> #include <fstream> #include <sstream> #include <string> #include <ctime> using namespace std; const int NX=200; //x方向 const int NY=200; //y方向 const int NZ=90; //z方向 const double PI=3.1415926; double z[NX+1][NY+1],h1[NX+1][NY+1],gauss[NX][NY]; int i,j,k,ii,jj,ip,jp,kp,n,q,lx,ly,M,k1,k2; double kk; double z1,z2,z3,z4,maxz,minz,sumz,minh,Ra,Ry,Rz,a,min,max,miu,sigma; void init(); void surface(); void surface_a(); void outputsurface1(); void outputsurface2(); void output(int m); double Normal(double x,double miu,double sigma); double AverageRandom(double min,double max); double gaussrand(double miu,double sigma,double min,double max); int main() { using namespace std; Ra=1; //平均粗糙度 for (i=0;i<=NX;i++) for (j=0;j<=NY;j++) { gauss[i][j]=gaussrand(0,1,-1,1); } surface(); for (i=0;i<=NX;i++) for (j=0;j<=NY;j++) { h1[i][j]=z[i][j]; } /* surface_a(); for (i=0;i<=NX;i++) for (j=0;j<=NY;j++) { h1[i][j]+=z[i][j]; }*/ outputsurface1(); outputsurface2(); return 0; } //==============================================生成表面形貌 void surface() { srand(unsigned(time(0))); maxz=0; minz=0; sumz=0; lx=10; //凸点距离 ly=10; for (i=0;i<=NX;i+=lx) for (j=0;j<=NY;j+=ly) { z[i][j]=gauss[i][j]; if ((i==0) || (i==NX)) z[i][j]=0; if (i==0 && j==0) z[i][j]=0; if (i==0 && j==NY) z[i][j]=0; if (i==NX && j==0) z[i][j]=0; if (i==NX && j==NY) z[i][j]=0; sumz+=fabs(z[i][j]); } sumz/=int(NX/lx)*int(NY/ly); a=Ra/sumz; for (i=0;i<=NX;i+=lx) for (j=0;j<=NY;j+=ly) { z[i][j]*=a; } /* for (i=0;i<NX;i+=lx) //双线性插值 for (j=0;j<NY;j+=ly) { z1=z[i][j],z2=z[i+lx][j],z3=z[i][j+ly],z4=z[i+lx][j+ly]; for (ii=i;ii<i+lx;ii++) for (jj=j;jj<j+ly;jj++) { z[ii][jj]=(z1*(i+lx-ii)*(j+ly-jj)+z2*(ii-i)*(j+ly-jj)+z3*(i+lx-ii)*(jj-j)+z4*(ii-i)*(jj-j))/lx/ly; } } */ for (i=0;i<NX;i+=lx) //双三次样条插值 for (j=0;j<NY;j+=ly) { z1=z[i][j],z2=z[i+lx][j],z3=z[i][j+ly];z4=z[i+lx][j+ly]; for (ii=i;ii<=i+lx;ii++) for (jj=j;jj<=j+ly;jj++) { z[ii][jj]=z1+3*(z2-z1)*(ii-i)*(ii-i)/lx/lx+3*(z3-z1)*(jj-j)*(jj-j)/ly/ly +9*(z4+z1-z2-z3)*(ii-i)*(ii-i)*(jj-j)*(jj-j)/lx/lx/ly/ly +2*(z1-z3)*(jj-j)*(jj-j)*(jj-j)/ly/ly/ly+2*(z1-z2)*(ii-i)*(ii-i)*(ii-i)/lx/lx/lx +6*(z2+z3-z4-z1)*(ii-i)*(ii-i)*(ii-i)*(jj-j)*(jj-j)/lx/lx/lx/ly/ly +6*(z2+z3-z4-z1)*(ii-i)*(ii-i)*(jj-j)*(jj-j)*(jj-j)/lx/lx/ly/ly/ly +4*(z4+z1-z2-z3)*(ii-i)*(ii-i)*(ii-i)*(jj-j)*(jj-j)*(jj-j)/lx/lx/lx/ly/ly/ly; } }//插值结束 } /*void surface_a() { maxz=0; minz=0; sumz=0; lx=40; //凸点距离 ly=6; for (i=0;i<=NX;i+=lx) for (j=0;j<=NY;j+=ly) { z[i][j]=gauss[i][j]; if (i==0 && j==0) z[i][j]=0; if (i==0 && j==NY) z[i][j]=0; if (i==NX && j==0) z[i][j]=0; if (i==NX && j==NY) z[i][j]=0; sumz+=fabs(z[i][j]); } sumz/=int(NX/lx)*int(NY/ly); a=Ra/sumz; for (i=0;i<=NX;i+=lx) for (j=0;j<=NY;j+=ly) { z[i][j]*=a; } for (i=0;i<NX;i+=lx) //双线性插值 for (j=0;j<NY;j+=ly) { z1=z[i][j],z2=z[i+lx][j],z3=z[i][j+ly],z4=z[i+lx][j+ly]; for (ii=i;ii<i+lx;ii++) for (jj=j;jj<j+ly;jj++) { z[ii][jj]=(z1*(i+lx-ii)*(j+ly-jj)+z2*(ii-i)*(j+ly-jj)+z3*(i+lx-ii)*(jj-j)+z4*(ii-i)*(jj-j))/lx/ly; } } for (i=0;i<NX;i+=lx) //双三次样条插值 for (j=0;j<NY;j+=ly) { z1=z[i][j],z2=z[i+lx][j],z3=z[i][j+ly];z4=z[i+lx][j+ly]; for (ii=i;ii<=i+lx;ii++) for (jj=j;jj<=j+ly;jj++) { z[ii][jj]=z1+3*(z2-z1)*(ii-i)*(ii-i)/lx/lx+3*(z3-z1)*(jj-j)*(jj-j)/ly/ly +9*(z4+z1-z2-z3)*(ii-i)*(ii-i)*(jj-j)*(jj-j)/lx/lx/ly/ly +2*(z1-z3)*(jj-j)*(jj-j)*(jj-j)/ly/ly/ly+2*(z1-z2)*(ii-i)*(ii-i)*(ii-i)/lx/lx/lx +6*(z2+z3-z4-z1)*(ii-i)*(ii-i)*(ii-i)*(jj-j)*(jj-j)/lx/lx/lx/ly/ly +6*(z2+z3-z4-z1)*(ii-i)*(ii-i)*(jj-j)*(jj-j)*(jj-j)/lx/lx/ly/ly/ly +4*(z4+z1-z2-z3)*(ii-i)*(ii-i)*(ii-i)*(jj-j)*(jj-j)*(jj-j)/lx/lx/lx/ly/ly/ly; } }//插值结束 } */ //=======================================================高斯分布 /*double gaussrand() { static double V1, V2, S; static int phase = 0; double guassrand; if ( phase == 0 ) { do { double U1 = (double)rand() / RAND_MAX; double U2 = (double)rand() / RAND_MAX; V1 = 2 * U1 - 1; V2 = 2 * U2 - 1; S = V1 * V1 + V2 * V2; } while(S >= 1 || S == 0); guassrand = V1 * sqrt(-2 * log(S) / S); } else guassrand = V2 * sqrt(-2 * log(S) / S); phase = 1 - phase; return guassrand; } */ double Normal(double x,double miu,double sigma) //概率密度函数 { return 1.0/sqrt(2*PI*sigma)*exp(-1*(x-miu)*(x-miu)/(2*sigma*sigma)); } double AverageRandom(double min,double max) { int minInteger = (int)(min*10000); int maxInteger = (int)(max*10000); int randInteger = rand()*rand(); int diffInteger = maxInteger - minInteger; int resultInteger = randInteger % diffInteger + minInteger; return resultInteger/10000.0; } double gaussrand(double miu,double sigma,double min,double max)//产生正态分布随机数 { double x; double dScope; double y; do { x = AverageRandom(min,max); y = Normal(x, miu, sigma); dScope = AverageRandom(0, Normal(x,miu,sigma)); }while( dScope > y); return x; } //==================================================求一个积分，详情见书 //===================================================输出表面形貌 void outputsurface1() { ostringstream name; name<<"surface1.dat"; ofstream out(name.str().c_str()); out<<"Title=\"SURFACE1\"\n"<<"VARIABLES=\"X\",\"Y\",\"Z\"\n"<<"ZONE T=\"BOX\",I="<<NX+1<<",J=" <<NY+1<<",F=POINT"<<endl; for (j=0;j<=NY;j++) for(i=0;i<=NX;i++) { out<<double(i)<<" "<<double(j)<<" "<<double (h1[i][j])<<endl; } } void outputsurface2() { ostringstream name; name<<"surface2.dat"; ofstream out(name.str().c_str()); out<<"Title=\"SURFACE1\"\n"<<"VARIABLES=\"X\",\"Y\",\"Z\"\n"<<"ZONE T=\"BOX\",I="<<NX+1<<",J=" <<NY+1<<",F=POINT"<<endl; for (j=0;j<=NY;j++) for(i=0;i<=NX;i++) { if(h1[i][j]<=0.0) { out<<double(i)<<" "<<double(j)<<" "<<double (0)<<endl; } else { out<<double(i)<<" "<<double(j)<<" "<<double (h1[i][j])<<endl; } } }