% m-file for Exercise 4. % Loads in wind profile data to calculate shear velcoity and test % sensitivity of input vairables close all; clear all; %% For the first question we need to load in the data and declare the function % To simplify the dataset (and make it easier to plot) we can group the % variables by surface type when loading in the data. The data will have % two columns; the first column is wind speed and the second column Zo grass = [10, 0.01; 3, 0.01]; desert = [20, 0.001; 10, 0.001]; urban = [10,0.01;5,0.01]; airfield = [10,0.0001;20,0.0001]; % Next we need to code in the formula to calculate shear velocity kappa = 0.4; Z = 10; % U* = U*k / ln(z/z0) ustar_g = (grass(:,1).*kappa)./log(Z./grass(:,2)); ustar_d = (desert(:,1).*kappa)./log(Z./desert(:,2)); ustar_u = (urban(:,1).*kappa)./log(Z./urban(:,2)); ustar_a = (airfield(:,1).*kappa)./log(Z./airfield(:,2)); % Another way to do this is to group all variables into a 3d matrix and run % a loop 4 times to select each one %data{1} = [grass,desert, urban, airfield]; %for n = 1:2:7 % ustar{n} = (data{1}(:,n).*kappa)./log(Z./data{1}(:,n+1)); %end % Now we can plot the data as requested. Ustar as a function of windspeed figure1 = figure; plot(grass(:,1),ustar_g,'ko'); hold on plot(desert(:,1),ustar_d,'r*'); plot(urban(:,1),ustar_u,'b.'); plot(airfield(:,1),ustar_a,'mo'); xlabel('Wind Speed (ms^-^1)'); ylabel('Shear Velocity (ms^-^1)'); legend('Grass','Desert','Urban','Airfield','Location','northwest'); xlim([0 25]); % For the second part we want to perform an univariate analysis of U and Z0 % for U*. This is done by inputing the full expected range of U while % holding Zo (and the other variables) constant to solve U*. This is then % done for Zo, while holding U & the other variables constant % So for U we can make it equal to a range from 0 to 20 m/s Uz = 0:0.2:20; % Then we solve for U8 giving the following constants z0 = 0.001; % And then solve for U* u_star = (Uz.*kappa)./log(10/z0); % We can now plot this univariate analysis figure2 = figure; plot(Uz,u_star,'m'); xlabel('Wind Speed (ms^-^1)'); ylabel('Shear Velocity (ms^-^1)'); title('Univariate analysis of U_* by wind speed'); % We can now do the same for z0. However as this is a semi-logarthmic % relationship we can range z0 with logarthmically equal spacing through z0 = 10.^(-5:0.1:-1); % Then solve for U* as we did setting U constant U = 10; u_star = (U*kappa)./log(10./z0); % We can now plot this univariate analysis figure3 = figure; plot(z0,u_star,'md'); xlabel('Roughness Length(m)'); ylabel('Shear Velocity (ms^-^1)'); title('Univariate analysis of U_* by roughness length'); % However as you will see in figure 3 the values of U* increase dramtically % as z0 approaches zero in a regular arithmetic plot. We can allow those % values to be visualized better in a plot by creating a semi-logarthmic % plot. This plots U* as a regular axis and z0 as a logarthmic one figure4 = figure; semilogx(z0,u_star,'b'); xlabel('Roughness Length(m)'); ylabel('Shear Velocity (ms^-^1)'); title('Univariate analysis of U_* by roughness length'); %% We can also plot both as a surface plot so that we see the influence of % both Zo and U on U* in one figure - this is called a bivariate % sensitivity analysis % First we can create matrices of the intersections (or nodes) of the range % of both U and Zo using the values we assigned above using meshgrid wind = 0:0.1:20; z0 = 10.^(-7:0.1:-1); [x_wind,y_z0] = meshgrid(wind,z0); z_ustar = (x_wind.*kappa)./log(10./y_z0); % Now we can create our figure using the contourf plot function figure5 = figure('name','Bivariate Sensitivity of Shear Velocity'); % This is how contourf plot is performed. Its input variables are the two % matrix ranges that we calculated from using the meshgrid function and % then the calculated variable (also the same size matrix) from our % algebraic expression. We will also take the two ouput variables C and h % so that we can modulate them for a better looking graph % Because we entered logarthmic values for z0 we ideally would like a % semilog plot. However that is not a function option of contour plot. % Instead we can take the log of the just the z0 values and then change the % axes labels accordingly so as if it was a semi-logarthmic plot [C,h] = contourf(x_wind,log10(y_z0),z_ustar); % We can now specify the color of the filled contours colormap (flipud(autumn)); % There are many choices for a colormap be sure to look these up! b = colorbar; % This displays the colorbar on the plot ylabel(b,'Shear Velocity (ms^-^1)'); % This adds a label to the colorbar % Regular plot labeling of the axes xlabel('Wind Speed (ms^-^1)'); ylabel('Roughness Length (m)'); % Here we are changing the labels for the z0 axis so to create a semi-log % contour plot. ylim([-7 -1]); % This just specifies the range of values z0_labels = num2str(exp(-7:1:-1)','%1.4f'); % This sets the values for the % labels and forces them to be in an exponential format set(gca,'YTickLabel',z0_labels); % This then assigns those vlabels to the plot % Using the clabel function we can add contour labels to each of the % contour lines within the plot itself. The finction just needs the two % output variables from when contourf was called clabel(C,h);