%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%                    LATEXtools                         %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  This  file  is  free software ;  it  is distributed  in  %
%%  the hope that  it  will  be  useful,  but  without  any  %
%%  warranty.  You  can  redistribute  it and/or modify  it  %
%%  under  the  terms  of  the  GNU  General Public License  %
%%  as  published  by  the   Free  Software  Foundation  at  %
%%       http://www.gnu.org/copyleft/gpl.html                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MVP_L0toL1 (MATLAB SCRIPT)
%%
%% Converts .raw files into .mat, performs initial user
%% QC of the MVP profiles.
%%
%% 07_05_16: Altered line 138 to read L0 .raw files in ESSMVP, instead
%% of looking for #, looks for first empty line.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 22/04/16 : louise.rousselet@mio.osupytheas.fr :
%%            andrea.doglioli@mio.osupytheas.fr  : 
%%            alain.de-verneil@mio.osupytheas.fr : creation 
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;close all;

MVP_cruise_params; % Load up the relevant variables for the cruise.

% Get the .raw files in L0
filesL0 = dir([pathL0 filesep '*.raw']);
nFilesL0 = length(filesL0);

% Extract the filenames at the L0 level.
fileNamesL0 = [];
for iFiles = 1:nFilesL0
    fileNamesL0 = [fileNamesL0 ;filesL0(iFiles).name(1:end-4)]; % Don't include .raw at the end.
end


%dirsL0 = dir(pathL0);
%nDirs = length(dirsL0);
%nFilesL0 = 0;
%fileNamesL0 = [];
%for iDirs = 3:nDirs
%    filesL0=dir([pathL0 filesep dirsL0(iDirs).name filesep '*.raw'])
%    nFilesL0 = nFilesL0 + length(filesL0)
%    fileNamesL0 = [fileNamesL0 , filesL0.name] % Don't include .raw at the end.
%end


% Get the pre-existing L1 files (if any)
filesL1 = dir([pathL1 filesep '*.mat']);
nFilesL1 = length(filesL1);

% Extract the filenames at the L1 level, and their QC status.
fileNamesL1 = [];
QCstatL1 = [];
for iFiles = 1:nFilesL1
    fileNamesL1 = [fileNamesL1 ;filesL1(iFiles).name(1:end-6)]; % Don't include .mat or qc flag.
    QCstatL1 = [QCstatL1 str2num(filesL1(iFiles).name(end-4))];
end

hasL1 = zeros(1,nFilesL0); % 1/0 True/False is .raw file previously processed...
maxQC = zeros(1,nFilesL0); % List of max QC for each .raw file already processed...

for iL0Files = 1:nFilesL0
    for iL1Files = 1:nFilesL1
        if strcmp(fileNamesL0(iL0Files,:),fileNamesL1(iL1Files,:))
            hasL1(iL0Files) = 1;
            if QCstatL1(iL1Files) > maxQC(iL0Files)
                maxQC(iL0Files) = QCstatL1(iL1Files);
            end
        end
    end
end

% Display the current QC status of the L1 files.
for iL0Files = 1:nFilesL0
    if hasL1(iL0Files)
        disp(sprintf('Files %s has been processed and has QC flag %d.',fileNamesL0(iL0Files,:),maxQC(iL0Files)));
    end
end
disp(char(10));
% Get the user's QC option.
cnt = -1;
disp('Please choose a quality control option. All files below this QC flag will be processed...');
prompt = 'Quality Control Option: 0 = No QC, 1 = Pressure QC, 2 = T/C QC, 3 = Pressure and T/C: ';
while ~ismember(cnt,[0 1 2 3])
    cnt = input(prompt);
    if isempty(cnt)
        cnt = -1;
    end
end

clc;
% Calculates how many profiles to process, based on user option.
if cumsum(maxQC)==0
    num2Proc = nFilesL0;
else
    num2Proc = length(find(maxQC<cnt));
end

disp(sprintf('Based on your choice, %d files will be processed.',num2Proc));
pause(1); disp(char(10));

% If it's only pressure QC, ask if they want it automated.
auto_flag = -1;
if cnt == 1
    disp('For Pressure QC, do you want the automated selection?');
    disp('Keep in mind, any profiles flagged as bad (QC=9) but are not will have to be manually reassigned.');
    while ~ismember(auto_flag,[0 1])
        auto_flag = input('Input (1 for yes, 0 for no): ');
        if isempty(auto_flag)
            auto_flag = -1;
        end
    end
end


%%
%%%%%%%% CONVERSION TO .MAT FILES %%%%%%%%%%%%%%%%%%%%%
% Now convert all the .raw's into .mat in the L1 folder...
% Work profile by profile

