首页 > 其他 > 详细

一维粗糙面的模拟

时间:2014-04-19 10:43:04      阅读:445      评论:0      收藏:0      [点我收藏+]

   对于一个给定的随机粗糙面,对光波来说可能呈现很粗糙,而对微波来说却可能呈现的很光滑,这主要是因为随机表面的粗糙度是以波长为度量单位的统计参数来表征的。而描述粗糙面的统计量除功率谱密度外,还有高度起伏的概率密度函数和均方高,相关函数和相关长度,结构函数,特征函数,均方斜度与曲率半径等。下面将逐一介绍它们统计描述,并说明它们和功率谱密度间的相互关系。

bubuko.com,布布扣图1 一维粗糙表面的模拟图(其中x轴和y轴都是以入射波长λbubuko.com,布布扣 为单位的)

1 粗糙面中各个基本参量

图1为一个一维粗糙表面,其高度起伏z=f(x)bubuko.com,布布扣 ,概率密度函数p(f)bubuko.com,布布扣 表征高度起伏的分布情况;均方根高度δbubuko.com,布布扣 又称为表面高度标准离差,对于无限长度的粗糙面,其均方根高度表示为δ=E[fbubuko.com,布布扣2bubuko.com,布布扣(x)]?{E[f(x)]}bubuko.com,布布扣2bubuko.com,布布扣bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣=bubuko.com,布布扣bubuko.com,布布扣?bubuko.com,布布扣 fbubuko.com,布布扣2bubuko.com,布布扣p(f)df?[bubuko.com,布布扣bubuko.com,布布扣?bubuko.com,布布扣 fp(f)df]bubuko.com,布布扣2bubuko.com,布布扣bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣       

Ebubuko.com,布布扣 表示沿整个粗糙面取均值,一般以z=0bubuko.com,布布扣 的平面作为参考面。对于有限长度Lbubuko.com,布布扣 的粗糙面,需要将上式积分上下限改为?L/2 bubuko.com,布布扣~bubuko.com,布布扣L/2bubuko.com,布布扣 ,对该有限表面进行离散,设离散采样点数Nbubuko.com,布布扣 个,离散采样间隔Δxbubuko.com,布布扣 一般在1/10波长范围内,然后通过下式计算δ=1bubuko.com,布布扣N?1bubuko.com,布布扣bubuko.com,布布扣 [bubuko.com,布布扣i=1bubuko.com,布布扣Nbubuko.com,布布扣(fbubuko.com,布布扣ibubuko.com,布布扣)bubuko.com,布布扣2bubuko.com,布布扣?N?(E(f))bubuko.com,布布扣2bubuko.com,布布扣]bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣                    

