%Calculate the mean temperature for every cycle
for i = 1:width(cycle)
    cycle(1,i).meanT_degC = mean([cycle(1,i).tempDegCDown; cycle(1,i).tempDegCDown]);
end

%

%Correct Vrs for change in HCl pH due to pressure and temperature
% calculate k2
TC = csvData{1}.ThermistorTemp;
pres = csvData{1}.PressuredBar;
Vrs = csvData{1}.VRSVOLTS;

meangamma = -0.000113 .* TC + 0.9067;
delVo = 16.38 + 0.0896 * TC - 0.001264 .*TC .^2;
TK = TC + 273.15;

% activity coefficients
Vct = -8.314462.* TK./96487 .* log((meangamma.*0.01).^2);
Vcp = -delVo./23061.*0.0242.*(pres-10)/10*1/1.01;

% Ktp is the calibration coefficient after chemical component has been removed
csvData{1}.VrsCorrected = Vrs - Vct - Vcp;
csvData{1}.VrsDiff = Vrs - csvData{1}.VrsCorrected;

figure
hold on
yyaxis left
plot(Vct)
ylabel('Vct Temperature Correction')
yyaxis right
plot(Vcp)
ylabel('Vcp Pressure Correction')
title('Tempreature and Pressure Correction')

%

%plot all data
cmap = colormap(jet(500));

x = csvData{1}.LabviewDateTime';
y = csvData{1}.PressuredBar' ;
z = csvData{1}.VrsCorrected';

lastIndex = 17045
x = x(1:lastIndex);
y = y(1:lastIndex);
z = z(1:lastIndex);
colorIndex = colorIndex(1:lastIndex);

figure
hold on
plot3(x, y, z, 'Color', [.75 .75 .71]);

scatter3(x, y, z, 1, cmap(colorIndex,:));

zLimitIncr = 0.002
set(gca, 'ydir', 'reverse');
zLimits = zlim;
set(gca, 'ztick', [zLimits(1):zLimitIncr:zLimits(2)]);
view(3);
datetick('x', 'mmmdd HH:MM');
grid on
colorbar('Ticks', [0, 1], 'TickLabels', {'0 degC', '50 degC'});
title(" All Data From File " + i + ": " + csvFileName, 'Interpreter', 'none')
xlabel('Date Time')
ylabel('Pressure (dBar)')
zlabel('Vrs - (Vcp - Vct) (volts)')
view(3)

%

%Closer look at low pressure data
indexes = find(csvData{1}.PressuredBar < 200);
figure
yst
old on
plot(csvData{1}.PressuredBar(indexes), csvData{1}.VRSVOLTS(indexes))
plot(csvData{1}.PressuredBar(indexes), csvData{1}.VrsCorrected(indexes))
legend('Raw Vrs', 'Vrs Corrected for Chemistry')

%

%Plot check to look for data that shouldn't be part of the ramps (hold data
%or temperature transitions or ...
figure
hold on
for i = 1:width(cycle)
    plot(cycle(1,i).pressureUp, cycle(1,i).tempDegCUp)
    plot(cycle(1,i).pressureDown, cycle(1,i).tempDegCDown)
end



%

%Plot Vrs v time and TempDegC v time for hold periods
maxPHoldIndexes = find(csvData{1}.PressuredBar > maxPressure_psi);
maxHoldVrsValues= csvData{1}.VRSVOLTS(maxPHoldIndexes);
maxHoldTempValues_degC = csvData{1}.ThermistorTemp(maxPHoldIndexes);
maxHoldTimes = csvData{1}.LabviewDateTime(maxPHoldIndexes);
figure
yyaxis left
plot(maxHoldTimes, maxHoldVrsValues, '-+');
yyaxis right
plot(maxHoldTimes, maxHoldTempValues_degC, '-+');
datetick('x', 'mmmdd HH:MM');

%

%Caclulate dTempDegC/dSeconds for pressure hold periods
temps_degC = [20 10 5 2 20];
diffs = diff(maxHoldTempValues_degC);
endIndexes = find(abs(diffs) > 1);
startIndex = 1;
figure
hold on
for i = 1:length(temps_degC)
    if( i == length(temps_degC))
        endIndex = length(maxHoldTempValues_degC);
    else
        endIndex = endIndexes(i) - 1;
    end
    meanTemp_degC(i) = mean(maxHoldTempValues_degC(startIndex:endIndex));
    y = maxHoldVrsValues(startIndex:endIndex);
    meanVrs(i) = mean(y);
    if(i > 1)
        K2(i-1) = (meanVrs(i) - meanVrs(i-1)) / (meanTemp_degC(i) - meanTemp_degC(i-1));
    end

    [x, origUnit] = time2num(maxHoldTimes(startIndex:endIndex), "seconds");
    xDateTimes = maxHoldTimes(startIndex:endIndex);
    polys(i,:)= polyfit(x, y, 1);
    yFit = polyval(polys(i,:), x);

    plot(xDateTimes, y)
    plot(xDateTimes, yFit)
    startIndex = endIndex + 2;
