%PROGRAM WITH TWO MOVING BOUNDARY CONDITIONS TO COMPUTE ZIRCON
%DISSOLUTION AND GROWTH WITHIN USER_SPECIFIED CELLS OF MELTS
%See Bindeman I.N. and Melnik O.E.(2015), "Zircon survival, rebirth…" Journal of Petrology

function  ZirconDif_movBC( varargin )
close all
global  n R A B D F alpha1 beta Xs UCR ZrPl
global Dplag Temp tyear XH2O Tsolidus Csupsat V Dscale tscale length Rd W t
fig_setup; %setting up figure for the plots. SEPARATE FILE, KEEP IN THE SAME DIRECTORY AS THIS PROGRAM!


%Parameters_____________________________________________________
ZirconRadius=40e-4; %zircon radius in cm, can be changed by the user
Temp=TemperatureHistory(0); % Initial temperature, function T(in K) =TemperatureHistory(time) defines temperature evolutuon with time  
XH2O=2.; %initial water content in melt, required for diffusion coefficient simulations.
U=0; %arbitrary undersaturation set up by user, - growth, + dissol
Tsolidus=400+273; %arbitrary solidus temperature for phase diagram used
Csupsat=3; %ppm supersaturation to cause nucleation of a new crystal upon cooling
tfin=tyear*10000; %total time in years of the process
nt=600;  %number of timesteps, increase for larger tfin if instbility occurs
Dplag=0.1; %Distribution coefficient for main minerals (0 -fully incompatoble)
UCR=1; %Critical concentration for microZircon nucleation on major minerals
Csat=ZrSaturation(Temp); %saturation equation defined below
Cmelt=Csat-U; %initial melt Zr concentration
%CbulkZr = 300; %bulk rock ppm, set up by user
%length = 100*ZirconRadius/(CbulkZr-Cmelt)^0.333; % computed cell size given size of zircon and CbulkZr
length=0.1 ; %initial user defined radius of the cell (cm) regardless of Cbulk and zircon radius inside; comment out if cell size to be computed

%END:Parameters_____________________________________________________

%SCALING-----------------
Dscale=DiffusionCoefficient(Temp,XH2O);
n=round(length*2500*8); %number of points in radial mesh. Can be chaged by user depends on desired accuracy
tscale=length^2/Dscale; %dimentionless scale for the time

tfin_d=tfin/tscale;
dt=tfin_d/(nt-1); %Timestep calculation
t=0;
Xs=ZirconRadius/length; %initial crystal boundary coordinate
Rd=1;
W=0; %W growth of major minerals rate
V=0; %V - dissolution rate in cm/s
%END:SCALING-----------------

matrixes; % storage matrixes, separate file, KEEP IN THE SAME DIRECTORY AS THIS PROGRAM!

%Initial Conditions
C0(1:n)=Cmelt; %U - Zr undersaturation in ppm
pause(1e-5)
%MAIn LOOP in time _______________________
for i=1:nt
  C=progonka(C0,dt);
  C0=C;
  CC(i,:)=C(:);
  xpl=Xplag(Temp);
  Rd=((1-xpl)+xpl*Xs^3)^(1/3);
  rr=R*(Rd-Xs)+Xs;
  set(f,'CurrentAxes',ha1)
  hold on
  plot (rr*1e4*length,C,'r-') %Linear R scale in um
  pause(1e-10);
  t=t+dt;
  Xs=max(0.,Xs-V*dt);
  XXs(i)=Xs*1e4*length; %zircon radius in um
  RRd(i)=Rd*1e4*length; %melt cell radius in um
  VV(i)=-V*length/tscale; %*length/tscale % array of dissolution rate
  tt(i)=t/tyear*tscale;
  UU(i)=C(1); %-C(n)
  Tsave(i)=Temp-273;
  ZrPls(i)=ZrPl;
  Xp_sav(i)=Xplag(Temp);
  Temp=TemperatureHistory(t);

end
%END:MAIn LOOP in time _______________________
%Saving results
Colheader={'Time','Velocity','Zr Radius','Cell Radius','Temperature','Saturation','Plag content','Zr in Plag'};
csvwrite_with_headers('time-evolution.csv',[tt,VV,XXs,RRd,Tsave,UU,Xp_sav(1:nt),ZrPls],Colheader)

fid=fopen('cons.bin','w+');
fwrite(fid,n,'int');
fwrite(fid,nt,'int');
fwrite(fid,R,'single');
fwrite(fid,tt,'single');
fwrite(fid,CC,'single');
fclose(fid);
%END:Saving results

%Plotting data ____________________________