粗糙面上任意两点间的关联程度,称为相关函数Cbubuko.com,布布扣 ,定义为C(xbubuko.com,布布扣1bubuko.com,布布扣,xbubuko.com,布布扣2bubuko.com,布布扣)=(1/δbubuko.com,布布扣2bubuko.com,布布扣)E[f(xbubuko.com,布布扣1bubuko.com,布布扣)f(xbubuko.com,布布扣2bubuko.com,布布扣)]bubuko.com,布布扣 ,其中δbubuko.com,布布扣 为均方根高度,xbubuko.com,布布扣1bubuko.com,布布扣bubuko.com,布布扣 xbubuko.com,布布扣2bubuko.com,布布扣bubuko.com,布布扣 为表面上任意两点横坐标,常用的相关函数有高斯相关函数和指数相关函数,分别如下
C(xbubuko.com,布布扣1bubuko.com,布布扣,xbubuko.com,布布扣2bubuko.com,布布扣)=exp(?(xbubuko.com,布布扣1bubuko.com,布布扣?xbubuko.com,布布扣2bubuko.com,布布扣)bubuko.com,布布扣2bubuko.com,布布扣/lbubuko.com,布布扣2bubuko.com,布布扣)bubuko.com,布布扣 (高斯),C(xbubuko.com,布布扣1bubuko.com,布布扣,xbubuko.com,布布扣2bubuko.com,布布扣)=exp(?|xbubuko.com,布布扣1bubuko.com,布布扣?xbubuko.com,布布扣2bubuko.com,布布扣|/l)bubuko.com,布布扣 (指数) 
其中lbubuko.com,布布扣 为相关长度,D=xbubuko.com,布布扣1bubuko.com,布布扣?xbubuko.com,布布扣2bubuko.com,布布扣bubuko.com,布布扣 为表面上两点间水平距离,当C(xbubuko.com,布布扣1bubuko.com,布布扣,xbubuko.com,布布扣2bubuko.com,布布扣)=1/ebubuko.com,布布扣 时需要(xbubuko.com,布布扣1bubuko.com,布布扣?xbubuko.com,布布扣2bubuko.com,布布扣)bubuko.com,布布扣2bubuko.com,布布扣/lbubuko.com,布布扣2bubuko.com,布布扣=1bubuko.com,布布扣 ,即此时D=xbubuko.com,布布扣1bubuko.com,布布扣?xbubuko.com,布布扣2bubuko.com,布布扣=lbubuko.com,布布扣 ,即为相关长度lbubuko.com,布布扣
对表面高度相关函数作傅立叶变换,即得到其功率谱密度函数,表示为W(k)=1bubuko.com,布布扣2πbubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣?bubuko.com,布布扣C(D)exp(?ikD)dDbubuko.com,布布扣
对高斯相关函数和指数相关函数分别作如上傅立叶变换,得到高斯谱密度函数和指数谱密度函数,分别如下W(k)=δbubuko.com,布布扣2bubuko.com,布布扣lbubuko.com,布布扣2πbubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣exp(?kbubuko.com,布布扣2bubuko.com,布布扣lbubuko.com,布布扣2bubuko.com,布布扣/4)bubuko.com,布布扣 (高斯),W(k)=δbubuko.com,布布扣2bubuko.com,布布扣lbubuko.com,布布扣π(1+kbubuko.com,布布扣2bubuko.com,布布扣lbubuko.com,布布扣2bubuko.com,布布扣)bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣 (指数)      
均方根斜率sbubuko.com,布布扣 表示为s=E[(df/dx)bubuko.com,布布扣2bubuko.com,布布扣]bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣?bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣 ,均方根斜率sbubuko.com,布布扣 与功率谱密度函数间有如下关系s=[kbubuko.com,布布扣2bubuko.com,布布扣S(k)dk]bubuko.com,布布扣1/2bubuko.com,布布扣bubuko.com,布布扣 ,对于高斯谱密度粗糙面,s=2bubuko.com,布布扣bubuko.com,布布扣bubuko.com,布布扣h/lbubuko.com,布布扣 ;对于指数谱密度粗糙面,因通常表面边缘形状较尖锐,积分不可积,均方根斜率sbubuko.com,布布扣 通常不存在。
结构函数定义为表面上两点高度差的均值,即
D(xbubuko.com,布布扣1bubuko.com,布布扣,xbubuko.com,布布扣2bubuko.com,布布扣)=E[[f(xbubuko.com,布布扣1bubuko.com,布布扣)?f(xbubuko.com,布布扣2bubuko.com,布布扣)]bubuko.com,布布扣2bubuko.com,布布扣]bubuko.com,布布扣            

上面费了好长时间,公式神马的实在是太麻烦了,真心想了解的同学可以下载下面的附件,我在这里就直接贴出代码了   

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

 

第一种方法 不使用傅里叶变换,直接叠加法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
// OneDimensionRoughSurface.cpp : 定义控制台应用程序的入口点。
//
 
#include "stdafx.h"
#include <iostream>
#include <complex>
#include <stdlib.h>
#include <fstream>
#include <time.h>
using namespace std;
 
 
 
 
const complex<double> unit_i(0.0,1.0);
const complex<double> unit_r(1.0,0.0);
const double wavespeed=3.0*pow(10.0,8);  //波速
const double frequence=1.0*pow(10.0,9);  //频率
const double lamda=wavespeed/frequence;  //波长
const double delta=0.2*lamda;           //均方根高度
const double lcor=lamda;             //相关长度
const int    M=1024;                     //总采样点数
const double m=8;                        //一个波长内的采样点数
const double dx=lamda/m;                 //x方向的间隔
const double L=(M-1)*dx;
//const double L=102.4*lamda;
const double pi=3.1415926;
 
 
 
