风。波高频率玫瑰图

工程水文学第三次作业第一题风向频率玫瑰图解：原始数据及整理后如下表：风向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-1l(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); 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-1l(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);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-1l(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); endendif 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-1l(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);endendend %以上计算阴影线在内维折线上的端点的r值，即用r2存储for i=1:16for j=1:y-1if j-(y-1)<0polar([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.00mP(%)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 0WNW 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.721.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 00.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 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值r(:,1)=Pc;for i=2:2:length(r(1,:))-2r(:,i)=r(:,i-1)+wave(:,i/2);r(:,i+1)=r(:,i);endr(:,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)) ;endtheta=zeros(16,12); %各角点的θt=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)==0r(i,j)=0;r(i,j+1)=0;endendendr(:,length(r(1,:))-1)=r(:,length(r(1,:))-2);for i=1:16if wave(i,length(r(1,:))/2)==0r(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;%以上计算得到θ1和θ2，为同一方位上下或左右角点的θr1=r./cos(dettheta);r2=r1; %以上计算得到r1和r2，为同一方位上下或左右角点的rfor 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;endendfor 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; %绘制无浪频率圆。