set(f,'CurrentAxes',ha)
hold on
xlim([0 max(RRd)])
ylim([0 tfin/tyear])
zlim([0 max(max(CC))])
for i=1:5:nt
    XXX(1:n)=R(1:n)*(RRd(i)-XXs(i))+XXs(i);
    YYY(1:n)=i*dt*tscale/tyear;
    hold on
    plot3(XXX,YYY,CC(i,:),'b-');
end
plot3(XXs,tt,CC(:,1),'g-');
aaa=zeros(nt,1);
plot3(RRd,tt,aaa,'r-');
plot3(XXs,tt,aaa,'r-');
plot3(RRd,tt,CC(:,n),'g-');
grid on
view(-45,30)
hold off

ylabel('time in years')
xlabel('radial distance, um')
zlabel('Zr concentration')

set(f,'CurrentAxes',ha2)
hold on
plot(tt(1:i-1), XXs(1:i-1),'b-')
set(f,'CurrentAxes',ha3)
hold on
plot(tt(1:i-1), VV(1:i-1),'b-')
set(f,'CurrentAxes',ha4)
hold on
plot(tt(1:i-1), Tsave(1:i-1),'b-')
set(f,'CurrentAxes',ha5)
hold on
plot(tt(1:i-1), Xplag(Tsave(1:i-1)+273),'b-')
%END:Plotting data __Various Figure options are commented out__________________________

%figure()
%plot (R(1:n),CC(i,:),'b-') %Linear R scale in um
% figure()
% fid=fopen('cons.bin','r');
% R=fread(fid,n+1,'single');
% tt=fread(fid,nt,'single');
% CC=reshape(fread(fid,n*nt,'single'),nt,n);
% fclose(fid);
% mesh(R(1:n),tt,CC)
end
function C=progonka(C0,dt)
global  n R  Rd Dplag ZrPl
global A B D F alpha1 beta Xs Temp XH2O Tsolidus V W Csupsat Dscale t UCR

Dif=DiffusionCoefficient(Temp,XH2O)/Dscale; %see below Diff Coeff dependednt on water and T in cm2/s
Csat=ZrSaturation(Temp);
Dflux=Dif*(C0(2) - C0(1))/(R(2)-R(1))/(Rd-Xs);
xpl=Xplag(Temp+1e-5);
dxp=(xpl-Xplag(Temp-1e-5))/2e-5;
dTdt=(TemperatureHistory(t+dt)-TemperatureHistory(t))/dt;
W=(-(1-Xs^3)*dxp*dTdt+3*Xs^2*xpl*V)/3/Rd^2;

V=Dflux/(C0(1)-490000.); %497871 ppm Zr in pure Zr
if(Xs<1e-5) && (C0(1)-Csat > Csupsat) %nucleation conditions if Zr radius decreased to 0 during dissolution
    Xs=1.005e-4;
    C0(1)=Csat;
    V=Dif*Csupsat/Xs/(C0(1)-490000.);
end
if(Temp<Tsolidus)
Dif=0;
end

%boundary conditions for i=1 and n
if(Xs>1e-5)
D(1)=1;
B(1)=0;
F(1)=Csat;
else
D(1)=1;
B(1)=-1;
F(1)=0;
V=0;
end
kpl=0.; %partition coefficient of Zr into major phases; change to >0 if they take up Zr or generate zircon microinclusions
if (C0(n)-Csat >UCR) && (W<0.)
    kpl=Dplag;
end
ZrPl=kpl*C0(n);
D(n)=-Dif-W*(R(n)-R(n-1))*(Rd-Xs)*(1-kpl);
A(n)=Dif;
F(n)=0;

%Coefficients for Thomas method
s=Xs;
for i=2:n-1
psi1=R(i-1);
psi2=R(i);
psi3=R(i+1);
t1 = (Dif * dt);
t5 = (psi1 * Rd - psi1 * s + s) ^ 2;
t6 = psi2 - psi1;
t8 = (t5 / t6);
t12 = (Rd * psi2);
t14 = ((-psi2 + 1) * s + t12) ^ 2;
t15 = Rd - s;
t20 = (-W + V) * psi2 - V;
A(i) = -t14 * t15 * dt * psi2 * t20 - t1 * t8;
t25 = (-psi2 * s + s + t12) ^ 2;
t28 = t25 / (psi3 - psi2);
B(i) = -t1 * t28;
t32 = -t15;
t33 = t32 ^ 2;
t34 = -t6;
t38 = (t32 * psi2 - s) ^ 2;
D(i) = -t1 * (-t28 - t8) - t33 * t34 * t38 - t20 * psi2 * dt * t38 * t32;
t44 = t15 ^ 2;
t48 = (t15 * psi2 + s) ^ 2;
F(i) = -t34 * t44 * t48 * C0(i);
end

