工程水文学第三次作业第一题风向频率玫瑰图解：原始数据及整理后如下表：风向1--3 级0.3~5.4m/s 4--5 级5.5~10.7m/s>=6级>=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级频率(％)、总频率(％)。

win d(:,1)=[364,201,389,238,429,255,416,315,541,332,525,263,331,148,292,242,0,0];wi nd(:,2)=[154,139,221,269,390,181,150,181,223,168,304,87,109,111,189,128,0,0];wi nd(:,3)=[12,16,61,103,104,7,1,8,3,17,32,1,5,22,19,15,0,0];win d(:,4)=wi nd(:,1)+w in d(:,2)+wi nd(:,3);win d(le ngth(wi nd(:,1))-1,4)=50;win d(le ngth(w in d(:,1)),4)=sum(w in d(:,4));win d(le ngth(wi nd(:,1)),3)=sum(wi nd(:,3));win d(:,5)=wi nd(:,3)./wi nd(le ngth(wi nd(:,1)),4).*100;win d(:,6)=wi nd(:,4)./wi nd(le ngth(wi nd(:,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)];定义十六个方向0 O计算无风频率，并将其分布在各方向上矩阵存储个方位总频率值。

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)p=wind1(i+1);q=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)A2+a1(i,j)A2);bate(i,j)=atan(a1(i,j)./l(i,j));theta1(i,j)=theta(i)-dettheta+bate(i,j);endendif 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)A2+a1(i,j)A2);bate(i,j)=atan(a1(i,j)./l(i,j));theta1(i,j)=theta(i)-bate(i,j);endendend % 以上为计算阴影线在外围线上的端点的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)A2+a1(i,j)A2); bate(i,j)=atan(a1(i,j)./l(i,j)); theta2(i,j)=theta(i)-dettheta+bate(i,j);endendif wind2(i+1)<wind2(i)q=wind2(i+1);p=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)A2+a1(i,j)A2);bate(i,j)=atan(a1(i,j)./l(i,j));theta2(i,j)=theta(i)-bate(i,j);endendend % 以上计算阴影线在内维折线上的端点的r 值，即用r2 存储for i=1:16for j=1:y-1if j-(y-1)<0 polar([theta1(i,j),theta2(i,j+1)],[r1(i,j),r2(i,j+1)],'-');hold on;endendend % 绘制阴影线条第二题海浪频率玫瑰图解：原始数据及整理后如下表：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 0WNW5.15 2.06 0.07 0.14 0.070 NW 2.951.92 0.62 0.07 0 0 NNW3.094.59 1.58 0.41 0.69 0.14c 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 00.14 ...0.07 0.14 2.06 1.92 4.59];wave(3,:)=[0.42 0.35 0.68 0 0 0 00 0 0 00.07 ...0 0.07 0.62 1.58];0 5% 10%------ 2 ■ti -a.oomI ____ >=3,00mwave(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 00.14];n=365*4; % 统计总次数Pc=0;%Pc 为无浪频率，以下计算Pc 的值for i=1:length(wave(:,1))Pc=Pc+sum(wave(i,:));endPc=100-Pc;% 计算得到Pc 的值wave=wave';%wave 转置，16 行表示16 个方向r=zeros(16,2*length(wave(1,:))); %r 存储各个方向各波高级到原点的r(:,1)=Pc;for i=2:2:length(r(1,:))-2r(:,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); %*************************** deta=[0.5 0.5 1.0 1.0 1.5 1.5 2.0 2.0 2.5以上得到r 矩阵2.53.0 3.0 ];% 各波级的宽度，即各小矩形的宽dettheta=zeros(16,length(r(1,:))); % 角度△ 0for i=1:length(r(1,:))dettheta(:,i)=atan(deta(i)./r(:,i)) ;endtheta=zeros(16,12); % 各角点的0t=pi/2:-2*pi/16:-3*pi/2+2*pi/16;t=t'; for i=1:12theta(:,i)=t;endfor i=1:16for j=2:2:length(r(1,:))-2if wave(i,j/2)==0 r(i,j)=0; r(i,j+1)=0;endend endr(:,length(r(1,:))-1)=r(:,length(r(1,:))-2);for i=1:16if wave(i,length(r(1,:))/2)==0 r(i,length(r(1,:)))=0;endendfor i=1:16if wave(:,length(r(1,:))/2)~=0r(:,length(r(1,:)))=r(:,length(r(1,:))-1)+wave(:,length(r(1,:))/2);endendtheta1=theta-dettheta;theta2=theta+dettheta;% 以上计算得到0 1和e 2，为同一方位上下或左右角点的0 r1=r./cos(dettheta);r2=r1; 以上计算得到r1和r2for i=1:16;for j=1:2:11if r(i,j+1)==0break;endpolar([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;endendj=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:50r0(i,1)=Pc;endpolar(theta0,r0,'-');hold on; % 绘制无浪频率圆。