if cnt == 0
    hWait = waitbar(0,'Processing .raw files...');
end

if cnt == 1 && auto_flag == 1
    hWait = waitbar(0,'Processing and Automated pressure QC of .raw files...');
end

numProc = 0;
% Looping through each profile
for iProf = 1 : nFilesL0
    
    if ~hasL1(iProf) || maxQC(iProf)<cnt % Unconverted or unprocessed files only...
        
        numProc = numProc + 1; % Keep tally of how many files made...
        
        % Read the file
        fileName = [pathL0 filesep filesL0(iProf).name];
        metadata.fileName = filesL0(iProf).name(1:end-4);
        
        % Open the file
        fid = fopen(fileName,'rt');
        format long;
        
        %%%%%%% READING THE HEADER AND MVP METADATA %%%%%%%%%
        % Read the header
        hdr = [];
        ln  = [];
        ln = fgets(fid);
%        while ~feof(fid) && ~isempty(strfind(ln,'#')) % OLD CONDITION
        while ~feof(fid) && ~isempty(deblank(ln))
            ln = fgets(fid);
            hdr = [hdr ln];
        end
        
        % - profile index
        iind = strfind(hdr,'Index:');
        if ~isempty(iind)
            metadata.profileIndex = sscanf(hdr(iind:end),'Index: %d');
        else
            disp('Wrong format for File Index!');
        end
        
        % - ship's log
        ilog = strfind(hdr,'Ship''s Log (kt):');
        if ~isempty(ilog)
            metadata.shipLog.value = sscanf(hdr(ilog:end),'Ship''s Log (kt): %f');
            metadata.shipLog.unit  = 'kt';
        else
            disp('Wrong format for Ship''s log!');
        end
        
        % - longitude
        ilon = strfind(hdr,'LON (dddmm.mmmmmmm,E):');
        if ~isempty(ilon)
            lon = sscanf(hdr(ilon:end),'LON (dddmm.mmmmmmm,E): %f,%s[\n]');
            EW  = char(lon(2));
            lon = lon(1);
            londeg = floor(lon/100); lonmin = (lon- londeg*100)/60;
            lon = londeg + lonmin;
            if strcmp(EW,'W')
                lon = 360-lon;
            end
            metadata.lon.value = lon;
            metadata.lon.unit  = 'degE';
        else
            disp('Wrong format for Longitude!');
        end
        
        % - latitude
        ilat = strfind(hdr,'LAT ( ddmm.mmmmmmm,N):');
        if ~isempty(ilat)
            lat = sscanf(hdr(ilat:end),'LAT ( ddmm.mmmmmmm,N): %f,%s[\n]');
            NS  = char(lat(2));
            lat = lat(1);
            latdeg = floor(lat/100); latmin = (lat-latdeg*100)/60;
            lat = latdeg + latmin;
            if strcmp(NS,'S')
                lat = -lat;
            end
            metadata.lat.value = lat;
            metadata.lat.unit  = 'degN';
        else
            disp('Wrong format for Latitude!');
        end
        
        % - time & date
        itime = strfind(hdr, 'Time (hh|mm|ss.s):');
        if ~isempty(itime)
            time = sscanf(hdr(itime:end),'Time (hh|mm|ss.s): %d:%d:%f');
        else
            disp('Wrong format for Time!');
        end
        idate = strfind(hdr, 'Date (dd/mm/yyyy):');
        if ~isempty(itime)
            date = sscanf(hdr(idate:end),'Date (dd/mm/yyyy): %d/%d/%d');
        else
            disp('Wrong format for Date!');
        end
        date_vec = datestr([date(3) date(2) date(1) time(1) time(2) time(3)]);
        
        metadata.date.value = datenum(date_vec);   % days since 01-Jan-0000 00:00:00
        metadata.date.unit  = 'days since 01-Jan-0000 00:00:00';
                
        %%%%%% END OF READING HEADER METADATA %%%%%%%%%%%%%%%%%
        
        %%%%%% READ THE RAW DATA INTO MATLAB STRUCTURES %%%%%%%
        
        % - read the data
        while ~feof(fid) && isempty(deblank(ln))    % empty lines
            ln = fgets(fid);
        end
        
        if ~feof(fid)
            rawData = double(fscanf(fid,dataFormat,[length(columnNames) inf]))';
        end
        
        % Build the variables
        if isempty(rawData)     % if a file doesn't have any data, just skip it
            data = [];
        else
            for iVar = 1 : length(columnNames)
                switch columnNames{iVar}
