tgenerate_plots.m - cosmo - front and backend for Markov-Chain Monte Carlo inversion of cosmogenic nuclide concentrations
(HTM) git clone git://src.adamsgaard.dk/cosmo
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) 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)/M);
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,:),fs.m_true(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),' of ',fsamp.ProposerType],...
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.TargetAcceptanceRatio)];
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),' of ',fsamp.ProposerType],...
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.TargetAcceptanceRatio)];
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 coordinate
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 we do not sample outside
144 xbins{i1} = linspace(mminmax_i1(1)-dxbin,mminmax_i1(2)+2*dxbin,Nbin+4); %one extra to make pcolor work
145 case 'logunif'
146 dfacxbin = (mminmax_i1(2)/mminmax_i1(1))^(1/(Nbin)); %we want outside bins to check that we do not sample outside
147 xbins{i1} = logspace(log10(mminmax_i1(1)/dfacxbin),log10(mminmax_i1(2)*dfacxbin^2),Nbin+4);
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 that we do not sample outside
158 ybins = linspace(mminmax_i2(1)-dybin,mminmax_i2(2)+2*dybin,Nbin+4);
159 case 'logunif'
160 dfacybin = (mminmax_i2(2)/mminmax_i2(1))^(1/Nbin); %we want outside bins to check that we do not sample outside
161 ybins = logspace(log10(mminmax_i2(1)/dfacybin),log10(mminmax_i2(2)*dfacybin^2),Nbin+4);
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,:)*1.0e6)
305 case 'logunif', set(gca,'xscale','log','xlim',mminmax(i1,:)*1.0e6)
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 walkers
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,:)*1.0e6)
416 case 'logunif', set(gca,'xscale','log','xlim',mminmax(i1,:)*1.0e6)
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_triinterp_2p6Ma
435 elseif strcmp(record, 'rec_20kyr')
436 d18O_filename = 'lisiecki_triinterp_2p6Ma_20ky.mat'; % zachos_triinterp_2p6Ma
437 elseif strcmp(record, 'rec_30kyr')
438 d18O_filename = 'lisiecki_triinterp_2p6Ma_30ky.mat'; % zachos_triinterp_2p6Ma
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.5);
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.2*dy, ylims(2) + 0.2*dy])
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 = 20,
555 % % % % iz = 1;
556 % % % % DepthExposureTimeDens(zsss,tss,dtfine,tfinemax,dzfine,zfinemax,fixed_stuff,iz)
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 concentration or ...
566 dXxfine = 0.01; Xxfinemax = 50; iz=1; %Xx may be depth or nucleide concentration or ...
567 dzfine = dXxfine; zfinemax = Xxfinemax;
568 [smoothgrid,histgrid,tsfine,Xxfine]=XxTimeDens(Xxsss,tss,dtfine,tfinemax,dXxfine,Xxfinemax,fixed_stuff,iz);
569 %pcolor(-tsfine,Xxfine,sqrt(smoothgrid));
570 %pcolor(tsfine,Xxfine,sqrt(histgrid)); %<<<BHJ: v?lge mellem glattet og r? farvel?gning
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,zsfine);
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,zsfine);
587 quants2(:,:,iwalk) = GetHistgridQuantiles2(histgrid,N,fractions,tsfine,zsfine);
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: added 2014 dec 09
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>\n'...
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>\n'...
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</sub> [‰]</td>\n'...
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) ' ‰</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);