% ps07_1.m is the m-file answer for problem set 6 question 1 
% created by T.Kenna 19981205 1406
% Modif'd by D. Glover 2002:11:08
% Modif'd 2006:11:09 DMG to correct the name of the file in the header
% Modif'd 2008:11:13 DMG to change the name of the file in the header

clear all
close all

global G
global g
global n

load natrsf.dat			% load in surface water history

Hs=natrsf(:,2);
t=natrsf(:,1);	% time in col. 1, trit in col 2

H=zeros(50,1);He=zeros(50,1);		% initialize data
H(1)=0.321; He(1)=0.179;			% starting values in 1950

tau=10; lambda=1/18; n=length(natrsf);	% other factors

for i=2:n
   H(i)=H(i-1)+(Hs(i)-H(i-1))/tau - lambda*H(i-1);  % c is Tritium
   He(i)=He(i-1)+(0 -He(i-1))/tau + lambda*H(i-1);	    % h is Helium
end

Age=18.*log(1+(He./H));

%clc
disp('We start this problem by modifying Jenkins'' boxtrans.m in which he');
disp('made all the desired assumptions such as a well-mixed water mass, with');
disp('a true age of 10 years (tau = 10) and begins with steady state values');
disp('for both tritium ans helium-3. The change that I made was to insert');
disp('the given equation for age. You could have done this as part of the');
disp('for loop or as a separate operation after you calculated all the');
disp('values of tritium and helium-3; I chose the latter. You also needed');
disp('to change the plot commands to include an age vs time plot. Let''s');
disp('run the model and make some plots. <return>');
pause

figure(1);clf			% plot up time history
subplot(311)
p=plot(t,H,'r');
set(p,'LineWidth',2);
set(gca,'LineWidth',2,'FontSize',12)
title('PS07-1A: ^3H-^3He Box Model: Constant Renewal Rate');
ylabel('^3H Concentration');
axis([1950 2000 0 8])
grid

subplot(312)
p=plot(t,He,'g'); 
set(p,'LineWidth',2);
set(gca,'LineWidth',2,'FontSize',12);
ylabel('^3He Concentration');
axis([1950 2000 0 8])
grid

subplot(313)
p=plot(t,Age,'b'); 
set(p,'LineWidth',2);
set(gca,'LineWidth',2,'FontSize',12);
xlabel('Year');ylabel('^3H-^3He Age');
%axis([1950 2000 0 8])
grid
t1=t;
H1=H;
He1=He;
Age1=Age;
%clc
disp(' ');
disp('In general, what we are seeing is the effect on the tritium-helium age');
disp('of a variable surface water input for tritium which is not only');
disp('varying due to variable input events but also due to radioactive');
disp('decay of tritium to helium, combined with a constant renewal rate of');
disp('.1 y-1. <return>');
pause
%clc
disp(' ');
disp('In summary of figure 1, from 1950-1960, tritium-helium age is a');
disp('constant steady state value for the first half of the decade as both');
disp('tritium and 3He remain unchanged.  The age curve then follows a more');
disp('or less linear decline as tritium initially increases over the helium');
disp('in the second half of the decade due to the radioactive decay lag.'); 
disp(' ');
disp('From 1960-1970, we observe a leveling off of the age curve in the');
disp('beginning of the decade as tritium concentration remains constant and');
disp('3He grows in. This is followed by a drop in the age as tritium begins');
disp('to increase.  We then observe a short leveling off of the age as 3He');
disp('and continues to grow in. Towards the end of the decade, we observe');
disp('the tritium-helium age begin a more or less steady increase. The');
disp('tritium concentration reaches its maximum and begins to decline, while');
disp('3He begins a fairly steady increase. This increase corresponds not');
disp('only to the 3He that is growing in from the earlier tritium inputs but');
disp('also to the "small" percentage (actually larger than the existing 3He');
disp('conc.) of 3He that is decaying from the more recent and significantly');
disp('larger inputs of tritium. <return>');
pause

%clc
disp(' ');
disp('From 1970-1980, the tritium concentration continues to decline and the');
disp('3He attains its maximum value near the end of the decade.  The age');
disp('increase for this decade is nearly constant there is however, a slight');
disp('slope decrease in the age curve as the 3He levels off.  In the');
disp('remaining two decades, the concentrations of both tritium and 3He');
disp('decline although tritium is declining more rapidly.  Ultimately, the');
disp('3He exceeds the tritium concentration; I suspect that the age will');
disp('continue to increase but will ultimately return to its steady state');
disp('value as the tritium input returns to normal levels and the water mass');
disp('is ''flushed of 3He''. <return>');
pause

%clc
disp(' ');
disp('Now for the same problem, but with a varying renewal rate (tau) instead');
disp('of a constant tau=10. Since tau is a function of t only it can be');
disp('inserted into the program before the main for loop. Let''s look a plot');
disp('of the renewal function before we run the model. Remember that tau is');
disp('the percentage of the water mass that is replaced annually. <return>');
pause

