function plot_null_hypotheses(colnull)

% This plots the null hypothesis on the data plot created by
% pot_main_figure.m. Requires trial_forward_integration.m. 
%
% Greg Balco
% Berkeley Geochronology Center
% Feb, 2022
% 
% Accompanies a paper entitled 'REVERSIBLE ICE SHEET THINNING IN THE
% AMUNDSEN SEA EMBAYMENT DURING THE LATE HOLOCENE,' by Balco, Brown,
% Nichols, Venturelli, and others. 
%
% This code is provided solely for review of the accompanying paper. It
% is not licensed for other use or distribution. 

% Assemble input data. All these things have to be passed in a structure. 

% d.data10 - data structure for N10 (not necessary unless misfit needed)
% d.data14 - data structure for N14 (not necessary unless misfit needed)
% d.elv_now - meters present ice surface elevation
% d.t1 - output of generate_te_history
% d.l14 - C-14 decay const
% d.l10 - Be-10 decay const
% d.Dmass, d.Lmass - ice density approx parameters
% d.rhoice - density of glacier ice 

if nargin < 1; colnull = [1 0 0];end

load KP_consts
d.l14 = consts.l14;
d.l10 = consts.l10;
d.l26 = consts.l26;
d.rhoice = consts.rhoice;

load KP_calculated_P
d.P = P;

d.N10inh = 0;
d.N26inh = 0;



load icedensity_params
d.Lmass = Lmass;

% Null data 

d.plot_z_cm = 5:5:135;

d.data10.z = [1];
d.data10.Natomsg_bounds = [1 2 3 4 5];
d.data14.z = [1];
d.data14.Natomsg_bounds = [1 2 3 4 5];
d.data26.z = [1];
d.data26.Natomsg_bounds = [1 2 3 4 5];



%% Core-specific. Make two for each core to span exposure-age data cloud. 
% Note: this includes core elevations, etc., which is redundant with
% the random search wrapper. 

% The age-elevation data defined here are null hypotheses in which the ice
% surface elevation reaches the present elevation in the middle Holocene
% and does not change thereafter. 

times1 = [0 4000 6000 35000];
times2 = [0 7000 9000 35000];

dnull1_H1 = d;
dnull1_H1.Dmass = Dmasses.H1;
dnull1_H1.site_elv_now = 84.2;
dnull1_H1.t1.times = times1;
dnull1_H1.t1.surface_elvs = [120 120 350 350];

dnull2_H1 = dnull1_H1;
dnull2_H1.t1.times = times2;

dnull1_H4 = d;
dnull1_H4.Dmass = Dmasses.H4;
dnull1_H4.site_elv_now = 80;
dnull1_H4.t1.times = times1;
dnull1_H4.t1.surface_elvs = [119.4 119.4 350 350];

dnull2_H4 = dnull1_H4;
dnull2_H4.t1.times = times2;

dnull1_H5 = d;
dnull1_H5.Dmass = Dmasses.H5;
dnull1_H5.site_elv_now = 77.8;
dnull1_H5.t1.times = times1;
dnull1_H5.t1.surface_elvs = [118.7 118.7 350 350];

dnull2_H5 = dnull1_H5;
dnull2_H5.t1.times = times2;

%% Get forward data

% H1

rnull1_H1 = trial_forward_integration(dnull1_H1,1);
rnull2_H1 = trial_forward_integration(dnull2_H1,1);

axes(findobj('tag','ax1l'));
plot(rnull1_H1.plot_N10,dnull1_H1.plot_z_cm,'-','color',colnull,'tag','null_1_H1');
plot(rnull2_H1.plot_N10,dnull2_H1.plot_z_cm,'-','color',colnull,'tag','null_2_H1');
axes(findobj('tag','ax1m'));
plot(rnull1_H1.plot_N14,dnull1_H1.plot_z_cm,'-','color',colnull,'tag','null_1_H1');
plot(rnull2_H1.plot_N14,dnull2_H1.plot_z_cm,'-','color',colnull,'tag','null_2_H1');

