close all;
% load the data and assign variables
% NOTE THAT I INCORRECLTY USE lam = 2*pi*sqrt(4*pi*kappa/f).  THIS SHOULD
% JUST BE lam = sqrt(4*pi*kappa/f);
% HENCE MY DIFFUSIVITY IS TOO LOW COMPARED WITH REALITY.  THE PHYSICS IS
% THE SAME THOUGH.

load boreholeT.txt;
T = boreholeT(:,2);
z = boreholeT(:,1); % z = depth

% take a look at the data
plot(T,-z);
xlabel('Temperature (C)','fontsize',14);
ylabel('Depth (m)');

% take a closer look at the upper part
axis([5 15 -50 0]);
% Note that the seasonal cycle has essentially disappeared
% below ~30 m.  A pretty reasonable estimate of today's mean
% annual temperature is give by the point where the seasonal cycle
% intersects the projection of the temperature-depth line, which is
% pretty much straight up and down.  It looks like about 8° C to me.
% I tried that and changed my mind -- 8.8 looks about right.
hold; plot([8.8 8.8],[-50,0],'k'); % This makes a straight line on the graph.

% Set the mean annual temperature today
Tm = 8.8; %°C

% The equation we want is:
% T(z,t) = T0 + DT * exp(-2*pi*z/lam)*sin(2*pi*(f*t-z/lam+phi);
% where
t = 0; % set to January 1st
phi = 1/4; % phase shift.  This ensures that the seasonal sine
           % wave is at a maximum when t=0 (mid-summer).
f = 1; % f = 1/year because we are looking for the annual cycle

% We are looking for kappa, which is unknown, but let's try
kappa = 1 % m^2/year;
lam = 2*pi*sqrt(4*pi*kappa/f);

DT = max(T)-Tm % the amplitude of the temperature cycle
t = 0 % in years (mid-summer);
Ttry = Tm + DT * exp(-2*pi.*z/lam).*sin(2*pi*(f*t-z/lam+phi));
plot(Ttry,-z,'r');

% Note bad, but let's try a bigger diffusivity;
kappa = 3;
lam = 2*pi*sqrt(4*pi*kappa/f);
Ttry = Tm + DT * exp(-2*pi.*z/lam).*sin(2*pi*(f*t-z/lam+phi));
plot(Ttry,-z,'m');

% Ok, how about kappa = 2.5?;
kappa = 2.5;
lam = 2*pi*sqrt(4*pi*kappa/f);
Tseasonal = Tm + DT * exp(-2*pi.*z/lam).*sin(2*pi*(f*t-z/lam+phi));
plot(Tseasonal,-z,'k+');

% That looks good.  Let's replot.
clf;plot(T,-z,'b'); hold;
plot(Ttry,-z,'k+');
xlabel('Temperature (C)','fontsize',14);
ylabel('Depth (m)');
legend('T','T for kappa = 2.5')

% Now, let's guess that the temperature change occurred 100 years ago, 
% and was 5°C
dt = 100;
dT = 5;

% Here is the transient solution;
T_step_z = dT *erfc(z./(2*(kappa * dt).^0.5) );

% Estimate the geothermal gradient from the slope of the
% lower part of the borehole
index = find(z>100);
zlowerdepths = z(index);
Tlowerdepths = T(index);
p = polyfit(zlowerdepths,Tlowerdepths,1);
geotherm = polyval(p,z);

% Add on initial steady state profile and seasonal cycle
T_soln = T_step_z + geotherm + (Tseasonal - Tm);
figure;
plot(T,-z); hold on;
plot(T_soln,-z,'r');

% Whoa, that's obviously too big a change.
% Try it again.
dT = 2;
T_soln = dT *erfc(z./((2*kappa * dt).^0.5)) + geotherm + (Tseasonal - Tm);
plot(T_soln,-z,'g');

dT = 0.5;
T_soln = dT *erfc(z./((2*kappa * dt).^0.5)) + geotherm + (Tseasonal - Tm);
plot(T_soln,-z,'m');

legend('T','dT = 5','dT = 2','dT = 0.5')

% Not bad, but the change is too high up
% so the change must have occurred earlier in time.
dt = 400;
dT = 2;
T_soln = dT *erfc(z./((2*kappa * dt).^0.5)) + geotherm + (Tseasonal - Tm);
plot(T_soln,-z,'k+');
legend('T','dT = 5','dT = 2','dT = 0.5','dT = 2; dt = 400 years')

%Not terrible, but let's try again with a T change earlier on.
dt = 1000;
dT = 2;
T_soln = dT *erfc(z./((2*kappa * dt).^0.5)) + geotherm + (Tseasonal - Tm);
plot(T_soln,-z,'b+');
legend('T','dT = 5','dT = 2','dT = 0.5','dT = 2; dt = 400 years', 'DT = 2.5; dt = 1000 years')

% I rather like that, though it looks like perhaps I underestimated the magnitude of
% the slope of the geotherm.
% Also, I could have figured out dT ~ 2 at the beginning by noting
% where the geotherm intersects the surface.
% Make a final plot to illustrate this.
figure;
plot(T,-z,'b'); hold on;
plot(geotherm,-z,'g')
plot(T_soln,-z,'r')
legend('T data','geothermal gradient','Solution');
xlabel('Temperature (C)','fontsize',14);
ylabel('Depth (m)','fontsize',14);
title('Solution: kappa = 2.5m^2/year; dt = 1000 years; dT = %2°C','fontsize',14)


% It is useful to take a look at the individual components of the solution.
figure
subplot(2,2,1)
plot(Tseasonal-Tm,-z,'b'); hold on;
plot(Tm*ones(length(z),1),-z,'r');
title('seasonal cycle and mean annual temperature')
legend('seasonal cycle','mean temperature (prior to step change)','location','best')

subplot(2,2,2)
plot(geotherm-Tm,-z,'r'); title('geotherm')

subplot(2,2,3);
plot(T_soln-geotherm-Tm-Tseasonal,-z,'r'); title('transient response to temperature change at surface')

subplot(2,2,4);
plot(T,-z,'b'); hold on; 
plot(T_soln,-z,'r');
title('solution')
legend('original data','geotherm + seasonal + mean + transient','location','best')