%{
Change the file path and name for the Table1 = readtable(... then run the
program.
%}
clear;  % Clears the variables

Table1 = readtable('C:\Users\Makit\Documents\MIT\Senior\Thesis\Data\500mPaTesting4.xlsx', 'sheet', 'Sheet1' ,'range', 'A:B');    % Import the data from Excel
AverageStart = 1;   % Position to start averaging
AverageEnd = 0; % Number of positions from the end to stop averaging
StartTrim = 5;  % Seconds to trim off from start of each speed
EndTrim = 1;    % Seconds to trim off from end of each speed
TestLength = 15;    % Seconds, total time in seconds of each speed testing

% Viscometer Dimensions
VaneRadius = 10;    % mm
VaneLength = 50;    % mm
CupRadius = 15; % mm
LoadCellLever = 22.937; % mm, Distance from center of cup to loadcell contact point

% Unit Conversions
VaneRadius = VaneRadius/1000;   % m
VaneLength = VaneLength/1000;   % m
CupRadius = CupRadius/1000; % m

Array1 = table2array(Table1);   % Convert imported data into an array
RPM1 = Array1(: , 1);   % RPM, make the RPM column vector
LoadCell1 = (-1).*Array1(: , 2); % Make the LoadCell column vector

% Trim 0 values at the end of the dataset
c = length(RPM1);
while RPM1(c) == 0  % Cycles while the values are 0
    c = c - 1;  % Increments to next value
end

RPM1 = RPM1(1:c, 1);    % Trim RPM
LoadCell1 = LoadCell1(1:c, 1);  % Trim Loadcell

% Trim ramping at start of data set
c = 2;
Num = RPM1(1);
while abs(Num-RPM1(c)) >= 0.01 || RPM1(c+1) > 60.1  % Cycles while the speed is changing and larger than 60
    Num = RPM1(c);
    c = c + 1;  % Increments to next value
end

RPM1 = RPM1((c+1):end, 1);  % Trim RPM
LoadCell1 = LoadCell1((c+1):end, 1);    % Trim Loadcell

% Viscosity Caluculations 
Torque1 = (4.672*10^(-8)).*LoadCell1;  % Nm, converts loadcell values to torque
RotationRate1=RPM1.*(2*pi()/60);    % rad/sec, converts RPM to rotation rate
ShearRate1 = RotationRate1.*(2/(1-(VaneRadius/CupRadius)^2));   % 1/s, calculates the shear rate
Sigma1 = Torque1./((2*pi()*VaneRadius^2*VaneLength*(1+2*VaneRadius/(3*VaneLength))));   % Pa, calculates Sigma
Viscosity1 = Sigma1./ShearRate1;    % Pa/s, calculates the viscosity

% Plotting
figure(1)
plot(LoadCell1);
xlabel('Samples')
ylabel('Loadcell Values')

% Trim the start and end of each dataset
c=1;
Num = RPM1(1);
for i = 1:length(RPM1)  % Cycles through RPM
    if RPM1(i) ~= Num   % Finds where the speed changes
        Num = RPM1(i);  % Updates the comparision number
        Split(c) = i;   % Records where the splits are
        c = c + 1;  % Increments the number of splits
    end
end

Split = [1, Split, length(RPM1)];   % Adds the end points to Split
RPM1Split = zeros([],1);    % Zeros the RPMSplit vector
LoadCell1Split = zeros([],1);   % Zeros the LoadCellSplit vector
for i = 1:(length(Split)-1) % Cycles through Split
    SplitRPM = RPM1(Split(i):(Split(i+1))-1, 1);    % Sets SplitRPM to each speed
    SplitLoadCell = LoadCell1(Split(i):(Split(i+1))-1, 1);  % Splits loadcell the same as SplitRPM
    Length = length(SplitRPM);  % Length of split RPM
    SplitRPM = SplitRPM(round(Length*(StartTrim/TestLength)):(Length - round(Length*(EndTrim/TestLength))));    % Trims SplitRPM
    SplitLoadCell = SplitLoadCell(round(Length*(StartTrim/TestLength)):(Length - round(Length*(EndTrim/TestLength)))); % Trims SplitLoadCell
    RPM1Split = [RPM1Split; SplitRPM];  % Recombines each RPM
    LoadCell1Split = [LoadCell1Split; SplitLoadCell];   % Recombines LoadCell the same as RPM
end

% Viscosity Calculations
Torque1 = (4.672*10^(-8)).*LoadCell1Split;  % Nm, converts loadcell values to torque
RotationRate1=RPM1Split.*(2*pi()/60);    % rad/sec, converts RPM to rotation rate
ShearRate1 = RotationRate1.*(2/(1-(VaneRadius/CupRadius)^2));   % 1/s, calculates the shear rate
Sigma1 = Torque1./((2*pi()*VaneRadius^2*VaneLength*(1+2*VaneRadius/(3*VaneLength))));   % Pa, calculates Sigma
Viscosity1 = Sigma1./ShearRate1;    % Pa/s, calculates the viscosity

% Trim torques lower than 10^-5
c = length(RPM1Split);

%Trims low torques and RPMs
while Torque1(c) < 10^(-5) || RPM1Split(c)<0.25  % Cycles while the torque or speed is low
    c = c - 1;  % Increments to next value
end

RPM1Trim = RPM1Split(1:c, 1);  % Trim RPM
LoadCell1Trim = LoadCell1Split(1:c, 1);    % Trim Loadcell

% Viscosity Calculations
Torque1 = (4.672*10^(-8)).*LoadCell1Trim;  % Nm, converts loadcell values to torque
RotationRate1=RPM1Trim.*(2*pi()/60);    % rad/sec, converts RPM to rotation rate
ShearRate1 = RotationRate1.*(2/(1-(VaneRadius/CupRadius)^2));   % 1/s, calculates the shear rate
Sigma1 = Torque1./((2*pi()*VaneRadius^2*VaneLength*(1+2*VaneRadius/(3*VaneLength))));   % Pa, calculates Sigma
Viscosity1 = Sigma1./ShearRate1;    % Pa/s, calculates the viscosity
Viscosity1Trim = Viscosity1;
ShearRate1Trim = ShearRate1;

Split = zeros([],1);
c=1;
Num = RPM1Trim(1);
for i = 1:length(RPM1Trim)  % Cycles through RPM
    if RPM1Trim(i) ~= Num   % Finds where the speed changes
        Num = RPM1Trim(i);  % Updates the comparision number
        Split(c) = i;   % Records where the splits are
        c = c + 1;  % Increments the number of splits
    end
end

Split = [1, Split, length(RPM1Trim)];   % Adds the end points to Split
RPM1Avg = zeros([],1);    % Zeros the RPMSplit vector
LoadCell1Avg = zeros([],1);   % Zeros the LoadCellSplit vector

for i = 1:(length(Split)-1) % Cycles through Split
    RPM1Avg(i) = mean(RPM1Trim(Split(i):Split(i+1)));   % Averages each RPM step
    LoadCell1Avg(i) = mean(LoadCell1Trim(Split(i):Split(i+1))); % Averages each RPM step
    LoadCell1Max(i) = max(LoadCell1Trim(Split(i):Split(i+1)));   % Averages each RPM step
    LoadCell1Min(i) = min(LoadCell1Trim(Split(i):Split(i+1)));   % Averages each RPM step
end

% Viscosity Caluculations 
Torque1 = (4.672*10^(-8)).*LoadCell1Avg;  % Nm, converts loadcell values to torque
RotationRate1=RPM1Avg.*(2*pi()/60);    % rad/sec, converts RPM to rotation rate
ShearRate1 = RotationRate1.*(2/(1-(VaneRadius/CupRadius)^2));   % 1/s, calculates the shear rate
Sigma1 = Torque1./((2*pi()*VaneRadius^2*VaneLength*(1+2*VaneRadius/(3*VaneLength))));   % Pa, calculates Sigma
Viscosity1 = Sigma1./ShearRate1;    % Pa.s, calculates the viscosity
Viscosity1Avg = Viscosity1;
AverageViscosity=mean(Viscosity1)
ShearRate1Avg = ShearRate1;

TorqueMax = (4.672*10^(-8)).*LoadCell1Max;  % Nm, converts loadcell values to torque
RotationRateMax=RPM1Avg.*(2*pi()/60);    % rad/sec, converts RPM to rotation rate
ShearRateMax = RotationRateMax.*(2/(1-(VaneRadius/CupRadius)^2));   % 1/s, calculates the shear rate
SigmaMax = TorqueMax./((2*pi()*VaneRadius^2*VaneLength*(1+2*VaneRadius/(3*VaneLength))));   % Pa, calculates Sigma
ViscosityMax = SigmaMax./ShearRateMax;    % Pa.s, calculates the viscosity

TorqueMin = (4.672*10^(-8)).*LoadCell1Min;  % Nm, converts loadcell values to torque
RotationRateMin=RPM1Avg.*(2*pi()/60);    % rad/sec, converts RPM to rotation rate
ShearRateMin = RotationRateMin.*(2/(1-(VaneRadius/CupRadius)^2));   % 1/s, calculates the shear rate
SigmaMin = TorqueMin./((2*pi()*VaneRadius^2*VaneLength*(1+2*VaneRadius/(3*VaneLength))));   % Pa, calculates Sigma
ViscosityMin = SigmaMin./ShearRateMin;    % Pa.s, calculates the viscosity

Sigmalim = (7*10^-5)/((2*pi()*VaneRadius^2*VaneLength*(1+2*VaneRadius/(3*VaneLength))));   % Pa, calculates min Sigma
Viscositylim = Sigmalim./ShearRate1Avg;    % Pa/s, calculates the min viscosity

minx = 1;
maxx = 30;
lsx = linspace(minx,maxx);

figure(2)
loglog(ShearRate1Avg, Viscosity1Avg, 'rs');
xlabel('Shear Rate (1/s)')
ylabel('Viscosity (Pa.s)')
hold on
loglog(ShearRateMax, ViscosityMax, 'r--');
loglog(lsx, Sigmalim./lsx, 'k')
loglog(ShearRateMin, ViscosityMin, 'r--');
legend({'Our Data', 'Cylinder Data Range', 'Min Torque Limit'},'FontSize',18)
set(gca,'fontsize', 22);
set(gcf,'position',[10,10,700,500])
hold off

figure(3)
loglog(ShearRate1Avg, Sigma1, 'rs');
xlabel('Shear Rate (1/s)')
ylabel('Shear Stress (Pa)')
hold on
f=fit(ShearRate1Avg', Sigma1', 'a+b*x^c')
loglog(ShearRateMax, SigmaMax, 'r--');
loglog([min(ShearRate1Avg), max(ShearRate1Avg)], [Sigmalim, Sigmalim], 'k')
loglog(ShearRateMin, SigmaMin, 'r--');
legend({'Our Data', 'Cylinder Data Range', 'Min Torque Limit'},'FontSize',18)
set(gca,'fontsize', 22);
set(gcf,'position',[10,10,700,500])
hold off
