% call RDAAM and ZRDAAM:
% - does handle U,Th,Sm profiles derived from e.g. LA-ICP-MS measurements
% - uses all important stopping distances making up >99% of all alphas
% - does calculate profile for short axis of cylinder, assuming >100
% microns prism length

clear all
clc

Avg=6.022E23; % atoms/mol
upramp=4; % min
downramp=50; % min
duration=10; % min
temp_high=380; % °C
downrampT=1-downramp/(60*24*365*1E6);  % in Myr
durationT=1+duration/(60*24*365*1E6); % in Myr
uprampT=durationT+upramp/(60*24*365*1E6);  % in Myr

% time-temperature path
time=[0 downrampT 1 durationT uprampT 10]; % in Myr
%time=[0 1 1 1.000000000007610 10];
temp=[20 20 temp_high temp_high 20 20];

% grain radius
radius=100; % grain radius
% radionuclide concentrations
rdim=128;
% homogenous concentration
U(1:rdim)=10;
Th(1:rdim)=10;
Sm(1:rdim)=10;
% zonation, with rim beeing enriched 5x
% U(100:rdim)=50;
% Th(100:rdim)=50;
% Sm(100:rdim)=50;

% x vector required for plotting
x=linspace(0,radius,rdim);

% HeProfile is in nmol/g similar to HeFTy
[ap_corr_age,ap_age,ap_total_He,HeProfCyl]=RDAAM_cylinder_ProDiffMEX(time,temp,radius,U,Th,Sm);
ap_age
% He data originally in nmol/g, converted here to at/g
figure 
hold on
plot(x,HeProfCyl*Avg*1E-9,LineWidth=2,Color='r')
xlabel('Distance (μm)')
ylabel('He (at/g)')
set(gca,'FontSize',14)


