利用Matlab模拟年降雨序列！

[wpf8898]

（初次尝试，不一定正确！）

function y = raint(n)
%UNTITLED2 此处显示有关此函数的摘要
%   此处显示详细说明
y1=[-3.5000   -3.3000   -2.4000   -2.1000   -1.2000   -0.8000    0.8000    1.2000    2.9000    3.7000    4.5000    5.5000];
m=normpdf(y1,1,2.2);
k1=[];
k2=[];
y=[];
for h=1:12
    k1=[0,1];k2=[1-m(h)/n*5,m(h)*5/n];
    y=[y,randsrc(1,30,[k1;k2])];
end
y=[y,randsrc(1,5,[k1;k2])];
end

%%生成年际变化序列
d=1:365;
y1=[-3.5000   -3.3000   -2.4000   -2.1000   -1.2000   -0.8000    0.8000    1.2000    2.9000    3.7000    4.5000    5.5000];
m=normpdf(y1,1,2.2);
d1=linspace(1,365,12);
y2=interp1(d1,m,d,'linear');

%%按随机降雨重现周期生成雨量序列
p = exprnd(1,[1,365]);
r1=67734.72*(1+0.841*log10(p))/1123.43448;
r1=r1.*raint(1);

p = exprnd(5,[1,365]);
r2=67734.72*(1+0.841*log10(p))/1123.43448;
r2=r2.*raint(5);

p = exprnd(20,[1,365]);
r3=67734.72*(1+0.841*log10(p))/1123.43448;
r3=r3.*raint(20);

p = exprnd(50,[1,365]);
r4=67734.72*(1+0.841*log10(p))/1123.43448;
r4=r4.*raint(50);

p = exprnd(100,[1,365]);
r5=67734.72*(1+0.841*log10(p))/1123.43448;
r5=r5.*raint(100);

r100=r1+r2+r3+r4+r5;

%for n=1:1:365    
  % if r100(1,n)<0
        r100(1,n)=0;
  %  end
%end

        
k1=[0,1];k2=[0.78,0.22];
k4=randsrc(1,365,[k1;k2]);

t=k4;

r101=r100.*y2.*t.*2.2;
r1000=r101';

s1=[sum(t),sum(r1000)];
s=[mean(r1000),max(r1000)];

disp('年降雨天数、降雨量')
disp(s1)
disp('日平均、最大降雨量')
disp(s)