% plot inflection timescale data for all clasts %
% 1: Cooling timescale
% 2: Initial T
% 3: Spherical equivalent radius
% 4: Clast number
% 5: clast mass (g)
% 6: Spherical equivalent surface area
% 7: Surface area times factor (1.5)
% 8: Mean between surface area and factor
% 9: Difference between mean and spherical equivalent surface area
%10: clast aspect ratio
Data=[8.66	139	1.38	1	7.699	23.92	47.84	35.88	11.96	1.8
9.83	170	1.38	1	7.699	23.92	47.84	35.88	11.96	1.8
12.55	205	1.38	1	7.699	23.92	47.84	35.88	11.96	1.8
14.11	257	1.38	1	7.699	23.92	47.84	35.88	11.96	1.8
34.26	458	1.38	1	7.699	23.92	47.84	35.88	11.96	1.8
35.72	523	1.38	1	7.699	23.92	47.84	35.88	11.96	1.8
37.47	566	1.38	1	7.699	23.92	47.84	35.88	11.96	1.8
41.85	670	1.38	1	7.699	23.92	47.84	35.88	11.96	1.8
8.30	111	2.29	2	38.896	65.87	131.73	98.80	32.93	1.54
16.92	165	2.29	2	38.896	65.87	131.73	98.80	32.93	1.54
26.43	198	2.29	2	38.896	65.87	131.73	98.80	32.93	1.54
35.23	256	2.29	2	38.896	65.87	131.73	98.80	32.93	1.54
82.94	449	2.29	2	38.896	65.87	131.73	98.80	32.93	1.54
90.16	562	2.29	2	38.896	65.87	131.73	98.80	32.93	1.54
97.91	656	2.29	2	38.896	65.87	131.73	98.80	32.93	1.54
31.00	457	1.17	6	5.7	17.19	34.39	25.79	8.60	1.4
20.00	265	1.4	3	7.85	24.62	49.24	36.93	12.31	1.7
34.00	550	1.4	3	7.85	24.62	49.24	36.93	12.31	1.7
100.00	681	2.71	7	63.79	92.24	184.48	138.36	46.12	2.32
35.00	705	1.4	3	7.85	24.62	49.24	36.93	12.31	1.7
208	452	4.3	5	302.2	232.23	464.47	348.35	116.12	1.6
71.00	478	2.7	7	63.79	91.56	183.12	137.34	45.78	2.32
46.67	488.7	2.85	8	60.7	81.65	163.30	122.48	40.83	1.6]; 
%%
hold on;
% Fit a linear line to the data
%figure;
SAlow=Data(:,6); %spherical equivalent surface area
SAhigh=Data(:,6).*Data(:,10);%weight by aspect ratio
SA=(SAlow+SAhigh)./2; %mean surface area
cp=2;

x=(Data(:,2)-20).*Data(:,5).*cp./SA; % plot expected with mean surface area
xhigh=(Data(:,2)-20).*Data(:,5).*cp./SAhigh;
err=abs(xhigh-x);%estimate error 
y=Data(:,1); %time

errorbar(x,y,err,'horizontal','k.');%,'MarkerType', Data(:,6));
hold on

% perform a least squares regression on the data
q=x\y;
q_save=1/q %determine the 
% plot model fit
x2=1:1200;
ycalc=q*x2; %calculated values from regression
hold on 
plot(x2,ycalc,'k')

xlabel('(T_i-T_{\infty})c_pMass_p/surface area (J cm^-^2)');
ylabel('Low-to-high flux transition time, t^* (s)');
set(gca,'FontSize',20);
ylim([0 300]);

%%
hold on
figure
% plot model fit with markers
mark={'<','d','h','>','s','^','o','x','+'};
cp=2;%J/gK
for i=1:length(Data(:,1))
    hold on
    errorbar(x(i),y(i),err(i),'horizontal','Marker',mark{Data(i,4)},'Color','k','MarkerSize',10,'MarkerFaceColor','k');
    hold on
   
end
hold on 
plot(x2,ycalc,'k')
set(gca,'FontSize',20);
ylim([0 300]);
xlabel('(T_i-T_{\infty})c_pMass_p/surface area (J cm^-^2)');
ylabel('Low-to-high flux transition time, t^* (s)');
 
%%
figure;
SA=4*pi.*Data(:,3).^2;
SA = Data(:,8);
err= Data(:,9);
x=(Data(:,2)-100).*Data(:,5).*cp./SA;
y=Data(:,1);
for i=1:length(Data(:,4))
MARK{i}=mark{Data(i,4)};
end
errorbar(x,y,err,'horizontal','MarkerType',MARK);%,'MarkerType', Data(:,6));
hold on


xlabel('(T_i-T_{\infty})c_pMass_p/surface area (J cm^-^2)');
ylabel('Low-to-high flux transition time, t^* (s)');
set(gca,'FontSize',20);
ylim([0 300]);

%% Now validate liquid flux estimate

q=7.5;
mark2={'ko','ks','k.','k^','kd'};
% clast numbers
% 1:  ML01
% 2:  ML02
% 3:  ML03
% 4:  Pum06
% 5:  SM16
ClastNum=[ 1 2 3 6 9];% label clasts
Ql=[0.13	0.30	0.14	0.0785	0.09032];
phi=[0.69 0.67 0.70 0.63 0.73];
SA=[33.48 83.6549 33.2370 20.6280 20.6280]; %These are made to be consistent with the surface areas used above
v=Ql./(SA.*phi);%velocity into the clast in cm/s
L=2260;
Qpred=q./(phi*L).*SA;
Qpred=q.*SA./L; %modified by Kristen in February 2018 not to include porosity here
figure;
%subplot(1,2,2)
set(gca,'FontSize',20)
for i=1:length(phi)
    hold on
    plot(Ql(i),Qpred(i),mark{ClastNum(i)},'MarkerSize',14,'Color','k','MarkerFaceColor','k');
end
xlabel('observed average liquid flux (g/s)');
ylabel('predicted average liquid flux (g/s)');
hold on
plot([0 0.35],[0 0.35],'k')