% 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);