clear all
close all

dirName = 'X:\CPF\Data\CPF-P2\Paragon Tests\2017Mar20-24\Engr Logs\';

dirData = dir(dirName);
fileDatenums = dirData.date
[sortedTime, i] = sort([dirData.datenum]);
sortedDirData = dirData(i);

dateIndex = 2;
timeIndex = 3;
stateIndex = 4;
commentIndex = 5;
pressureIndex = 6;
bellowsPosIndex = 7;
batteryVoltsIndex = 8;
batteryAmpsIndex = 9;
ctdIndex = 11;
rawPressureIndex = 12;
gpsIndex = 38;
iridiumIndex = 39;

pressureFig = figure;
title('Pressure vs Time');
hold on;
% bellowsPosFig = figure;
% title('Bellows Position vs Time');
% hold on;
powerFig = figure;
title('Power vs Time');
hold on;

numSatsIndex = 1;
numSBDIXIndex = 1;
profileStartTimeIndex = 1;
glideOnIndex = 1;
glideOffIndex = 1;
errorIndex = 1;
bpLimitErrorIndex = 1;
processedFileNum = 0;
plotFilename = true;

for fileNum = 4:4
    if (length(sortedDirData(fileNum).name) >= 4)
        if ( (strcmp(sortedDirData(fileNum).name(end-3:end), '.eng')) && (sortedDirData(fileNum).bytes > 100000) )
            processedFileNum = processedFileNum + 1;
            fname = strtrim(sortedDirData(fileNum).name);
            fprintf('%s',[datestr(now, 'HH:MM:SS ') 'Importing file: ' [dirName sortedDirData(fileNum).name] ' ...']);
            fprintf(1, '\r');
            
            fid = fopen([dirName sortedDirData(fileNum).name]);
            
            headerLine = fgetl(fid);
            header  = textscan(headerLine,'%s','Delimiter','\t');
            header  = header{1,1};
            
            formatString = '';
            for i = 1:length(header)
                formatString = [formatString,' %s'];
            end
            
            cellData = textscan(fid, formatString, 'Delimiter','\t');
            fclose(fid);
            
            %fprintf([datestr(now, 'HH:MM:SS ') 'Processing date/time ...\r']);
            d = char(cellData{dateIndex});
            t = char(cellData{timeIndex});
            dt = [d t];
            dateTime = datetime(dt, 'InputFormat', 'dd-MMM-yyyyHH:mm:ss.SSS');
            profileStartTime(profileStartTimeIndex) = dateTime(1);
            profileStartTimeIndex = profileStartTimeIndex + 1;
            
            %fprintf([datestr(now, 'HH:MM:SS ') 'Processing pressure ...\r']);
            a = cellfun(@str2num, cellData{pressureIndex}, 'UniformOutput', false);
            emptyCellIndexes = cellfun('isempty', a);
            a(emptyCellIndexes) = {nan};
            pressure = cell2mat(a);

            glideOnFind = strfind(cellData{commentIndex}, 'glide on');
            glideOnMsgIndexes = find(cellfun('isempty', glideOnFind) == 0);
            for i = 1:length(glideOnMsgIndexes)
                for j = 0:5
                    glideOnDepth(glideOnIndex) = str2double(char(cellData{pressureIndex}(glideOnMsgIndexes(i)-j)));
                    if ~isnan(glideOnDepth(glideOnIndex)) 
                        break;
                    end
                end
                glideOnTime(glideOnIndex) = dateTime(glideOnMsgIndexes(i));
                glideOnIndex = glideOnIndex + 1;
            end
            
            glideOffFind = strfind(cellData{commentIndex}, 'glide off');
            glideOffMsgIndexes = find(cellfun('isempty', glideOffFind) == 0);
            for i = 1:length(glideOffMsgIndexes)
                for j = 0:5
                    glideOffDepth(glideOffIndex) = str2double(char(cellData{pressureIndex}(glideOffMsgIndexes(i)-j)));
                    if ~isnan(glideOffDepth(glideOffIndex)) 
                        break;
                    end
                end
                glideOffTime(glideOffIndex) = dateTime(glideOffMsgIndexes(i));
                glideOffIndex = glideOffIndex + 1;
            end

            figure(pressureFig);
            plot(dateTime(~isnan(pressure)), pressure(~isnan(pressure)));
            if plotFilename
                text(datenum(dateTime(1)), 0, sortedDirData(fileNum).name, 'rotation', 270, 'HorizontalAlignment', 'left');
            end
            title('Pressure (dBar) vs. Time)');
            
            if length(glideOnMsgIndexes) > 0
                plot(glideOnTime, glideOnDepth, 'gO')
            end
            glideOnMsgIndexes = [];
            if length(glideOffMsgIndexes) > 0
                plot(glideOffTime, glideOffDepth, 'rX')
            end
            glideOffMsgIndexes = [];
            
            %fprintf([datestr(now, 'HH:MM:SS ') 'Processing bellows position ...\n']);
            a = cellfun(@str2num, cellData{bellowsPosIndex}, 'UniformOutput', false);
            emptyCellIndexes = cellfun('isempty', a);
            a(emptyCellIndexes) = {nan};
            bellowsPos = cell2mat(a);
            
