% Script to calulate interface position in the clast by
% (1) using the flux calculated from pumice weight measurements (run
% AnalyzeSaturationML0X first)
% (2) using approximate solution from Pedroso and Domoto (1973)
% (3) using approximate solution from London and Seban (1948)
% created by Kristen in September 2017

%close all
%% Load flux data
%Q=ones(1,1000)*0.15;% [cm^3/s] (constant flux)
Q=Flux(2:300); % [g/s] Load Flux from AnalyzeSaturationML0X
%Q=Q./1000;% flux in kg
%% Clast specific parameters
% Change these depending on the clast %
R=1.38;% [cm] clast radius 
L=0.25;% [cm]initial interface distance from exterior of the clast (estimated) for i=7, L=0.22
I=R-L;% [cm] distance of interface from the center of the clast
phi=0.7; % [dimensionless] porosity
ho=200; %heat transfer coefficient (W/m2/K)
%% Preliminaries 
Tf=100;% [C] fusion temperature
Tw=20;% [C] background temperature
LH=2269;% [J/g] heat of condensation
rho=1;% [g/cm^3] density of liquid water 
%k=2/100;% [Watt/cm/K] thermal conductivity of rock-water mixture
kp = 1.05/100; % [W/mK]  %thermal conductivity of glass (there is a range) (perhaps I should use 2)
kw = 0.58/100; % [W/mK] thermal conductivity of water 
k=(phi*kw+(1-phi)*kp); % volume weighted thermal conductivity

%k=0.28;
kapp=30*k; %apparent thermal conductivity 30 times the effective value
cp=4.185;% [J/g/K] heat capacity of water
e=cp*(Tf-Tw)/LH;% [dimensionless] pertubation physical parameter (Stefan number)
%% Specifics for numerical calculation
dt=1;% timestep in seconds
t=1;% initial time in seconds
i=2;% counter
%%
% (1) Calculate interface location as a function of flux
i=1
C=2.6;
SA=(4*pi*R^2)*phi/C;
while i<200
    t(i+1)=t(i)+dt; %advance time
    I(i+1)=(R-L(i)); %distance of interface from center of the clast
    %SA=4*pi*I(i)^2; %surface area along the spherical interface
    %SA=4*pi*R^2;
    l(i+1)=Q(t(i+1))/SA/phi*dt; %length to advance along the interface
    L(i+1)=L(i)+l(i+1);% advance interface
    V(i+1)=L(i+1)./dt; %calculate interface velocity
    i=i+1;%counter
end

% Make results dimensionless %
tau=k*(Tf-Tw)*t/(rho*LH*R^2);%dimensionless time
RI=I./R; %dimensionless distance

% Plot results %
%figure;
% plot(RI,t,'x');%dimensionless interface distance vs. time.
% ylabel('time (s)');
% xlabel('normalized front position');
%hold on;

 hold on
plot(t,RI,'.')
xlabel('time (s)');
ylabel('normalized front position');
xlim([0 50])
%% (2) Estimate time vs. position from Pedroso and Domoto (1973) %%
%k=kapp;
e=3;
phi=0.12:0.001:5; % made-up variable
% Dimensionless front position %
RSTAR=phi.*(60.*phi.^2+e.*(32.*phi-13))./(60*phi.^2+e.*(22.*phi-3)); 
% Dimensionless time %
TSTAR=((60.*phi.^2.*(1+2.*phi)+e.*(1+2.*phi+120.*phi.^2)).*(1-phi).^2)./(6.*(60.*phi.^2+e));
% Time in seconds% 
TT=TSTAR*0.0138^2/(2.5*10^-7);
TT=TSTAR/(Tf-Tw)*rho*LH*0.0138^2/(2.5*10^-7); %[s] time
%plot Pedroso and Domoto (1973) results
hold on;
plot(TT,RSTAR);
%xlim([0 1]);

% % % Calculate flux from Pedroso and Domoto Results %
% % for i=2:length(RSTAR)
% %     Qtest(i)=0.75*pi*RSTAR(i-1)^3-0.75*pi*RSTAR(i)^3;
% %     dt_test(i)=TSTAR(i)-TSTAR(i-1);
% %     Q2test(i)=Qtest(i)/dt_test(i);
% % end
% % 
% % figure;plot(TSTAR,Q2test,'x');

%% (3) Estimate time vs. position from London and Seban (1948)
r=0.02:0.02:R;
rstar=r./R; % dimensionless radius
Ri=(1./r-1/R)./(4*pi*k);% [K/W] resistance offered by the liquid water
Ro=1/(4*pi*R^2*ho);% [K/W] surface resistance
Rstar=Ri./Ro;%[dimensionless] resistance
tau=1/3.*(1./Rstar-1).*(1-rstar.^3)+1/2.*(1-rstar.^2); %dimensionless time
tt=tau*rho*LH*R^2/((Tf-Tw)*k);
hold on;

plot(tt,rstar);

%plot(tau,rstar);