function main
clear all
rand('state',0);
randn('state',0);
N_dot=1;
p=1000;ax=[0];ay=[0];az=[10];axi=[ax' ay' az'];
v10=10;alpha=0.16;k=0.03;delta_t=0.1;Uz=v10*(az/10).^alpha;
%产生功率谱
N=p*N_dot;
posi=zeros(N_dot,N_dot,N);
ch=zeros(N_dot,N_dot,N);
S=zeros(N_dot,N_dot,N);
H=zeros(N_dot,N_dot,N);
B=zeros(N_dot,N_dot,p);
G=zeros(N_dot,N_dot,p);
G1=zeros(size(G));
u=zeros(N_dot,p);
for m=1:N
mid1_S_H=zeros(N_dot,N_dot);
mid2_S_H=zeros(N_dot,N_dot);
n(m)=m/N;%频率
f(m)=1200*n(m)/v10;
Pn(m)=4*0.03*v10.^2*f(m).^2/(1+f(m).^2).^(4/3); %风速谱,如计算竖向 Panofsky 谱将其替换
为相应函数
Uf=1;%摩擦速度
Su(m)=Pn(m)/n(m);
for g=1:N_dot
for h=1:N_dot
if g==h
ch(g,h,m)=1;
else
delta_z=sqrt(36*(ax(g)-ax(h)).^2+256*(ay(g)-ay(h)).^2+100*(az(g)-az(h)).^2);
ch(g,h,m)=exp(-2*n(m)*delta_z/(Uz(g)+Uz(h)));
f_star=2*n(m)*abs(az(g)-az(h))/(Uz(g)+Uz(h));
if f_star<=0.1
posi(g,h,m)=pi*f_star/4;
else if f_star<=0.125
posi(g,h,m)=-10*pi+f_star+1.25;
else posi(g,h,m)=pi+2*pi*rand(1);
end
end
end
end
for g=1:N_dot