%             figure(bellowsPosFig);
%             plot(dateTime(~isnan(bellowsPos)), bellowsPos(~isnan(bellowsPos)));
%             title('Bellows Position (mm) vs. Time'); 

            %fprintf([datestr(now, 'HH:MM:SS ') 'Calculating Power ...\n']);
            clear a
            a = cellfun(@str2num, cellData{batteryVoltsIndex}, 'UniformOutput', false);
            emptyCellIndexes = cellfun('isempty', a);
            a(emptyCellIndexes) = {nan};
            batteryVolts = cell2mat(a);
            
            clear a
            a = cellfun(@str2num, cellData{batteryAmpsIndex}, 'UniformOutput', false);
            emptyCellIndexes = cellfun('isempty', a);
            a(emptyCellIndexes) = {nan};
            batteryAmps = cell2mat(a);
            power = batteryVolts .* batteryAmps;
            
            figure(powerFig);
            plot(dateTime(~isnan(power)), power(~isnan(power)));
            title('Power (watts) vs. Time');
            
            %fprintf([datestr(now, 'HH:MM:SS ') 'Calculating Energy ... \r']);
            powerValuesIndex = ~isnan(power);
            tempPower = power(powerValuesIndex);
            energySumDateTime = dateTime(powerValuesIndex);
            energyIndexes = powerValuesIndex;
            energy = tempPower(1:end-1) .* seconds((energySumDateTime(2:end) - energySumDateTime(1:end-1)));
            
            energySum = cumsum(energy);
            energySum(length(energySum) + 1) = energySum(end);

            totalEnergy = sum(energy);
            energyHistoryTime(processedFileNum) = dateTime(end);
            energyHistoryEnergy(processedFileNum) = totalEnergy;
            fprintf('Total Energy this profile = %.0f Joules\r\r', totalEnergy);
            
            ubxFind = strfind(cellData{gpsIndex}, 'UBX');
            ubxMsgIndexes = find(cellfun('isempty', ubxFind) == 0);
            try
            for i = 1:length(ubxMsgIndexes)
                ubxMsg = char(cellData{gpsIndex}(ubxMsgIndexes(i)));
                ubxSplit = strsplit(ubxMsg, ',');
                numSats(numSatsIndex) = str2double(char(ubxSplit(18)));
                numSatsTime(numSatsIndex) = dateTime(ubxMsgIndexes(i));
                numSatsIndex = numSatsIndex + 1;
            end
            catch me
                fprintf(1, 'Error parsing UBX message\r');
            end
                
            
            sbdixFind = strfind(cellData{iridiumIndex}, '+SBDIX:');
            sbdixMsgIndexes = find(cellfun('isempty', sbdixFind) == 0);
            for i = 1:length(sbdixMsgIndexes)
                sbdixMsg = char(cellData{iridiumIndex}(sbdixMsgIndexes(i)));
                sbdixSplit = strsplit(sbdixMsg, ',');
                numSBDIX(numSBDIXIndex) = i;
                numSBDIXTime(numSBDIXIndex) = dateTime(sbdixMsgIndexes(i));
                numSBDIXIndex = numSBDIXIndex + 1;
            end
            
            errorFind = strfind(cellData{commentIndex}, 'Error');
            errorMsgIndexes = find(cellfun('isempty', errorFind) == 0);
            for i = 1:length(errorMsgIndexes)
                errorMsg(errorIndex) = cellData{commentIndex}(errorMsgIndexes(i));
                errorMsgTime(errorIndex) = dateTime(errorMsgIndexes(i));
                errorIndex = errorIndex + 1;
            end

            bpLimitErrorFind = strfind(cellData{commentIndex}, 'exceeds limits');
            bpLimitErrorMsgIndexes = find(cellfun('isempty', bpLimitErrorFind) == 0);
            for i = 1:length(bpLimitErrorMsgIndexes)
                bpLimitErrorMsg(bpLimitErrorIndex) = cellData{commentIndex}(bpLimitErrorMsgIndexes(i));
                bpLimitErrorMsgTime(bpLimitErrorIndex) = dateTime(bpLimitErrorMsgIndexes(i));
                bpLimitErrorIndex = bpLimitErrorIndex + 1;
            end
        end
    end
