风。波高频率玫瑰图

工程水文学第三次作业
第一题风向频率玫瑰图
解：
原始数据及整理后如下表：
风向1--3级4--5级>=6级总次数
0.3~5.4m/s 5.5~10.7m/s >=10.8m/s
>=6级频率（%）总频率N 364 154 12 530 0.136970665 6.049538 NNE 201 139 16 356 0.182627554 4.063463 NE 389 221 61 671 0.696267549 7.658943 ENE 238 269 103 610 1.175664878 6.962675 ENE 429 390 104 923 1.187079101 10.53533 ESE 255 181 7 443 0.079899555 5.0565 SE 416 150 1 567 0.011414222 6.471864 SSE 315 181 8 504 0.091313777 5.752768 S 541 223 3 767 0.034242666 8.754708 SSW 332 168 17 517 0.194041776 5.901153 SW 525 304 32 861 0.365255108 9.827645 WSW 263 87 1 351 0.011414222 4.006392 WSW 331 109 5 445 0.057071111 5.079329 WNW 148 111 22 281 0.251112887 3.207396 NW 292 189 19 500 0.21687022 5.707111 NNW 242 128 15 385 0.171213332 4.394476 C 0 0 0 50 0 0.570711 sum 0 0 426 8761 4.862458623 100
用Matlab编程绘图如下：
附1（Matlab程序代码）：
wind=zeros(18,6);%wind矩阵，前三列分别存储着1-3级、4-5级、》=6级风向频率原始
%资料，后三列分别存储总次数、>=6级频率（%）、总频率（%）。

wind(:,1)=[364,201,389,238,429,255,416,315,541,332,525,263,331,148,29 2,242,0,0];
wind(:,2)=[154,139,221,269,390,181,150,181,223,168,304,87,109,111,189 ,128,0,0];
wind(:,3)=[12,16,61,103,104,7,1,8,3,17,32,1,5,22,19,15,0,0];
wind(:,4)=wind(:,1)+wind(:,2)+wind(:,3);
wind(length(wind(:,1))-1,4)=50;
wind(length(wind(:,1)),4)=sum(wind(:,4));
wind(length(wind(:,1)),3)=sum(wind(:,3));
wind(:,5)=wind(:,3)./wind(length(wind(:,1)),4).*100;
wind(:,6)=wind(:,4)./wind(length(wind(:,1)),4).*100;
%*********************************************以上为由原始数据求的后几列数据
theta=1*pi/2:-1*2*pi/16:-1*3/2*pi;%定义十六个方向θ。

nowind=zeros(1,length(theta));
nowind(:,:)=wind(length(wind(:,1))-1,6);%计算无风频率，并将其分布在各方向上a=wind(:,6)'; %a矩阵为中间变量，wind1矩阵存储个方位总频率值。