axes(findobj('tag','ax1r'));
plot(dnull1_H1.t1.times,dnull1_H1.t1.surface_elvs,'-','color',colnull,'tag','null_1_H1');
plot(dnull2_H1.t1.times,dnull2_H1.t1.surface_elvs,'-','color',colnull,'tag','null_2_H1');

% H4

rnull1_H4 = trial_forward_integration(dnull1_H4,1);
rnull2_H4 = trial_forward_integration(dnull2_H4,1);

axes(findobj('tag','ax4l'));
plot(rnull1_H4.plot_N10,dnull1_H4.plot_z_cm,'-','color',colnull,'tag','null_1_H4');
plot(rnull2_H4.plot_N10,dnull2_H4.plot_z_cm,'-','color',colnull,'tag','null_2_H4');
axes(findobj('tag','ax4m'));
plot(rnull1_H4.plot_N14,dnull1_H4.plot_z_cm,'-','color',colnull,'tag','null_1_H4');
plot(rnull2_H4.plot_N14,dnull2_H4.plot_z_cm,'-','color',colnull,'tag','null_2_H4');

axes(findobj('tag','ax4r'));
plot(dnull1_H4.t1.times,dnull1_H4.t1.surface_elvs,'-','color',colnull,'tag','null_1_H4');
plot(dnull2_H4.t1.times,dnull2_H4.t1.surface_elvs,'-','color',colnull,'tag','null_2_H4');

% H5

rnull1_H5 = trial_forward_integration(dnull1_H5,1);
rnull2_H5 = trial_forward_integration(dnull2_H5,1);

axes(findobj('tag','ax5l'));
plot(rnull1_H5.plot_N10,dnull1_H5.plot_z_cm,'-','color',colnull,'tag','null_1_H5');
plot(rnull2_H5.plot_N10,dnull2_H5.plot_z_cm,'-','color',colnull,'tag','null_2_H5');
axes(findobj('tag','ax5m'));
plot(rnull1_H5.plot_N14,dnull1_H5.plot_z_cm,'-','color',colnull,'tag','null_1_H5');
plot(rnull2_H5.plot_N14,dnull2_H5.plot_z_cm,'-','color',colnull,'tag','null_2_H5');

axes(findobj('tag','ax5r'));
plot(dnull1_H5.t1.times,dnull1_H5.t1.surface_elvs,'-','color',colnull,'tag','null_1_H5');
plot(dnull2_H5.t1.times,dnull2_H5.t1.surface_elvs,'-','color',colnull,'tag','null_2_H5');



%% Spit out max-min

disp('Null hypothesis:');
disp(['H1 - Be-10: ' sprintf('%0.0f',min(rnull1_H1.plot_N10)) ' - ' sprintf('%0.0f',max(rnull2_H1.plot_N10)) ' atoms/g']);
disp(['H1 - C-14: ' sprintf('%0.0f',min(rnull1_H1.plot_N14)) ' - ' sprintf('%0.0f',max(rnull2_H1.plot_N14)) ' atoms/g']);
disp(['H1 - Al-26: ' sprintf('%0.0f',min(rnull1_H1.plot_N26)) ' - ' sprintf('%0.0f',max(rnull2_H1.plot_N26)) ' atoms/g']);


disp(['H4 - Be-10: ' sprintf('%0.0f',min(rnull1_H4.plot_N10)) ' - ' sprintf('%0.0f',max(rnull2_H4.plot_N10)) ' atoms/g']);
disp(['H4 - C-14: ' sprintf('%0.0f',min(rnull1_H4.plot_N14)) ' - ' sprintf('%0.0f',max(rnull2_H4.plot_N14)) ' atoms/g']);


disp(['H5 - Be-10: ' sprintf('%0.0f',min(rnull1_H5.plot_N10)) ' - ' sprintf('%0.0f',max(rnull2_H5.plot_N10)) ' atoms/g']);
disp(['H5 - C-14: ' sprintf('%0.0f',min(rnull1_H5.plot_N14)) ' - ' sprintf('%0.0f',max(rnull2_H5.plot_N14)) ' atoms/g']);










