CMOD5模型的matlab实现

%% 数据准备
clear
clc
%% 模型系数准备
global c
c=zeros(1,28);
c(1:13)=[-0.318004,-1.416264,0.006174,-0.058233,0.001851,0.023727,0.071528,0.016858,4.889050,-0.541012,-1.762694,0.487292,-0.058692];
c(14:18)=[0.128558,0.010419,0.08296,0.032500,16.730787];
c(19:28)=[1.466281,7.358231,13.328972,-9.012029,0.999416,2.52457,-0.524396,0.407859,2.227098,-0.306554];
v=0:0.1:20;
theta=30;psai=0;p=1.6;
sigma0=zeros(2,length(v));
for i=1:length(v)
    sigma0(1,i)=sigma0_hh(v(i),psai,theta,p);
end
c(1:13)=[-0.688,-0.793,0.338,-0.173,0,0.004,0.111,0.0162,6.34,2.57,-2.18,0.4,-0.6];
c(14:18)=[0.045,0.007,0.33,0.012,22];
c(19:28)=[1.95,3,8.39,-3.44,1.36,5.35,1.99,0.29,3.8,1.53];
%网站上Fortran程序的参数
%cmod5n(5.0,0.0,45.0) =   0.00872497 这是程序给出的一个参考结果
%c=[-0.6878, -0.7957,  0.3380, -0.1728, 0.0000,  0.0040, 0.1103, 0.0159,  6.7329,  2.7713, -2.2885, 0.4971, -0.7250, 0.0450, 0.0066,  0.3222,  0.0120, 22.7000, 2.0813,  3.0000, 8.3659, -3.3428,  1.3236,  6.2437,  2.3893, 0.3249,  4.1590, 1.6930];
for i=1:length(v)
    sigma0(2,i)=sigma0_vv(v(i),psai,theta,p);
end
res=sigma0_vv(5,0,45,p);
%% 修饰区
plot(v,sigma0(1,:),'b',v,sigma0(2,:),'r','LineWidth',1.5);
xlabel('wind speed(m/s)');
ylabel('\sigma_0');
legend('CMOD5-HH','CMOD5-VV')
figure
plot(v,std2dB(sigma0(1,:)),'b',v,std2dB(sigma0(2,:)),'r','LineWidth',1.5);
xlabel('wind speed(m/s)');
ylabel('\sigma_0(dB)');
legend('CMOD5-HH','CMOD5-VV')
figure
plot(v,std2dB(sigma0(2,:)./sigma0(1,:)),'LineWidth',1.5);
xlabel('wind speed(m/s)');
ylabel('PR(dB)');
%% 核心函数
function res=sigma0_vv(v,psai,theta,p)
res=B0(v,theta)*(1+B1(v,theta)*cosd(psai)+B2(v,theta)*cosd(2*psai))^p;
end

function res=sigma0_hh(v,psai,theta,p)
res=(B0(v,theta)*(1+B1(v,theta)*cosd(psai)+B2(v,theta)*cosd(2*psai)))^p;
end

function res=std2dB(sigma0)
res=10*log10(sigma0);
end
%% B0函数群
function res=B0(v,theta)
x=(theta-40)/25;
res=(10^(a0(x)+a1(x)*v))*(f(v,theta)^gamma(x));
end

function res=f(v,theta)
x=(theta-40)/25;
s=a2(x)*v;
s_0=s0(x);
if s<s_0
    res=((s/s_0)^alpha(s_0))*g(s_0);
else
    res=g(s);
end
end

function res=g(s)
res=1/(1+exp(-s));
end

function res=alpha(s0)
res=s0*(1-g(s0));
end

%% 只和x有关的函数
function res=a0(x)
global c
temp1=c(1:4);
temp2=ones(4,1);
for i=2:4
    temp2(i)=temp2(i-1)*x;
end
res=temp1*temp2;
end

function res=a1(x)
global c
res=c(5)+c(6)*x;
end

function res=a2(x)
global c
res=c(7)+c(8)*x;
end

function res=gamma(x)
global c
temp1=c(9:11);temp2=ones(3,1);
for i=2:3
    temp2(i)=temp2(i-1)*x;
end
res=temp1*temp2;
end

function res=s0(x)
global c
res=c(12)+c(13)*x;
end

function res=v0(x)
global c
temp1=c(21:23);temp2=ones(3,1);
for i=2:3
    temp2(i)=temp2(i-1)*x;
end
res=temp1*temp2;
end

function res=d1(x)
global c
temp1=c(24:26);temp2=ones(3,1);
for i=2:3
    temp2(i)=temp2(i-1)*x;
end
res=temp1*temp2;
end

function res=d2(x)
global c
res=c(27)+c(28)*x;
end
%% B1
function res=B1(v,theta)
global c
x=(theta-40)/25;
temp1=c(14)*(1+x)-c(15)*v*(0.5+x-tanh(4*(x+c(16)+c(17)*v)));
temp2=1+exp(0.34*(v-c(18)));
res=temp1/temp2;
end

function res=B2(v,theta)
x=(theta-40)/25;
res=(-d1(x)+d2(x)*v2(v,theta))*exp(-v2(v,theta));
end

function res=v2(v,theta)
global c
x=(theta-40)/25;
y=(v+v0(x))/v0(x);
y0=c(19);n=c(20);
a=y0-(y0-1)/n;
b=1/(n*((y0-1)^(n-1)));
if y<y0
    res=a+b*((y-1)^n);
else
    res=y;
end
end

Donate