Age versus length of horizontal fission tracks represented by 2D-Histograms

The code calculates the age distribution of horizontal fission tracks versus their lengths. The code is for tracks measured in the old way where the angles to the c-axis were not considered.  Equations relating 2D track length histograms and associated surface track densities to track ages are given in an early version in Jensen et al. (1992) and illustrated in a poster (Jensen and Hansen, 2018). The present code is based on further development (Jensen, 2021) and inversion is now performed using a probabilistic theory (Tarantola, 2005). The code is developed for the examples given in Jensen and Hansen (2021). A code for tracks measurements including the angle to the c-axis is under preparation for upload.

Main programme, running in Octave:

%
%   **********  LsqTrack - angle 2D  ********************** 17.02.2021
%
%     Peter Klint Jensen, Technical University of Denmark, 
%     Dept. Civil Engineering, Geotechnichs and Geology.
%
% This programme calculates fission track ages as a function of their lengths 
% based on 2D-track length histograms.
%
% The equations relating track length to age is given in 
% Jensen, P.K., Hansen, K. (2021) Age distribution of horizontal fission tracks.
% submitted to Geochronology, feb. 2021. The probailistic least least squares 
% technique applied is described in Tarantola (2005).
%
% An early set of equations and invertion by the use of mathematical simulated
% annealing is described in:  
%
% Jensen, P.K., Hansen, K., and Kunzendorf, H. (1992). A numerical model for 
% the thermal history of rocks based on confined horizontal fission tracks. 
% Nucl. Tracks. Radiat. Meas., Vol. 20, No. 2, 349-359.	
%
% Jensen, P. K. and Hansen, K.: Identifying the post-sedimentary part of 
% fission track length histograms with inherited tracks, 
% Thermo 2018, 16th International Conference on Thermochronology, 
% Quedlinburg, Germany; 16-21 September 2018, 
% https://doi.org/10.1002/essoar.10500031.1, 2018. 
%   
% The computer code is written in GNU Octave, version 6.1.0.
% Please report bugs to: peterklintjensen@gmail.com.
%
%  Preparations of files:
%   
%   All scripts should be in the same folder:
%   
%   Scripts:
%
%   LsqTrack.M: This main programme.
%   Functions:
%	  Agedec.M : Surface track density of post-depositional tracks.
%   MeanL.M : Mean length.
%	  MisFit.M: Difference between measured and calculated histograms.
%
%   Input data:
%
%   InSampleName.txt: General Input data.
%   PriorSampleName.txt: Prior histogram.
%   PriorUnan.txt: Prior of inverted track lengths of an unannealed track length
%   distribution.
%   LSampleName.txt: Measured track lengths.
%   LGreenUnan.txt: Unannealed track length distribution.
%   zir1.txt: Interpolated filters, output from AgeVar.m
%   Var.txt: Variance of interpolated filters, output from AgeVar.m
%   xi.txt and yi.txt: Center of mesh, output from AgeVar.m
%
% Launch 'Octave' and load the packages 'optim': Example:
% > pkg load statistics optim
%
%

clear all ; clf
pkg load optim

global fac1 fac2 fac3 Binmtil

%   ***** Choose sample input data *****

%   East Greenland, Jameson Land, GGU103113.
load('InGGU103113.txt') ;
InheritIn = InGGU103113 ;
load('LGGU103113.txt') ; % Individual track lengths.
Lengths = LGGU103113 ;
load('PriorGGU103113.txt') ; % Prior information.
mpriorU = PriorGGU103113 ;


%  **********************************************************************

load('zir1.txt') ; % Interpolated filter generated by FilterSmooth.m 
load('Var.txt') ;  % Variance of Interpolated filter, zir1.  
load('xi.txt') ;   % Gridpoints corresponding to zir1.     
load('yi.txt') ;   % 

disp('                          ')

yesAngle = InheritIn(1) ; % 1 for 
useprior = InheritIn(2) ; % 1 for user values for prior.
mPl = InheritIn(3) ;      % Lower limit for model.
mPu = InheritIn(4) ;      % Upper.
unitm = InheritIn(5) ;    % Unit of data model bin. 
facVdobs = InheritIn(6) ; % Factor increasing data variance.                       
facVmp = InheritIn(7) ;   % Factor increasing prior variance.                         
Smoothd = InheritIn(8) ;  % Smooth input histogram. 
SmoothCD = InheritIn(9) ; % Increase data correlation distance.
AgeSed = InheritIn(10) ;  % Deposition age /Ma.
RHOs = InheritIn(11) ;    % Surface density.
RHOsN = InheritIn(12) ;   % Number counts of RHOs.
ppm238 = InheritIn(13) ;  % Uranium U238 koncentration, ppm .
ppmD = InheritIn(14)  ;   % Standard diviation, ppm238.

%   ---- Input data usually not to be changed ------

unitd = InheritIn(15) ;     % Unit of data bin.
L0 = InheritIn(16) ;        % Mean track length of induced tracks.
L0D = InheritIn(17) ;       % Std. dev. L0
ksipp = InheritIn(18) ;     % ksi prime, calibration constant.
ksippD = InheritIn(19) ;    % Standard div. ksipp
bpar = InheritIn(20) ;      % Calibration constant L/L0(rho/rho0).
bparD = InheritIn(21) ;     % Std. bpar.

MxBin = floor (24/unitd) ; % Max number of bins.

LambF = 8.46E-17 ; % Fission decay constant, 1/yr, +- 0.06 
LambFD = 0.06E-17 ; % Total decay constant
LambD = 1.55125e-10 ;% Total decay constant, 1/yr, approx. = alfa-decay.
RHOMIN = 3.2E3 ; % Density of glasstandards (kg/m**3)
UMASS = 3.952E-25 ; % Mass of uranium atom (kg)

%  Calculation of uranium concentration and std.
c238 = ppm238*1.0E-6*RHOMIN/UMASS ;
c238D = ppmD*1.0E-6*RHOMIN/UMASS ;

