function [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)
%-----------------------------------------------
% ZirTiDiS - Zirconium-in-Titanite Diffusion Software
% S. Schorn & E. Moulas - Mainz, 13.05.2024
%for more details see Zenodo documentation
%doi: 10.5281/zenodo.11184086
%-----------------------------------------------
draw_now = plot_now;                          %[0,1] - 1 for plotting on the fly
activity_from_pseudosection = lookuptable;    %[0,1] - 1 to load precompiled data
%---------------------------------
%diffusion parameters from Cherniak
%Contrib Mineral Petrol (2006) 152:639–647
%-----------
%Zr-in-titanite equation from Hayden et al
%Contrib Mineral Petrol (2008) 155:529–540
%---------------------------------
%+++++++++++++++++++++++++++++++++++++++++++++++++++++
%define your system here
GR       = grsz;             %grain radius        [µm]
Tinitial = Tin;              %initial temperature [°C]
Tfinal   = Tfin;             %final temperature   [°C]
CR       = cool_r;           %cooling rate        [°C/myr]
%if isobaric: set PGastop = PGPa0
PGPa0    = Pin;              %initial pressure [GPa]
PGPastop = Pfin;             %final pressure   [GPa]
%variable activities - preset values for linear change
aTiO20   = aTi_in;           %starting activity TiO2 [0..1]
aTiO2end = aTi_fin;          %ending activity TiO2 [0..1]
aSiO20   = aSi_in;           %starting activity SiO2 [0..1]
aSiO2end = aSi_fin;          %ending activity SiO2 [0..1]
%+++++++++++++++++++++++++++++++++++++++++++++++++++++
%Physics - DO NOT CHANGE
D0      = 5.33*1e-7;         %pre-exponent [m2/s]
Q       = 325000;            %activation energy [J/mol]
R       = 8.314;             %gas constant [J/mol K]
C2K     = 273;               %convert °C to K
myr2sec = 60*60*24*365*1e6;  %milion year [s]
L       = GR*1e-6;           %domain size (grain radius) [m]
TK0     = Tinitial + C2K;    %starting T [K]
TKstop  = Tfinal   + C2K;    %ending T [K]
%-----------------------------
%Numerics
dt      = tstep;             %time step (~1e10 - 1e11)
nr      = 200;               %spatial resolution
dr      = L/(nr-1);          %spatial step
r       = 0:dr:L;            %discretize radius; n elements
rc      = 0.5*r(1:end-1) + 0.5*r(2:end);    %centerpoints between nodes; n-1 elements
rcL     = rc(1:end-1);       %centerpoints for left fluxes; n-2 elements
rcR     = rc(2:end);         %centerpoints for right fluxes; n-2 elements
nout    = 1000;              %plot each nout times
%Preporcessing - creating a linear cooling/decompression path
CR      = CR/myr2sec;               %CR [°C/s]
t_tot   = (TK0-TKstop)/CR;          %total time given by T interval and CR [s]
DPR     = (PGPa0-PGPastop)/t_tot;   %decompression rate [GPa/s]
DaTiO2  = (aTiO20-aTiO2end)/t_tot;  %rate of change in aTiO2
DaSiO2  = (aSiO20-aSiO2end)/t_tot;  %rate of change in aTiO2
%initialize
A       = sparse(nr,nr);                %create a sparse nr*nr matrix
TK      = TK0;                          %starting value of T
PGPa    = PGPa0;                        %starting value for P
aTiO2   = aTiO20;                       %starting value for aTiO2
aSiO2   = aSiO20;                       %starting value for aSiO2
t_ar    = linspace(0,t_tot);            %time array
T_ar    = linspace(TK0,TKstop);         %Temperature array
P_ar    = linspace(PGPa0,PGPastop);     %Pressure array
%define Zr-in-titanite temperature [K] from concentration - Hayden equation
TZK  = @(PGPa,aTiO2,aSiO2,C) (7708 + 960*PGPa) ./ (10.52 - log10(aTiO2) - log10(aSiO2) - log10(C));    %[K]
%define equilibrium Zr concentrartion @ starting PT (= BC) - Hayden equation
ZiT  = @(TK,PGPa,aTiO2,aSiO2) 10.^(10.52 - 7708./TK - 960 * PGPa./TK - log10(aTiO2) - log10(aSiO2));   %[ppm]
%-----------------------------------------------------------------
%executed if activities are set by user (NOT from lookup table)
%plot activities & Zr-concentrations
if activity_from_pseudosection == 0
    %1) activities are set by user and constant
    if aTiO2 == aTiO2end && aSiO2 == aSiO2end
        %calculate equilibrium Zr in titanite [ppm] over PT range
        ns = 10;
        ts = linspace(500,1000,ns)+C2K;
        ps = linspace(0,3.0,ns);
        [Ts,Ps] = meshgrid(ts,ps);
        Zr_conc = ZiT(Ts,Ps,aTiO2,aSiO2);
        %plot
        figure(1);
        subplot(121)
        contourf(Ts-C2K,Ps,log10(Zr_conc),100,'edgecolor','none'),shading flat
        colorbar
        hold on
        [Cpl,hpl] = contour(Ts-C2K,Ps,log10(Zr_conc),'k'); %Extract plot properties
        h_line = plot(T_ar-273,P_ar,'LineWidth',2,'Color','r');
        legend([h_line],'PT path')
        legend boxoff
        hold off
        title(['Zr in titanite log10(ppm) | aTiO_2 = ', num2str(aTiO20), ' | aSiO_2 = ', num2str(aSiO20), ' (user fixed values)'])
        xlabel('Temperature (°C)')
        ylabel('Pressure (GPa)')
        axis square
        clabel(Cpl,hpl,'LabelSpacing',1200,'FontSize',12)     %Add values on contours
        set(gca, 'FontSize', 14)
        subplot(122)
        plot(T_ar - C2K,P_ar)
        title('Temperature vs. Pressure')
        xlabel('Temperature (°C)')
        ylabel('Pressure (GPa)')
        grid on
        set(gca, 'FontSize', 14)
        %2) activities are set by user and NOT constant
    else
        aTi_ar = linspace(aTiO20,aTiO2end);
        aSi_ar = linspace(aSiO20,aSiO2end);
        figure(1);
        subplot(131)
        plot(T_ar-C2K,aTi_ar, T_ar-C2K,aSi_ar)
        xlabel('Temperature (°C)')
        ylabel('Activities')
        title('Activity vs. Temperature')
        legend('aTiO2', 'aSiO2', Location='southeast')
        legend boxoff
        grid on
        set(gca, 'FontSize', 14)
        subplot(132)
        plot(P_ar,aTi_ar, P_ar,aSi_ar)
        xlabel('Pressure (GPa)')
        ylabel('Activities')
        title('Activity vs. Pressure')
        legend('aTiO2', 'aSiO2', Location='southeast')
        legend boxoff
        grid on
        set(gca, 'FontSize', 14)
        subplot(133)
        plot(T_ar-C2K,P_ar)
        xlabel('Temperature (°C)')
        ylabel('Pressure (GPa)')
        title('Temperature vs. Pressure')
        grid on
        set(gca, 'FontSize', 14)
    end
    %initialize
    CBC  = ZiT(TK,PGPa,aTiO2,aSiO2);%[ppm]
    C0   = ones(1,nr)*CBC;          %starting with equilibrium Zr-concentration from starting T
    C    = C0;                      %initialize C
    t    = 0;                       %time starts @ 0
    it   = 0;                       %counter starts @ 0
    itp  = 0;                       %counter starts @ 0
    tic
    %solve diffusion equation
    while t < t_tot                 %until time is out
        it = it + 1;                %increase counter
        t  = t + dt;                %track time
        %if we overshoot
        if t > t_tot
            dt    = dt - (t_tot-t); %updated time step
            t     = t_tot;          %final time
            TK    = TKstop;         %final T
            PGpa  = PGPastop;       %final P
            aTiO2 = aTiO2end;       %final aTiO2
            aSiO2 = aSiO2end;       %final aSiO2
        end
        D = D0*exp(-Q/(R*TK));      %update T-dependent diffusivity
        %set up matrix for solver
        for i = 2:nr-1                              %loop through matrix
            j = i-1;                                %differnt index for rcR & rcL
            S = (D*dt)/(r(i)^2*dr^2);               %constants
            A(i,i-1) = -S*(rcL(j)^2);               %below diagonal
            A(i,i)   = 1 + S*(rcL(j)^2 + rcR(j)^2); %@ diagonal
            A(i,i+1) = -S*(rcR(j)^2);               %above diagonal
        end
        A(1,1)   = 1;                               %apply BC
        A(nr,nr) = 1;                               %apply BC
        b        = C(:);                            %define concentration vector
        X        = A\b;                             %apply backslash magic
        C        = X(:)';                           %transpose vector
        C(1)     = C(2);                            %zero-flux Neumann BC at center for symmetry
        C(end)   = CBC;                             %set target concentration
        TK       = TK    - CR*dt;                   %update T @ each step
        PGPa     = PGPa  - DPR*dt;                  %update P @ each step
        aTiO2    = aTiO2 - DaTiO2*dt;               %update aTiO2 @ each step
        aSiO2    = aSiO2 - DaSiO2*dt;               %update aSiO2 @ each step
        CBC      = ZiT(TK,PGPa,aTiO2,aSiO2);        %update BC - Hayden equation
        if mod(it,nout) == 1
            itp        = itp + 1;  %plotting counter
            tp(itp)    = t;        %actual time
            Tp(itp)    = TK;       %actual Temperature
            Pp(itp)    = PGPa;     %actual Pressure
            %get recorded core-T from Zr
            Tcore(itp) = TZK(PGPa,aTiO2,aSiO2,C(1)); %[K]
            %overall grain T from Zr
            TiT        = TZK(PGPa,aTiO2,aSiO2,C);    %[K]
            %plotting on the fly
            if draw_now == 1
                figure(2)
                %Temperature vs. Pressure
                subplot(141)
                plot(T_ar-C2K,P_ar, TK-C2K,PGPa,'o', Tcore-C2K,Pp,'r', Tcore(end)-C2K,Pp(end),'r *')
                %time vs. Temperature
                subplot(142)
                plot(t_ar/myr2sec,T_ar-C2K, t/myr2sec,TK-C2K,'o', tp/myr2sec,Tcore-C2K,'r', tp(end)/myr2sec,Tcore(end)-C2K,'r *')
                %radius vs. concentration
                subplot(143)
                plot(r/L,C)
                %radius vs. Zr-Temperature
                subplot(144)
                plot(r/L,TiT-C2K)
                drawnow
            end
        end
    end
    %update to plot final figures
    itp        = itp + 1;  %plotting counter
    tp(itp)    = t;        %actual time
    Tp(itp)    = TK;       %actual Temperature
    Pp(itp)    = PGPa;     %actual Pressure
    Tcore(itp) = TZK(PGPa,aTiO2,aSiO2,C(1)); %get recorded core-T from Zr [K]
    TiT        = TZK(PGPa,aTiO2,aSiO2,C);    %[K]
    toc
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%3) loading data from precompiled lookup table
%figure plotting is done via activity_data_function
if activity_from_pseudosection == 1
    %nearest neighbour interpolation to get activities from werami data
    [interp_aTi_obj, interp_aSi_obj, interp_aTi, interp_aSi] = activity_data_function(R, T_ar, P_ar, C2K);
    %-----------------------------
    %initialize
    aTiO2 = interp_aTi(1);           %initial aTiO2
    aSiO2 = interp_aSi(1);           %initial aSiO2
    CBC   = ZiT(TK,PGPa,aTiO2,aSiO2);%[ppm]
    C0    = ones(1,nr)*CBC;          %starting with equilibrium Zr-concentration from starting T
    C     = C0;                      %initialize C
    t    = 0;                        %time starts @ 0
    it   = 0;                        %counter starts @ 0
    itp  = 0;                        %counter starts @ 0
    tic
    %solve diffusion equation
    while t < t_tot                  %until time is out
        it = it + 1;                 %increase counter
        t  = t + dt;                 %track time
        %if we overshoot
        if t > t_tot
            dt    = dt - (t_tot-t);  %updated time step
            t     = t_tot;           %final time
            TK    = TKstop;          %final T
            PGpa  = PGPastop;        %final P
            aTiO2 = interp_aTi(end); %final aTiO2
            aSiO2 = interp_aSi(end); %final aSiO2
        end
        D = D0*exp(-Q/(R*TK));       %update T-dependent diffusivity
        %set up matrix for solver
        for i = 2:nr-1                              %loop through matrix
            j = i-1;                                %differnt index for rcR & rcL
            S = (D*dt)/(r(i)^2*dr^2);               %constants
            A(i,i-1) = -S*(rcL(j)^2);               %below diagonal
            A(i,i)   = 1 + S*(rcL(j)^2 + rcR(j)^2); %@ diagonal
            A(i,i+1) = -S*(rcR(j)^2);               %above diagonal
        end
        A(1,1)   = 1;                               %apply BC
        A(nr,nr) = 1;                               %apply BC
        b        = C(:);                            %define concentration vector
        X        = A\b;                             %apply backslash magic
        C        = X(:)';                           %transpose vector
        C(1)     = C(2);                            %zero-flux Neumann BC at center for symmetry
        C(end)   = CBC;                             %set target concentration
        TK       = TK    - CR*dt;                   %update T @ each step
        PGPa     = PGPa  - DPR*dt;                  %update P @ each step
        aTiO2    = interp_aTi_obj(TK,PGPa*1e4);     %interp aTiO2 on the fly
        aSiO2    = interp_aSi_obj(TK,PGPa*1e4);     %interp aSiO2 on the fly
        CBC      = ZiT(TK,PGPa,aTiO2,aSiO2);        %update BC - Hayden equation
        if mod(it,nout) == 1
            itp        = itp + 1;  %plotting counter
            tp(itp)    = t;        %actual time
            Tp(itp)    = TK;       %actual Temperature
            Pp(itp)    = PGPa;     %actual Pressure
            %get recorded core-T from Zr [K]
            Tcore(itp) = TZK(PGPa,aTiO2,aSiO2,C(1)); %[K]
            %overall grain T from Zr
            TiT        = TZK(PGPa,aTiO2,aSiO2,C);    %[K]
            %plotting on the fly
            if draw_now == 1
                figure(2)
                %Temperature vs. Pressure
                subplot(141)
                plot(T_ar-C2K,P_ar, TK-C2K,PGPa,'o', Tcore-C2K,Pp,'r', Tcore(end)-C2K,Pp(end),'r *')
                %time vs. Temperature
                subplot(142)
                plot(t_ar/myr2sec,T_ar-C2K, t/myr2sec,TK-C2K,'o', tp/myr2sec,Tcore-C2K,'r', tp(end)/myr2sec,Tcore(end)-C2K,'r *')
                %radius vs. concentration
                subplot(143)
                plot(r/L,C)
                %radius vs. Zr-Temperature
                subplot(144)
                plot(r/L,TiT-C2K)
                drawnow
            end
        end
    end
    %update to plot final figures
    itp        = itp + 1;  %plotting counter
    tp(itp)    = t;        %actual time
    Tp(itp)    = TK;       %actual Temperature
    Pp(itp)    = PGPa;     %actual Pressure
    Tcore(itp) = TZK(PGPa,aTiO2,aSiO2,C(1)); %final core ZiT Temperature [K]
    TiT        = TZK(PGPa,aTiO2,aSiO2,C);    %final grain ZiT Temperature [K]
    toc
