clear all
close all

%

gdfFloatNumbers = readlines("gdfFloatNumbers.txt");
numGDFFloats=length(gdfFloatNumbers) - 1;
%Note: 5907149, 6990688 and 7902263 were not QC'd
badGDF_SN = [30 31 32 25];
badGDF_UWID = [21441 21479 21485 21546];
badGDF_wmoNumbers = string([2903473 2903455 2903454 2903471]);
assembler = ['BD '; 'MBB'; 'TW '; 'PW?']
failureCycle= [17 2 1 1];

%% Set up the Import Options and import the data

opts = delimitedTextImportOptions("NumVariables", 31);

% Specify range and delimiter
opts.DataLines = [36, Inf];
opts.Delimiter = "\t";

% Specify column names and types
opts.VariableNames = ["Cruise", "Station", "Type", "mondayyr", "hhmm", "LonE", "LatN", "QF", "IceEvRec", "Pressuredbar", "QF1", "TemperatureC", "QF2", "Salinitypss", "QF3", "Sigma_thetakgm3", "QF4", "Depthm", "QF5", "Oxygenmolkg", "QF6", "OxygenSat", "QF7", "Nitratemolkg", "QF8", "Chl_amgm3", "QF9", "b_bp7001m", "QF10", "pHinsituTotal", "QFpHinsituTotal"];
opts.VariableTypes = ["double", "double", "categorical", "datetime", "datetime", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];

% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";

% Specify variable properties
%opts = setvaropts(opts, ["Type", "hhmm"], "EmptyFieldRule", "auto");
opts = setvaropts(opts, ["Type", "hhmm"], "EmptyFieldRule", "missing");
opts = setvaropts(opts, "mondayyr", "InputFormat", "MM/dd/yyyy");
opts = setvaropts(opts, "hhmm", "InputFormat", "HH:mm");
opts = setvaropts(opts, ["Cruise", "Station", "LonE", "LatN", "QF", "IceEvRec", "Pressuredbar", "QF1", "TemperatureC", "QF2", "Salinitypss", "QF3", "Sigma_thetakgm3", "QF4", "Depthm", "QF5", "Oxygenmolkg", "QF6", "OxygenSat", "QF7", "Nitratemolkg", "QF8", "Chl_amgm3", "QF9", "b_bp7001m", "QF10", "pHinsituTotal", "QFpHinsituTotal"], "ThousandsSeparator", ",");

% Import the data
for i=1:numGDFFloats
    fileDir = "DataFiles\";
    fileName(i) = fileDir + gdfFloatNumbers(i) + ".txt";
    dataTable{i} = readtable(fileName(i), opts);

    %get rid of the NaN rows between profiles
    dataTable{i} = rmmissing(dataTable{i});

    %Clean up the time of day column
    dataTable{i}.hhmm=timeofday(dataTable{i}.hhmm);

    %Synthesize time stamps for each sample in the profile
    ascendRate_dbarPerSec = 0.1;
    sampleDateTimes(1) = dataTable{i}.mondayyr(1) + dataTable{i}.hhmm(1);
    sampleSeconds = (dataTable{i}.Pressuredbar(1:end-1)-dataTable{i}.Pressuredbar(2:end))/ascendRate_dbarPerSec;
    sampleSeconds = [0;cumsum(sampleSeconds)];
    dataTable{i}.runningDT =  dataTable{i}.mondayyr + dataTable{i}.hhmm + seconds(sampleSeconds);
end


% Clear temporary variables
clear opts

%

for i=1:length(dataTable)
    numGoodQFs(i) = length(find(dataTable{i}.QFpHinsituTotal == 0 ));
    numMissingQFs(i) = length(find(dataTable{i}.QFpHinsituTotal == 1 ));
    numQuesQFs(i) = length(find(dataTable{i}.QFpHinsituTotal == 4 ));
    numBadQFs(i) = length(find(dataTable{i}.QFpHinsituTotal == 8 ));
    deployedHours(i) = hours(dataTable{i}.runningDT(end) - dataTable{i}.runningDT(1));
    deploymentDate(i) = dataTable{i}.runningDT(1);
    numProfiles(i) = dataTable{i}.Station(end);
    wmoNumbers{i} = string(dataTable{i}.Cruise(1));
    numRows(i) = height(dataTable{i});
end
    
numGDFsDeployed = length(dataTable)
totalDeployedHours = sum(deployedHours)
numFailures = length(badGDF_UWID)
FITRate= (numFailures/totalDeployedHours)*10^9
MTFF = totalDeployedHours/numFailures
totaProfiles = sum(numProfiles)
maxProfiles = max(numProfiles)
maxProfileIndex = find(numProfiles == maxProfiles)
maxProfileWMONumber = wmoNumbers(maxProfileIndex)

projectedNumGDFs = 10;
FITProjectionNoFailures = (numFailures/(totalDeployedHours +(projectedNumGDFs*24*365)))*10^9
FITProjection1Failures = ((numFailures + 1)/(totalDeployedHours+(projectedNumGDFs*24*365)))*10^9
FITRate - FITProjectionNoFailures
FITRate - FITProjection1Failures

projectedNumGDFs = 10;
numFailures = 1
FITRate = (numFailures/totalDeployedHours) * 10^9
FITProjectionNoFailures = (numFailures/(totalDeployedHours + (projectedNumGDFs*24*365)))*10^9
FITProjection1Failures = ((numFailures + 1)/(totalDeployedHours + (projectedNumGDFs*24*365)))*10^9
FITRate - FITProjectionNoFailures
FITRate - FITProjection1Failures


