% Written by Varun Venoor and David Kazmer, April 1, 2021
 
clear all;close all; clc;
bSubplot=false; 
iFile=4;
 
% Define column titles
sAries{1}='Time (seconds)';
sAries{2}='Oscillation strain';
sAries{3}='Oscillation stress';
sAries{4}='Tan(delta)';
sAries{5}='Storage modulus';
sAries{6}='Loss modulus';
sAries{7}='{\eta}* (Pa.s)';
sAries{8}='Oscillation torque, uN.mm';
sAries{9}='Gap';
sAries{10}='Gel content, wt pct';
sAries{11}='Gel inner radius, mm';
sAries{12}='Modeled Torque, uN.mm';
sAries{13}='Modeled Bulk Torque, uN.mm';
sAries{14}='Modeled Edge Torque, uN.mm';
sAries{15}='Bulk complex viscosity';
sAries{16}='Ring complex viscosity';
sGel{1}='Temperature';
sGel{2}='Time, s';
sGel{3}='Gel content, wt pct';
 
% Read File
sFileAries{3}='TPA_255C_Dry_Air_10rad.csv';
sFileAries{4}='TPA_270C_Dry_AIr_10rad.csv';
sFileGel{3}='insoluble gel in (wt.-%)-255C_dry_Air.csv';
sFileGel{4}='insoluble gel in (wt.-%)-270C_dry_Air.csv';
Aries=csvread(sFileAries{iFile},2);
Gel=csvread(sFileGel{iFile},1);
 
Gel(:,2)=60*Gel(:,2);
Gel(:,3)=0.01*Gel(:,3);
 
% Plot complex viscosity and define initial evolution of viscosity
figure(1);clf;
if bSubplot
    subplot(4,1,1);hold on;
else
    hold on;
end
plot(Aries(:,1),Aries(:,7),'k');
ylabel(sAries{7});
set(gca,'yScale','Log')
grid on
 
% Interpolate gel weight percent and plot radial thickness vs time
Aries(:,10)=interp1(Gel(:,2),Gel(:,3),Aries(:,1),'spline');
if bSubplot
    subplot(4,1,2);hold on;
else
    figure(2);clf;
end
plot(Aries(:,1),Aries(:,10),'k.');
ylabel(sAries{10});
grid on
 
% Calculate cross-link inner radius: gel%=(pi*(RO^2-RI^2))/(pi*RO^2), so...
RO=12.5;
Aries(:,11)=(RO^2-Aries(:,10)*RO^2).^0.5;
if bSubplot
    subplot(4,1,3);hold on;
else
    figure(3);clf;
end
plot(Aries(:,1),Aries(:,11),'k.');

% Calculate torque & torque contributions
% First, estimate gel < KGEL%
KGEL= 2.99/100;  % KGEL for 255 0.6799
iEvolve=min(find(Aries(:,10)>KGEL)); 
if bSubplot
    subplot(4,1,1);set(gca,'yScale','Log')
else
    figure(1);
end
EtaBulk=Aries(iEvolve,7);
iBig=find(Aries(:,7)>EtaBulk);
Aries(:,15)=Aries(:,7); Aries(:,7)=EtaBulk; Aries(iBig,15)=EtaBulk;
plot(Aries(:,1),Aries(:,15),'g');
grid on
 
% Now, model the torques removing the strain limiter
Aries(:,12)=max(Aries(:,2))*Aries(:,8)./Aries(:,2);
% Now, compute the torque:viscosity ratio of early data
kT=mean(Aries(1:iEvolve,12))/mean(Aries(1:iEvolve,15))/RO^2;
% Now, compute the torque applied by the bulk viscosity
Aries(:,13)=kT*Aries(:,15).*Aries(:,11).^2;
% Now, compute the torque by subtracting coulmn 13 from 12
Aries(:,14)=Aries(:,12)-Aries(:,13);
if bSubplot
    subplot(4,1,4);hold on;
else
    figure(4);clf;hold on;
end
plot(Aries(:,1),Aries(:,8),'k');
plot(Aries(:,1),Aries(:,12),'b');
plot(Aries(:,1),Aries(:,13),'g');    
plot(Aries(:,1),Aries(:,14),'r');
ylim([0.01*Aries(1,8) Aries(end,12)])
set(gca,'yScale','Log')
ylabel('Torque, uN.mm');
xlabel(sAries{1});
legend('Measured','Unlimited','Bulk','Gel')
grid on
 
% Finally, compute the viscosity of the cross-linked material or gel 
Aries(:,16)=Aries(:,14)./(RO^2-Aries(:,11).^2)/kT;
if bSubplot
    subplot(4,1,1);
else
    figure(1);
end
Aries(1:iEvolve,16)=nan;
plot(Aries(:,1),Aries(:,16),'r');
ylabel(sAries{7});
xlabel(sAries{1});
legend('Experimental Data','Apparent Bulk','Apparent Gel')
opt.Markers = {'o','o','o'};
% title(sTitle(1:end-4))
