clear,clc,close all
%-----------------------------------------------
% ZirTiDiS - Zirconium-in-Titanite Diffusion Software
% S. Schorn & E. Moulas - Mainz, 18.05.2024
%for more details see Zenodo documentation
%doi: 10.5281/zenodo.11184086
%%+++++++++++++++++++++++++++++++++++++++++++++++++++
%CALCULATE COOLING HISTORIES
%%+++++++++++++++++++++++++++++++++++++++++++++++++++
%----------------
%IMPLICIT SOLVER
%----------------
plot_now    = 1;              %[0,1] - 1 for plotting on the fly
lookuptable = 0;              %[0,1] - 1 to load precompiled data
safe_data   = 0;              %[0,1] - 1 to safe data
%+++++++++++++++++++++++++++++++++++++++++++++++++++++
%define your system here
grsz     = 50;                %grain radius        [µm]
Tin      = 800;               %initial temperature [°C]
Tfin     = 600;               %final temperature   [°C]
cool_r   = 20;                %cooling rate        [°C/myr]
tstep    = 1e10;              %time step (~1e10 - 1e11)
C2K      = 273;               %convert °C to K
myr2sec  = 60*60*24*365*1e6;  %milion year [s]
%if isobaric: set Pin = Pfin
Pin      = 1.0;               %initial pressure [GPa]
Pfin     = 0.4;               %final pressure   [GPa]
%set activities - used if lookup tabel == 0
aTi_in   = 1.0;               %starting activity TiO2
aTi_fin  = 1.0;               %ending activity TiO2
aSi_in   = 1.0;               %starting activity SiO2
aSi_fin  = 1.0;               %ending activity SiO2
%+++++++++++++++++++++++++++++++++++++++++++++++++++++
%call ZirTiDiS
[r,C,T_ar,P_ar,t_ar,TiT,Tcore,Tp,Pp] = ZirTiDiS_fun(plot_now,lookuptable,grsz,Tin,Tfin,cool_r,Pin,Pfin,aTi_in,aTi_fin,aSi_in,aSi_fin,tstep);

%saving data
if safe_data == 1
    radius  = r;
    Zr_conc = C;
    Tpath   = T_ar;
    Ppath   = P_ar;
    tpath   = t_ar/myr2sec;
    CRate   = cool_r;
    ZiTT    = TiT;
    ZiTcore = Tcore;
    Tpoint  = Tp;
    Ppoint  = Pp;
    clearvars -except radius Zr_conc Tpath Ppath tpath CRate ZiTT ZiTcore Tpoin Ppoint
    save ZirTiDiS_data
end

%output key:
%r     = grain radius [µm]
%C     = final Zr concentration [ppm]
%T_ar  = input Temperature array [K]
%P_ar  = input Pressure array [GPa]
%t_ar  = resulting timre array from cooling rate [s]
%TiT   = final Zr-in-Titanite Temperature for entire grain [K]
%Tcore = Zr-in-Titanite Temperature at the core at each time step [K]
%Tp    = actual Temperature at each time step [K]
%Pp    = Pressure at each time step [GPa]
