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

• Main Author: Jensen, Peter Klint
• Format: Article in Journal/Newspaper
• Language: English
• Published: Zenodo 2021
• Subjects: fission tracks; age length distribution; probabilistic inversion; Greenland; Rho; Eta; Jameson Land; East Greenland
• Online Access: https://dx.doi.org/10.5281/zenodo.4553144; https://zenodo.org/record/4553144

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: [e-mail]. % % 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 : {"references": ["Jensen, P.K., Hansen, K.: Age distribution of horizontal fission tracks, submitted to Geochronology.", "Jensen, P. K., Hansen, K., and Kunzendorf, H.: A numerical model for the thermal history of rocks based on confined horizontal fission tracks. International Journal of Radiation Applications and Instrumentation, Part D. Nucl. Tracks Radiat. Meas., 20, 2, 349-359, doi: 10.1016/1359-0189(92)90064-3, 1992.", "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.", "Tarantola, A.: Inverse problem theory, SIAM, Philadelphia, doi:0.1137/1.9780898717921, 2005."]}