end
datetick('x', "mm/dd HH:MM");



%

lineColors = ["red"; "green"; "blue"; "yellow"; "magenta"; "cyan"; "black"];
yLowLim = min(csvData{1}.VRSVOLTS);
yHighLim = max(csvData{1}.VRSVOLTS);
colorIndex = 1;
polyOrder = 3;
xUp = [];
yUp = [];
xDown = [];
yDown = [];
figure
hold on
for i = 2:5
    xUp = [xUp; cycle(1,i).pressureUp];
    yUp = [yUp; cycle(1,i).VrsUp];
    xDown = [xDown; cycle(1,i).pressureDown];
    yDown = [yDown; cycle(1,i).VrsDown];
    plot(cycle(1,i).pressureUp, cycle(1,i).VrsUp, '-+', 'Color', lineColors(colorIndex))
    plot(cycle(1,i).pressureDown, cycle(1,i).VrsDown, '-o', 'Color', lineColors(colorIndex))
    %ylim([yLowLim, yHighLim]);
    colorIndex = colorIndex + 1;
end
[polyPUp, polySUp, muUp] = polyfit(xUp, yUp, polyOrder);
[VrsFitUp, polyDeltaUp] = polyval(polyPUp, xUp, polySUp, muUp);
[polyPDown, polySDown, muDown] = polyfit(xDown, yDown, polyOrder);
[VrsFitDown, polyDeltaDown] = polyval(polyPDown, xDown, polySDown, muDown);
residualsUp = yUp - VrsFitUp;
residualsUpMean = mean(abs(residualsUp))
residualsUpMax = max(abs(residualsUp))
residualsDown = yDown - VrsFitDown;
residualsDownMean = mean(abs(residualsDown))
residualsDownMax = max(abs(residualsDown))

plot(xUp, VrsFitUp, 'LineWidth', 2, 'Color', 'black');
plot(xDown, VrsFitDown, 'LineWidth', 2, 'Color', 'black');

%hyst = VrsFitUp - VrsFitDown;

%

xUp = [];
yUp = [];
xDown = [];
yDown = [];

figure
hold on
colorIndex = 1;
for i = 6:10
    xUp = [xUp; cycle(1,i).pressureUp];
    yUp = [yUp; cycle(1,i).VrsUp];
    xDown = [xDown; cycle(1,i).pressureDown];
    yDown = [yDown; cycle(1,i).VrsDown];
    plot(cycle(1,i).pressureUp, cycle(1,i).VrsUp, '-+', 'Color', lineColors(colorIndex))
    plot(cycle(1,i).pressureDown, cycle(1,i).VrsDown, '-o', 'Color', lineColors(colorIndex))
    %ylim([yLowLim, yHighLim]);
    colorIndex = colorIndex + 1;
end
[polyPUp, polySUp, muUp] = polyfit(xUp, yUp, polyOrder);
[VrsFitUp, polyDeltaUp] = polyval(polyPUp, xUp, polySUp, muUp);
[polyPDown, polySDown, muDown] = polyfit(xDown, yDown, polyOrder);
[VrsFitDown, polyDeltaDown] = polyval(polyPDown, xDown, polySDown, muDown);
residualsUp = yUp - VrsFitUp;
residualsUpMean = mean(abs(residualsUp))
residualsUpMax = max(abs(residualsUp))
residualsDown = yDown - VrsFitDown;
residualsDownMean = mean(abs(residualsDown))
residualsDownMax = max(abs(residualsDown))

plot(xUp, VrsFitUp, 'LineWidth', 2, 'Color', 'black');
plot(xDown, VrsFitDown, 'LineWidth', 2, 'Color', 'black');

xUp = [];
yUp = [];
xDown = [];
yDown = [];