void random(double start,double end,double seed,double size,double *myrand);
double gaussrand();
 
 
int _tmain(int argc, _TCHAR* argv[])
{  
    double *x=new double [M];
    for(int i=0;i<M;++i)
        x[i]=(-0.5*M+i)*dx;
    double *kj=new double [M];
    complex<double> *Fkj=new complex<double> [M];
    double *S=new double [M];
    for(int i=0;i<M;++i)
    {
        kj[i]=(-0.5*M+i)*2*pi/L;
        S[i]=delta*delta*lcor/(2*sqrt(pi))*exp(-kj[i]*kj[i]*lcor*lcor/4);
    }
    double *myrand=new double [M]; 
    random(0.0,1.0,time(0),M,myrand);
 
    Fkj[0]=sqrt(2*pi*L*S[0])*myrand[0];
    Fkj[M/2]=sqrt(2*pi*L*S[M/2])*myrand[1];
 
    for(int i=1;i<M/2;++i)
    {
        Fkj[i]=sqrt(2*pi*L*S[i])*(myrand[2*i]+unit_i*myrand[2*i+1])/sqrt(2.0); 
        Fkj[M-i]=conj(Fkj[i]);
    }  
    complex<double> *res=new complex<double> [M];
    for(int i=0;i<M;++i)
    {
        double xn=(-0.5*M+i)*dx;
        //double xn=(i)*dx;
        for(int j=0;j<M;++j)
        {  
            res[i]+=Fkj[j]*exp(unit_i*kj[j]*xn);
        }
        res[i]/=L;
    }
    ofstream sout("diejia.dat",ios::out);
    for(int i=0;i<M;++i)
        sout<<res[i].real()/lamda<<" "<<res[i].imag()/lamda<<endl;
    sout.close();
 
    delete [] res; 
    delete [] x;
    delete [] S;
    delete [] myrand;
    delete [] kj;
    delete [] Fkj;
    return 0;
}
void random(double start,double end,double seed,double size,double *myrand)  //在用二维数组作为形参时,必须指明列数
{
    double s=65536.0;
    double w=2053.0;
    double v=13849.0;
    double T=0.0;
    int m;
    for(int i=0;i<size;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;
        }
        myrand[i]=start+end*(T-6.0);
    }
}

  第二种,使用傅里叶变换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
// OneDimensionRoughSurface.cpp : 定义控制台应用程序的入口点。
//
 
#include "stdafx.h"
#include <iostream>
#include <complex>
#include <stdlib.h>
#include <fstream>
#include <time.h>
using namespace std;
 
 
 
const complex<double> unit_i(0.0,1.0);
const complex<double> unit_r(1.0,0.0);
const double wavespeed=3.0*pow(10.0,8);  //波速
const double frequence=1.0*pow(10.0,9);  //频率
const double lamda=wavespeed/frequence;  //波长
const double delta=0.2*lamda;           //均方根高度
const double lcor=lamda;             //相关长度
const int    M=1024;                     //总采样点数
const double m=8;                        //一个波长内的采样点数
const double dx=lamda/m;                 //x方向的间隔
const double L=(M-1)*dx;
const double pi=3.1415926;                         
 
 
void random(double start,double end,double seed,double size,double *myrand);
double gaussrand();
inline void swap (float &a, float &b);
void bitrp (float xreal [], float ximag [], int n);
void   IFFT (float xreal [], float ximag [], int n);
 
