% calc spot ages from exported he profile for an infinite cylinder

clear
clc

addpath '/Users/cglotzbach/Documents/src/MATLAB/thermo_models/RDAAM/RDAAM_pro_input_SD_cylinder'

%% input: tT-path
%time=[0 10 200]; % in Ma
%temp=[10 50 210]; % in °C
%time=[0 60 60.5]; % in Ma
%temp=[10 10 210]; % in °C
time=[0 200]; % in Ma
temp=[10 210]; % in °C

%% input: grain and spot data
radius=100; % grain radius
zpos=100; % vertical position in grain
U=10; % U concentration ppm
Th=10; % Th concentration ppm
Sm=10; % Sm concentration ppm
spot_size=15; % spot radius in microns
spot_depth=5; % spot depth in microns

%% fixed parameters
Avg=6.022E23; % atoms/mol
ApRho=3.2; % density of apatite in g/cm^3
% molar masses in g/mol
mm238=238.0508;
mm235=235.0439;
mm232=232.0380;
mm147=146.9149;
R235238=1/136.1384; % natural U isotope ratio

%% input vectors for U,Th,Sm
%Uvec(1:100)=U;
%Thvec(1:100)=Th;
%Smvec(1:100)=Sm;
% zoned grain with 5x higher in the outer 10%
Uvec(1:90)=U; Uvec(91:100)=U*5;
Thvec(1:90)=Th; Thvec(91:100)=Th*5;
Smvec(1:90)=Sm; Smvec(91:100)=Sm*5;
% interpolate to rdim vector length
rdim=128;
Uvec=interp1(1:length(Uvec),Uvec,linspace(1,length(Uvec),rdim));
Thvec=interp1(1:length(Thvec),Thvec,linspace(1,length(Thvec),rdim));
Smvec=interp1(1:length(Smvec),Smvec,linspace(1,length(Smvec),rdim));


%% calculate He profile of a cylinder
[ap_corr_age,ap_age,ap_total_He,HeProfile]=RDAAM_cylinder_ProDiffMEX(time,temp,radius,Uvec,Thvec,Smvec);


%% simplify HeProfile to length=100
sizeHeProfile=length(HeProfile);
if length(HeProfile)>100
    HeProfile=interp1(linspace(0,1,sizeHeProfile),HeProfile,linspace(0,1,101));
end


%% transform raw He profile data
dx=radius/(length(HeProfile)-1);
HeGrain=HeProfile;
DistanceGrain=0:dx:radius;


%% calculate area and volume of spot, and volume of column of pixels
spot_area=pi()*spot_size^2; % in micrometer
spot_volume=spot_area*spot_depth; % in micrometer
pixel_mass=dx*dx*dx*ApRho*1E-12; % in gram
spot_mass=spot_volume*ApRho*1E-12; % in gram


%% transform 1D profile to 3D sphere
M=prof1Dto3Dcyl(length(HeProfile)-1);
%M=squeeze(M(round(zpos/dx),:,:));
M=round(M*(length(HeProfile)-1))+1;
HeGrainAtGr=HeGrain*1E-9*Avg; % first He in at/g
HeGrainAtPixel=HeGrainAtGr*pixel_mass; % second He in each column in at
% map He to circle plane
HeGrain3D=zeros(size(M));
for i=1:size(HeGrain3D,1)
    for j=1:size(HeGrain3D,2)
        for k=1:size(HeGrain3D,3)
            if M(i,j,k)>0
                HeGrain3D(i,j,k)=HeGrainAtPixel(M(i,j,k));
            end
        end
    end
end

% make pixels in the center different than 1
M(radius/dx+1,radius/dx+1,:)=1.0001;
% make values outside the grain to zero, assuming no He outside the grain
I=find(M==1);
I(I==sub2ind(size(M),101,101,101))=[];
HeGrain3D(I)=0;
% create disk filter
disk_filter=fspecial("disk",spot_size/dx);
disk_filter=disk_filter./max(disk_filter(:));
% calculate cumulative He for specific spot sizes throughout the grain
HeGrainSpot=zeros(size(M,1),size(M,2));
for i=1:floor(spot_depth/dx)
    HeGrainSpot=HeGrainSpot+filter2(disk_filter,squeeze(HeGrain3D(round(zpos/dx)+i-1,:,:)));
end
rest=rem(spot_depth,dx);
% add He in-between pixels
if rest>0
    HeGrainSpot=HeGrainSpot+rest/dx*filter2(disk_filter,squeeze(HeGrain3D(round(zpos/dx)+i-1,:,:)));
end
% get values along a profile through the grain
HeGrainSpotProfile=HeGrainSpot(:,round(radius/dx)+1);
% make values to NaNs at borders (spots that overlap with grain border)
HeGrainSpotProfile(1:round(spot_size/dx))=NaN;
HeGrainSpotProfile(end-round(spot_size/dx)+1:end)=NaN;


%% calculate uncertainty of modelled He atoms 
LOD=0.000079; % limit of detection in ncc (see Calculate LOD.xlsx for details)
HeMolMass=4.46E-14; % mol/ncc
LODAt=LOD*Avg*HeMolMass; % limit of detection in atoms
HeSDGrainSpotProfile=sqrt(LODAt./HeGrainSpotProfile).^2*100;


%% initial parent nuclide concentration in atoms
U238spot=U*1E-6*spot_mass*Avg/mm238;
U235spot=U*1E-6*R235238*spot_mass*Avg/mm235;
Th232spot=Th*1E-6*spot_mass*Avg/mm232;
Sm147spot=Sm*1E-6*spot_mass*Avg/mm147;


%% calculate ages along transect
transect_ages=zeros(size(HeGrainSpotProfile));
for i=1:length(HeGrainSpotProfile)
    if HeGrainSpotProfile(i)>0
        [transect_ages(i),~] = calc_age_uncorr('Ap',HeGrainSpotProfile(i),0,U238spot,U235spot,0,Th232spot,0,Sm147spot,0);
    end
end
transect_ages(transect_ages<=0)=NaN;


%% plot He along transect
cmp=colormap('parula');
%cmp = flipud(cmp);
colormap(cmp);
x=[DistanceGrain radius+DistanceGrain(1:end-1)+1];
%figure
yyaxis left
patch(x-radius,transect_ages,HeSDGrainSpotProfile,'LineWidth',2,'EdgeColor','interp')
clim([5 50])
xlabel('Distance from center (\mum)')
ylabel('Apatite (U-Th)/He age (Ma)')
hold on
yline(ap_corr_age,'--','LineWidth',2)
yyaxis right
plot(x-radius,[fliplr(HeGrainAtGr) HeGrainAtGr(2:end)],'LineWidth',2)
ylabel('He concentration (At/g)')
set(gca,'FontSize',14)
colorbar

