% 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 parameters
radius=100; % grain radius
% radionuclides of rim in ppm
U=100;
Th=100;
Sm=100;
U(1:100)=10;
Th(1:100)=10;
Sm(1:100)=10;

% interpolate to rdim vector length
rdim=128;
U=interp1(1:length(U),U,linspace(1,length(U),rdim));
Th=interp1(1:length(Th),Th,linspace(1,length(Th),rdim));
Sm=interp1(1:length(Sm),Sm,linspace(1,length(Sm),rdim));
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)