end

fprintf(1, 'Summary of errors in %d files processed: \r', processedFileNum );

for i = 1:errorIndex-1
    fprintf(1, 'Error %d at time %s : %s \r', i, datestr(errorMsgTime(i)), errorMsg{i});
end

fprintf(1, '\r\nTotal number of Bellow Position Limit Errors = %d \r', bpLimitErrorIndex - 1);

figure
plot(energySumDateTime, energySum);
title('Energy vs Time')

figure
plotyy(dateTime(~isnan(pressure)), pressure(~isnan(pressure)), energySumDateTime, energySum);

energySumSurfaceOps = 0;
energySumABRetractFast = 0;
energySumABRetractSlow = 0;
energySumDescend = 0;
energySumPark = 0;
energySumAnchor = 0;
energySumStartCPMode = 0;
energySumDeAnchor = 0;
energySumAscend = 0;
energySumInitProfile = 0;
energySumSurfaceOpsSetBellows = 0;
energySumDumpCPData = 0;

energyStates = cellData{1,4}(energyIndexes);
for i=1:length(energyStates)-1
    switch lower(energyStates{i})
        case 'surface ops'
            energySumSurfaceOps = energySumSurfaceOps + energy(i);
        case 'ab retract fast'
            energySumABRetractFast = energySumABRetractFast + energy(i);
        case 'ab retract slow'
            energySumABRetractSlow = energySumABRetractSlow + energy(i);
        case 'descend'
            energySumDescend = energySumDescend + energy(i);
        case 'park'
            energySumPark = energySumPark + energy(i);
        case 'anchor'
            energySumAnchor = energySumAnchor + energy(i);
        case 'start cp mode'
            energySumStartCPMode = energySumStartCPMode + energy(i);
        case 'deanchor'
            energySumDeAnchor = energySumDeAnchor + energy(i);
        case 'ascend'
            energySumAscend = energySumAscend + energy(i);
        case 'init profile'
            energySumInitProfile = energySumInitProfile + energy(i);
        case 'surface ops set bellows'
            energySumSurfaceOpsSetBellows = energySumSurfaceOpsSetBellows + energy(i);
        case 'dump cp data'
            energySumDumpCPData = energySumDumpCPData + energy(i);
    end
end
    fprintf(1, '\r\r');
    fprintf(1, 'Energy sum surface ops %.0f\r', energySumSurfaceOps);
    fprintf(1, 'Energy sum ab retract fast %.0f\r', energySumABRetractFast);
    fprintf(1, 'Energy sum ab retract slow %0.f\r', energySumABRetractSlow);
    fprintf(1, 'Energy sum descend %0.f\r', energySumDescend);
    fprintf(1, 'Energy sum park %0.f\r', energySumPark);
    fprintf(1, 'Energy sum anchor %0.f\r', energySumAnchor);
    fprintf(1, 'Energy sum start cp mode %0.f\r', energySumStartCPMode);
    fprintf(1, 'Energy sum deanchor %0.f\r', energySumDeAnchor);
    fprintf(1, 'Energy sum ascend %0.f\r', energySumAscend);
    fprintf(1, 'Energy sum init profile %0.f\r', energySumInitProfile);
    fprintf(1, 'Energy sum surface ops set bellows %0.f\r', energySumSurfaceOpsSetBellows);
    fprintf(1, 'Energy sum dumpcpdata %0.f\r', energySumDumpCPData);
    
    fprintf(1, 'Energy sum double check = %.0f\r', energySumSurfaceOps + ...
        energySumABRetractFast + energySumABRetractSlow + ...
        energySumDescend + energySumPark + ...
        energySumAnchor + energySumStartCPMode + ...
        energySumDeAnchor + energySumAscend + ...
        energySumInitProfile + energySumSurfaceOpsSetBellows + ...
        energySumDumpCPData);
    
    % figure
    % plot(numSatsTime, numSats, '+-');
% title('Number of GPS Satellites vs Time')
% 
% figure
% plot(numSBDIXTime, numSBDIX, '+-');
% title('Number of SBDIX Attempts vs Time')

