tadd new plots of exhumation history - 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
---
(DIR) commit 7406efb705d7485353262c46115774eb3da73cf9
(DIR) parent f4b59b8b4e63bb3a55c7abbf6d770d6ca90c00b6
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Tue, 10 Nov 2015 16:33:40 +0100
add new plots of exhumation history
Diffstat:
M foot.html | 2 +-
M matlab/generate_plots.m | 169 +++++++++++++++++++++++++++++++
2 files changed, 170 insertions(+), 1 deletion(-)
---
(DIR) diff --git a/foot.html b/foot.html
t@@ -48,7 +48,7 @@
href="mailto:anders.damsgaard@geo.au.dk">Anders
Damsgaard</a> and designed using
<a class="white-text text-lighten-3"
- href="http://materializecss.com">Materialize</a>
+ href="http://materializecss.com">Materialize</a>.
</div>
</div>
</footer>
(DIR) diff --git a/matlab/generate_plots.m b/matlab/generate_plots.m
t@@ -163,6 +163,7 @@ for i1=1:M
end
%% Histogram plots for all four parameters, one subplot per walker, fh(4)
+disp('Histogram plots for all four parameters, one subplot per walker')
fh = [fh;figure('visible', show_figures)];
for i1 = 1:M
for iwalk=1:min(4,Nwalkers)
t@@ -199,6 +200,7 @@ if titles
end
%% Plot one histogram per model parameter per walker
+disp('One histogram per model parameter per walker')
fh = [fh;figure('visible', show_figures)];
for i1 = 1:M % for each model parameter
for iwalk=1:min(4,Nwalkers) % for each walker
t@@ -251,6 +253,7 @@ end
%% Plot one histogram per model parameter, including data from all walkers
+disp('One histogram per model parameter')
fh = [fh;figure('visible', show_figures)];
for i1 = 1:M % for each model parameter
t@@ -335,6 +338,7 @@ end
%% Plot d18O curve, from PlotSmoothZachos.m
+disp('d18O curve, exposure, and erosion history')
fh = [fh;figure('visible', show_figures)];
if strcmp(record, 'rec_5kyr')
t@@ -440,6 +444,170 @@ linkaxes(axh,'x')
% stairs(-tStarts,relExpos+2,'r')
+%% Exhumation history from InspectDepthConcTracks_medTrueDepthIndsat.m
+disp('Exhumation history')
+%[fn,pn]=uigetfile('*.mat')
+%load([pn,fn])
+DoConfLevelTrim = 1;
+
+fh = [fh;figure('visible', show_figures)];
+
+Nwalks = length(Ss);
+for iwalk = 1:Nwalks
+ % iwalk=input(['What iwalk?[1..',num2str(length(Ss)),']']),
+ subplot(Nwalks,1,iwalk)
+ lump_MetHas = Ss{iwalk}.lump_MetHas;
+ fixed_stuff = Ss{iwalk}.fs;
+ % if ~isempty(lump_MetHas.zsss{1})
+ % % % % zsss = lump_MetHas.zsss;
+ % % % % tss = lump_MetHas.tss;
+ % % % % dtfine = -1000; tfinemax = -20e5; dzfine = 0.1; zfinemax = 20,
+ % % % % iz = 1;
+ % % % % DepthExposureTimeDens(zsss,tss,dtfine,tfinemax,dzfine,zfinemax,fixed_stuff,iz)
+ % % % % title('z(time)')
+ % % % % axis ij
+ % % % % hold on
+
+ %making quantiles to plot ontop of exhumation densities:
+ tss = lump_MetHas.tss;
+ dtfine = -1000; tfinemax = -20e5;
+ Xxsss = lump_MetHas.zsss;
+ dXxfine = 0.01; Xxfinemax = 10; iz=1; %Xx may be depth or nucleide concentration or ...
+ dzfine = dXxfine; zfinemax = Xxfinemax;
+ [smoothgrid,histgrid,tsfine,Xxfine]=XxTimeDens(Xxsss,tss,dtfine,tfinemax,dXxfine,Xxfinemax,fixed_stuff,iz);
+ pcolor(tsfine,Xxfine,sqrt(smoothgrid));
+ hold on
+ %plot a selection of depth histories
+ for i=1:ceil(length(tss)/50):length(tss)
+ plot(tss{i},Xxsss{i},'.-w')
+ end
+ %Compute and plot quantiles
+ N = sum(histgrid(:,1));
+ fractions = [0.25,0.5,0.75]; %quartiles
+ tsfine = 0:dtfine:tfinemax; %bin boundaries
+ Ntfine = length(tsfine);
+ zsfine = 0:dzfine:zfinemax; %bin boundaries
+ % quants{iwalk} = GetHistgridQuantiles(histgrid,N,fractions,tsfine,zsfine);
+ % plot(tsfine,quants{iwalk}(1,:),'r','linewidth',3)
+ % plot(tsfine,quants{iwalk}(2,:),'k','linewidth',3)
+ % plot(tsfine,quants{iwalk}(3,:),'b','linewidth',3)
+ quants(:,:,iwalk) = GetHistgridQuantiles(histgrid,N,fractions,tsfine,zsfine);
+
+ grid on; shading flat; axis tight; set(gca,'fontsize',8); hold on
+ lh(1)=plot(tsfine,quants(1,:,iwalk),'r','linewidth',2);
+ lh(2)=plot(tsfine,quants(2,:,iwalk),'k','linewidth',2);
+ lh(3)=plot(tsfine,quants(3,:,iwalk),'g','linewidth',2);
+ legend(lh,'25%','median','75%','location','northwest')
+ axis ij
+ %track_handle=AddTrueModelDepthTrack(Ss{iwalk}.fs,'*-m'); %<<< BHJ: added 2014 dec 09
+end
+colormap(1-copper)
+% subplot(2,2,1);
+
+% the figure below is really computationally intensive
+% fh = [fh;figure('visible', show_figures)];
+% for iwalk = 1:Nwalks
+% plot(tsfine,quants(1,:,iwalk),'r','linewidth',3); hold on
+% plot(tsfine,quants(2,:,iwalk),'k','linewidth',3)
+% plot(tsfine,quants(3,:,iwalk),'g','linewidth',3)
+% end
+% %axis([-2e6 0 0 4])
+% axis ij
+
+% load handel; sound(y(1:5000))
+% iwalk=input(['What iwalk?[1..',num2str(length(Ss)),']']);
+% lump_MetHas = Ss{iwalk}.lump_MetHas;
+% fixed_stuff = Ss{iwalk}.fs;
+% % if ~isempty(lump_MetHas.zsss{1})
+% fh = [fh;figure('visible', show_figures)];
+% zsss = lump_MetHas.zsss;
+% tss = lump_MetHas.tss;
+% dtfine = -1000; tfinemax = -20e5; dzfine = 0.1; zfinemax = 20,
+% iz = 1;
+% DepthExposureTimeDens(zsss,tss,dtfine,tfinemax,dzfine,zfinemax,fixed_stuff,iz)
+% title('z(time)')
+% axis ij
+% hold on
+%
+% %making quantiles to plot ontop of exhumation densities:
+% Xxsss = lump_MetHas.zsss;
+% dXxfine = 0.1; Xxfinemax = 10;
+% [smoothgrid,histgrid,tsfine,Xxfine]=XxTimeDens(Xxsss,tss,dtfine,tfinemax,dXxfine,Xxfinemax,fixed_stuff,iz);
+% %plot a selection of depth histories
+% for i=1:ceil(length(tss)/50):length(tss)
+% plot(tss{i},zsss{i},'.-w')
+% end
+%
+% %Compute and plot quantiles
+% N = sum(histgrid(:,1));
+% fractions = [0.25,0.5,0.75]; %quartiles
+% tsfine = 0:dtfine:tfinemax; %bin boundaries
+% Ntfine = length(tsfine);
+% zsfine = 0:dzfine:zfinemax; %bin boundaries
+% quants = GetHistgridQuantiles(histgrid,N,fractions,tsfine,zsfine);
+% plot(tsfine,quants(1,:),'r','linewidth',3)
+% plot(tsfine,quants(2,:),'k','linewidth',3)
+% plot(tsfine,quants(3,:),'g','linewidth',3)
+
+% fh = [fh;figure('visible', show_figures)];
+% zsss = lump_MetHas.zsss;
+% tss = lump_MetHas.ExposureTimeSinceNows;
+% dtExposurefine = 1000; tExposurefinemax = 20e5; dzfine = 0.1; zfinemax = 20,
+% iz = 1;
+% DepthExposureTimeDens(zsss,tss,dtExposurefine,tExposurefinemax,dzfine,zfinemax,fixed_stuff,iz)
+% title('z(ExposureTime)')
+% % end
+% axis ij
+% hold on
+% for i=1:ceil(length(tss)/50):length(tss)
+% plot(tss{i},zsss{i},'.-w')
+% end
+%
+% figure %five subplots:
+% tss = lump_MetHas.tss;
+% dtfine = -1000; tfinemax = -20e5;
+% iz = 1;
+%
+% axh(1)=subplot(5,1,1) %depth track
+% Xxsss = lump_MetHas.zsss;
+% dXxfine = 0.1; Xxfinemax = 10;
+% [smoothgrid,histgrid,tsfine,Xxfine]=XxTimeDens(Xxsss,tss,dtfine,tfinemax,dXxfine,Xxfinemax,fixed_stuff,iz);
+% pcolor(tsfine,Xxfine,sqrt(smoothgrid));
+% % axis(sqrt(3)*[-1,1,-1,1])
+% % set(gca,'xtick',[-1.5:0.5:1.5],'ytick',[-1.5:0.5:1.5])
+% grid on; shading flat; axis tight; set(gca,'fontsize',8); hold on
+% lh(1)=plot(tsfine,quants(1,:),'r','linewidth',2)
+% lh(2)=plot(tsfine,quants(2,:),'k','linewidth',2)
+% lh(3)=plot(tsfine,quants(3,:),'g','linewidth',2)
+% legend(lh,'25%','median','75%','location','northwest')
+% axis ij
+%
+% Xxsss{2}=lump_MetHas.c14Csss;
+% Xxsss{3}=lump_MetHas.c10Besss;
+% Xxsss{4}=lump_MetHas.c26Alsss;
+% Xxsss{5}=lump_MetHas.c21Nesss;
+%
+% for ic=2:5
+% axh(ic)=subplot(5,1,ic) %10Be track
+% dXxfine = 0.01*Xxsss{ic}{end}(1,end); Xxfinemax = 1.3*Xxsss{ic}{end}(1,end);
+% [smoothgrid,histgrid,tsfine,Xxfine]=XxTimeDens(Xxsss{ic},tss,dtfine,tfinemax,dXxfine,Xxfinemax,fixed_stuff,iz);
+% pcolor(tsfine,Xxfine,sqrt(smoothgrid));
+% % axis(sqrt(3)*[-1,1,-1,1])
+% % set(gca,'xtick',[-1.5:0.5:1.5],'ytick',[-1.5:0.5:1.5])
+% grid on; shading flat; axis tight; set(gca,'fontsize',8,'tickdir','out');
+% hold on
+% end
+% for ic=1:4
+% set(axh(ic),'xticklabel',[])
+% end
+% subplot(5,1,1); title(['Depth track:',fn],'interpreter','none'); set(gca,'layer','top')
+% subplot(5,1,2); title('^{14}C'); set(gca,'layer','top')
+% subplot(5,1,3); title('^{10}Be'); set(gca,'layer','top')
+% subplot(5,1,4); title('^{26}Al'); set(gca,'layer','top')
+% subplot(5,1,5); title('^{21}Ne'); set(gca,'layer','top')
+% linkaxes(axh,'x')
+% colormap(flipud(gray))
+
%% position figure windows at certain coordinates on the screen
if strcmp(show_figures, 'on')
t@@ -455,6 +623,7 @@ if strcmp(show_figures, 'on')
end
%% save all figures
+disp('Saving figures')
for i=1:length(fh)
figure_save_multiformat(fh(i), ...
strcat(save_file, '-', num2str(i)), ...