int _tmain(int argc, _TCHAR* argv[])
{  
    double *x=new double [M];
    for(int i=0;i<M;++i)
        x[i]=-L/2+i*dx;
    double *kj=new double [M];
    complex<double> *Fkj=new complex<double> [M];
    double *S=new double [M];
    for(int i=0;i<=M/2;++i)
    {
        kj[i]=i*2*pi/L;
        S[i]=delta*delta*lcor/(2*sqrt(pi))*exp(-kj[i]*kj[i]*lcor*lcor/4);
    }
    double *myrand=new double [M]; 
    random(0.0,1.0,time(0),M,myrand);
 
    Fkj[0]=sqrt(2*pi*L*S[0])*myrand[0];
    Fkj[M/2]=sqrt(2*pi*L*S[M/2])*myrand[1];
 
    for(int i=1;i<M/2;++i)
    {
        Fkj[i]=sqrt(2*pi*L*S[i])*(myrand[2*i]+unit_i*myrand[2*i+1])/sqrt(2.0); 
        Fkj[M-i]=conj(Fkj[i]);
    }
     
    float xreal [M] = {}, ximag [M] = {};
    for(int i=0;i<M;++i)
    {
        xreal[i]=Fkj[i].real();
        ximag[i]=Fkj[i].imag();        
    }
    IFFT (xreal, ximag, M);
    ofstream sout("fft.dat",ios::out);
    for(int i=0;i<M;++i)
    {
        sout<<xreal[i]*M/L/lamda<<" "<<ximag[i]*M/L/lamda<<endl;
    }
    sout.close();  
     
    delete [] x;
    delete [] S;
    delete [] myrand;
    delete [] kj;
    delete [] Fkj;
    return 0;
}
void random(double start,double end,double seed,double size,double *myrand)  //在用二维数组作为形参时,必须指明列数
{
    double s=65536.0;
    double w=2053.0;
    double v=13849.0;
    double T=0.0;
    int m;
    for(int i=0;i<size;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;
        }
        myrand[i]=start+end*(T-6.0);
    }
}
inline void swap (float &a, float &b)
{
    float t;
    t = a;
    a = b;
    b = t;
}
 
void bitrp (float xreal [], float ximag [], int n)
{
    // 位反转置换 Bit-reversal Permutation
    int i, j, a, b, p;
    for (i = 1, p = 0; i < n; i *= 2)
    {
        p ++;
    }
    for (i = 0; i < n; i ++)
    {
        a = i;
        b = 0;
        for (j = 0; j < p; j ++)
        {
            b = (b << 1) + (a & 1);     // b = b * 2 + a % 2;
            a >>= 1;         // a = a / 2;
        }
        if ( b > i)
        {
            swap (xreal [i], xreal [b]);
            swap (ximag [i], ximag [b]);
        }
    }
}
void   IFFT (float xreal [], float ximag [], int n)
{
    // 快速傅立叶逆变换
    float wreal [M / 2], wimag [M / 2], treal, timag, ureal, uimag, arg;
    int m, k, j, t, index1, index2;
 
    bitrp (xreal, ximag, n);
 
    // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1,   , n / 2 - 1
    arg = 2 * pi / n;
    treal = cos (arg);
    timag = sin (arg);
    wreal [0] = 1.0;
    wimag [0] = 0.0;
    for (j = 1; j < n / 2; j ++)
    {
        wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;
        wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;
    }
    for (m = 2; m <= n; m *= 2)
    {
        for (k = 0; k < n; k += m)
        {
            for (j = 0; j < m / 2; j ++)
            {
                index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;     // 旋转因子 w 的实部在 wreal [] 中的下标为 t
                treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];
                timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2];
                ureal = xreal [index1];
                uimag = ximag [index1];
                xreal [index1] = ureal + treal;
                ximag [index1] = uimag + timag;
                xreal [index2] = ureal - treal;
                ximag [index2] = uimag - timag;
            }
        }
    }
    for (j=0; j < n; j ++)
    {
        xreal [j] /= n;
        ximag [j] /= n;
    }
}

  结果如上图所示

一维粗糙面的模拟,布布扣,bubuko.com

一维粗糙面的模拟

原文:http://www.cnblogs.com/xd-jinjian/p/3674149.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!