function be_box_model(action)

% Implements hokey box model for transport of Be-10 from the preglacial regolith
% to tills in Minnesota and South Dakota. 
%
% Greg Balco -- 

if nargin < 1
	action = 'start';
end;

if strcmp(action,'start');
	
	close all;
	
	% basic definitions
	
	thefig = figure('resize','off','color',[0.735 0.735 0.735],...
		'pos',[16 375 560 420]);
	theplot = axes('position',[0.13 0.4 0.775 0.5]);
	
	topText = uicontrol(thefig,'style','text','pos',[72 385 433 20],...
		'string','Simplistic box model for Be-10 in tills',...
		'fontsize',14);	
	
	% text boxes
	
	textE = uicontrol(thefig,'style','text','Pos',[30 130 190 15],...
	'string','Erosion amount E (cm)'...
	);	
	textfK = uicontrol(thefig,'style','text','Pos',[30 108 190 15],...
	'string','Rock fraction fK'...
	);
	textfR = uicontrol(thefig,'style','text','Pos',[30 86 190 15],...
	'string','Recycle fraction fR'...
	);
	textID = uicontrol(thefig,'style','text','Pos',[30 64 190 15],...
	'string','Interglacial inventory ID (at/cm2)'...
	);
	textn = uicontrol(thefig,'style','text','Pos',[30 42 190 15],...
	'string','Number of cycles n'...
	);
	textD = uicontrol(thefig,'style','text','Pos',[30 20 190 15],...
	'string','Length of each cycle D (yr)'...
	);
	
	% editable text boxes
	
	showE = uicontrol(thefig,'style','edit','pos',[225 128 60 19],'string','800');
	showfK = uicontrol(thefig,'style','edit','pos',[225 106 60 19],'string','0.1');
	showfR = uicontrol(thefig,'style','edit','pos',[225 84 60 19],'string','0.1');
	showID = uicontrol(thefig,'style','edit','pos',[225 62 60 19],'string','2e10');
	shown = uicontrol(thefig,'style','edit','pos',[225 40 60 19],'string','20');
	showD = uicontrol(thefig,'style','edit','pos',[225 18 60 19],'string','1e5');
	
	% run button
	
	runbutton = uicontrol(thefig,'style','push','pos',[375 88 80 30],...
		'string','Run','callback','be_box_model(''run'')'...
		);
		
	% clear button
		
	clearbutton = uicontrol(thefig,'style','push','pos',[375 53 80 30],...
		'string','Clear','callback','be_box_model(''clear'')'...
		);	
		
	% reset button
	
	defaultsbutton = uicontrol(thefig,'style','push','pos',[375 18 80 30],...
		'string','Defaults','callback','be_box_model');
		
	% soilprofile menu
	
	pptext = uicontrol(thefig,'style','text','pos',[300 130 120 15],...
		'string','Regolith Be-10 profile');
	
	ppicker = uicontrol(thefig,'style','popup','pos',[420 127 90 19],...
	'string','VA Gneiss|VA Granite');
	
	% stash important handles in root object userdata

	h.thefig = thefig;
	h.theplot = theplot;
	h.showE = showE;
	h.showfK = showfK;
	h.showfR = showfR;
	h.shown = shown;
	h.showD = showD;
	h.showID = showID;
	h.ppicker = ppicker;
	
	set(0,'userdata',h);
	
end;


if strcmp(action,'clear');


	h = get(0,'userdata');
	delete(get(h.theplot,'children'));

	
end;



if strcmp(action,'run');

	
% parameters


h = get(0,'userdata');

% assign initial-inventory function I(z)

if get(h.ppicker,'value') == 1;
	% initial inventory function from Pavich paper
	data = [0	0
	5.60E+10	30
	9.40E+10	50
	1.87E+11	90
	4.37E+11	180
	6.21E+11	280
	7.92E+11	420
	8.31E+11	565
	8.45E+11	655
	8.56E+11	745	
	8.63E+11	837.5
	8.72E+11	942.5
	8.75E+11	1032.5
	8.78E+11	1202.5
	8.80E+11	1450
	8.81E+11    5000
	8.81E+11	100000];