load natrsf.dat			% load in surface water history

Hs=natrsf(:,2);
t=natrsf(:,1);	% time in col. 1, trit in col 2

H=zeros(50,1);He=zeros(50,1);		% initialize data
H(1)=0.321; He(1)=0.179;			% starting values in 1950

lambda=1/18; n=length(natrsf);	% other factors
tau=10-5*sin(2*pi*(t-1950)/46);
oldtau=10;

figure(2);clf
subplot(211)
p=plot(t,natrsf(:,2),'r'); 
set(p,'LineWidth',2);
set(gca,'LineWidth',2,'FontSize',12);
xlabel('Year');ylabel('Surface Water Tritium');
title('Surface Water Tritium Concentration and Varying Renewal Rate')

subplot(212)
p=plot(t,tau,'b'); 
set(p,'LineWidth',2);
set(gca,'LineWidth',2,'FontSize',12);
xlabel('Year');ylabel('Renewal Rate (tau)');

%clc
disp(' ');
disp('The renewal function is a sinusoid curve that is lower than the');
disp('constant rate of 10 in the first half of the plot and greater than');
disp('the constant rate in the second half. If we compare the surface tritium');
disp('concentration to this renewal function; we observe that the renewal');
disp('rate is near its lowest value when the surface concentration is at its');
disp('highest value.  This is important in understanding the next plot');
disp('which we can do now. <return>');
pause

for i=2:n
   H(i)=H(i-1)+(Hs(i)-H(i-1))/tau(i) - lambda*H(i-1);  % c is Tritium
   He(i)=He(i-1)+(0 -He(i-1))/tau(i) + lambda*H(i-1);	    % h is Helium
 end
Age=18.*log(1+He./H);

figure(3);clf			% plot up time history
subplot(311)
p=plot(t,H,'r',t,H1,'r--');
set(p,'LineWidth',2);
set(gca,'LineWidth',2,'FontSize',12)
title('PS07-1B: ^3H-^3He Box Model: Climatically Varying Renewal Rate');
ylabel('^3H Concentration');
axis([1950 2000 0 10])
text(1971,8.5,'Note: Dashed lines are for constant renewal rate') 
grid

subplot(312)
p=plot(t,He,'g',t,He1,'g--'); 
set(p,'LineWidth',2);
set(gca,'LineWidth',2,'FontSize',12);
ylabel('^3He Concentration');
axis([1950 2000 0 8])
grid

subplot(313)
p=plot(t,Age,'b',t,Age1,'b--'); 
set(p,'LineWidth',2);
set(gca,'LineWidth',2,'FontSize',12);
xlabel('Year');ylabel('^3H-^3He Age');
%axis([1950 2000 0 8])
grid

%clc
disp(' ');
disp('Now we can observe the effect of a varying renewal rate on the');
disp('concentrations of tritium and helium and the resulting tritium-helium');
disp('age . Although the curves in both figure 1 and figure 2 are similar in');
disp('shape, the effect of a changing renewal rate is evident.');
disp(' ');
disp('From 1950-1960 renewal rate is increasing, but the effect on the curves');
disp('is small.  As tritium increases towards the end of the decade, renewal');
disp('rate is also reaching its maximum; this has the effect of mixing more');
disp('tritium into the water mass while 3He input is still controlled by');
disp('radioactive decay. With a relatively higher tritium concentration than');
disp('in 1a  and a similar 3He concentration, the tritium-helium age is');
disp('driven lower.');
disp(' ');
disp('From 1960-1970, the renewal rate reaches its maximum and begins to');
disp('decline although renewal rate for the entire decade is still greater');
disp('than .1/y. The effect is similar to the previous decade only more');
disp('pronounced.  Renewal rate is faster when tritium concentration is');
disp('high; and although we see a corresponding increase in 3He for this');
disp('decade, the age is still lower than in 1a due to the higher tritium');
disp('concentration.');
pause
%clc
disp(' ');
disp('From 1970-1980, the renewal rate transitions from >.1/yr to <.1/yr.');
disp('Now we begin to see a reverse in the process.  Not only is tritium');
disp('declining due to radioactive decay, but the amount we are mixing in');
disp('is also declining; in effect we are putting less live tritium and');
disp('keeping more dead tritium which results in an older tritium-helium');
disp('age than in 1a by the end of the decade. This effect continues as the');
disp('combination of a reduced renewal and tritium decay limits the amount');
disp('of tritium available for input similarly the 3He is remaining longer');
disp('as evidenced by the relatively broader maximum in the concentration curve.');  
disp('As a result, the age curve is increasing with a larger slope and the 3He');
disp('concentration exceeds that of tritium almost a decade earlier.  Because the'); 
disp('renewal rate begins increasing mid 1980s and continues to do so, I imagine');
disp('the age curve would attain a higher maximum but would probably return to');
disp('steady state faster than the 1a age curve. <return>');
pause

%clc
disp(' ');
disp('End of Problem 1');