end
%---------------------------------
%%FINAL PLOT
figure(2);
%Temperature vs. Pressure
if activity_from_pseudosection == 0
    sgtitle(['Grain radius = ', num2str(GR), ' µm | Cooling rate = ', ...
        num2str(CR*myr2sec), ' °C/m.y. | Decompression rate = ', num2str((PGPa0-PGPastop)/t_tot*myr2sec), ' GPa/m.y. |', ...
        ' aTiO_2 = ', num2str(aTiO2), ' | aSiO_2 = ', num2str(aSiO2)])
else
    sgtitle(['Grain radius = ', num2str(GR), ' µm | Cooling rate = ', ...
        num2str(CR*myr2sec), ' °C/m.y. | Decompression rate = ', num2str((PGPa0-PGPastop)/t_tot*myr2sec), ' GPa/m.y.'])
end
subplot(141)
plot(T_ar-C2K,P_ar,Tcore-C2K,Pp,'r')
title('Temperature vs. Pressure')
xlabel('Temperature (°C)')
ylabel('Pressure (GPa)')
legend([{'real T'}, {'T in titanite'}],'Location','NW')
grid on
legend boxoff
set(gca, 'FontSize', 14)
%time vs. Temperature
subplot(142)
plot(t_ar/myr2sec,T_ar-C2K, tp/myr2sec,Tcore-C2K)
title('Time vs. Temperature')
xlabel('time (m.y.)')
ylabel('Temperature (°C)')
legend([{'real T'}, {'T in titanite'}],'Location','NE')
grid on
legend boxoff
set(gca, 'FontSize', 14)
%radius vs. concentration
subplot(143)
plot(r/L,C)
title('Norm. radius vs. Zr conc')
xlabel('normalized radius')
ylabel('Zr in titanite (ppm)')
grid on
set(gca, 'FontSize', 14)
%radius vs. Zr-Temperature
subplot(144)
plot(r/L, TiT-C2K)
title('Norm. radius vs. ZiT temperature')
xlabel('normalized radius')
ylabel('Zr-in-titanite Temperature (°C)')
grid on
set(gca, 'FontSize', 14)
%---------------------------------------------------------
%some output in the console
disp(['Initial Zr concentration in the grain is ' num2str(C0(1) ) ' ppm.'])
disp(['Final rim Zr concentration is ' num2str(C(end)) ' ppm.'])
disp(['Final core Zr concentration is ' num2str(C(1)) ' ppm.'])
disp(['Zr-in-titanite temperature at the core is ' num2str(Tcore(itp)-C2K) ' °C.'])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ACTIVITY DATA FUNCTION - GRABS µ FROM DATATABLE AND CALCULATES ACTIVITIES
    function [interp_aTi_obj, interp_aSi_obj, interp_aTi, interp_aSi] = activity_data_function(R, T_ar, P_ar, C2kel)
        werami   = load('data.dat');                                   %load datafile
        nx       = sqrt(length(werami(:,1)));                          %WERAMI makes nx*nx* datagrid
        T        = werami(:,1); T         = reshape(T,nx,nx);          %grab T and reshape as vector
        P        = werami(:,2); P         = reshape(P,nx,nx);          %grab P and reshape
        mu_Ti_sm = werami(:,3); mu_Ti_sm  = reshape(mu_Ti_sm,nx,nx);   %grab µTiO2 from SM89 and reshape
        mu_Ti_ru = werami(:,4); mu_Ti_ru  = reshape(mu_Ti_ru,nx,nx);   %grab µTiO2 from rutile and reshape
        mu_Si_sm = werami(:,5); mu_Si_sm  = reshape(mu_Si_sm,nx,nx);   %grab µSiO2 from SM89 and reshape
        mu_Si_qz = werami(:,6); mu_Si_qz  = reshape(mu_Si_qz,nx,nx);   %grab µSiO2 from quartz and reshape
        %calculate activities
        a_TiO2   = exp((mu_Ti_sm - mu_Ti_ru) ./ (R*T));                %TiO2 activity in SM89
        a_SiO2   = exp((mu_Si_sm - mu_Si_qz) ./ (R*T));                %SiO2 activity in SM89
        %calculate equilibrium Zr in titanite [ppm]
        Zr_cont = 10.^(10.52 - 7708./T - 960 * (P./1e4)./T - log10(a_TiO2) - log10(a_SiO2));
        %WARNINGS ------------------------------------
        if max(T_ar) > max(T)
            disp('The chosen maximum temperature is outside range of activity data. Fix and try again.')
            return
        elseif min(T_ar) < min(T)
            disp('The chosen minimum temperature is outside range of activity data. Fix and try again.')
            return
        elseif max(P_ar*1e4) > max(P)
            disp('The chosen maximum pressure is outside range of activity data. Fix and try again.')
            return
        elseif min(P_ar*1e4) < min(P)
            disp('The chosen minimum pressure is outside range of activity data. Fix and try again.')
            return
        end
        %------------------------------------------------
        %reshape data to use with scatteredInterpolant
        Terp = T(:);
        Perp = P(:);
        aTi  = a_TiO2(:);
        aSi  = a_SiO2(:);
        %create scatteredInterpolant object with 'nearest' interpolation method
        interp_aTi_obj = scatteredInterpolant(Terp,Perp,aTi,'nearest');
        interp_aSi_obj = scatteredInterpolant(Terp,Perp,aSi,'nearest');
        %get the data
        interp_aTi = interp_aTi_obj(T_ar,P_ar*1e4);
        interp_aSi = interp_aSi_obj(T_ar,P_ar*1e4);
        %------------------------------------------------
        %FIGURE PLOTTING
        figure(1)
        subplot(231)
        contourf(T-273,P./1e4,a_TiO2,200,'edgecolor','none'),shading flat
        colorbar
        title('SM89 - aTiO_2')
        xlabel('Temperature (°C)')
        ylabel('Pressure (GPa)')
        axis square
        hold on
        [Cpl,hpl] = contour(T-273,P./1e4,a_TiO2,'k'); %Extract plot properties
        h_line = plot(T_ar-273,P_ar,'LineWidth',2,'Color','r');
        legend([h_line],'PT path')
        legend boxoff
        hold off
        clabel(Cpl,hpl,'LabelSpacing',1200,'FontSize',12)           %Add values on contours
        set(gca, 'FontSize', 14)
        subplot(232)
        contourf(T-273,P./1e4,a_SiO2,200,'edgecolor','none'),shading flat
        colorbar
        title('SM89 - aSiO_2')
        xlabel('Temperature (°C)')
        ylabel('Pressure (GPa)')
        axis square
        hold on
        [Cpl,hpl] = contour(T-273,P./1e4,a_SiO2,'k'); %Extract plot properties
        h_line = plot(T_ar-273,P_ar,'LineWidth',2,'Color','r');
        legend([h_line],'PT path')
        legend boxoff
        hold off
        clabel(Cpl,hpl,'LabelSpacing',1200,'FontSize',12)           %Add values on contours
        set(gca, 'FontSize', 14)
        subplot(233)
        contourf(T-273,P./1e4,log10(Zr_cont),500,'edgecolor','none'),shading flat
        colorbar
        title('SM89 - Zr in titanite log10(ppm)')
        xlabel('Temperature (°C)')
        ylabel('Pressure (GPa)')
        axis square
        hold on
        [Cpl,hpl] = contour(T-273,P./1e4,log10(Zr_cont),'k'); %Extract plot properties
        h_line = plot(T_ar-273,P_ar,'LineWidth',2,'Color','r');
        legend([h_line],'PT path')
        legend boxoff
        hold off
        clabel(Cpl,hpl,'LabelSpacing',1200,'FontSize',12)     %Add values on contours
        set(gca, 'FontSize', 14)
        subplot(234)
        plot(T_ar-C2kel,P_ar)
        title('Chosen P-T path')
        xlabel('Temperature (°C)')
        ylabel('Pressure (GPa)')
        axis square
        grid on
        set(gca, 'FontSize', 14)
        subplot(235)
        plot(T_ar-C2kel, interp_aTi, T_ar-C2kel, interp_aSi)
        title('Activities vs. Temperature')
        xlabel('Temperature (°C)')
        ylabel('Activity')
        axis square
        legend('aTiO_2', 'aSiO_2', Location='southeast')
        legend boxoff
        grid on
        set(gca, 'FontSize', 14)
        subplot(236)
        plot(P_ar, interp_aTi, P_ar, interp_aSi)
        title('Activities vs. Pressure')
        xlabel('Pressure (GPa)')
        ylabel('Activity')
        axis square
        legend('aTiO_2', 'aSiO_2',  Location='southeast')
        legend boxoff
        grid on
        set(gca, 'FontSize', 14)
    end
end