%Forward Thomas path
alpha1(2)=-B(1)/D(1);
beta(2)=F(1)/D(1);
for i=2:n-1
    alpha1(i+1)=-B(i)/(A(i)*alpha1(i)+D(i));
    beta(i+1)=(F(i)-A(i)*beta(i))/(A(i)*alpha1(i)+D(i));
end
%Backward Thomas path
C(n)=(F(n)-A(n)*beta(n))/(D(n)+A(n)*alpha1(n));
for i=n-1:-1:1
    C(i)=C(i+1)*alpha1(i+1)+beta(i+1);
end

end

function Dif=DiffusionCoefficient(T,x)  %defining Zr diffusion coefficients in melts as f(T,X2O)
theta=1000/T;
lnD=-(11.4*x+3.13)/(0.84*x+1)-(21.4*x+47)/(1.06*x+1)*theta;%best fit of Zr Diff coefficients (several workers) and WH83 dependence on XH2O
Dif=exp(lnD)*1e4;% in cm2/s
%Dif=0.1*exp(-235980/8.31/T); %Watson 96, Eq 2 for 3wt% water in cm2/s to check his results; switch off for our zircon
% These diffusion coefficient parametrizations (this work) as f(T, H2O) can
% be used for monazite and apatite
%for Monazite:
%logDif=(-1.92422063 *x- 10.91566931)/(x*2.16346615610509+2.309291627128083)* (theta - 5.86)-8;% fitting Dphos as a function of water
%Dif=10^(logDif);% in cm2/s, last two rows are for MONAZITE, switch off for Apatite:
%lnDif=(x*9.46939482772897101-181.512270550850587)/(x+2.37368068088965156)*(theta-0.54)-26.7;% fitting Dphos as a function of water
%Dif=exp(lnDif)*1e4;% in cm2/s

end

function Csat=ZrSaturation(T) %defining Zr saturation conditions
%Csat=4.414e7/exp(13352/T)/2; %Watson 96, Eq 1, in ppm Zr for checking. (divide by 1),or mol Zr (divide by 2)
%Mfactor = 0.0000048*(T)^2-0.0083626*(T)+4.8484463; % empirical relations from magma Fig.
%differentiation calc (file  M_factorsforOleg.xlsx
Mfactor=1.3;
Csat=490000/exp(10108/T-1.16*(Mfactor-1)-1.48); % Boehnkeetal2013ChemGeol351,324 assuming 490,000 ppm in Zircon
%Csat=490000/(exp(12900/T-0.85*(Mfactor-1)-3.80));% Watson and Harrison 1983

% for Monazite:
% H2O = 1wt% use below expression (Table 4, Rapp Watson 86)
%Csat=0.0000190*exp(0.0143872*T);
%Csat=600000/(exp(-0.0144*T+24.177));
% H2O = 6wt% use below expression (Table 4, Rapp Watson 86)
%Csat=0.00012*exp(0.01352*T);
%Csat=600000/(exp(-0.0135*T+22.296));


%for Apatite:
%SiO2=0.68;
%Csat=430000 / exp((-4800 + 26400 * SiO2) / T + 3.10 - 12.4 * SiO2);%Harrison Watson 84

end

function Th=TemperatureHistory(t)
%Th=950+273; % constnant temperature, remove for proper T history
global  tscale tyear
tt=0;
if(t>0)
    tt=t*tscale;
end
Td=950+273; %T of the dike, set up by user
T0=650+273; %T of country rocks
%Th=max(T0,Td-0.1*tt/tyear); % temperature decrease from Td to T0 by 0.1 deg per year
%Th=min(Td,T0+0.1*tt/tyear); % temperature increase from T0 to Td by 0.1 deg per year

%Dyke conductive heating/cooling, coordinate x is from the center of the
%dyke, after x>h, it is in the country rocks
h=250; %half thickness of the dike
Dist=10; %distance from the dyke-rock interface, negative values - inside the dyke
x=h+Dist;% coordinate for T evolution
k=0.6e-6; %temeparature conduction coeff in m^2/s
Th=(Td - T0) * (erf((h - x) * (k * tt)^(-0.5)/2)/2 + erf((h + x) * (k *tt)^(-0.5)/2)/2) + T0; 
%conductive slab analytical solution from Crank (1975)

end

function Xpl=Xplag(T); %phase diagram
x=T-273;
Y = (1.13656E-07*x^4-3.8926E-04*x^3+0.49848*x^2-283.186*x+60299.573)/100;
%liquidus -solidus relation with eutectic based on experiments from Piwinski (1968)
Xpl=min(0.99,max(0.,Y)); %there is always 1%melt remains around zircon
Xpl=0; %- no crystalization, no motion of outer boundary
end