elseif get(h.ppicker,'value') == 2;
	% initial inventory function from Brown granite profile
	data = [0 0
	2.02E+11	107
	3.05E+11	230
	3.33E+11	380
	4.20E+11	565
	4.95E+11	750
	5.77E+11	958
	7.56E+11	1340
	9.23E+11	1787
	9.62E+11	2046
	1.02E+12	100000];
end;


IR = data(:,1); % atoms/cm2
Z = data(:,2); % cm

% -------- ADJUSTABLE PARAMETERS HERE... -------------------

% unpack adjustable parameters from text boxes 

E = str2num(get(h.showE,'string'));
fK = str2num(get(h.showfK,'string'));
fR = str2num(get(h.showfR,'string'));
ID = str2num(get(h.showID,'string'));
n = str2num(get(h.shown,'string'));
D = str2num(get(h.showD,'string'));

% other parameters

rhoT = 2; % g/cm3
l10 = 4.62e-7; % /yr

% ----------------------------------------------------------

% simple calculation without decay...

C(1) = (1-fK)*interp1(Z,IR,E) / (E * rhoT);
for j = 2:n;
	C(j) = ( (1-fR)/(E*rhoT) )*((interp1(Z,IR,(E*j)) - interp1(Z,IR,(E*j-E)))*(1-fK) + ID) + ...
	(fR / (rhoT*E*(j-1))) * ( interp1(Z,IR,(E*j-E))*(1-fK) + ID*(j-2) );
end;

% more complicated calculation with decay...

CD(1) = (1-fK)*interp1(Z,IR,E)*exp(-l10*D*n) / (E * rhoT);
for j = 2:n;
	% transported inventory corrected to the end of period j
	inv1 = ( (1-fR)/(E*rhoT) )*((interp1(Z,IR,(E*j)) - interp1(Z,IR,(E*j-E)))*(1-fK)*exp(-l10*D*j) + ID*exp(-l10*D));
	% average of previous deposits
	if j == 2; 
		inv2 = interp1(Z,IR,(j*E-E))*(1-fK)*exp(-l10*2*D);
	else
		inv2 = (interp1(Z,IR,(j*E-E))*(1-fK)*exp(-l10*2*D)) + ((sum(exp(-l10.*D.*(2:j-1))))*ID);
	end;
	CD(j) = exp(-l10*D*(n-j)) * (inv1 + inv2*fR/(rhoT*(j*E-E)));
end;

% plotting...

% include sequence of lowest tills from 
% holes that seem to go to the bottom of the 
% glacial sequence - lines

sdcodata = [1.07E+07 4
	1.24E+07 3
	1.11E+07 2
	1.88E+08 1];
swra3data = [1.67E+07 3
	2.86E+07 2
	1.21E+08 1];
umrb3data = [1.36E+07 5
	1.70E+07 4
	2.37E+07 3
	1.51E+07 2
	4.68E+07 1];

% throw in some data from 
% Wisconsinan tills - blue dots

wiscdata = [7.73E+06
	1.38E+07
	1.43E+07
	9.92E+06
	9.12E+06
	9.73E+06
	9.46E+06
	8.99E+06
	1.23E+07
	4.61E+06];


% plot the observations....	

axes(h.theplot);
semilogy(sdcodata(:,2),sdcodata(:,1),'b','linewidth',2);
hold on;
plot(swra3data(:,2),swra3data(:,1),'g','linewidth',2);
plot(umrb3data(:,2),umrb3data(:,1),'c','linewidth',2);
plot((zeros(size(wiscdata))+n),wiscdata,'o','markerfacecolor','b','markersize',5);
ylabel('atoms/g','fontsize',12);axis([0 n+1 1e6 1e10]);grid on;

% plot the model results...

%plot(C,'o','MarkerFaceColor','g'); % without decay
plot(CD,'o','MarkerFaceColor','y','markeredgecolor','r');

end;