%                     % - conversion from ADC counts to analog voltage
%                     case 'ADCcounts'
%                         name                = 'ADCcounts';
%                         data.(name).value   = rawData(:,iVar) * 5 / 2^12;
%                         data.(name).unit    = 'V';
                    case 'LOPCcounts'
                        % - counts/0.5s -> counts/s (cf. MVP MSFFF LOPC Integration Guide, p. 15)
                        name                = 'LOPCcounts';
                        data.(name).value   = rawData(:,iVar) * 2;
                        data.(name).unit    = 'V';
                    otherwise
                        name                = columnNames{iVar};
                        data.(name).value   = rawData(:,iVar);
                        data.(name).unit    = columnUnits{iVar};
                end
            end
        end
        
        % Build the parameter structure
        params.instrType = instrType;
        params.fSampl.value = fSample;
        params.fSampl.unit  = 'Hz';
        params.QCcast       = 0;      % <- Default; Signifies the cast hasn't even been looked at!
        % Close the file
        fclose(fid);
        
        if isempty(data)    % At least flag the empty data...
            disp(sprintf('Cast%03d, %s, is empty! ',iProf,fileNames(iProf).name));
            params.QCcast = 9;
            pause(1)
        end
        
        %%%%%%%% END OF READING RAW DATA %%%%%%%%%%%
        
        %%%%%%%% PERFORMING INITIAL QC MEASURES %%%%
        
        % Various QC choices put into action...
        switch cnt
            
            case 1   %%%%% QC OPTION 1 (only Pressure) %%%%
                
                if ~isempty(data)                % Bad cast is empty
                    
                    [maxpres,maxind] = max(data.pres.value);
                    
                    if auto_flag % automated QC
                        
                        if max(data.pres.value) < 100
                            params.QCcast = 9;
                        else
                            if maxind>= round(length(data.pres.value)/2)
                                params.QCcast = 9;
                            end
                            params.QCcast = 1;
                        end
                        
                    else
                        % Check if short cast...
                        if max(data.pres.value)<100
                            plot(data.pres.value);
                            xlabel('Time series (interval depends on instrument...','fontsize',15);
                            ylabel('Pressure (dBar) ','fontsize',15);
                            title([metadata.fileName ' Pressure (short cast?) '],'fontsize',18);
                            axis ij
                            drawnow
                            temp = -1;
                            prompt = 'Is this short profile good or bad? (1 good, 0 bad): ';
                            while ~ismember(temp,[0 1])
                                temp = input(prompt);
                                if isempty(temp)
                                    temp = -1;
                                end
                            end
                            disp(char(10))
                            if temp == 1
                                params.QCcast = 1;
                            else
                                params.QCcast = 9;
                            end
                        % Check if strange shape cast... (max is more than half-way through timeseries)    
                        elseif maxind>= round(length(data.pres.value)/2)
                            plot(data.pres.value);
                            xlabel('Time series (interval depends on instrument...','fontsize',15);
                            ylabel('Pressure (dBar) ','fontsize',15);
                            title([metadata.fileName ' Pressure (short cast?) '],'fontsize',18);
                            axis ij
                            drawnow
                            prompt = 'Is this profile good or bad? Station mode (symmetric) is ok (1 good, 0 bad): ';
                            temp = -1;
                            while ~ismember(temp,[0 1])
                                temp = input(prompt);
                                if isempty(temp)
                                    temp = -1;
                                end
                            end
                            disp(char(10));
                            if temp == 1
                                params.QCcast = 1;
                            else
                                params.QCcast = 9;
                            end
                        end
                        params.QCcast = 1;
                    end
                end
                
            case 2   %%%% QC OPTION 2 (only Temp/Conduc/SoundVel Data %%%%%%%
                
                if ~isempty(data)                % Bad cast is empty
                    
                    subplot(1,3,1)
                    plot(data.pres.value);
                    title([metadata.fileName sprintf(' cast #: %03d',iProf)],'fontsize',18);
                    xlabel('Time'); ylabel('P (dB)');
                    axis ij
                    subplot(1,3,2)
                    plot(data.temp.value,data.pres.value); axis ij;
                    title('Temperature ','fontsize',18);
                    xlabel('T (deg C)'); ylabel('P (dB)');
                    subplot(1,3,3)
                    if isfield(data,'cond')
                        plot(data.cond.value,data.pres.value); axis ij;
                        title('Conductivity ','fontsize',18);
                        xlabel('C (mS/cm)'); ylabel('P (dB)');
                    else
                        plot(data.soundVel.value,data.pres.value); axis ij;
                        title('Sound Speed ','fontsize',18);
                        xlabel('SV (m/s) '); ylabel('P (dB)');
                    end
                    drawnow
                    prompt = 'Is this a good or bad cast? (1 good, 0 bad)';
                    temp = -1;
                    while ~ismember(temp,[0 1])
                        temp = input(prompt);
                        if isempty(temp)
                            temp = -1;
                        end
                    end
                    if temp ==1
                        params.QCcast = 2;
                    else
                        params.QCcast = 9;
                    end
                end
                
            case 3   %%%%% QC OPTION 3 (both Pressure AND Temp/Cond/SoundVel) %%%%%%
                
                if ~isempty(data)                % Bad cast is empty
                    
                    [maxpres,maxind] = max(data.pres.value);
                    % Check if short cast...
                    if max(data.pres.value)<100
                        %                     figure(1)
                        plot(data.pres.value);
                        xlabel('Time series (interval depends on instrument...','fontsize',15);
                        ylabel('Pressure (dBar) ','fontsize',15);
                        title([metadata.fileName ' Pressure (short cast?) '],'fontsize',18);
                        axis ij
                        drawnow
                        prompt = 'Is this short profile good or bad? (1 good, 0 bad): ';
                        temp = -1;
                        while ~ismember(temp,[0 1])
                            temp = input(prompt);
                            if isempty(temp)
                                temp = -1;
                            end
                        end
                        disp(char(10))
                        if temp == 1
                            params.QCcast = 1;
                        else
                            params.QCcast = 9;
                        end
                    elseif maxind>= round(length(data.pres.value)/2) % Check if weird shape...
                        %                     figure(1)
                        plot(data.pres.value);
                        xlabel('Time series (interval depends on instrument...','fontsize',15);
                        ylabel('Pressure (dBar) ','fontsize',15);
                        title([metadata.fileName ' Pressure (short cast?) '],'fontsize',18);
                        axis ij
                        drawnow
                        prompt = 'Is this profile good or bad? Station mode (symmetric) is ok (1 good, 0 bad): ';
                        temp = -1;
                        while ~ismember(temp,[0 1])
                            temp = input(prompt);
                            if isempty(temp)
                                temp = -1;
                            end
                        end
                        disp(char(10));
                        if temp == 1
                            params.QCcast = 1;
                        else
                            params.QCcast = 9;
                        end
                    end
                    
                    if params.QCcast ~= 9
                        disp('Now to check the data profiles...');
                        disp(char(10));
                        subplot(1,3,1)
                        plot(data.pres.value);
                        title([metadata.fileName sprintf(' cast #: %03d',iProf)],'fontsize',18);
                        xlabel('Time'); ylabel('P (dB)');
                        axis ij
                        subplot(1,3,2)
                        plot(data.temp.value,data.pres.value); axis ij;
                        title('Temperature ','fontsize',18);
                        xlabel('T (deg C)'); ylabel('P (dB)');
                        subplot(1,3,3)
                        if isfield(data,'cond')
                            plot(data.cond.value,data.pres.value); axis ij;
                            title('Conductivity ','fontsize',18);
                            xlabel('C (mS/cm)'); ylabel('P (dB)');
                        else
                            plot(data.soundVel.value,data.pres.value); axis ij;
                            title('Sound Speed ','fontsize',18);
                            xlabel('SV (m/s) '); ylabel('P (dB)');
                        end
                        drawnow
                        prompt = 'Is this a good or bad cast? (1 good, 0 bad)';
                        temp = -1;
                        while ~ismember(temp,[0 1])
                            temp = input(prompt);
                            if isempty(temp)
                                temp = -1;
                            end
                        end
                        if temp ==1
                            params.QCcast = 3;
                        else
                            params.QCcast = 9;
                        end
                    end
                end
        end
        %%%%%%% END OF INITIAL QC %%%%%%%
        
        %%%% SAVING PROFILE DATA %%%%%%
        % Save the profile data into a .mat file
        [~, matFileName, ~] = fileparts(fileName);
        save([pathL1 filesep matFileName '_' num2str(params.QCcast)],'params','metadata','data');
        
        % Clear the structures
        clear params metadata data
        clc;
    end
    % Waitbar (only useful is QC is set to 0)
    if cnt ==0
        waitbar(iProf/nFilesL0);
    end
    if cnt ==1 && auto_flag == 1
        waitbar(iProf/nFilesL0);
    end
    if cnt >1
        clf;
    end
    
end

if cnt ==0
    close(hWait);
end
if cnt ==1 && auto_flag == 1
    close(hWait);
end

% Log message
disp(blanks(1)');
disp(['*** Done - ' num2str(numProc) ' .mat files created ***']);
disp(blanks(1)');


disp('You have successfully completed L0 to L1.');
disp('You should be able to now run L1toL2');
disp('Press any key to exit.');
pause

clear all; clc; close all;

%%%%%%%%% END of L0toL1 %%%%%%%%%%%%%%%%%%%%%%%%
