同一维粗糙面模拟一样，公式和基本理论神马的我就不敲了 有兴趣的同学可以参考
http://files.cnblogs.com/xd-jinjian/%E7%AC%AC%E4%BA%8C%E7%AB%A0%EF%BC%88%E5%BB%BA%E6%A8%A1%EF%BC%89.pdf
二维粗糙面的效果图如下所示
具体代码如下
// 2D_GS_RS.cpp : 定义控制台应用程序的入口点。
//*************本程序基于Monte Carlo方法生成二维高斯粗糙面**************
//by changwei
//**********************************************************************
#include "stdafx.h"
#include<iostream>
#include<iomanip> //格式化输出控制符必包含的头文件
#include<fstream>
#include<cmath>
#include<stdlib.h>
#include<complex>
#define M 64 //x方向采样点数
#define N 64 //y方向采样点数
#define Mdx 8 //x方向单位波长内的采样点数
#define Ndy 8 //y方向单位波长内的采样点数
#define pi 3.141592627 //定义常数pi
using namespace std;
complex<double> unit_i(0.0,1.0);
complex<double> unit_r(1.0,0.0);
//产生随机数的子函数
void random2(double start,double end,double seed,double *rand2) //在用二维数组作为形参时，必须指明列数
{
double s=65536.0;
double w=2053.0;
double v=13849.0;
double T=0.0;
int m;
for(int i=0;i<N*M;i++)
{
T=0.0;
for(int k=0;k<12;k++)
{
seed=w*seed+v;
m=seed/s;
seed=seed-m*s;
T=T+seed/s;
}
rand2[i]=start+end*(T-6.0);
}
return;
}
//2D傅里叶变换的子程序
void fft2D(double *data1,int *nn,int ndim,int isign)
{
int i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1,k2,n,nprev,nrem,ntot;
double tempi,tempr;
double theta,wi,wpi,wpr,wr,wtemp;
ntot=1;
for(idim=0;idim<ndim;idim++)
{
ntot=ntot*nn[idim];
}
nprev=1;
for(idim=0;idim<ndim;idim++)
{
n=nn[idim];
nrem=ntot/(n*nprev);
ip1=2*nprev;
ip2=ip1*n;
ip3=ip2*nrem;
i2rev=1;
for(i2=1;i2<=ip2;i2=i2+ip1)
{
if(i2<i2rev)
{
for(i1=i2;i1<=i2+ip1-2;i1=i1+2)
{
for(i3=i1;i3<=ip3;i3=i3+ip2)
{
i3rev=i2rev+i3-i2;
tempr=data1[i3-1];
tempi=data1[i3];
data1[i3-1]=data1[i3rev-1];
data1[i3]=data1[i3rev];
data1[i3rev-1]=tempr;
data1[i3rev]=tempi;
}
}
}
ibit=ip2/2;
while((ibit>=ip1)&&(i2rev>ibit))
{
i2rev=i2rev-ibit;
ibit=ibit/2;
}
i2rev=i2rev+ibit;
}
ifp1=ip1;
while(ifp1<ip2)
{
ifp2=2*ifp1;
theta=isign*2.0*pi/(ifp2/ip1);
wpr=-2.0*sin(0.5*theta)*sin(0.5*theta);
wpi=sin(theta);
wr=1.0;
wi=0.0;
for(i3=1;i3<=ifp1;i3=i3+ip1)
{
for(i1=i3;i1<=i3+ip1-2;i1=i1+2)
{
for(i2=i1;i2<=ip3;i2=i2+ifp2)
{
k1=i2;
k2=k1+ifp1;
tempr=wr*data1[k2-1]-wi*data1[k2];
tempi=wr*data1[k2]+wi*data1[k2-1];
data1[k2-1]=data1[k1-1]-tempr;
data1[k2]=data1[k1]-tempi;
data1[k1-1]=data1[k1-1]+tempr;
data1[k1]=data1[k1]+tempi;
}
}
wtemp=wr;
wr=wr*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
ifp1=ifp2;
}
nprev=n*nprev;
}
return;
}
void Gauss_roughsurface(double hrms,double lcor,double seed,double dx,double dy,double x[][N],double y[][N],double z[][N])
{
double LX=M*dx; //粗糙面x方向的长度
double LY=N*dy; //粗糙面y方向的长度
double lx,ly; //x，y方向的相关长度
lx=lcor;
ly=lcor;
double kx,ky; //x,y方向的离散空间波数
double rand2[M*N]; //随机数
double s[M][N]; //功率谱
double dataF[2*M*N],dataFR[2*M*N],dataFG[2*M*N];
double tr1,tr2,ti1,ti2;
int Nij;
int NN[2];
NN[0]=N;
NN[1]=M;
//相关长度的变换
lcor=lcor/sqrt(dx*dx+dy*dy);
//计算坐标轴
for(int i=0;i<M;i++)
{
for(int j=0;j<N;j++)
{
x[i][j]=-LX/2.0+j*dx; //x方向坐标
y[i][j]=-LY/2.0+i*dy; //y方向坐标
}
}
//计算功率谱
/*for(int i=0;i<M;i++)
{
kx=2*pi*i/LX; //求x方向的空间波数
for(int j=0;j<N;j++)
{
ky=2*pi*j/LY; //求y方向的空间离散波数
s[i][j]=hrms*hrms*lx*ly/4.0/pi*exp(-kx*kx*lx*lx/4.0-ky*ky*ly*ly/4.0); //求高斯功率谱
//s[i][j]=2.0/sqrt(pi)/lcor*exp(-2.0*(x[i][j]*x[i][j]+y[i][j]*y[i][j])/lcor/lcor);
}
}*/
double tl1,tl2,temp1,temp2;
double xx,yy;
tl1=lcor*lcor/2.0;
tl2=sqrt(pi)*lcor/2.0;
for(int i=0;i<M;i++)
{
xx=-(M-1)/2+i-0.5;
for(int j=0;j<N;j++)
{
yy=-(N-1)/2+j-0.5;
temp1=xx*xx+yy*yy;
temp2=exp(-1.0*temp1/tl1);
s[i][j]=temp2/tl2;
}
}
//滤波函数的二维傅里叶变换
for(int i=0;i<M;i++)
{
for(int j=0;j<N;j++)
{
Nij=i*M+j;
dataFG[2*Nij]=s[i][j];
dataFG[2*Nij+1]=0.0;
}
}
fft2D(dataFG,NN,2,1);
//正态分布随机数并转换为傅里叶变换数组
random2(0.0,1.0,seed,rand2);
for(int i=0;i<M;i++)
{
for(int j=0;j<N;j++)
{
Nij=i*M+j;
dataFR[2*Nij]=rand2[Nij];
dataFR[2*Nij+1]=0.0;
}
}
fft2D(dataFR,NN,2,1);
//滤波函数与随机数的乘积及其二维逆傅里叶变换
for(int i=0;i<M*N;i++)
{
tr1=dataFR[2*i]*dataFG[2*i];
tr2=dataFR[2*i+1]*dataFG[2*i+1];
dataF[2*i]=tr1-tr2;
ti1=dataFR[2*i]*dataFG[2*i+1];
ti2=dataFR[2*i+1]*dataFG[2*i];
dataF[2*i+1]=ti1+ti2;
}
fft2D(dataF,NN,2,-1);
//高度函数的生成
for(int i=0;i<M;i++)
{
for(int j=0;j<N;j++)
{
Nij=i*M+j;
z[i][j]=hrms*dataF[2*Nij]/(M*N);
}
}
return;
}
int _tmain(int argc, _TCHAR* argv[])
{
double wavespeed=3.0*pow(10.0,8); //波速
double frequence=1.0*pow(10.0,9); //频率
double lamd=wavespeed/frequence; //波长
double hrms=0.1*lamd; //均方根高度
double lcor=0.5*lamd; //相关长度
double seed=123.456; //随机数种子
double dx=lamd/Mdx; //x方向采样间距
double dy=lamd/Ndy; //y方向采样间距
double x[M][N]; //二维粗糙面x方向坐标
double y[M][N]; //二维粗糙面y方向坐标
double z[M][N]; //二维粗糙面高度起伏值
Gauss_roughsurface(hrms,lcor,seed,dx,dy,x,y,z);
ofstream foutrs("2DGRS.dat",ios::out);
for(int i=0;i<M;i++)
for(int j=0;j<N;j++)
foutrs<<setiosflags(ios::left)<<setw(8)<<x[i][j]<<setiosflags(ios::left)<<setw(8)<<y[i][j]<<setiosflags(ios::left)<<setw(8)<<z[i][j]<<endl;
foutrs.close();
return 0;
}