wind1=[a(1:(length(wind(:,6))-2)),wind(1,6)];
a=wind(:,5)';%wind2矩阵存储<6级风的个方位频率值
wind2=[a(1:(length(wind(:,5))-2)),wind(1,5)];
wind2=wind1-wind2;
%******************************************************************** *****
polar(theta,wind2,'-');
hold on; %绘制<6级风的频率折线图，即内层折线
polar(theta,wind1,'-');
hold on; %绘制总频率折线图，即图像最外围折线
polar(theta,nowind,'-');
hold on; %绘制无风频率圆
for k=1:length(theta)
polar([theta(k),theta(k)],[wind1(k),wind(length(wind(:,1))-1,6)],'-') ;
hold on;
end %绘制连接圆与总频率折线的各个方向的放射线
%******************************************************************** ****
%以下绘制>6级与总频率折线图之间的阴影线
y=10; %设置每个方位阴影线的条数
dettheta=2*pi/(length(wind1)-1);
p=zeros(1,(length(wind1)-1));
q=zeros(1,(length(wind1)-1));
a=zeros(1,(length(wind1)-1));
b=zeros(1,(length(wind1)-1));
alpha=zeros(1,(length(wind1)-1));
det=zeros(1,(length(wind1)-1));
r1=zeros((length(wind1)-1),y-1);
r2=r1;
a1=r1;
l=r1;
bate=r1;
theta1=r1;
theta2=r1; %以上定义数据处理过程中变量
for i=1:(length(wind1)-1)
if wind1(i+1)>wind1(i)
b(i)=p-q;
a(i)=q.*sin(dettheta);
det(i)=b(i)/y;
alpha(i)=atan(a(i)/b(i));
for j=1:y-1
l(i,j)=q*cos(dettheta)+j*det(i);
a1(i,j)=(p-l(i,j))*tan(alpha(i));
r1(i,j)=sqrt(l(i,j).^2+a1(i,j).^2);
bate(i,j)=atan(a1(i,j)./l(i,j));
theta1(i,j)=theta(i)-dettheta+bate(i,j); end
end
if wind1(i+1)<wind1(i)
q=wind1(i+1);p=wind1(i);
b(i)=p-q;
a(i)=q.*sin(dettheta);
det(i)=b(i)/y;
alpha(i)=atan(a(i)/b(i));
for j=1:y-1
l(i,j)=q*cos(dettheta)+(y-j)*det(i);
a1(i,j)=(p-l(i,j))*tan(alpha(i));
r1(i,j)=sqrt(l(i,j).^2+a1(i,j).^2);
bate(i,j)=atan(a1(i,j)./l(i,j));
theta1(i,j)=theta(i)-bate(i,j);
end
end
end %以上为计算阴影线在外围线上的端点的r值，即用r1存储
for i=1:(length(wind1)-1)
if wind2(i+1)>wind2(i)
p=wind2(i+1);q=wind2(i);
b(i)=p-q;
a(i)=q.*sin(dettheta);
det(i)=b(i)/y;
alpha(i)=atan(a(i)/b(i));
for j=1:y-1
l(i,j)=q*cos(dettheta)+j*det(i);
a1(i,j)=(p-l(i,j))*tan(alpha(i));
r2(i,j)=sqrt(l(i,j).^2+a1(i,j).^2);
bate(i,j)=atan(a1(i,j)./l(i,j));
theta2(i,j)=theta(i)-dettheta+bate(i,j); end
end
if wind2(i+1)<wind2(i)
b(i)=p-q;
a(i)=q.*sin(dettheta);
det(i)=b(i)/y;
alpha(i)=atan(a(i)/b(i));
for j=1:y-1
l(i,j)=q*cos(dettheta)+(y-j)*det(i);
a1(i,j)=(p-l(i,j))*tan(alpha(i));
r2(i,j)=sqrt(l(i,j).^2+a1(i,j).^2);
bate(i,j)=atan(a1(i,j)./l(i,j));
theta2(i,j)=theta(i)-bate(i,j);
end
end
end %以上计算阴影线在内维折线上的端点的r值，即用r2存储
for i=1:16
for j=1:y-1
if j-(y-1)<0
polar([theta1(i,j),theta2(i,j+1)],[r1(i,j),r2(i,j+1)],'-'); hold on;
end
end
end %绘制阴影线条
第二题海浪频率玫瑰图
解：
原始数据及整理后如下表：
波向0--0.60m 0.61--1.20m 1.21--1.80m 1.81--2.40m 2.41--3.00m >=3.00m
P(%)
N 6.45 3.15 0.42 0.55 0.49 0.44 NNE 5.63 2.26 0.35 0.55 0.07 0 NE 6.25 1.99 0.68 0.14 0.14 0 ENE 10.64 0.48 0 0.14 0 0 E 2.33 0.21 0 0 0 0 ESE 2.2 0 0 0 0 0 SE 2.82 0 0 0 0 0 SSE 1.72 0 0 0 0 0 S 1.72 0 0 0 0 0 SSW 0.56 0 0 0 0 0 SW 8.24 0.14 0 0 0 0 WSW 9.61 0.07 0.07 0 0 0 WSW 0.69 0.14 0 0 0 0
WNW 5.15 2.06 0.07 0.14 0.07 0 NW 2.95 1.92 0.62 0.07 0 0 NNW 3.09 4.59 1.58 0.41 0.69 0.14 c 0.511
总次数365*4= 1460
用Matlab编程绘制波高玫瑰图如下：
附2（Matlab编程代码如下）：
wave=zeros(6,16); %wave矩阵用于存储6个波级在个方向上的原始数据
wave(1,:)=[6.45 5.63 6.25 10.64 2.33 2.2 2.82 1.72
1.72 ...
0.56 8.24 9.61 0.69 5.15 2.95 3.09];
wave(2,:)=[3.15 2.26 1.99 0.48 0.21 0 0 0 0 0
0.14 ...
0.07 0.14 2.06 1.92 4.59];
wave(3,:)=[0.42 0.35 0.68 0 0 0 0 0 0 0 0 0.07 ...
0 0.07 0.62 1.58];
wave(4,:)=[0.55 0.55 0.14 0.14 0 0 0 0 0 0 0 0 ...
0 0.14 0.07 0.41];
wave(5,:)=[0.49 0.07 0.14 0 0 0 0 0 0 0 0 ...
0 0 0.07 0 0.69];
wave(6,:)=[0.44 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0.14];
n=365*4; %统计总次数
Pc=0;%Pc为无浪频率，以下计算Pc的值
for i=1:length(wave(:,1))
Pc=Pc+sum(wave(i,:));
end
Pc=100-Pc;%计算得到Pc的值
wave=wave';%wave转置，16行表示16个方向
r=zeros(16,2*length(wave(1,:))); %r存储各个方向各波高级到原点的r值
r(:,1)=Pc;
for i=2:2:length(r(1,:))-2
r(:,i)=r(:,i-1)+wave(:,i/2);
r(:,i+1)=r(:,i);
end
r(:,length(r(1,:))-1)=r(:,length(r(1,:))-2);
r(:,length(r(1,:)))=r(:,length(r(1,:))-1)+wave(:,length(r(1,:))/2); %***************************以上得到r矩阵
deta=[0.5 0.5 1.0 1.0 1.5 1.5 2.0 2.0 2.5 2.5 3.0 3.0 ];
%各波级的宽度，即各小矩形的宽
dettheta=zeros(16,length(r(1,:))); %角度Δθ
for i=1:length(r(1,:))
dettheta(:,i)=atan(deta(i)./r(:,i)) ;
end
theta=zeros(16,12); %各角点的θ
t=pi/2:-2*pi/16:-3*pi/2+2*pi/16;
t=t';
for i=1:12
theta(:,i)=t;
end
for i=1:16
for j=2:2:length(r(1,:))-2
if wave(i,j/2)==0
r(i,j)=0;
r(i,j+1)=0;
end
end
end
r(:,length(r(1,:))-1)=r(:,length(r(1,:))-2);
for i=1:16
if wave(i,length(r(1,:))/2)==0
r(i,length(r(1,:)))=0;
end
end
for i=1:16
if wave(:,length(r(1,:))/2)~=0
r(:,length(r(1,:)))=r(:,length(r(1,:))-1)+wave(:,length(r(1,:))/2); end
end
theta1=theta-dettheta;
theta2=theta+dettheta;%以上计算得到θ1和θ2，为同一方位上下或左右角点的θ
r1=r./cos(dettheta);
r2=r1; %以上计算得到r1和r2，为同一方位上下或左右角点的r
for i=1:16;
for j=1:2:11
if r(i,j+1)==0
break;
end
polar([theta1(i,j),theta1(i,j+1)],[r1(i,j),r1(i,j+1)],'-'); hold on;
polar([theta2(i,j),theta2(i,j+1)],[r2(i,j),r2(i,j+1)],'-'); hold on;
end
end
for i=1:16;
for j=3:2:length(theta1(1,:));
polar([theta1(i,j),theta2(i,j)],[r1(i,j),r2(i,j)],'-');
hold on;
end
end
j=length(theta1(1,:));
for i=1:16;
polar([theta1(i,j),theta2(i,j)],[r1(i,j),r2(i,j)],'-');
hold on;
end
%***************************************以上绘制图像上的各小矩形
theta0=pi/2:-2*pi/49:-3*pi/2;
theta0=theta0';
r0=zeros(49,1);
for i=1:50
r0(i,1)=Pc;
end
polar(theta0,r0,'-');
hold on; %绘制无浪频率圆。