GopherProxy

	tgenerate_plots.m - cosmo - front and backend for Markov-Chain Monte Carlo inve…
	git clone git://src.adamsgaard.dk/cosmo
	Log
	Files
	Refs
	README
	LICENSE
	---
	tgenerate_plots.m (37503B)
	---
	1 function generate_plots(Ss, matlab_scripts_folder,...
	2 save_file, formats, show_figures, ...
	3 sample_id, name, email, lat, long, ...
	4 be_conc, al_conc, c_conc, ne_conc, ...
	5 be_uncer, al_uncer, c_uncer, ne_uncer, ...
	6 zobs, ...
	7 be_prod_spall, al_prod_spall, c_prod_spall, ne_prod_spall, ...
	8 be_prod_muons, al_prod_muons, c_prod_muons, ne_prod_muons, ...
	9 rock_density, ...
	10 epsilon_gla_min, epsilon_gla_max, ...
	11 epsilon_int_min, epsilon_int_max, ...
	12 t_degla_min, t_degla_max, ...
	13 record, ...
	14 record_threshold_min, record_threshold_max, ...
	15 nwalkers)
	16
	17 %% Copied from m_pakke2014maj11/CompareWalks2.m
	18 % Generates and saves all relevant figures
	19
	20 S = Ss{1};
	21 fixed_stuff = S.fs;
	22 Nwalkers = fixed_stuff.Nwalkers;
	23 M = size(fixed_stuff.mminmax,1);
	24
	25 % put titles on figures?
	26 titles = 0;
	27
	28 % Burn-in and convergence overview MCMC plot, fh(1)
	29 fh(1) = figure('visible', show_figures);
	30
	31 % generate matrices for percentiles
	32 epsilon_int_25 = zeros(Nwalkers,1);
	33 epsilon_int_50 = zeros(Nwalkers,1);
	34 epsilon_int_75 = zeros(Nwalkers,1);
	35 epsilon_gla_25 = zeros(Nwalkers,1);
	36 epsilon_gla_50 = zeros(Nwalkers,1);
	37 epsilon_gla_75 = zeros(Nwalkers,1);
	38 record_threshold_25 = zeros(Nwalkers,1);
	39 record_threshold_50 = zeros(Nwalkers,1);
	40 record_threshold_75 = zeros(Nwalkers,1);
	41 E_25 = zeros(Nwalkers,1);
	42 E_50 = zeros(Nwalkers,1);
	43 E_75 = zeros(Nwalkers,1);
	44
	45 for iwalk=1:min(4,Nwalkers) %only up to the first four walks
	46 QsBurnIns(:,iwalk)=Ss{iwalk}.QsBurnIn;
	47 Qss(:,iwalk) =Ss{iwalk}.Qs;
	48 acceptss(:,iwalk)=Ss{iwalk}.accepts;
	49 normdmpropspr(:,iwalk) = sqrt(sum(Ss{iwalk}.lump_MetHas.dmpropspr.^2)/…
	50 normdmcurspr(:,iwalk) = sqrt(sum(Ss{iwalk}.lump_MetHas.dmcurspr.^2)/M);
	51 end
	52 subplot(2,4,1)
	53 semilogy((1:size(QsBurnIns,1)),QsBurnIns,'-');
	54 ylim([0.1,1e4])
	55 ylabel('St. sum of sq.')
	56 title('Burn in')
	57 grid on
	58
	59 subplot(2,4,2:4)
	60 semilogy(size(QsBurnIns,1)+(1:size(Qss,1)),Qss,'.')
	61 legend('1','2','3','4','5','6','7','8','9','10','location','NorthEast')
	62 ylim([0.1,1e4])
	63 % ylabel('St. sum of sq.')
	64 title(['Sampling. Result file =',save_file],'interp','none')
	65 grid on
	66
	67 subplot(2,4,6)
	68 hist(normdmcurspr)
	69 xlabel('\|\|dm_{cur}\|\|')
	70 ylimit = ylim;
	71 title('Cur2cur step lengths (Stdz)')
	72
	73 subplot(2,4,5)
	74 hist(normdmpropspr)
	75 xlabel('\|\|dm_{prop}\|\|')
	76 ylim(ylimit)
	77 title('Proposed step lengths (Stdz)')
	78
	79 subplot(2,4,7)
	80 textfontsize=8;
	81 fs = Ss{1}.fs;
	82 title('Data and Model')
	83 %text out nucleides, RelError, depths
	84 axis([-0.03,1,-1,9])
	85 str = ['Nucleides'];
	86 for i=1:length(fs.Nucleides)
	87 str = [str,' - ',fs.Nucleides{i}];
	88 end
	89 str = {str,['>> Rel.Error=',num2str(fs.RelErrorObs)],...
	90 ['>> zobs = ',sprintf('%g, ',fs.zobs)]};
	91 text(0,1,str,'fontsize',textfontsize,'interp','none')
	92 %text out prior intervals
	93 for im=1:size(fs.mminmax,1)
	94 %str=[fs.mname{im},'_minmax=[',sprintf('%g,%g]. True=%g',fs.mminmax(im…
	95 text(0,im+2,str,'fontsize',textfontsize,'interp','none')
	96 end
	97
	98 axis ij
	99 box on
	100
	101 subplot(2,4,8)
	102 title('MCMC-analysis')
	103 str = ['Nwalkers:',num2str(fs.Nwalkers),'. ',fs.WalkerStartMode];
	104 text(0,1,str,'fontsize',textfontsize,'interp','none')
	105
	106 fsamp = fs.BurnIn;
	107 str = {['BurnIn: ',num2str(fsamp.Nsamp),' skip ',num2str(fsamp.Nskip),' …
	108 ['>> StepFactor ',fsamp.StepFactorMode]},
	109 switch fsamp.StepFactorMode
	110 case 'Fixed', str{2} = [str{2},'=',num2str(fsamp.StepFactor)];
	111 case 'Adaptive', str{2} = [str{2},', AccRat=',num2str(fsamp.TargetAcce…
	112 end
	113 text(0,2,str,'fontsize',textfontsize,'interp','none')
	114
	115 fsamp = fs.Sampling;
	116 str = {['Sampling: ',num2str(fsamp.Nsamp),' skip ',num2str(fsamp.Nskip),…
	117 ['>> StepFactor ',fsamp.StepFactorMode]},
	118 switch fsamp.StepFactorMode
	119 case 'Fixed', str{2} = [str{2},'=',num2str(fsamp.StepFactor)];
	120 case 'Adaptive', str{2} = [str{2},', AccRat=',num2str(fsamp.TargetAcce…
	121 end
	122 text(0,3.5,str,'fontsize',textfontsize,'interp','none')
	123
	124 axis([-0.03,1,0,10])
	125 axis ij
	126 box on
	127
	128 % Compare sampling cross-plots, fh(2) and fh(3)
	129 fh = [fh;figure('visible', show_figures)];
	130 Nbin = 50;
	131
	132 mminmax = fixed_stuff.mminmax; %bounds of uniform prior intervals
	133 mDistr = fixed_stuff.mDistr; %'uniform' or 'logunif' for each coordina…
	134 % M = size(ms,1);
	135
	136 isub1 = 0;
	137 for i1=1:M
	138 mminmax_i1 = mminmax(i1,:);
	139 % mDistr_i1 = mDistr(i1,:);
	140 % xi = ms(i1,:);
	141 switch mDistr(i1,:)
	142 case 'uniform',
	143 dxbin = diff(mminmax_i1)/Nbin; %we want outside bins to check that…
	144 xbins{i1} = linspace(mminmax_i1(1)-dxbin,mminmax_i1(2)+2*dxbin,Nbi…
	145 case 'logunif'
	146 dfacxbin = (mminmax_i1(2)/mminmax_i1(1))^(1/(Nbin)); %we want outs…
	147 xbins{i1} = logspace(log10(mminmax_i1(1)/dfacxbin),log10(mminmax_i…
	148 end
	149
	150 % for i2 = 1:M
	151 for i2 = (i1+1):M
	152 % yi = ms(i2,:);
	153
	154 mminmax_i2 = mminmax(i2,:);
	155 switch mDistr(i2,:)
	156 case 'uniform',
	157 dybin = diff(mminmax_i2)/Nbin; %we want outside bins to check th…
	158 ybins = linspace(mminmax_i2(1)-dybin,mminmax_i2(2)+2*dybin,Nbin+…
	159 case 'logunif'
	160 dfacybin = (mminmax_i2(2)/mminmax_i2(1))^(1/Nbin); %we want outs…
	161 ybins = logspace(log10(mminmax_i2(1)/dfacybin),log10(mminmax_i2(…
	162 end
	163 W = ones(2); W = conv2(W,W); W = W/sum(W(:));
	164 %loop over walks
	165 isub1 = isub1+1;
	166 for iwalk = 1:min(4,Nwalkers)
	167 xi = Ss{iwalk}.ms(i1,:);yi = Ss{iwalk}.ms(i2,:);
	168 [smoothgrid,histgrid] = smoothdens(xi,yi,xbins{i1},ybins,W);
	169 isub2 = iwalk; isub = isub2 + (isub1-1)*4;
	170 if isub>4*5
	171 isub1 = 1; isub=1;
	172 fh = [fh;figure('visible', show_figures)];
	173 end
	174 subplot(5,4,isub)
	175 pcolor(xbins{i1},ybins,smoothgrid);
	176 xlabel(fixed_stuff.mname{i1})
	177 ylabel(fixed_stuff.mname{i2})
	178 grid on; shading flat; axis tight; set(gca,'fontsize',8);
	179 switch mDistr(i1,:)
	180 case 'uniform', set(gca,'xscale','lin')
	181 case 'logunif', set(gca,'xscale','log')
	182 end
	183 switch mDistr(i2,:)
	184 case 'uniform', set(gca,'yscale','lin')
	185 case 'logunif', set(gca,'yscale','log')
	186 end
	187 axis([mminmax_i1,mminmax_i2])
	188 end
	189 end
	190 end
	191
	192 %% Histogram plots for all four parameters, one subplot per walker, fh(4)
	193 disp('Histogram plots for all four parameters, one subplot per walker')
	194 fh = [fh;figure('visible', show_figures)];
	195 for i1 = 1:M
	196 for iwalk=1:min(4,Nwalkers)
	197 isub = (i1-1)*4 + iwalk;
	198 %subplot(5,4,isub)
	199 subplot(5,4,isub)
	200 Nhistc=histc(Ss{iwalk}.ms(i1,:),xbins{i1});
	201 bar(xbins{i1},Nhistc,'histc')
	202 xlabel(fixed_stuff.mname{i1})
	203 %set (gca,'xtick',[1e-7:1e-3]);
	204
	205 switch mDistr(i1,:)
	206 case 'uniform', set(gca,'xscale','lin','xlim',mminmax(i1,:))
	207 case 'logunif', set(gca,'xscale','log','xlim',mminmax(i1,:))
	208 end
	209 end
	210 end
	211
	212 %% Titles for Mads' plots
	213 if titles
	214 %Putting in titles over figure 2:4
	215 figure(fh(2)); set(fh(2), 'Visible', show_figures)
	216 subplot(5,4,1)
	217 title(['Density cross-plots A. Result file = ',save_file])
	218
	219 figure(fh(3)); set(fh(3), 'Visible', show_figures)
	220 subplot(5,4,1)
	221 title(['Density cross-plots B. Result file = ',save_file])
	222
	223 figure(fh(4)); set(fh(4), 'Visible', show_figures)
	224 subplot(M,Nwalkers,1)
	225 %title(['Histograms. Result file =',save_file],'interp','none')
	226 title('Distribution of model parameter values')
	227 end
	228
	229 %% Plot one histogram per model parameter per walker
	230 epsilon_int_data_w = cell(1, Nwalkers);
	231 epsilon_gla_data_w = cell(1, Nwalkers);
	232 t_degla_data_w = cell(1, Nwalkers);
	233 record_threshold_data_w = cell(1, Nwalkers);
	234
	235 disp('One histogram per model parameter per walker')
	236 fh = [fh;figure('visible', show_figures)];
	237 for i1 = 1:M % for each model parameter
	238 for iwalk=1:min(4,Nwalkers) % for each walker
	239 %isub = (i1-1)*4 + iwalk;
	240 %isub = (iwalk-1)*M + i1;
	241 isub = (i1-1)*Nwalkers + iwalk;
	242 %i1, M, iwalk, Nwalkers, isub
	243 %subplot(5,2,isub)
	244 subplot(M,Nwalkers,isub)
	245 %Nhistc=histc(Ss{iwalk}.ms(i1,:),xbins{i1});
	246 %bar(xbins{i1},Nhistc,'histc')
	247
	248 if i1 == 1 \|\| i1 == 2
	249 % change units from m/yr to m/Myr
	250 histogram(Ss{iwalk}.ms(i1,:)1.0e6, xbins{i1}1.0e6);
	251 else
	252 histogram(Ss{iwalk}.ms(i1,:), xbins{i1});
	253 end
	254
	255 if i1 == 1
	256 title(['MCMC walker ' num2str(iwalk)])
	257 %xlabel('Interglacial erosion rate [m/yr]')
	258 %xlabel('Interglacial erosion rate [m/Myr]')
	259 xlabel('\epsilon_{int} [m/Myr]')
	260 text(0.02,0.98,'a', 'Units', ...
	261 'Normalized', 'VerticalAlignment', 'Top')
	262 epsilon_int_25(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 25);
	263 epsilon_int_50(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 50);
	264 epsilon_int_75(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 75);
	265 epsilon_int_data_w{iwalk} = Ss{iwalk}.ms(i1,:)*1.0e6;
	266
	267 elseif i1 == 2
	268 %xlabel('Glacial erosion rate [m/yr]')
	269 %xlabel('Glacial erosion rate [m/Myr]')
	270 xlabel('\epsilon_{gla} [m/Myr]')
	271 text(0.02,0.98,'b', 'Units', ...
	272 'Normalized', 'VerticalAlignment', 'Top')
	273 epsilon_gla_25(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 25);
	274 epsilon_gla_50(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 50);
	275 epsilon_gla_75(iwalk) = prctile(Ss{iwalk}.ms(i1,:)*1.0e6, 75);
	276 epsilon_gla_data_w{iwalk} = Ss{iwalk}.ms(i1,:)*1.0e6;
	277
	278 elseif i1 == 3
	279 %xlabel('Timing of last deglaciation [yr]')
	280 xlabel('t_{degla} [yr]')
	281 text(0.02,0.98,'c', 'Units', ...
	282 'Normalized', 'VerticalAlignment', 'Top')
	283 t_degla_data_w{iwalk} = Ss{iwalk}.ms(i1,:);
	284
	285 elseif i1 == 4
	286 %xlabel('$\delta^{18}$O$_\mathrm{threshold}$ [$^o/_{oo}$]', ...
	287 %'Interpreter', 'LaTeX')
	288 xlabel(['\delta^{18}O_{threshold} [' char(8240) ']'])
	289 text(0.02,0.98,'d','Units', ...
	290 'Normalized', 'VerticalAlignment', 'Top')
	291 record_threshold_25(iwalk) = prctile(Ss{iwalk}.ms(i1,:), 25);
	292 record_threshold_50(iwalk) = prctile(Ss{iwalk}.ms(i1,:), 50);
	293 record_threshold_75(iwalk) = prctile(Ss{iwalk}.ms(i1,:), 75);
	294 record_threshold_data_w{iwalk} = Ss{iwalk}.ms(i1,:);
	295
	296 else
	297 disp(['Using mname for i1 = ' i1])
	298 xlabel(fixed_stuff.mname{i1})
	299 end
	300 %set (gca,'xtick',[1e-7:1e-3]);
	301
	302 if i1 == 1 \|\| i1 == 2 % shift axes for new units
	303 switch mDistr(i1,:)
	304 case 'uniform', set(gca,'xscale','lin','xlim',mminmax(i1,:)*…
	305 case 'logunif', set(gca,'xscale','log','xlim',mminmax(i1,:)*…
	306 end
	307 else
	308 switch mDistr(i1,:)
	309 case 'uniform', set(gca,'xscale','lin','xlim',mminmax(i1,:))
	310 case 'logunif', set(gca,'xscale','log','xlim',mminmax(i1,:))
	311 end
	312 end
	313 end
	314 end
	315
	316
	317
	318
	319 %% Plot one histogram per model parameter, including data from all walke…
	320 disp('One histogram per model parameter')
	321 fh = [fh;figure('visible', show_figures)];
	322 for i1 = 1:M % for each model parameter
	323
	324 subplot(M,1,i1)
	325
	326 data = [];
	327 for iwalker=1:Nwalkers
	328 if i1 == 1 \|\| i1 == 2
	329 data = [data, Ss{iwalker}.ms(i1,:)*1.0e6];
	330 else
	331 data = [data, Ss{iwalker}.ms(i1,:)];
	332 end
	333 end
	334
	335 hold on
	336 %Nhistc=histc(data, xbins{i1});
	337 %bar(xbins{i1},Nhistc,'histc')
	338 if i1 == 1 \|\| i1 == 2
	339 xbins{i1} = xbins{i1}*1.0e6; % change to m/Myr
	340 end
	341 histogram(data, xbins{i1});
	342
	343 % 2nd quartile = median = 50th percentile
	344 med = median(data);
	345 plot([med, med], get(gca,'YLim'), 'm-')
	346
	347 % save median values for later
	348 if i1 == 1
	349 epsilon_int_med = med;
	350 epsilon_int_data = data;
	351 elseif i1 == 2
	352 epsilon_gla_med = med;
	353 epsilon_gla_data = data;
	354 elseif i1 == 3
	355 t_degla_med = med;
	356 t_degla_data = data;
	357 elseif i1 == 4
	358 record_threshold_med = med;
	359 record_threshold_data = data;
	360 else
	361 error('Unknown parameter');
	362 end
	363 %ylims = get(gca,'YLim');
	364 %text(med, ylims(1) + (ylims(2) - ylims(1))*0.9, ...
	365 %['\leftarrow' num2str(med)]);
	366
	367 % 1st quartile = 25th percentile
	368 prctile25 = prctile(data, 25);
	369 plot([prctile25, prctile25], get(gca,'YLim'), 'm--')
	370
	371 % 3rd quartile = 75th percentile
	372 prctile75 = prctile(data, 75);
	373 plot([prctile75, prctile75], get(gca,'YLim'), 'm--')
	374
	375 hold off
	376
	377 if i1 == 1
	378 %xlabel(['Interglacial erosion rate [m/yr], median = ' ...
	379 %xlabel(['Interglacial erosion rate [m/Myr], median = ' ...
	380 xlabel(['\epsilon_{int} [m/Myr], median = ' ...
	381 num2str(med, 4) ' m/Myr'])
	382 %num2str(med, 4) ' m/yr'])
	383 text(0.02,0.98,'a', 'Units', ...
	384 'Normalized', 'VerticalAlignment', 'Top')
	385 elseif i1 == 2
	386 %xlabel(['Glacial erosion rate [m/yr], median = ' ...
	387 %xlabel(['Glacial erosion rate [m/Myr], median = ' ...
	388 xlabel(['\epsilon_{gla} [m/Myr], median = ' ...
	389 num2str(med, 4) ' m/Myr'])
	390 %num2str(med, 4) ' m/yr'])
	391 text(0.02,0.98,'b', 'Units', ...
	392 'Normalized', 'VerticalAlignment', 'Top')
	393 elseif i1 == 3
	394 %xlabel(['Timing of last deglaciation [yr], median = ' ...
	395 xlabel(['t_{degla} [yr], median = ' ...
	396 num2str(med, 4) ' yr'])
	397 text(0.02,0.98,'c', 'Units', ...
	398 'Normalized', 'VerticalAlignment', 'Top')
	399 elseif i1 == 4
	400 %xlabel(['$\delta^{18}$O$_\mathrm{threshold}$ [$^o/_{oo}$]'...
	401 %', median = ' num2str(med) ' $^o/_{oo}$'], ...
	402 %'Interpreter', 'LaTeX')
	403 xlabel(['\delta^{18}O_{threshold} [' char(8240) '], median = '...
	404 num2str(med, 4) ' ' char(8240)])
	405 text(0.02,0.98,'d','Units', ...
	406 'Normalized', 'VerticalAlignment', 'Top')
	407 else
	408 disp(['Using mname for i1 = ' i1])
	409 xlabel(fixed_stuff.mname{i1})
	410 end
	411 %set (gca,'xtick',[1e-7:1e-3]);
	412
	413 if i1 == 1 \|\| i1 == 2 % shift axes for new units
	414 switch mDistr(i1,:)
	415 case 'uniform', set(gca,'xscale','lin','xlim',mminmax(i1,:)*…
	416 case 'logunif', set(gca,'xscale','log','xlim',mminmax(i1,:)*…
	417 end
	418 else
	419 switch mDistr(i1,:)
	420 case 'uniform', set(gca,'xscale','lin','xlim',mminmax(i1,:))
	421 case 'logunif', set(gca,'xscale','log','xlim',mminmax(i1,:))
	422 end
	423 end
	424 end
	425
	426
	427
	428
	429 %% Plot d18O curve, from PlotSmoothZachos.m
	430 disp('d18O curve, exposure, and erosion history')
	431 fh = [fh;figure('visible', show_figures)];
	432
	433 if strcmp(record, 'rec_5kyr')
	434 d18O_filename = 'lisiecki_triinterp_2p6Ma_5ky.mat'; % zachos_triint…
	435 elseif strcmp(record, 'rec_20kyr')
	436 d18O_filename = 'lisiecki_triinterp_2p6Ma_20ky.mat'; % zachos_triin…
	437 elseif strcmp(record, 'rec_30kyr')
	438 d18O_filename = 'lisiecki_triinterp_2p6Ma_30ky.mat'; % zachos_triin…
	439 else
	440 error(['record ' record ' not understood']);
	441 end
	442
	443 load([matlab_scripts_folder, d18O_filename])
	444 xs = ti;
	445 dt = diff(ti(1:2));
	446 ys = d18O_triang;
	447 %colpos = [0.5,0.5,1];
	448 colpos = [0,0,1];
	449 %colneg = [0.5,1,0.5];
	450 colneg = [1,0,0];
	451 %midvalue = 3.75; %%% THIS IS THE THRESHOLD
	452 midvalue = record_threshold_med;
	453 axh(1)=subplot(3,1,1);
	454 [~,~,~,i_cross]=fill_red_blue(xs,ys,colpos,colneg,midvalue);
	455 text(0.02,0.98,'a', 'Units', 'Normalized', 'VerticalAlignment', 'Top')
	456 xlabel('Age BP [Ma]')
	457 ylabel('\delta^{18}O')
	458 axis tight
	459 xlim([min(xs), max(xs)])
	460 %axis([-0.1,2.7,3.0,5.2])
	461 %axis([-0.05,0.3,3.0,5.2])
	462 %xlim([-0.1,2.7])
	463 %xlim([-0.05,0.3])
	464
	465 hold on
	466 %plot(A(:,1)/1000,A(:,2),'.-m')
	467 axis ij
	468 xs_cross = (i_cross-1)*dt;
	469 xs_cross = [0;xs_cross];
	470 %title([fn,'. minrad = ',num2str(minrad),' Ma'],'interp','none')
	471
	472 % deglaciation timing
	473 plot([t_degla_med, t_degla_med], get(gca,'YLim'), 'k--');
	474
	475 xs_cross(2) = 0.011;
	476
	477 % Exposure
	478 axh(2)=subplot(3,1,2);
	479 %stairs(xs_cross,(1+-1*(-1).^(1:length(xs_cross)))/2,'b','linewidth',1.5…
	480 exposure = (1+-1(-1).^(1:length(xs_cross)))/2 100;
	481 stairs(xs_cross, exposure, ...
	482 'b','linewidth',1.5);
	483 hold on
	484 start1 = [xs_cross(end), xs(end)];
	485 start2 = [exposure(end), exposure(end)];
	486 plot(start1,start2,'b','linewidth',1.5);
	487 %plot(start1,start2,'b','linewidth',1.5);
	488 %title('Exposure. 0 = glaciated, 1 = not glaciated')
	489 %axis([-0.1,2.7,-0.5,1.5])
	490 %axis([-0.05,0.3,-0.5,1.5])
	491 text(0.02,0.98,'b', 'Units', 'Normalized', 'VerticalAlignment', 'Top')
	492 xlabel('Age BP [Ma]')
	493 ylabel('Exposure [%]')
	494 axis tight
	495 xlim([min(xs), max(xs)])
	496 ylim([-30, 130])
	497 % deglaciation timing
	498 plot([t_degla_med, t_degla_med], get(gca,'YLim'), 'k--');
	499
	500 % Erosion rate
	501 axh(2)=subplot(3,1,3);
	502 %stairs(xs_cross,(1+-1*(-1).^(1:length(xs_cross)))/2 ,'r','linewidth',1.…
	503 erosionrate = epsilon_gla_med + ...
	504 -epsilon_gla_med(1+-1(-1).^(1:length(xs_cross)))/2 + ...
	505 epsilon_int_med(1+-1(-1).^(1:length(xs_cross)))/2;
	506 stairs(xs_cross, erosionrate, 'r','linewidth',1.5);
	507 hold on
	508 start1 = [xs_cross(end), xs(end)];
	509 start2 = [erosionrate(end), erosionrate(end)];
	510 plot(start1,start2,'r','linewidth',1.5);
	511 %title('Exposure. 0 = glaciated, 1 = not glaciated')
	512 %axis([-0.1,2.7,-1,2])
	513 %axis([-0.05,0.3,-1,2])
	514 text(0.02,0.98,'c', 'Units', 'Normalized', 'VerticalAlignment', 'Top')
	515 xlabel('Age BP [Ma]')
	516 %ylabel('Erosion rate [m/yr]')
	517 ylabel('Erosion rate [m/Myr]')
	518 axis tight
	519 xlim([min(xs), max(xs)])
	520 ylims = get(gca,'YLim');
	521 dy = ylims(2)-ylims(1);
	522 ylim([ylims(1) - 0.2dy, ylims(2) + 0.2dy])
	523
	524 % deglaciation timing
	525 plot([t_degla_med, t_degla_med], get(gca,'YLim'), 'k--');
	526
	527 hold on
	528 d18Oth = midvalue;
	529 %ErateInt=1e-6; ErateGla=1e-7;
	530 ErateInt = epsilon_int_med;
	531 ErateGla = epsilon_gla_med;
	532 [~,~] = ExtractHistory2(ti,d18O_triang,d18Oth,ErateInt,ErateGla);
	533
	534 linkaxes(axh,'x')
	535 % stairs(-tStarts,relExpos+2,'r')
	536
	537
	538 %% Exhumation history from InspectDepthConcTracks_True_plot.m
	539 disp('Exhumation history');
	540 fh = [fh;figure('visible', show_figures)];
	541 for iwalk = 1:Nwalkers
	542 % iwalk=input(['What iwalk?[1..',num2str(length(Ss)),']']),
	543
	544 % if all walker results are written to a single figure, it exceeds
	545 % system memory limits
	546 subplot(Nwalkers,1,iwalk)
	547 %fh = [fh;figure('visible', show_figures)];
	548
	549 lump_MetHas = Ss{iwalk}.lump_MetHas;
	550 fixed_stuff = Ss{iwalk}.fs;
	551 % if ~isempty(lump_MetHas.zsss{1})
	552 % % % % zsss = lump_MetHas.zsss;
	553 % % % % tss = lump_MetHas.tss;
	554 % % % % dtfine = -1000; tfinemax = -20e5; dzfine = 0.1; zfinemax = 2…
	555 % % % % iz = 1;
	556 % % % % DepthExposureTimeDens(zsss,tss,dtfine,tfinemax,dzfine,zfinem…
	557 % % % % title('z(time)')
	558 % % % % axis ij
	559 % % % % hold on
	560
	561 %making quantiles to plot ontop of exhumation densities:
	562 tss = lump_MetHas.tss;
	563 dtfine = -500; tfinemax = -20e5;
	564 Xxsss = lump_MetHas.zsss;
	565 % dXxfine = 0.05; Xxfinemax = 50; iz=1; %Xx may be depth or nucleide…
	566 dXxfine = 0.01; Xxfinemax = 50; iz=1; %Xx may be depth or nucleide c…
	567 dzfine = dXxfine; zfinemax = Xxfinemax;
	568 [smoothgrid,histgrid,tsfine,Xxfine]=XxTimeDens(Xxsss,tss,dtfine,tfin…
	569 %pcolor(-tsfine,Xxfine,sqrt(smoothgrid));
	570 %pcolor(tsfine,Xxfine,sqrt(histgrid)); %<<<BHJ: v?lge mellem glattet…
	571 hold on
	572 %plot a selection of depth histories
	573 for i=1:ceil(length(tss)/50):length(tss)
	574 % plot(tss{i},Xxsss{i},'.-c')
	575 end
	576 %Compute and plot quantiles
	577 N = sum(histgrid(:,1));
	578 fractions = [0.25,0.5,0.75]; %quartiles
	579 tsfine = 0:dtfine:tfinemax; %bin boundaries
	580 Ntfine = length(tsfine);
	581 zsfine = 0:dzfine:zfinemax; %bin boundaries
	582 % quants{iwalk} = GetHistgridQuantiles(histgrid,N,fractions,tsfine,z…
	583 % plot(tsfine,quants{iwalk}(1,:),'r','linewidth',3)
	584 % plot(tsfine,quants{iwalk}(2,:),'k','linewidth',3)
	585 % plot(tsfine,quants{iwalk}(3,:),'b','linewidth',3)
	586 quants(:,:,iwalk) = GetHistgridQuantiles(histgrid,N,fractions,tsfine…
	587 quants2(:,:,iwalk) = GetHistgridQuantiles2(histgrid,N,fractions,tsfi…
	588
	589 grid on; shading flat; axis tight;
	590 %set(gca,'fontsize',8);
	591 hold on
	592 % lh(1)=plot(tsfine,quants(1,:,iwalk),'r','linewidth',2)
	593 % lh(2)=plot(tsfine,quants(2,:,iwalk),'k','linewidth',2)
	594 % lh(3)=plot(tsfine,quants(3,:,iwalk),'g','linewidth',2)
	595
	596 lh2(1)=plot(-tsfine,quants2(1,:,iwalk),'r','linewidth',1); % 25%
	597 lh2(2)=plot(-tsfine,quants2(2,:,iwalk),'k','linewidth',1); % 50%
	598 lh2(3)=plot(-tsfine,quants2(3,:,iwalk),'r','linewidth',1); % 75%
	599
	600 %legend(lh2,'25%','median','75%','location','northwest')
	601 axis ij
	602 %track_handle=AddTrueModelDepthTrack(Ss{iwalk}.fs,'-m'); %<<< BHJ: a…
	603 %set(track_handle,'linewidth',2);
	604
	605 % Save erosion magnitudes over the last 1 Myr
	606 E_25(iwalk) = quants2(1, 2001, iwalk);
	607 E_50(iwalk) = quants2(2, 2001, iwalk);
	608 E_75(iwalk) = quants2(3, 2001, iwalk);
	609
	610 %axis([0 1e6 0 25])
	611 %caxis([0 25])
	612 %set (gca,'xtick',[0:0.1e5:1e6]);
	613 xlim([0, 1e6]);
	614 %set (gca,'ytick',[0:3:12]);
	615
	616 title(['MCMC walker ' num2str(iwalk)])
	617 xlabel('Time BP [yr]')
	618 ylabel('Depth [m]')
	619 end
	620 colormap(1-copper(15))
	621 %subplot(2,2,1);
	622 %title(fn,'interpreter','none')
	623 %axis([-2e6 0 0 40])
	624
	625 %colorbar
	626 %set (gca,'xtick',[-2e6:0.2e6:0]);
	627
	628 %XTicksAt = ([-2e6:0.2e6:0]);
	629
	630
	631
	632 %% position figure windows at certain coordinates on the screen
	633 if strcmp(show_figures, 'on')
	634 figpos1 = [6 474 1910 504];
	635 figpos2 =[ 12 94 645 887];
	636 figpos3 =[ 610 94 645 887];
	637 figpos4 =[ 1207 94 740 887];
	638 set(fh(2),'pos',figpos2)
	639 set(fh(3),'pos',figpos3)
	640 set(fh(4),'pos',figpos4)
	641 set(fh(1),'pos',figpos1)
	642 figure(fh(1))
	643 end
	644
	645 %% save all figures
	646 disp('Saving figures')
	647 for i=1:length(fh)
	648 disp(['Saving figure ' num2str(i) ' of ' num2str(length(fh))])
	649 figure_save_multiformat(fh(i), ...
	650 strcat(save_file, '-', num2str(i)), ...
	651 formats);
	652 end
	653
	654 % for i1 = 1:M
	655 % hold off
	656 % subplot(M,M,(M-1)*M+i1)
	657 % xlabel(fixed_stuff.mname{i1})
	658 % subplot(M,M,(i1-1)*M+1)
	659 % ylabel(fixed_stuff.mname{i1})
	660 % end
	661 % i=1;
	662
	663 %% generate html table of results
	664 filename = [save_file, '.html'];
	665 disp('Saving html table of results')
	666
	667 % header
	668 html = [...
	669 '\n'...
	670 '<table class="highlight">\n'...
	671 ' <thead>\n'...
	672 ' <tr>\n'...
	673 ' <th data-field="param">Parameter</th>\n'...
	674 ' <th data-field="param">Percentile</th>\n'];
	675
	676 for i=1:Nwalkers
	677 html = [html, ...
	678 ' <th data-field="w1">Walker ', num2str(i), '</th>\n'…
	679 end
	680
	681 % epsilon_int
	682 html = [html, ...
	683 ' <th data-field="avg">Average</th>\n'...
	684 ' </tr>\n'...
	685 ' </thead>\n'...
	686 ' <tbody>\n'...
	687 ' <tr>\n'...
	688 ' <td> </td>\n'...
	689 ' <td align="center">25%%</td>\n'];
	690 for i=1:Nwalkers
	691 html = [html, ' <td>',...
	692 num2str(epsilon_int_25(i),3),'</td>\n'];
	693 end
	694
	695 html = [html, ' <td>',...
	696 num2str(sum(epsilon_int_25)/Nwalkers,3),'</td>\n'...
	697 ' </tr>\n'...
	698 ' <tr>\n'...
	699 ' <td align="center">ε<sub>int</sub> [m/Myr]</td>…
	700 ' <td align="center">50%%</td>\n'];
	701
	702 for i=1:Nwalkers
	703 html = [html, ' <td>',...
	704 num2str(epsilon_int_50(i),3),'</td>\n'];
	705 end
	706
	707
	708 html = [html, ' <td>',...
	709 num2str(sum(epsilon_int_50)/Nwalkers,3),'</td>\n'...
	710 ' </tr>\n'...
	711 ' <tr style="border-bottom:1px solid #D0D0D0">\n'...
	712 ' <td> </td>\n'...
	713 ' <td align="center">75%%</td>\n'];
	714
	715 for i=1:Nwalkers
	716 html = [html, ' <td>',...
	717 num2str(epsilon_int_75(i),3),'</td>\n'];
	718 end
	719
	720 html = [html, ' <td>',...
	721 num2str(sum(epsilon_int_75)/Nwalkers,3),'</td>\n'...
	722 ' </tr>\n'];
	723
	724
	725 % epsilon_gla
	726 html = [html, ...
	727 ' <tr>\n'...
	728 ' <td> </td>\n'...
	729 ' <td align="center">25%%</td>\n'];
	730 for i=1:Nwalkers
	731 html = [html, ' <td>',...
	732 num2str(epsilon_gla_25(i),3),'</td>\n'];
	733 end
	734
	735 html = [html, ' <td>',...
	736 num2str(sum(epsilon_gla_25)/Nwalkers,3),'</td>\n'...
	737 ' </tr>\n'...
	738 ' <tr>\n'...
	739 ' <td align="center">ε<sub>gla</sub> [m/Myr]</td>…
	740 ' <td align="center">50%%</td>\n'];
	741
	742 for i=1:Nwalkers
	743 html = [html, ' <td>',...
	744 num2str(epsilon_gla_50(i),3),'</td>\n'];
	745 end
	746
	747
	748 html = [html, ' <td>',...
	749 num2str(sum(epsilon_gla_50)/Nwalkers,3),'</td>\n'...
	750 ' </tr>\n'...
	751 ' <tr style="border-bottom:1px solid #D0D0D0">\n'...
	752 ' <td> </td>\n'...
	753 ' <td align="center">75%%</td>\n'];
	754
	755 for i=1:Nwalkers
	756 html = [html, ' <td>',...
	757 num2str(epsilon_gla_75(i),3),'</td>\n'];
	758 end
	759
	760 html = [html, ' <td>',...
	761 num2str(sum(epsilon_gla_75)/Nwalkers,3),'</td>\n'...
	762 ' </tr>\n'];
	763
	764
	765 % record_threshold
	766 html = [html, ...
	767 ' <tr>\n'...
	768 ' <td> </td>\n'...
	769 ' <td align="center">25%%</td>\n'];
	770 for i=1:Nwalkers
	771 html = [html, ' <td>',...
	772 num2str(record_threshold_25(i),3),'</td>\n'];
	773 end
	774
	775 html = [html, ' <td>',...
	776 num2str(sum(record_threshold_25)/Nwalkers,3),'</td>\n'...
	777 ' </tr>\n'...
	778 ' <tr>\n'...
	779 ' <td align="center">δ<sup>18</sup>O<sub>threshold<…
	780 ' <td align="center">50%%</td>\n'];
	781
	782 for i=1:Nwalkers
	783 html = [html, ' <td>',...
	784 num2str(record_threshold_50(i),3),'</td>\n'];
	785 end
	786
	787
	788 html = [html, ' <td>',...
	789 num2str(sum(record_threshold_50)/Nwalkers,3),'</td>\n'...
	790 ' </tr>\n'...
	791 ' <tr style="border-bottom:1px solid #D0D0D0">\n'...
	792 ' <td> </td>\n'...
	793 ' <td align="center">75%%</td>\n'];
	794
	795 for i=1:Nwalkers
	796 html = [html, ' <td>',...
	797 num2str(record_threshold_75(i),3),'</td>\n'];
	798 end
	799
	800 html = [html, ' <td>',...
	801 num2str(sum(record_threshold_75)/Nwalkers,3),'</td>\n'...
	802 ' </tr>\n'];
	803
	804 % E
	805 html = [html, ...
	806 ' <tr>\n'...
	807 ' <td> </td>\n'...
	808 ' <td align="center">25%%</td>\n'];
	809 for i=1:Nwalkers
	810 html = [html, ' <td>',...
	811 num2str(E_25(i),3),'</td>\n'];
	812 end
	813
	814 html = [html, ' <td>',...
	815 num2str(sum(E_25)/Nwalkers,3),'</td>\n'...
	816 ' </tr>\n'...
	817 ' <tr>\n'...
	818 ' <td align="center">E [m/Myr]</td>\n'...
	819 ' <td align="center">50%%</td>\n'];
	820
	821 for i=1:Nwalkers
	822 html = [html, ' <td>',...
	823 num2str(E_50(i),3),'</td>\n'];
	824 end
	825
	826
	827 html = [html, ' <td>',...
	828 num2str(sum(E_50)/Nwalkers,3),'</td>\n'...
	829 ' </tr>\n'...
	830 ' <tr style="border-bottom:1px solid #D0D0D0">\n'...
	831 ' <td> </td>\n'...
	832 ' <td align="center">75%%</td>\n'];
	833
	834 for i=1:Nwalkers
	835 html = [html, ' <td>',...
	836 num2str(E_75(i),3),'</td>\n'];
	837 end
	838
	839 html = [html, ' <td>',...
	840 num2str(sum(E_75)/Nwalkers,3),'</td>\n'...
	841 ' </tr>\n'];
	842
	843
	844 % footer
	845 html = [html, ' </tbody>\n'...
	846 '</table>\n'...
	847 ];
	848 fileID = fopen(filename,'w');
	849 fprintf(fileID, html);
	850 fclose(fileID);
	851
	852 %% generate csv table of results
	853 filename = [save_file, '.csv'];
	854 disp('Saving CSV table of results')
	855 % header
	856 csv = [...
	857 'Parameter;Percentile;'];
	858
	859 for i=1:Nwalkers
	860 csv = [csv, ...
	861 'Walker ', num2str(i), ';'];
	862 end
	863
	864 % epsilon_int
	865 csv = [csv, ...
	866 'Average\n'...
	867 ';25%%;'];
	868 for i=1:Nwalkers
	869 csv = [csv, num2str(epsilon_int_25(i),3),';'];
	870 end
	871
	872 csv = [csv, num2str(sum(epsilon_int_25)/Nwalkers,3),'\n'...
	873 'epsilon_int [m/Myr];50%%;'];
	874
	875 for i=1:Nwalkers
	876 csv = [csv, num2str(epsilon_int_50(i),3),';'];
	877 end
	878
	879
	880 csv = [csv, num2str(sum(epsilon_int_50)/Nwalkers,3),'\n'...
	881 ';75%%;'];
	882
	883 for i=1:Nwalkers
	884 csv = [csv, num2str(epsilon_int_75(i),3),';'];
	885 end
	886
	887 csv = [csv, num2str(sum(epsilon_int_75)/Nwalkers,3),'\n'];
	888
	889
	890 % epsilon_gla
	891 csv = [csv, ...
	892 ';25%%;'];
	893 for i=1:Nwalkers
	894 csv = [csv, num2str(epsilon_gla_25(i),3),';'];
	895 end
	896
	897 csv = [csv, num2str(sum(epsilon_gla_25)/Nwalkers,3),'\n'...
	898 'epsilon_gla [m/Myr];50%%;'];
	899
	900 for i=1:Nwalkers
	901 csv = [csv, num2str(epsilon_gla_50(i),3),';'];
	902 end
	903
	904
	905 csv = [csv, num2str(sum(epsilon_gla_50)/Nwalkers,3),'\n'...
	906 ';75%%;'];
	907
	908 for i=1:Nwalkers
	909 csv = [csv, num2str(epsilon_gla_75(i),3),';'];
	910 end
	911
	912 csv = [csv, num2str(sum(epsilon_gla_75)/Nwalkers,3),'\n'];
	913
	914 % record_threshold
	915 csv = [csv, ...
	916 ';25%%;'];
	917 for i=1:Nwalkers
	918 csv = [csv, num2str(record_threshold_25(i),3),';'];
	919 end
	920
	921 csv = [csv, num2str(sum(record_threshold_25)/Nwalkers,3),'\n'...
	922 'd18O_threshold [permille];50%%;'];
	923
	924 for i=1:Nwalkers
	925 csv = [csv, num2str(record_threshold_50(i),3),';'];
	926 end
	927
	928
	929 csv = [csv, num2str(sum(record_threshold_50)/Nwalkers,3),'\n'...
	930 ';75%%;'];
	931
	932 for i=1:Nwalkers
	933 csv = [csv, num2str(record_threshold_75(i),3),';'];
	934 end
	935
	936 csv = [csv, num2str(sum(record_threshold_75)/Nwalkers,3),'\n'];
	937
	938
	939 % E
	940 csv = [csv, ...
	941 ';25%%;'];
	942 for i=1:Nwalkers
	943 csv = [csv, num2str(E_25(i),3),';'];
	944 end
	945
	946 csv = [csv, num2str(sum(E_25)/Nwalkers,3),'\n'...
	947 'E [m/Myr];50%%;'];
	948
	949 for i=1:Nwalkers
	950 csv = [csv, num2str(E_50(i),3),';'];
	951 end
	952
	953
	954 csv = [csv, num2str(sum(E_50)/Nwalkers,3),'\n'...
	955 ';75%%;'];
	956
	957 for i=1:Nwalkers
	958 csv = [csv, num2str(E_75(i),3),';'];
	959 end
	960
	961 csv = [csv, num2str(sum(E_75)/Nwalkers,3),'\n'];
	962
	963
	964 fileID = fopen(filename,'w');
	965 fprintf(fileID, csv);
	966 fclose(fileID);
	967
	968
	969 %% generate html table of input parameters
	970 filename = [save_file, '-input.html'];
	971 disp('Saving html table of input values')
	972
	973 % header
	974 html = ['\n' ...
	975 '<table class="highlight">\n'...
	976 ' <thead>\n'...
	977 ' <tr>\n'...
	978 ' <th data-field="param">Parameter</th>\n'...
	979 ' <th data-field="val">Value</th>\n'...
	980 ' </tr>\n'...
	981 ' </thead>\n'...
	982 ' <tbody>\n'...
	983 ' <tr>\n'...
	984 ' <td>Sample ID</td>\n'...
	985 ' <td>' sample_id{1} '</td>\n'...
	986 ' </tr>\n'...
	987 ' <td>Name</td>\n'...
	988 ' <td>' name{1} '</td>\n'...
	989 ' </tr>\n'...
	990 ' <td>Email</td>\n'...
	991 ' <td>' email{1} '</td>\n'...
	992 ' </tr>\n'...
	993 ' <tr>\n'...
	994 ' <td>Latitude</td>\n'...
	995 ' <td>' lat{1} '</td>\n'...
	996 ' </tr>\n'...
	997 ' <tr>\n'...
	998 ' <td>Longitude</td>\n'...
	999 ' <td>' long{1} '</td>\n'...
	1000 ' </tr>\n'...
	1001 ' <tr>\n'...
	1002 ' <td><sup>10</sup>Be concentration</td>\n'...
	1003 ' <td>' num2str(be_conc/1000.) ' atoms/g</td>\n'...
	1004 ' </tr>\n'...
	1005 ' <tr>\n'...
	1006 ' <td><sup>26</sup>Al concentration</td>\n'...
	1007 ' <td>' num2str(al_conc/1000.) ' atoms/g</td>\n'...
	1008 ' </tr>\n'...
	1009 ' <tr>\n'...
	1010 ' <td><sup>14</sup>C concentration</td>\n'...
	1011 ' <td>' num2str(c_conc/1000.) ' atoms/g</td>\n'...
	1012 ' </tr>\n'...
	1013 ' <tr>\n'...
	1014 ' <td><sup>21</sup>Ne concentration</td>\n'...
	1015 ' <td>' num2str(ne_conc/1000.) ' atoms/g</td>\n'...
	1016 ' </tr>\n'...
	1017 ' <tr>\n'...
	1018 ' <td><sup>10</sup>Be conc. uncertainty</td>\n'...
	1019 ' <td>' num2str(be_uncer*100.) ' %%</td>\n'...
	1020 ' </tr>\n'...
	1021 ' <tr>\n'...
	1022 ' <td><sup>26</sup>Al conc. uncertainty</td>\n'...
	1023 ' <td>' num2str(al_uncer*100.) ' %%</td>\n'...
	1024 ' </tr>\n'...
	1025 ' <tr>\n'...
	1026 ' <td><sup>14</sup>C conc. uncertainty</td>\n'...
	1027 ' <td>' num2str(c_uncer*100.) ' %%</td>\n'...
	1028 ' </tr>\n'...
	1029 ' <tr>\n'...
	1030 ' <td><sup>21</sup>Ne conc. uncertainty</td>\n'...
	1031 ' <td>' num2str(ne_uncer*100.) ' %%</td>\n'...
	1032 ' </tr>\n'...
	1033 ' <tr>\n'...
	1034 ' <td>Observation depth</td>\n'...
	1035 ' <td>' num2str(zobs) ' m</td>\n'...
	1036 ' </tr>\n'...
	1037 ' <tr>\n'...
	1038 ' <td><sup>10</sup>Be production (spallation)</td>\n'...
	1039 ' <td>' num2str(be_prod_spall/1000.) ' atoms/g/yr</td>\n'…
	1040 ' </tr>\n'...
	1041 ' <tr>\n'...
	1042 ' <td><sup>26</sup>Al production (spallation)</td>\n'...
	1043 ' <td>' num2str(al_prod_spall/1000.) ' atoms/g/yr</td>\n'…
	1044 ' </tr>\n'...
	1045 ' <tr>\n'...
	1046 ' <td><sup>14</sup>C production (spallation)</td>\n'...
	1047 ' <td>' num2str(c_prod_spall/1000.) ' atoms/g/yr</td>\n'.…
	1048 ' </tr>\n'...
	1049 ' <tr>\n'...
	1050 ' <td><sup>21</sup>Ne production (spallation)</td>\n'...
	1051 ' <td>' num2str(ne_prod_spall/1000.) ' atoms/g/yr</td>\n'…
	1052 ' </tr>\n'...
	1053 ' <tr>\n'...
	1054 ' <td><sup>10</sup>Be production (muons)</td>\n'...
	1055 ' <td>' num2str(be_prod_muons/1000.) ' atoms/g/yr</td>\n'…
	1056 ' </tr>\n'...
	1057 ' <tr>\n'...
	1058 ' <td><sup>26</sup>Al production (muons)</td>\n'...
	1059 ' <td>' num2str(al_prod_muons/1000.) ' atoms/g/yr</td>\n'…
	1060 ' </tr>\n'...
	1061 ' <tr>\n'...
	1062 ' <td><sup>14</sup>C production (muons)</td>\n'...
	1063 ' <td>' num2str(c_prod_muons/1000.) ' atoms/g/yr</td>\n'.…
	1064 ' </tr>\n'...
	1065 ' <tr>\n'...
	1066 ' <td><sup>21</sup>Ne production (muons)</td>\n'...
	1067 ' <td>' num2str(ne_prod_muons/1000.) ' atoms/g/yr</td>\n'…
	1068 ' </tr>\n'...
	1069 ' <tr>\n'...
	1070 ' <td>Rock density</td>\n'...
	1071 ' <td>' num2str(rock_density) ' kg/m<sup>3</sup></td>\n'.…
	1072 ' </tr>\n'...
	1073 ' <tr>\n'...
	1074 ' <td>ε<sub>gla</sub></td>\n'...
	1075 ' <td>' num2str(epsilon_gla_min*1.0e6) ...
	1076 ' to ' num2str(epsilon_gla_max*1.0e6) ' m/Myr</td>\n'...
	1077 ' </tr>\n'...
	1078 ' <tr>\n'...
	1079 ' <td>ε<sub>int</sub></td>\n'...
	1080 ' <td>' num2str(epsilon_int_min*1.0e6) ...
	1081 ' to ' num2str(epsilon_int_max*1.0e6) ' m/Myr</td>\n'...
	1082 ' </tr>\n'...
	1083 ' <tr>\n'...
	1084 ' <td><i>t</i><sub>degla</sub></td>\n'...
	1085 ' <td>' num2str(t_degla_min) ...
	1086 ' to ' num2str(t_degla_max) ' yr</td>\n'...
	1087 ' </tr>\n'...
	1088 ' <tr>\n'...
	1089 ' <td>Climate record</td>\n'...
	1090 ' <td>' record{1} '</td>\n'...
	1091 ' </tr>\n'...
	1092 ' <tr>\n'...
	1093 ' <td>δ<sup>18</sup>O<sub>threshold</sub></td>\n'...
	1094 ' <td>' num2str(record_threshold_min) ...
	1095 ' to ' num2str(record_threshold_max) ' &permil;</td>\n'...
	1096 ' </tr>\n'...
	1097 ' <tr>\n'...
	1098 ' <td>MCMC walkers</td>\n'...
	1099 ' <td>' num2str(nwalkers) '</td>\n'...
	1100 ' </tr>\n'...
	1101 ' </tbody>\n'...
	1102 '</table>\n'];
	1103 fileID = fopen(filename,'w');
	1104 fprintf(fileID, html);
	1105 fclose(fileID);
	1106
	1107 %% Save general data
	1108 disp('Saving general data files');
	1109 dlmwrite([save_file, '-eps_int.txt'], epsilon_int_data);
	1110 dlmwrite([save_file, '-eps_gla.txt'], epsilon_gla_data);
	1111 dlmwrite([save_file, '-t_degla.txt'], t_degla_data);
	1112 dlmwrite([save_file, '-d18O_threshold.txt'], record_threshold_data);
	1113
	1114 %% Save per-walker data
	1115 for i=1:Nwalkers
	1116 dlmwrite([save_file, '-eps_int-w' num2str(i) '.txt'],...
	1117 epsilon_int_data_w{i});
	1118 dlmwrite([save_file, '-eps_gla-w' num2str(i) '.txt'],...
	1119 epsilon_gla_data_w{i});
	1120 dlmwrite([save_file, '-t_degla-w' num2str(i) '.txt'],...
	1121 t_degla_data_w{i});
	1122 dlmwrite([save_file, '-d18O_threshold-w' num2str(i) '.txt'],...
	1123 record_threshold_data_w{i});
	1124 end
	1125
	1126 % create HTML snippet with results
	1127 disp('Creating html snippet for per-walker data file links')
	1128 idstring = strsplit(save_file, '/');
	1129 id = idstring(end);
	1130
	1131 html = '\n';
	1132 for i=1:Nwalkers
	1133 html = [html, ...
	1134 ' <br><a href="output/', id{1}, '-eps_int-w', .…
	1135 num2str(i), '.txt"\n', ...
	1136 ' target="_blank">Walker ', num2str(i), ...
	1137 ' ε<sub>int</sub> ',...
	1138 'data</a>\n'];
	1139 html = [html, ...
	1140 ' <a href="output/', id{1}, '-eps_gla-w', ...
	1141 num2str(i), '.txt"\n', ...
	1142 ' target="_blank">Walker ', num2str(i), ...
	1143 ' ε<sub>gla</sub> ',...
	1144 'data</a>\n'];
	1145 html = [html, ...
	1146 ' <a href="output/', id{1}, '-t_degla-w', ...
	1147 num2str(i), '.txt"\n', ...
	1148 ' target="_blank">Walker ', num2str(i), ...
	1149 ' <i>t</i><sub>degla</sub> ',...
	1150 'data</a>\n'];
	1151 html = [html, ...
	1152 ' <a href="output/', id{1}, '-eps_int-w', ...
	1153 num2str(i), '.txt"\n', ...
	1154 ' target="_blank">Walker ', num2str(i), ...
	1155 ' δ<sup>18</sup>O', ...
	1156 '<sub>threshold</sub> data</a>\n'];
	1157 end
	1158 fileID = fopen([save_file, '-walker-data.html'],'w');
	1159 fprintf(fileID, html);
	1160 fclose(fileID);