disp('LsqTrack 18.01.2021')
disp('Peter Klint Jensen')
disp('   ')
printf('Uranium concentration: %2.1e +/- %2.1e nr./m^3  \n', c238,c238D );

%     Data: histogram ; 

Bin = [1:MxBin]' ; % Track length intervals.
Hmeas = hist( Lengths,(Bin-0.5)*unitd,'facecolor','gr' )' ;
[MaxH, iMaxH] = max(Hmeas) ;
dobs = Hmeas ; length (Hmeas) ; % Observed data.

disp('   ')
disp(' *********  Measured histogram data **********')

NHmeas = sum(Hmeas) ; 
printf('Number Measured hist.: %2.1i   \n', NHmeas );

LHmeas = MeanL(Hmeas,Bin,unitd) ;
printf('Mean length measured hist: %2.1f micro. m \n', LHmeas );
printf('Measured surface density measured hist: %3.2f 10^6 cm^-2 \n', RHOs/1e6 );

figure(1); % Measured histogram.   
	bar( (Bin-0.5)*unitd,Hmeas,1.0, 'facecolor','gr');
   title ('Measured histogram','fontsize',14)
	 xlabel ('Track length / \mum','fontweight', 'bold','fontsize',14);
	 ylabel ('Number','fontweight', 'bold','fontsize',14);
   set(gca (),'fontsize',14,'tickdir','out'); % sets font of numbers on axes 		
	 axis ([0 20 0 1.1*MaxH]); % min and max limits
   text (1,0.9*MaxH,['Tracks ',num2str(NHmeas);num2str( (RHOs)/1000000,2), ...
   'cm^-^2';num2str( LHmeas,3 ), ' \mum'  ],'fontweight', 'bold','fontsize',14)	
 	print figure(1) -mono fig1Meas.pdf ; % jpg.,.pdf,.svg.,.eps.
 
Kernel = [1 4 6 4 1] / 16 ; % Smoothing of measured histogram. 
if ( Smoothd == 1 )
Hmeas = conv (Hmeas,Kernel,"same") ; 
else
endif 
 
figure(2); % Smoothed histogram.   
	bar( (Bin-0.5)*unitd,Hmeas,1.0, 'facecolor','gr');
   title ('Smoothed histogram','fontsize',14)
	 xlabel ('Track length / \mum','fontweight', 'bold','fontsize',14);
	 ylabel ('Number','fontweight', 'bold','fontsize',14);
   set(gca (),'fontsize',10,'tickdir','out'); % sets font of numbers on axes 		
	 axis ([0 20 0 1.1*MaxH]); % min and max limits

 
% ---- Extracting variance of filter, examples -----------------

  y1 = [0:0.2:25] ; % micr. m. as in filterVar.m.                   