LambdaRWS_BPowerSupplyMTBFrcr9102=444089;
LambdaRWS_BPowerSupplyMTBFrcr9102B=218172;
LambdaRWS_BPowerSupplyMTBFtelcordia25degC=2235743;
LambdaRWS_BPowerSupplyMTBFtelcordia40degC=1063230;
FITLambda = (1 / LambdaRWS_BPowerSupplyMTBFrcr9102) * 10^9

%

close all

 for i=1:length(dataTable)
     xLabels(i) = strcat(string(i), ' wmo', string(wmoNumbers(i)));
     if(ismember(string(wmoNumbers(i)), badGDF_wmoNumbers))
         xLabels(i) = strcat("Bad ", xLabels(i));
     end

%         xLabels(i) = 'BAD ' + string(i) + '  wmo' + string(sortedData{i}.Cruise(1));
%     else
%         xLabels(i) = string(i) + '  wmo' + string(sortedData{i}.Cruise(1));
%     end
 end

bc = bar(xLabels, [numGoodQFs; numBadQFs; numQuesQFs; numMissingQFs]', 'stacked')
legend('Good QFs', 'Missing QFs', 'Ques QFs', 'Bad QFs')
ax = gca
ax.YAxis.Scale = 'log';
ylimits = ax.YLim;
ylimits(1) = 1;
ylim(ylimits);
xtickangle(90)
title 'Quality Flags Per Float';

hold off

%

%Sort the data table by date of first profile
for i=1:length(dataTable)
    firstProfileDates(i) = dataTable{i}.mondayyr(1);
end
firstProfileDates
[sortedDates sortedIndexes] = sort(firstProfileDates);
for i = 1:length(dataTable)
    sortedData{i} = dataTable{sortedIndexes(i)};
end

%

for i=1:length(sortedData)
    numGoodQFs(i) = length(find(sortedData{i}.QFpHinsituTotal == 0 ));
    numMissingQFs(i) = length(find(sortedData{i}.QFpHinsituTotal == 1 ));
    numQuesQFs(i) = length(find(sortedData{i}.QFpHinsituTotal == 4 ));
    numBadQFs(i) = length(find(sortedData{i}.QFpHinsituTotal == 8 ));
    deployedHours(i) = hours(sortedData{i}.runningDT(end) - sortedData{i}.runningDT(1));
    deploymentDate(i) = sortedData{i}.runningDT(1);
    numProfiles(i) = sortedData{i}.Station(end);
    %wmoNumbers{i} = string(gdfFloatNumbers{i}(:,1:7));
    wmoNumbers{i} = string(sortedData{i}.Cruise(1));
    numRows(i) = height(sortedData{i});
end
    
numGDFsDeployed = length(dataTable)
totalDeployedHours = sum(deployedHours)
numFailures = length(badGDF_UWID)
FITRate= (numFailures/totalDeployedHours)*10^9
MTFF = totalDeployedHours/numFailures
totaProfiles = sum(numProfiles)
maxProfiles = max(numProfiles)
maxProfileIndex = find(numProfiles == maxProfiles)
maxProfileWMONumber = wmoNumbers(maxProfileIndex)

projectedNumGDFs = 10;
FITProjectionNoFailures = (numFailures/(totalDeployedHours +(projectedNumGDFs*24*365)))*10^9
FITProjection1Failures = ((numFailures + 1)/(totalDeployedHours+(projectedNumGDFs*24*365)))*10^9
FITRate - FITProjectionNoFailures
FITRate - FITProjection1Failures

numFailures = 1
FITProjectionNoFailures = (numFailures/(totalDeployedHours +(projectedNumGDFs*24*365)))*10^9
FITProjection1Failures = ((numFailures + 1)/(totalDeployedHours+(projectedNumGDFs*24*365)))*10^9
FITRate - FITProjectionNoFailures
FITRate - FITProjection1Failures

LambdaRWS_BPowerSupplyMTBFrcr9102=444089;
LambdaRWS_BPowerSupplyMTBFrcr9102B=218172;
LambdaRWS_BPowerSupplyMTBFtelcordia25degC=2235743;
LambdaRWS_BPowerSupplyMTBFtelcordia40degC=1063230;
FIT1 = (1 / LambdaRWS_BPowerSupplyMTBFrcr9102) * 10^9



%

close all

 for i=1:length(sortedData)
     xLabels(i) = strcat(string(i), ' wmo', string(wmoNumbers(i)));
     if(ismember(string(wmoNumbers(i)), badGDF_wmoNumbers))
         xLabels(i) = strcat("Bad ", xLabels(i));
     end

%         xLabels(i) = 'BAD ' + string(i) + '  wmo' + string(sortedData{i}.Cruise(1));
%     else
%         xLabels(i) = string(i) + '  wmo' + string(sortedData{i}.Cruise(1));
%     end
 end

bc = bar(xLabels, [numGoodQFs; numBadQFs; numQuesQFs; numMissingQFs]', 'stacked')
legend('Good QFs', 'Missing QFs', 'Ques QFs', 'Bad QFs')
ax = gca
ax.YAxis.Scale = 'log';
ylimits = ax.YLim;
ylimits(1) = 1;
ylim(ylimits);
xtickangle(90)
title 'Quality Flags Per Float';

hold off

%

suspectFloats = [13 14 15 16 17 18 26 41 44 45 46];
for i = 1:length(suspectFloats)
    yyaxis right
    plot(sortedData{suspectFloats(i)}.runningDT, sortedData{suspectFloats(i)}.QFpHinsituTotal)
    yyaxis left
    plot(sortedData{suspectFloats(i)}.runningDT, sortedData{suspectFloats(i)}.pHinsituTotal)
    
    title('WMO' + string(wmoNumbers(i)));
end