figure
hold on
colorIndex = 1;
for i = 11:15
    xUp = [xUp; cycle(1,i).pressureUp];
    yUp = [yUp; cycle(1,i).VrsUp];
    xDown = [xDown; cycle(1,i).pressureDown];
    yDown = [yDown; cycle(1,i).VrsDown];
    plot(cycle(1,i).pressureUp, cycle(1,i).VrsUp, '-+', 'Color', lineColors(colorIndex))
    plot(cycle(1,i).pressureDown, cycle(1,i).VrsDown, '-o', 'Color', lineColors(colorIndex))
    %ylim([yLowLim, yHighLim]);
    colorIndex = colorIndex + 1;
end
[polyPUp, polySUp, muUp] = polyfit(xUp, yUp, polyOrder);
[VrsFitUp, polyDeltaUp] = polyval(polyPUp, xUp, polySUp, muUp);
[polyPDown, polySDown, muDown] = polyfit(xDown, yDown, polyOrder);
[VrsFitDown, polyDeltaDown] = polyval(polyPDown, xDown, polySDown, muDown);
residualsUp = yUp - VrsFitUp;
residualsUpMean = mean(abs(residualsUp))
residualsUpMax = max(abs(residualsUp))
residualsDown = yDown - VrsFitDown;
residualsDownMean = mean(abs(residualsDown))
residualsDownMax = max(abs(residualsDown))

plot(xUp, VrsFitUp, 'LineWidth', 2, 'Color', 'black');
plot(xDown, VrsFitDown, 'LineWidth', 2, 'Color', 'black');

xUp = [];
yUp = [];
xDown = [];
yDown = [];

figure
hold on
colorIndex = 1;
for i = 16:20
    xUp = [xUp; cycle(1,i).pressureUp];
    yUp = [yUp; cycle(1,i).VrsUp];
    xDown = [xDown; cycle(1,i).pressureDown];
    yDown = [yDown; cycle(1,i).VrsDown];
    plot(cycle(1,i).pressureUp, cycle(1,i).VrsUp, '-+', 'Color', lineColors(colorIndex))
    plot(cycle(1,i).pressureDown, cycle(1,i).VrsDown, '-o', 'Color', lineColors(colorIndex))
    %ylim([yLowLim, yHighLim]);
    colorIndex = colorIndex + 1;
end
[polyPUp, polySUp, muUp] = polyfit(xUp, yUp, polyOrder);
[VrsFitUp, polyDeltaUp] = polyval(polyPUp, xUp, polySUp, muUp);
[polyPDown, polySDown, muDown] = polyfit(xDown, yDown, polyOrder);
[VrsFitDown, polyDeltaDown] = polyval(polyPDown, xDown, polySDown, muDown);
residualsUp = yUp - VrsFitUp;
residualsUpMean = mean(abs(residualsUp))
residualsUpMax = max(abs(residualsUp))
residualsDown = yDown - VrsFitDown;
residualsDownMean = mean(abs(residualsDown))
residualsDownMax = max(abs(residualsDown))

plot(xUp, VrsFitUp, 'LineWidth', 2, 'Color', 'black');
plot(xDown, VrsFitDown, 'LineWidth', 2, 'Color', 'black');

xUp = [];
yUp = [];
xDown = [];
yDown = [];

figure
hold on
colorIndex = 1;
for i = 21:23
    xUp = [xUp; cycle(1,i).pressureUp];
    yUp = [yUp; cycle(1,i).VrsUp];
    xDown = [xDown; cycle(1,i).pressureDown];
    yDown = [yDown; cycle(1,i).VrsDown];
    plot(cycle(1,i).pressureUp, cycle(1,i).VrsUp, '-+', 'Color', lineColors(colorIndex))
    plot(cycle(1,i).pressureDown, cycle(1,i).VrsDown, '-o', 'Color', lineColors(colorIndex))
    %ylim([yLowLim, yHighLim]);
    colorIndex = colorIndex + 1;
end
[polyPUp, polySUp, muUp] = polyfit(xUp, yUp, polyOrder);
[VrsFitUp, polyDeltaUp] = polyval(polyPUp, xUp, polySUp, muUp);
[polyPDown, polySDown, muDown] = polyfit(xDown, yDown, polyOrder);
[VrsFitDown, polyDeltaDown] = polyval(polyPDown, xDown, polySDown, muDown);
residualsUp = yUp - VrsFitUp;
residualsUpMean = mean(abs(residualsUp))
residualsUpMax = max(abs(residualsUp))
residualsDown = yDown - VrsFitDown;
residualsDownMean = mean(abs(residualsDown))
residualsDownMax = max(abs(residualsDown))

plot(xUp, VrsFitUp, 'LineWidth', 2, 'Color', 'black');
plot(xDown, VrsFitDown, 'LineWidth', 2, 'Color', 'black');