% Soil Physics Theory 6583 
%
%DIFFUSIVITYHW calculates the soil thermal diffusivity by applying the sinusoidal
%boundary condition model.
%
% INSTRUCTIONS: There are four missing lines in this code. The blank lines
% are 36, 37, 52, and 53. Use your notes and reference book to write the
% missing code and restore the functionality of the script. See that the
% variable names are already assigned. SUGGESTION: Prior
% writing any missing code, study the script and its variables. Try using
% the cell mode CTRL+ENTER OR CTRL+SHIFT+ENTER to run the code cell by cell.
%

% Andres Patrignani and Tyson Ochsner                  24-Oct-2013 13:02:30

%% Load data
num=xlsread('Heat_flux_Hw_spreadsheet.xlsx');
time = num(:,3); % Time 
T0 = num(:,4); % Measured temperature at the surface [°C].
T = num(:,5); % Measured temperature at depth z [°C].
z = -0.025; % Soil Depth [m]
L = 24; % 24 hours

%% Calculate Tmean, Amplitude (A), and constant phase (phi)
T0mean = mean(T0); % Mean temperature of the surface [°C]
A = (max(T0)-min(T0))/2; % Temperature amplitude at the surface [°C].
omega = 2*pi/L; % Angular frequency (t-1).
[~,TmaxIdx] = max(T0); % Find index at which Tmax happens.
timeTmax = time(TmaxIdx); % Time at which Tmax occurs [t]
phi = pi/2-omega*timeTmax; % Phase constant for synchronization [unitless].

%% Calculate thermal diffusivity (alpha) and damping depth (d).
fun = @(d) (sqrt(sum( (T- (T0mean+A*exp(z/d)*sin(omega*time+z/d+phi)) ).^2 ))/length(T)); %RMSE eq between T and Tsim at z.
[d, fval] = fminsearch(fun,0.1);

alpha = % BLANK #1: write the code that calculate the thermal diffusivity [m2 s-1] --> Make sure the result of your code is in these units!
Tsim = % BLANK #2; Write the code that calculate the soil temperature at depth z.

%% Calculate RMSD and plot measured and predicted T at z
RMSE = sqrt(sum((T-Tsim).^2))/length(T);
plot(time,T,'v',time,Tsim,'-k')
legend('T2.5 Measured','T2.5 Predicted','Location','SouthEast')
legend boxoff
display(['RMSE of the soil temperature fit is ',num2str(RMSE),' °C'])
display(['Soil thermal diffusivity is ',num2str(alpha),'m2 s-1'])
ylabel('Soil temperature in °C','FontSize',13)
xlabel('Time in Hours','FontSize',13) 

%% Calculate thermal conductivity (\Lambda) and soil heat flux (G)
lambda = 2.16*10.^6*alpha; % Thermal conductivity (W m-1 K-1)

G0 =  % BLANK #3 Write the code that calculate the soil heat flux at the surface.
G25 = % BLANK #4 Write the code that calculate the soil heat flux at depth z.

figure,plot(time,G0,'-k',time,G25,'--k')
legend('G0','G25','Location','SouthEast')
legend boxoff
ylabel('Soil heat flux W m^-^2','FontSize',13)
xlabel('Time in Hours','FontSize',13)