figure(3) ; 
  plot( y1'+0.1 , zir1 (:,31 ) ) ; hold on
  plot( y1'+0.1 , zir1 (:,51 ) ) ; hold on
  plot( y1'+0.1 , zir1 (:,76 ) ) ;
   title ('Selected interpolated filters','fontsize',14)  
   legend ('6  \mum','10  \mum','15 \mu m')
   xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',14);
   ylabel ('Frequency','fontweight', 'bold','fontsize',14);   
hold off


%   ----- Deconvolution of measured histogram using  ---------------
%   ----- Least squares probabilistic inversion ---- 

% Covariance of data (measured filters) using multinomial variance.
% Variance is diagonal elements of covariance matrix.

eta(1:length (Hmeas)) = 1.*10^-3 ; % To avoid inf in inverse matrix.
Kernel = [1 4 6 4 1] / 16 ; % Smoothing measured histogram. 
if ( SmoothCD == 1 )
Hblur = conv (Hmeas,Kernel,"same") ;
else
Hblur = Hmeas ;
endif 

% Diagonal elements. 

prb = Hblur /sum ( Hblur ) ; Ey = eye( length(dobs) ) ; 
Diag = facVdobs*sum ( Hblur ) * prb .* ( 1- prb ) + eta' ;
CDdiag = eye ( length(dobs) ) .* Diag ;

% Off diagonal elements.
ons0 = ones ( length(dobs) ) - eye ( length(dobs) ) ; % ones with 0 diagonal.
AA = ons0 .* prb ; CDoff = -facVdobs*sum ( Hblur )*( AA .* AA' ) ; 
       
% Data covariance matrix.
CD = CDoff + CDdiag ; % CD is recalculated below to include modelization var.
  
figure(4) 
 hold on;
 imagesc (CD) ; 
  title ('Covariance of data ','fontsize',14) ; colorbar
  xlabel ('Columns','fontweight', 'bold','fontsize',14);
  ylabel ('Rows','fontweight', 'bold','fontsize',14); 
 hold off
  
% ------------  Prior model  -------------------------------

NBinm = floor( (mPu - mPl)/unitm ) ; % Number of prior bins. 
if ( useprior == 1 )
mprior = mpriorU ; % User prior model.
else
mprior(1:NBinm ) = NHmeas/NBinm ; % Default prior model, equal distribution.
endif

% Diagonal elements. 
etaP(1:length (mprior)) = 1.*10^-3 ; % To avoid inf in inverse matrix.
prbP = mprior /sum ( mprior ) ; ; EyP = eye( length(mprior) ) ;
DiagP = facVmp*sum ( mprior )* prbP .* ( 1- prbP )+ etaP' ;
CMdiag = EyP .* DiagP ;

% Off diagonal elements.
ons0P = ones ( length(mprior) ) - EyP ; % ones with 0 diagonal.
AAP = ons0P .* prbP ; 
CMoff = -facVmp*sum ( mprior )*( AAP .* AAP' ) ;
% Model covariance.
CM = CMoff + CMdiag ; % J = imsmooth(CM, "Gaussian", 0.5)


% **************  Modelization matrix G  *****************************

% Near middle of 1 micr. m intervals selected.

Gdum = zir1( 1 :(1/0.2): length(dobs)/0.2 , ...
(( (mPl+1) /0.2)) : (1/0.2) : (mPu/0.2) ) ; 
Gsum = sum ( Gdum(:,:),1 ) ;
G = ( Gdum' ./ Gsum' )' ; sum ( G(:,:),1 ) ; % Normalization.

%  Center of posterior Gaussian, mtil model transposed.

mtil = mprior' + ( G'*CD^-1*G + CM^-1 )^-1 * G'*CD^-1*( dobs - G*mprior' );
dt = G*mtil ; % mtil posterior.
%eps = 1.5*10^-1 ; mtil = ( ( G'*G + eps^2 * eye(,) )^-1 ) * G'*dobs ; %Tikonov
sum(dobs) ; sum(mtil) ; % Check sums are equal.


%-------------------------- Modelization variance --------------------
%                 That is the error caused by the model itself.

GVar = Var( 1 :(1/0.2): length(dobs)/0.2 , ...
(( (mPl+1) /0.2)) : (1/0.2) : (mPu/0.2) ) ;

dGVar = GVar*mtil ; % Variance of prediction due to modelization variance.
[dobs sqrt(abs(dGVar))] ;

figure(5) % 
 plot ( sqrt(abs(dGVar)) ) ;
  title ('Forward modelization std. deviation','fontsize',14) 
  xlabel ('Mean track length /  \mum','fontweight', 'bold','fontsize',14);
  ylabel ('Frequency','fontweight', 'bold','fontsize',14);    
hold off


%  Recalculating posterior model including modelization uncertainty. 

for ij = 1:length (Hmeas) ;  
  CD(ij,ij) =   CD(ij,ij) + dGVar(ij) ; 
endfor 
 
figure(6) 
 hold on;
 imagesc ( CD ) ; 
  title ('Covariance of data with modelization uncertainty ','fontsize',14) ; colorbar
  xlabel ('Columns','fontweight', 'bold','fontsize',14);
  ylabel ('Rows','fontweight', 'bold','fontsize',14); 
hold off
  
%  New coordinate system for model.

mtil = mprior' + ( G'*CD^-1*G + CM^-1 )^-1 * G'*CD^-1*( dobs - G*mprior' );
dt = G*mtil ; % mtil posterior. 
sum(dobs) ; sum(mtil) ; % Check summes are equal. 
 
%  mtil element number converted to data bin number.
dumC1 = zeros (1,mPl)' ; mtilCat1 = cat (1,dumC1,mtil) ; 
dumC2 = zeros ( 1,max(Bin)-mPu )' ; mtilCat2 = cat (1,mtilCat1,dumC2) ;
LOpt = MeanL(mtilCat2,Bin,unitd) ;

figure(7) 
 bar( (Bin-0.5)*unitd,dobs,1.0, 'facecolor','gr'); hold on
 plot ( (Bin-0.5)*unitd,dt ) ;
  title ('Measured and simulated','fontsize',14) 
  legend ( 'Observed', 'Simulated' )
  axis ([0 20 0 1.1*MaxH]);
  xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',14);
  ylabel ('Frequency','fontweight', 'bold','fontsize',14); 
  print figure(7) -mono figMeasSim.pdf 
hold off 
 
xm = mPl+0.5 : mPu-0.5 ;
 
figure(8) 
 plot ( xm,mtil ) ; hold on; plot ( xm,mprior' ) ; %plot ( Bin-0.5,dobs )
  title ('Posterior and prior ','fontsize',14) 
  legend ( 'Posterior' , 'Prior' )
  xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',14);
  ylabel ('Frequency','fontweight', 'bold','fontsize',14);    
hold off
 
figure(9) % Deconvolved of measured histogram.  
	bar( (Bin-0.5)*unitd,mtilCat2,1.0, 'facecolor','gr'); hold on ; 
   title ('Posterior','fontsize',14) 
	 xlabel ('Track length/  \mum','fontweight', 'bold','fontsize',14);
	 ylabel ('Number','fontweight', 'bold','fontsize',14);
   set(gca (),'fontsize',10,'tickdir','out'); % sets font of numbers on axes 		
	 axis ([0 20 0 1.1*max(mtil)]); % min and max limits.
   text (1,0.7*max(mtilCat2),['Tracks ',num2str( sum(mtil),2 );num2str( (RHOs)/1000000,2), ...
   'cm^-^2';num2str( LOpt,2 ), ' \mum'  ],'fontweight', 'bold','fontsize',14)	 
	 print figure(9) -mono fig1Dec.pdf
hold off ;

cummtilC2 = flipud (cumsum ( flipud (mtilCat2) )) ;

figure(10) 
 plot ( (Bin-0.5)*unitd,cummtilC2' )
  title ('Cummulated posterior models ','fontsize',14) 
  xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',14);
  ylabel ('Frequency','fontweight', 'bold','fontsize',14);    
hold off
  
%  -------------- Covariance of model ----------

CMt = ( G'*CD^-1*G + CM^-1 )^-1 ; % Covariance of best estimate model.

figure(11) 
 hold on;
 imagesc ( CMt ) ;  
  title ('Covariance of posterior model ','fontsize',14) ; colorbar
  xlabel ('Columns','fontweight', 'bold','fontsize',14);
  ylabel ('Rows','fontweight', 'bold','fontsize',14); 
hold off

Vmtil = diag(CMt)  ; % Variance. 
Dmtil = sqrt(Vmtil) ; % Vmtil in accordance truncated Bin.
VmtilCat1 = cat (1,dumC1,Vmtil) ; % Matrix extended.
VmtilCat2 = cat (1,VmtilCat1,dumC2) ; 
DmtilC2 = sqrt(VmtilCat2) ;
 
figure(12) 
 errorbar((Bin-0.5)*unitd,mtilCat2,DmtilC2) 
  title ('Posterior model and Std. dev.','fontsize',14) 
  xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',14);
  ylabel ('Frequency','fontweight', 'bold','fontsize',14);    
hold off

% *********   Variance-Covariance of cumulated posterior model  ***************

mirCMt = flip (flip(CMt,2),1) ; 

figure(13) %  rotated CM (preparation for cumulation).
 hold on;
 imagesc ( mirCMt ) ;  
  title ('Rotated Covariance of posterior ','fontsize',14) ; colorbar
  xlabel ('Column','fontweight', 'bold','fontsize',14);
  ylabel ('Row','fontweight', 'bold','fontsize',14); 
hold off

 
% Summation, backwards from long to small tracks.
    dumCMtr = zeros ( length(mprior),1 ) ;
for kk = 1:length(mprior) ;
    dumres = resize( mirCMt,kk ) ;
    dumCMtr(kk) = sum( dumres(:) ) ;
endfor

flpCMtr = ( flipud ( dumCMtr ) )' ; 
cumVmtilC2 = (zeros ( 1,length(mtilCat2) ) )' ;
cumVmtilC2(mPu-length(flpCMtr)+1 : mPu) = flpCMtr(1: length(mtil)  ); 
cumDmtilC2 = sqrt(cumVmtilC2) ; % Cummulated variance considering correlation.

figure(14) 
 errorbar ( (Bin-1)*unitd,cummtilC2, cumDmtilC2 )
  title ('Cummulated posterior','fontsize',14) 
  xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',14);
  ylabel ('No. tracks','fontweight', 'bold','fontsize',14);    
hold off
   
   
%  Resolution matrix. High values = bad resolution.
Res = eye( length(mtil),length(mtil)) - CMt*(CM^-1) ;

figure(15) 
 hold on;
  imagesc ( Res ) ; colorbar
  title ('Resolution','fontsize',14) 
  xlabel ('Column','fontweight', 'bold','fontsize',14);
  ylabel ('Row','fontweight', 'bold','fontsize',14);    
 hold off

 
% ************************************************************************
%	     Calculation of track ages using deconvolved histogram 
% ************************************************************************

% New binning prepared used to relate tr. density to mean track length.
a = 1  ;
Binp = L0*( 1 + log ( xm' ./ (a*L0) )/bpar ); % Valid only for L/L0 > a*exp(-b)
Llim = L0*a*e^-bpar ; % Binp should not be used when Binp >= Llim
printf('Lower limit for use of length, micrometer: %2.1f micr. m  \n', Llim );
Binp0 = max(Binp,0) ; % Negative values zeroed.


% *** Vector of time intervals and surface densities ***
AgeRho = Agedec(RHOs,ksipp,LambD,LambF,c238,mPl,mPu,mtil,G,Binp0,L0) ;
OldAge = AgeRho (1,1) ; 
Nmtil = sum(mtil) ; disp('   ')
disp('******* Deconvolution of measured histogram **********  ') 
printf('Number deconv. hist. corrected: %3.1f  \n', Nmtil )
printf('Mean Length decon hist.: %3.1f micr. m \n', LOpt )
RHOmtil = sum ( AgeRho (:,3) ) ;
printf('Surface density decon. hist: %3.2f 10^6 cm^-2 \n',...
RHOmtil/1e6 );
printf('Old. tr. age decon.: %4.1f Ma  \n', OldAge );


%    ***************** Variance of ages ******************************

tiV = AgeVar(RHOs,RHOsN,ksipp,ksippD,LambF,LambFD,c238,c238D,...
  mtil,Vmtil,CMt,mirCMt,mtilCat2,flpCMtr,bpar,bparD,mPu,Binp,Binp0,L0,L0D ) ;
tiD = sqrt ( abs(tiV) ) ; 
printf('Std. dev. old. tr. age decon.: %4.1f Ma  \n', tiD(1) );

cumAge = AgeRho(:,1) ;    
%  cumAge and tiD are moved to original binning, Bin.
dumCA1 = zeros (1,mPl)' ; cumAgeC1 = cat (1,dumCA1,cumAge) ; 
dumCA2 = zeros ( 1,max(Bin)-mPu )' ; cumAgeC2 = cat (1,cumAgeC1,dumCA2) ;
dumCt1 = zeros (1,mPl)' ; tiDC1 = cat (1,dumCt1,tiD) ; 
dumCt2 = zeros ( 1,max(Bin)-mPu )' ; tiDC2 = cat (1,tiDC1,dumCt2) ;

figure(16) 
 errorbar ( (Bin-1)*unitd,cumAgeC2,tiDC2 )
  title ('Cummulated track age','fontsize',14) 
  xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',14);
  ylabel ('No. tracks','fontweight', 'bold','fontsize',14);    
 hold off

%   Selecting post-depositional histogram, oldest column idep.

idepdum = find(cumAgeC2 >= 170) ; idep = max(idepdum) ;
dum0(1:idep) = 0 ;
cumAgepost = [ dum0 , [cumAgeC2(idep+1 : length(cumAgeC2) )]' ]' ;
  
mtilpost = [ dum0 , [mtilCat2(idep+1 : length(cumAgeC2) )]' ]' ;
Dmtilp = [ dum0 , [DmtilC2(idep+1 : length(cumAgeC2) )]' ]' ; 
OldAgep = max (cumAgepost) ; 
 
figure(17) 
 bar( (Bin-0.5)*unitd,cumAgeC2,1.0, 'facecolor','gr'); hold on ; 
  title ('Cummulated age ','fontsize',14) 
	xlabel ('Track length/  \ mum','fontweight', 'bold','fontsize',14);
	ylabel ('Age/Ma','fontweight', 'bold','fontsize',14);
  set(gca (),'fontsize',10,'tickdir','out'); % sets font of numbers on axes 		
	axis ([0 20 0 1.1*( max(cumAgeC2)+max(tiDC2) )]); % min and max limits
  text (0.5,1.15*AgeSed,'Deposition','fontweight', 'bold','fontsize',14)	
	text (0.5,0.85*AgeSed,'Age','fontweight', 'bold','fontsize',14)	
 plot( [0 20'],[AgeSed AgeSed]' ); hold on
 bar( (Bin-0.5)*unitd,cumAgepost,1.0, 'facecolor','r'); hold on
% errorbar( (Bin-0.5)*unitd,cumAgeC2,tiDC2,'.') ;
 h = errorbar( (Bin-0.5)*unitd,cumAgeC2,tiDC2,'.') ; % Handle.
  hc = get(h,'Children') ;
  set(hc(1),'color','r')
  set(hc(2),'color','black')  
  ylabel ('Age / Ma','fontweight', 'bold','fontsize',14);  
	print figure(17) -mono fig12Age.pdf  
  hold off

Hist = mtilpost ;  
Lmpost = MeanL(Hist,Bin,unitd) ;
  
   
%  *******  Forward convolution of deconvolved post- sed histogram.  
  
mtilp = mtilpost (mPl+1:mPu) ; mtilcon = G*mtilp ;

% Vector of surface densities are calculated:

cumRho = cumsum( AgeRho(:,3) ) ;
dumRho1 = zeros (1,mPl)' ; cumRho1 = cat (1,dumRho1,cumRho) ; 
dumCR2 = zeros ( 1,max(Bin)-mPu )' ; cumRho2 = cat (1,cumRho1,dumCR2) ;
cumRhopost = [ dum0 , [cumRho2(idep+1 : length(cumRho2) )]' ]' ;
RHOsp = sum ( cumRhopost ) ; 

cumAgepost = zeros( length(cumRho2) ,1) ;
cumAgepost = [ dum0 , [cumAgeC2(idep+1 : length(cumRho2) )]' ]' ;

 figure(18) % Posterior.  
	bar( (Bin-0.5)*unitd,mtilpost,1.0, 'facecolor','gr'); hold on ; 
   title ('Deconvolved post-depositional ','fontsize',14) 
	 xlabel ('Track length/  \mum','fontweight', 'bold','fontsize',14);
	 ylabel ('Number','fontweight', 'bold','fontsize',14);
   set(gca (),'fontsize',10,'tickdir','out'); % sets font of numbers on axes 		
	 axis ([0 20 0 1.1*max(mtilpost)]); % min and max limits
   text (1,0.6*max(mtilpost+Dmtilp),['Tracks ',num2str( sum(mtilpost),2 ); ...
   num2str( (RHOsp)/1000000,2),' cm^-^2';num2str( Lmpost,3 ), ' \mum'  ], ...
   'fontweight', 'bold','fontsize',14)	
	 print figure(18) -mono fig1Dec.pdf
hold off ;  

 figure(19) 
	bar( (Bin-0.5)*unitd,mtilcon,1.0, 'facecolor','gr'); hold on ; 
   title ('Convolved of posterior post-dep. ','fontsize',14) 
	 xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',14);
	 ylabel ('Number','fontweight', 'bold','fontsize',14);
   set(gca (),'fontsize',10,'tickdir','out'); % sets font of numbers on axes 		
	 axis ([ 0 20 0 1.1*max(mtilcon) ]); % min and max limits
   text (1,0.9*max(mtilcon),['Tracks ',num2str( sum(mtilp),2 ); ...
   num2str( (RHOsp)/1000000,2),' cm^-^2';num2str( Lmpost,3 ), ...
   ' \mum'  ],'fontweight', 'bold','fontsize',14)	
	 print figure(19) -mono fig1Dec.pdf
hold off ;  
  
%  Variance of forward model of data.
Vmtilp = Dmtilp .^2 ;
Devfd = sqrt ( G*Vmtilp(mPl+1:mPu) + GVar*mtilp ) ;
  
  sqDia = sqrt(Diag) ;  Fsqdia = find( sqDia >= 0.5 ) ; 
figure(20) ; % Plot all. 
 hold off ;
 subplot(3,2,1); % Measured hist.
	bar( (Bin-0.5)*unitd,Hmeas,1.0, 'facecolor','gr'); hold on;
  errorbar( (Bin(Fsqdia)-0.5)*unitd,Hmeas(Fsqdia),sqDia(Fsqdia),"."); hold on;
  plot ( (Bin-0.5)*unitd,dt ) ; hold on ; 
   title ('Measured and simulated','fontsize',11) ;   
	 xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
	 ylabel ('Number','fontweight', 'bold','fontsize',11);		
	 axis ([0 20 0 1.1*MaxH]); % min and max limits
   set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes    
   text (17,0.9*MaxH,'a','fontweight', 'bold','fontsize',11)	
   text (1,0.7*MaxH,['Tracks ',num2str(NHmeas);num2str( (RHOs)/1000000,2), ...
   'cm^-^2';num2str( LHmeas,3 ), ' \mum'  ],'fontweight', 'bold','fontsize',11)	
   hold off;

   sqDia = DmtilC2 ;  Fsqdia = find( sqDia >= 0.5 ) ;
 subplot(3,2,2); % Decon. hist.
 hold off ;
 bar( (Bin-0.5)*unitd,mtilCat2,1.0, 'facecolor','gr'); hold on ; 
 errorbar( (Bin(Fsqdia)-0.5)*unitd,mtilCat2(Fsqdia),DmtilC2(Fsqdia),".") ;
  title ('Deconvolved of measured ','fontsize',11) 
	xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
	ylabel ('Number','fontweight', 'bold','fontsize',11);
  set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes 		
	axis ([0 20 0 1.1*max(mtilCat2+DmtilC2) ]); % min and max limits
  text (17,0.9*max(mtilCat2+DmtilC2),'b','fontweight', 'bold','fontsize',11)  
  text (1,0.7*max(mtilCat2+DmtilC2),['Tracks ',num2str( sum(mtil),2 );num2str( (RHOs)/1000000,2), ...
  'cm^-^2';num2str( LOpt,2 ), ' \mum'  ],'fontweight', 'bold','fontsize',11)	  
  hold off ;
  sqDia = tiDC2 ;  Fsqdia = find( sqDia >= 0.5 ) ;
subplot(3,2,3); % Cummulated age.
  bar( (Bin-0.5)*unitd,cumAgeC2,1.0, 'facecolor','gr'); hold on ; 
	 xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
   ylabel ('Age / Ma','fontweight', 'bold','fontsize',11);     
   set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes 		
	 axis ([0 20 0 1.1*max(cumAgeC2+tiDC2)]); % min and max limits
   title ('Cummulated age ','fontsize',11) 
 plot( [0 20'],[AgeSed AgeSed]'); hold on 
 bar( (Bin-0.5)*unitd,cumAgepost,1.0, 'facecolor','r'); hold on;
  text (17,0.9*max(cumAgeC2+tiDC2),'c','fontweight', 'bold','fontsize',11)  
	text (0.5,1.3*AgeSed,['Dep.';'age';num2str( AgeSed,3 ),' Ma'],'fontweight', ... 
  'bold','fontsize',11)  
 h = errorbar( (Bin(Fsqdia)-0.5)*unitd,cumAgeC2(Fsqdia),tiDC2(Fsqdia),".") ;hold on; 
  hc = get(h,'Children') ;
  set(hc(1),'color','r')
  set(hc(2),'color','black')  
  hold off
  sqDia = Dmtilp ;  Fsqdia = find( sqDia >= 0.5 ) ; 
subplot(3,2,4); % Decon. post. hist.
 bar( (Bin-0.5)*unitd,mtilpost,1.0, 'facecolor','r'); hold on ; 
 errorbar( (Bin(Fsqdia)-0.5)*unitd,mtilpost(Fsqdia),Dmtilp(Fsqdia),".") ; hold on;  
  title ('Deconvolved post-depositional ','fontsize',11) 
	xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
	ylabel ('Number','fontweight', 'bold','fontsize',11);
  set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes 		
	axis ([0 20 0 1.1*max(mtilpost+Dmtilp)]); % min and max limits
  text (17,0.9*max(mtilpost+Dmtilp),'d','fontweight', 'bold','fontsize',11)  
  text (1,0.6*max(mtilpost+Dmtilp),['Tracks ',num2str( sum(mtilpost),2 ); ...
  num2str( (RHOsp)/1000000,2),' cm^-^2';num2str( Lmpost,3 ), ' \mum'  ], ...
  'fontweight', 'bold','fontsize',11)	
  hold off;
   sqDia = Devfd ;  Fsqdia = find( sqDia >= 0.5 ) ;
	subplot(3,2,5); % Colvolved post-dep. 
	 bar( (Bin-0.5)*unitd,mtilcon,1.0, 'facecolor','r'); hold on ; 
   errorbar( (Bin(Fsqdia)-0.5)*unitd,mtilcon(Fsqdia),Devfd(Fsqdia),".") ; hold on;   
    title ('Convolved of deconvolved post-dep. ','fontsize',11) 
	  xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
	  ylabel ('Number','fontweight', 'bold','fontsize',11);
    set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes 		
	  axis ([0 20 0 1.1*max(mtilcon+Devfd)]); % min and max limits
    text (17,0.9*max(mtilcon+Devfd),'e','fontweight', 'bold','fontsize',11)  
    text (1,0.7*max(mtilcon+Devfd),['Tracks ',num2str( sum(mtilp),2 ); ...
    num2str( (RHOsp)/1000000,2),' cm^-^2';num2str( Lmpost,3 ), ...
    ' \mum'  ],'fontweight', 'bold','fontsize',11)	   
%	  print -mono fig4all20.jpg ; 
	  print figure(20) -mono fig4all20.jpg ; % pdf ; %png ; jpg ;
hold off ;  


 %  Figure 21. Plot all on large paper.
    figRx = 1200 ; figRy = 1200  ; % 1500.
figure('Position',[7 7 figRx figRy]) ; % 
  hold off ;
 subplot(3,2,1); % Measured hist. 
	bar( (Bin-0.5)*unitd,Hmeas,1.0, 'facecolor','gr'); hold on;
  errorbar( (Bin(Fsqdia)-0.5)*unitd,Hmeas(Fsqdia),sqDia(Fsqdia),"."); hold on;
  plot ( (Bin-0.5)*unitd,dt ) ; hold on ; 
   title ('Measured and simulated','fontsize',11) ;   
%	 xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
	 ylabel ('Number','fontweight', 'bold','fontsize',11);		
	 axis ([0 20 0 30]); % min and max limits, 1.1*max(mtilCat2+DmtilC2)
   set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes    
   text (17,0.9*max(mtilCat2+DmtilC2),'a','fontweight', 'bold','fontsize',11)	
   text (1,0.7*max(mtilCat2+DmtilC2),['Tracks ',num2str(NHmeas);num2str((RHOs)/1000000,2), ...
   'cm^-^2';num2str( LHmeas,3 ), ' \mum'],'fontweight', 'bold','fontsize',11)	
   hold off;
   
   sqDia = DmtilC2 ;  Fsqdia = find( sqDia >= 0.5 ) ;
 subplot(3,2,2); % Decon. hist.
 hold off ;
 bar( (Bin-0.5)*unitd,mtilCat2,1.0, 'facecolor','gr'); hold on ; 
 errorbar( (Bin(Fsqdia)-0.5)*unitd,mtilCat2(Fsqdia),DmtilC2(Fsqdia),".") ;
  title ('Deconvolved of measured ','fontsize',11) 
%	xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
	ylabel ('Number','fontweight', 'bold','fontsize',11);
  set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes 		
	axis ([0 20 0 30 ]); % min and max limits, 1.1*max(mtilCat2+DmtilC2)
  text (17,0.9*max(mtilCat2+DmtilC2),'b','fontweight', 'bold','fontsize',11)  
  text (1,0.7*max(mtilCat2+DmtilC2),['Tracks ',num2str( sum(mtil),2 );num2str( (RHOs)/1000000,2), ...
  'cm^-^2';num2str( LOpt,2 ), ' \mum'  ],'fontweight', 'bold','fontsize',11)	  
  hold off ;

  sqDia = tiDC2 ;  Fsqdia = find( sqDia >= 0.5 ) ;
subplot(3,2,3); % Cummulated age.
  bar( (Bin-0.5)*unitd,cumAgeC2,1.0, 'facecolor','gr'); hold on ; 
%	 xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
   ylabel ('Age / Ma','fontweight', 'bold','fontsize',11);     
   set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes 		
	 axis ([0 20 0 1.1*max(cumAgeC2+tiDC2)]); % min and max limits
   title ('Cummulated age ','fontsize',11) 
 plot( [0 20'],[AgeSed AgeSed]'); hold on 
 bar( (Bin-0.5)*unitd,cumAgepost,1.0, 'facecolor','r'); hold on;
  text (17,0.9*max(cumAgeC2+tiDC2),'c','fontweight', 'bold','fontsize',11)  
	text (0.5,1.4*AgeSed,['Dep.';'age';num2str( AgeSed,3 ),' Ma'],'fontweight', ... 
  'bold','fontsize',11)  
 h = errorbar( (Bin(Fsqdia)-0.5)*unitd,cumAgeC2(Fsqdia),tiDC2(Fsqdia),".") ;hold on; 
  hc = get(h,'Children') ;
  set(hc(1),'color','r')
  set(hc(2),'color','black')  
  hold off
  
  sqDia = Dmtilp ;  Fsqdia = find( sqDia >= 0.5 ) ; 
subplot(3,2,4); % Decon. post. hist.
 bar( (Bin-0.5)*unitd,mtilpost,1.0, 'facecolor','r'); hold on ; 
 errorbar( (Bin(Fsqdia)-0.5)*unitd,mtilpost(Fsqdia),Dmtilp(Fsqdia),".") ; hold on;  
  title ('Deconvolved post-depositional ','fontsize',11) 
	xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
	ylabel ('Number','fontweight', 'bold','fontsize',11);
  set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes 		
	axis ([0 20 0 1.1*max(mtilpost+Dmtilp)]); % min and max limits
  text (17,0.9*max(mtilpost+Dmtilp),'d','fontweight', 'bold','fontsize',11)
%  text (17,0.9*max(mtilCat2+DmtilC2),'b','fontweight', 'bold','fontsize',11)  
  text (1,0.6*max(mtilCat2+DmtilC2),['Tracks ',num2str( sum(mtilpost),2 );num2str( (RHOsp)/10^12,2), ...
  'cm^-^2';num2str( LOpt,2 ), ' \mum'  ],'fontweight', 'bold','fontsize',11)	  
%  text (1,0.6*max(mtilpost+Dmtilp),['Tracks ',num2str( sum(mtilpost),2 ); ...
%  num2str( RHOsp/1000000,2),' cm^-^2';num2str( Lmpost,3 ), ' \mum'  ], ...
%  'fontweight', 'bold','fontsize',11)	
  hold off;  
  
   sqDia = Devfd ;  Fsqdia = find( sqDia >= 0.5 ) ;
	subplot(3,2,5); % Colvolved post-dep. 
	 bar( (Bin-0.5)*unitd,mtilcon,1.0, 'facecolor','r'); hold on ; 
   errorbar( (Bin(Fsqdia)-0.5)*unitd,mtilcon(Fsqdia),Devfd(Fsqdia),".") ; hold on;   
    title ('Convolved of deconvolved post-dep. ','fontsize',11) 
	  xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
	  ylabel ('Number','fontweight', 'bold','fontsize',11);
    set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes 		
	  axis ([0 20 0 1.1*max(mtilcon+Devfd)]); % min and max limits
    text (17,0.9*max(mtilcon+Devfd),'e','fontweight', 'bold','fontsize',11)  
    text (1,0.7*max(mtilcon+Devfd),['Tracks ',num2str( sum(mtilp),2 ); ...
    num2str( (RHOsp)/10^12,2),' cm^-^2';num2str( Lmpost,3 ), ...
    ' \mum'  ],'fontweight', 'bold','fontsize',11)	
    
	print figure(21) -mono fig4all21.png ; % pdf ; %png ; jpg ;  
    hold off;
  
 
figure(22)
 bar( (Bin-0.5)*unitd,cumAgeC2,1.0, 'facecolor','gr'); hold on ; 
	xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',10);
  ylabel ('Age / Ma','fontweight', 'bold','fontsize',10);     
  set(gca (),'fontsize',10,'tickdir','out'); % sets font of numbers on axes 		
	axis ([0 20 0 1.1*max(cumAgeC2+tiDC2)]); % min and max limits
  title ('Cummulated age ','fontsize',10) 
 plot( [0 20'],[AgeSed AgeSed]'); % hold on 
 bar( (Bin-0.5)*unitd,cumAgepost,1.0, 'facecolor','r'); 
	text (0.5,1.3*AgeSed,['Dep.';'age';num2str( AgeSed,3 ),' Ma'],'fontweight', ... 
  'bold','fontsize',10) 
 h = errorbar( (Bin-0.5)*unitd,cumAgeC2,tiDC2,'.') ; % Handle.
  hc = get(h,'Children') ;
  set(hc(1),'color','r')
  set(hc(2),'color','black')  
hold off
  
figure(23)  
 bar( (Bin(Fsqdia)-0.5)*unitd,mtilCat2(Fsqdia),1.0, 'facecolor','gr'); hold on ; 
 errorbar( (Bin(Fsqdia)-0.5)*unitd,mtilCat2(Fsqdia),DmtilC2(Fsqdia),".") ;
  title ('Posterior model and Std. dev.','fontsize',14) 
  xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',14);
  ylabel ('Frequency','fontweight', 'bold','fontsize',14);    
 hold off ; 
  
figure(24)   
 hold off ;
 bar( (Bin-0.5)*unitd,mtilCat2,1.0, 'facecolor','gr'); hold on ; 
 errorbar( (Bin-0.5)*unitd,mtilCat2,DmtilC2,".") ;
% errorbar( (Bin(Fsqdia)-0.5)*unitd,mtilCat2(Fsqdia),DmtilC2(Fsqdia),".") ; 
  title ('Deconvolved of measured ','fontsize',11) 
	xlabel ('Track length /  \mum','fontweight', 'bold','fontsize',11);
	ylabel ('Number','fontweight', 'bold','fontsize',11);
  set(gca (),'fontsize',11,'tickdir','out'); % sets font of numbers on axes 		
	axis ([0 20 0 1.1*max(mtilCat2+DmtilC2) ]); % min and max limits 
  text (1,0.7*max(mtilCat2+DmtilC2),['Tracks ',num2str( sum(mtil),2 );num2str( (RHOs)/1000000,2), ...
  'cm^-^2';num2str( LOpt,2 ), ' \mum'  ],'fontweight', 'bold','fontsize',11)	  
  hold off ;  

Function:

%
%     Function calculating: Surface density, timeintervals, cummulated ages.
%     Program called by LsqTrack.
%

function AgeRho = Agedec (RHOs,ksipp,LambD,LambF,c238,mPl,mPu,mtil,G,Binp0,L0) ;

global fac1 fac2 fac3 Binmtil

Binmtil = [mPl:mPu-1]+0.5 ;


% --------- Calculate time intervals, Ma -----------------------

% Calculate time at start of column n, Ma.
% LambF, 1/yr.
fac1 = 2*RHOs*1e4/( ksipp*LambF*1e6*c238 ) ; 
fac2j = sum( (G .* mtil') ./ Binmtil ) ; 
fac2 = triu( ones(length(mtil)) )*fac2j' ;
fac3 = sum( fac2j .*Binp0' )*1e-6  ;

delt = fac1*fac2 / fac3 ;
ti = log(1 + (fac1*LambD*(fac2 ./ fac3)) ) / LambD ; %  Ma, cumulative age.
delti = [(-diff ( [ti;0] )) ] ;

% --------- Cummulated surface density ---------------------
%dum1 = mtil .* Binp0 ;
RHOsi = ( RHOs*Binp0 .* fac2j' )  / fac3 ; % Surface density, tr./cm^2 for each column.
AgeRho = [ ti delti RHOsi ] ; 

endfunction

Function:

%
%     Function calculating variance of age for each column. 
%   

function tiV = AgeVar (RHOs,RHOsN,ksipp,ksippD,LambF,LambFD,c238,c238D,...
  mtil,Vmtil,CMt,mirCMt,mtilCat2,flpCMtr,bpar,bparD,mPu,Binp,Binp0,L0,L0D )

% ********   Variance of age  ***************************
%            ignoring decreasing U- concentration.

% Time Ma., LambF, 1/yr.

global fac1 fac2 fac3 Binmtil
% Binmtil: variance of lambda length.

% **** Variance for factor fac1 *******

RHOsD = RHOs*sqrt(RHOsN)/RHOsN ; % St. d. surface density.
fac1D = fac1*sqrt( (c238D/c238)^2 + (1.0*10^-12*RHOsD/RHOs)^2 + (ksippD/ksipp)^2 + ...
(LambFD/LambF)^2 ) ;

fac1V = fac1D^2 ;


% ************ Variance of factor fac2 *********************
% Summation, backwards from long to small tracks.
    dumCMtr = zeros ( length(mtil),1 ) ;
    mirCMt1 = (mirCMt ./ Binmtil') ./ Binmtil ;     % check !!!!!!!!  
for kk = 1:length(mtil) ;
    dumres = resize( mirCMt1,kk ) ;
    dumCMtr(kk) = sum( dumres(:) ) ;
endfor

flpCMtr = ( flipud ( dumCMtr ) )' ; 

%cumVmtilC2 = (zeros ( 1,length(mtilCat2) ) )' ;
%cumVmtilC2(mPu-length(flpCMtr)+1 : mPu) = flpCMtr(1: length(mtil)  ); 
%cumDmtilC2 = sqrt(cumVmtilC2) ; % Cummulated variance considering correlation.

fac2V = flpCMtr  ;


%fac2jV = sum( (G .* flpCMtr) ./ Binmtil ) ; % See function Agedec.
%fac2V = triu( ones(length(flpCMtr)) )*(fac2jV)' ;


% **********  Variance of factor fac3  *********************
%variance of lambda length.
%fac3iV = 1.0*10^-12 *(Binp0.^2).*fac2V + (fac2.^2).*BinmtilV ;
%fac3V = sum(fac3iV)

% Fac31 as summation, backwards from long to small tracks.
    dumCMtr = zeros ( length(mtil),1 ) ;
    mirCMt1 = (mirCMt ./ (Binmtil'./Binp0' ) ) ./ (Binmtil./Binp0) ; % check !!!!!!!!  
for kk = 1:length(mtil) ;
    dumres = resize( mirCMt1,kk ) ;
    dumCMtr(kk) = sum( dumres(:) ) ;
endfor

fac31V = sum( ( flipud ( dumCMtr ) )' ) ; 
%flpCMtr = ( flipud ( dumCMtr ) )' ; 

VarBinp0 = L0^2*( log(Binmtil/L0) ).^2*bparD^2 / bpar^4 ;
fac32V = sum( ( mtil ./ Binmtil' ).^2 .*VarBinp0' ) ;

%cumVmtilC2 = (zeros ( 1,length(mtilCat2) ) )' ;
%cumVmtilC2(mPu-length(flpCMtr)+1 : mPu) = flpCMtr(1: length(mtil)  ); 
%cumDmtilC2 = sqrt(cumVmtilC2) ; % Cummulated variance considering correlation.
 
fac3V = 1.0*10^-12 * (fac31V + fac32V);

% ************   Variance of ages   *****************************

% fac2V har længde forskellig fra de andre, skal tiV være 7 eller 24?

tiV = [ ( (fac1 * fac2)/fac3^2 ) .^2 ] * fac3V + ... 
 ( (fac1 / fac3)^2 ) .* fac2V' + ( (fac2 ./ fac3).^2 ) * fac1V ;     
 
endfunction

Function:

%
%         Mean track length of histograms in coulombs.
%

function Lmean = MeanL(Hist,Bin,unit)
NHist = sum(Hist) ;
ML_Hist = sum( Hist .*(Bin-0.5)*unit ) ;
Lmean = ML_Hist ./ NHist ;
endfunction