%
%  gliderOptode_Generate_L1_Data.m
%
%  Purpose: Generate Level 1 data mat files that include optode DO data.
%
%  Requires:
%  MATLAB folder - contains util
%
%  Authors:  Adapted from William Stark by Harvey Seim and Chris Calloway
%            Marine Sciences Department
%            UNC-Chapel Hill
%
%  Created: March 2013
%
%///////////////////////////////////////////////////////////////////////////////

clear all;

% add paths for required files...
addpath('MATLAB/util/');

% populate arrays for the deployment start and end dates...
% ex. strStart(2, 3) is start date for Ramses, Deployment 3
strStart = {'26-Jan-2012', '16-Feb-2012', '16-Mar-2012'; ...
            '26-Jan-2012', '16-Feb-2012', '16-Mar-2012'};
strEnd   = {'14-Feb-2012', '08-Mar-2012', '04-Apr-2012'; ...
            '14-Feb-2012', '12-Mar-2012', '03-Apr-2012'};

% SET THE GLIDER INDEX (Pelagia = 1, Ramses = 2) ...
for gliderIndex=1:2
    
    % SET THE DEPLOYMENT NUMBER (1, 2 or 3) ...
    for deploymentNumber=1:3
        
        clearvars -except gliderIndex deploymentNumber strStart strEnd;

        % glider name string...
        if (gliderIndex==1)
            strGliderName = 'Pelagia';
        else
            strGliderName = 'Ramses';
        end

        % deployment number string...
        strDeploymentNumber = num2str(deploymentNumber);

        % deployment start date string...
        strStartDate = strStart(gliderIndex, deploymentNumber);

        % deployment end date string...
        strEndDate = strEnd(gliderIndex, deploymentNumber);

        % define the path to the glider ascii files...
        datadir = strcat('GLIDER_DATA_LEVEL0/', strGliderName, ...
                         '_Deployment', strDeploymentNumber, '/');

        % define optode foil calibration coefficients
        % by glider and deployment
        switch gliderIndex
            case 1 % Pelagia
                switch deploymentNumber
                    case {1, 2} % Deployment 1 or 2;
                                % before foil replacement;
                                % 21 June 2004 calibration date;
                                % batch 2204
                        C0coef=[3.12899E+03,-1.05764E+02,2.06640E+00,-1.68983E-02];
                        C1coef=[-1.67671E+02,4.80582E+00,-8.73323E-02,6.61507E-04];
                        C2coef=[3.73685E+00,-8.78300E-02,1.47560E-03,-9.96701E-06];
                        C3coef=[-3.96052E-02,7.46930E-04,-1,17804E-05,6.677619E-08];
                        C4coef=[1.61999E-04,-2.37870E-06,3.63223E-08,-1.62194E-10];
                    case 3 % Deployment 3;
                           % replaced foil on 12 March 2012;
                           % 18 August 2010 calibration date;
                           % batch 1023
                        C0coef=[4.27019336E+03,-1.32723585E+02,2.15629751E+00,-1.40275831E-02];
                        C1coef=[-2.29729690E+02,5.74242078E+00,-6.85357898E-02,1.88612346E-04];
                        C2coef=[5.06401550E+00,-9.62084932E-02,5.22180779E-04,7.70889717E-06];
                        C3coef=[-5.26332308E-02,7.15467419E-04,3.31185072E-06,-1.86124024E-07];
                        C4coef=[2.10916841E-04,-1.84087896E-06,-4.28645540E-08,1.11120317E-09];
                end
            case 2 % Ramses all deployments;
                   % 2 June 2010 calibration date;
                   % batch 5009
                C0coef=[4.53793e3 -1.62595e2 3.29574 -2.79285e-2];
                C1coef=[-2.50953e2 8.02322 -1.58398e-1 1.31141e-3];
                C2coef=[5.66417 -1.59647e-1 3.07910e-3 -2.46265e-5];
                C3coef=[-5.99449e-2 1.48326e-3 -2.82110e-5 2.15156e-7];
                C4coef=[2.43614e-4 -5.26759e-6 1.00064e-7 -7.14320e-10];
        end

        %*** READ IN EBD DATA ****
        % declare arrays for accumulating data
        oxyw_oxygen = [];
        oxyw_saturation = [];
        oxyw_temp = [];
        oxyw_dphase = [];
        % oxyw_bphase = [];
        % oxyw_rphase = [];
        % oxyw_bamp = [];
        % oxyw_bpot = [];
        % oxyw_ramp = [];
        % oxyw_rawtemp = [];
        % oxyw_time = [];
        % oxyw_installed = [];
        ptime_ebd = [];

        % try to load all *.ebdasc files at once...
        [files, Dstruct] = wilddir(datadir, '.ebdasc');
        nfile = size(files, 1);

        for i=1:nfile-1
            % protect against empty ebd file
            if(Dstruct(i).bytes>0)
                data = read_gliderasc2([datadir, files(i,:)]);

                % if the number of values (in data.data) is less than the number
                % of vars (in data.vars), this means that the data were not
                % completely read in. To correct this, pad data.data with NaNs
                % until its length equals that of data.vars...
                if (length(data.data) < length(data.vars))
                    data.data = padarray(data.data, ...
                        [0 length(data.vars)-length(data.data)], ...
                        NaN, 'post');
                end

                % concatenate variables with data...
                if(~isempty(data.data))
                    oxyw_oxygen = [oxyw_oxygen; ...
                        data.data(:,strmatch('sci_oxy3835_wphase_oxygen',...
                                             data.vars, 'exact'))];
                    oxyw_saturation = [oxyw_saturation; ...
                        data.data(:,strmatch('sci_oxy3835_wphase_saturation',...
                                             data.vars, 'exact'))];
                    oxyw_temp = [oxyw_temp; ...
                        data.data(:,strmatch('sci_oxy3835_wphase_temp',...
                                             data.vars, 'exact'))];
                    oxyw_dphase = [oxyw_dphase; ...
                        data.data(:,strmatch('sci_oxy3835_wphase_dphase',...
                                             data.vars, 'exact'))];
                    % oxyw_bphase = [oxyw_bphase; ...
                    %     data.data(:,strmatch('sci_oxy3835_wphase_bphase',...
                    %                          data.vars, 'exact'))];
                    % oxyw_rphase = [oxyw_rphase; ...
                    %     data.data(:,strmatch('sci_oxy3835_wphase_rphase',...
                    %                          data.vars, 'exact'))];
                    % oxyw_bamp = [oxyw_bamp; ...
                    %     data.data(:,strmatch('sci_oxy3835_wphase_bamp',...
                    %                          data.vars, 'exact'))];
                    % oxyw_bpot = [oxyw_bpot; ...
                    %     data.data(:,strmatch('sci_oxy3835_wphase_bpot',...
                    %                          data.vars, 'exact'))];
                    % oxyw_ramp = [oxyw_ramp; ...
                    %     data.data(:,strmatch('sci_oxy3835_wphase_ramp',...
                    %                          data.vars, 'exact'))];
                    % oxyw_rawtemp = [oxyw_rawtemp; ...
                    %     data.data(:,strmatch('sci_oxy3835_wphase_rawtemp',...
                    %                          data.vars, 'exact'))];
                    % oxyw_time = [oxyw_time; ...
                    %     data.data(:,strmatch('sci_oxy3835_wphase_timestamp',...
                    %                          data.vars, 'exact'))];
                    % oxyw_installed = [oxyw_installed; ...
                    %     data.data(:,strmatch('sci_oxy3835_wphase_is_installed',...
                    %                          data.vars, 'exact'))];
                    ptime_ebd = [ptime_ebd; ...
                        data.data(:,strmatch('sci_m_present_time',...
                                             data.vars, 'exact'))];
                end

                data = [];
            end  
        end
        %*** END READ IN EBD DATA ****

        % remove nans from EBD data...
        i = find(~isnan(oxyw_dphase));
        oxyw_dphase = oxyw_dphase(i);
        ptime_ebd = ptime_ebd(i);
        oxyw_oxygen = oxyw_oxygen(i);
        oxyw_saturation = oxyw_saturation(i);
        oxyw_temp = oxyw_temp(i);
        % oxyw_bphase = oxyw_bphase(i);
        % oxyw_rphase = oxyw_rphase(i);
        % oxyw_bamp = oxyw_bamp(i);
        % oxyw_bpot = oxyw_bpot(i);
        % oxyw_ramp = oxyw_ramp(i);
        % oxyw_rawtemp = oxyw_rawtemp(i);
        % oxyw_time = oxyw_time(i);
        % oxyw_installed = oxyw_installed(i);

        % apply the sort() function to ptime_oxy
        % to make sure it increases monotonically...
        [Y,I] = sort(ptime_ebd);
        ptime_ebd = Y;
        oxyw_oxygen = oxyw_oxygen(I);
        oxyw_saturation = oxyw_saturation(I);
        oxyw_temp = oxyw_temp(I);
        oxyw_dphase = oxyw_dphase(I);
        % oxyw_bphase = oxyw_bphase(I);
        % oxyw_rphase = oxyw_rphase(I);
        % oxyw_bamp = oxyw_bamp(I);
        % oxyw_bpot = oxyw_bpot(I);
        % oxyw_ramp = oxyw_ramp(I);
        % oxyw_rawtemp = oxyw_rawtemp(I);
        % oxyw_time = oxyw_time(I);
        % oxyw_installed = oxyw_installed(I);

        % load CTD data from existing mat file...
        strMatFilePath = strcat('GLIDER_DATA_LEVEL1/', strGliderName, ...
            '_Deployment', strDeploymentNumber, '_CTD_L1.mat');
        load(strMatFilePath);

        % interpolate CTD to times of oxyw_dphase 
        tempi = interp1(ptime,temp,ptime_ebd);
        salini = interp1(ptime,salinCorrected,ptime_ebd);
        depthi = interp1(ptime,depth,ptime_ebd);

        % first, implement temperature-dependent correction to DO concentration
        % that utilizes dphase as the input (manual page 30), coefficients from cal
        % sheets (a 5x4 matrix of values) and will be glider (and foil) dependent.
        C0 = C0coef(1) + (C0coef(2) .* tempi) + ...
            (C0coef(3) .* (tempi.^2)) + (C0coef(4) .* (tempi.^3));
        C1 = C1coef(1) + (C1coef(2) .* tempi) + ...
            (C1coef(3) .* (tempi.^2)) + (C1coef(4) .* (tempi.^3));
        C2 = C2coef(1) + (C2coef(2) .* tempi) + ...
            (C2coef(3) .* (tempi.^2)) + (C2coef(4) .* (tempi.^3));
        C3 = C3coef(1) + (C3coef(2) .* tempi) + ...
            (C3coef(3) .* (tempi.^2)) + (C3coef(4) .* (tempi.^3));
        C4 = C4coef(1) + (C4coef(2) .* tempi) + ...
            (C4coef(3) .* (tempi.^2)) + (C4coef(4) .* (tempi.^3));

        o2_tcorr = C0 + (C1 .* oxyw_dphase) + (C2 .* (oxyw_dphase.^2)) + ...
            (C3 .* (oxyw_dphase.^3)) + (C4 .* (oxyw_dphase.^4));

        % second, implement the salinity correction to DO concentration (page 31)
        % ASSUMES DEFAULT SALINITY IN SENSOR WAS SET TO ZERO, OTHERWISE ANOHTER
        % SCALING IS REQUIRED
        B0=-6.24097e-3;
        B1=-6.93498e-3;
        B2=-6.90358e-3;
        B3=-4.29155e-3;
        C0=-3.11680e-7;

        Ts = log((298.15 - tempi) ./ (273.15 + tempi));

        o2_tscorr = o2_tcorr .* ...
            exp((salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + ...
                            (B3 .* (Ts.^3))))  + ...
            (C0 .* (salini.^2)));

        % third, implement the pressure correction to DO concentration (page 32)
        o2_tspcorr = o2_tscorr .* (1 + (0.04 .* depthi ./ 1000));

        % use polynomial to calculate DO satuations using the measured temp and
        % sal (manual page 30)
        A0 = 2.00856;
        A1 = 3.22400;
        A2 = 3.99063;
        A3 = 4.80299;
        A4 = 9.78188e-1;
        A5 = 1.71069;
        B0 = -6.24097e-3;
        B1 = -6.93498e-3;
        B2 = -6.90358e-3;
        B3 = -4.29155e-3;
        C0 = -3.11680e-7;

        % need interpolated temp, salinity, pressure at times of DO obs
        rslt = (A0 + (A1 .* Ts) + (A2 .* (Ts.^2)) + (A3 .* (Ts.^3)) + ...
               (A4 .* (Ts.^4)) + (A5 .* (Ts.^5))) + ...
            (salini .* (B0 + (B1 .* Ts) + (B2 .* (Ts.^2)) + (B3 .*(Ts.^3)))) + ...
            (C0 .* (salini.^2));

        o2_sol = exp(rslt);

        o2_sat = (o2_tspcorr .* 2.2414) ./ o2_sol;

        % convert ptime into datenum style...
        ptime_ebd_datenum = (ptime_ebd/3600/24) + datenum(1970, 1, 1, 0, 0, 0);

        % create configuration struct...
        units = struct( ...
            'oxyw_oxygen', '10e-6 mol/dm3', ...
            'oxyw_saturation', 'percent', ...
            'oxyw_temp', 'degrees C', ...
            'oxyw_dphase', 'degrees', ...
            ... % 'oxyw_bphase', 'degrees', ...
            ... % 'oxyw_rphase', 'degrees', ...
            ... % 'oxyw_bamp', 'mA', ...
            ... % 'oxyw_bpot', 'mV', ...
            ... % 'oxyw_ramp', 'mA', ...
            ... % 'oxyw_rawtemp', 'degrees C', ...
            ... % 'oxyw_time', 'timestamp', ...
            ... % 'oxyw_installed', 'bool', ...
            'ptime_ebd', 'seconds since 0000-01-01T00:00', ...
            'ptime_ebd_datenum', 'days since 1970-01-01T00:00', ...
            'tempi', 'degrees C', ...
            'salini', 'psu', ...
            'depthi', 'meters', ...
            'o2_tcorr', '10e-6 mol/dm3', ...
            'o2_tscorr', '10e-6 mol/dm3', ...
            'o2_tspcorr', '10e-6 mol/dm3', ...
            'o2_sol', 'cm3/liter at 1031 hPa', ...
            'o2_sat', 'percent');
        
        variable_description = struct( ...
            'oxyw_oxygen', 'dissolved oxygen', ...
            'oxyw_saturation', 'dissolved oxygen saturation', ...
            'oxyw_temp', 'water temperature', ...
            'oxyw_dphase', 'phase difference', ...
            ... % 'oxyw_bphase', 'blue phase', ...
            ... % 'oxyw_rphase', 'red phase', ...
            ... % 'oxyw_bamp', 'blue current bias', ...
            ... % 'oxyw_bpot', 'blue voltage bias', ...
            ... % 'oxyw_ramp', 'red current bias', ...
            ... % 'oxyw_rawtemp', 'raw water temperature', ...
            ... % 'oxyw_time', 'optode timestampe', ...
            ... % 'oxyw_installed', 'bool', ...
            'ptime_ebd', 'science computer time', ...
            'ptime_ebd_datenum', 'science computer date', ...
            'tempi', 'interpolated CTD water temperature', ...
            'salini', 'interpoloted CTD salinity', ...
            'depthi', 'interpolated CTD depth', ...
            'o2_tcorr', 'temperature corrected dissolved oxygen', ...
            'o2_tscorr', 'temperature and salinity corrected dissolved oxygen', ...
            'o2_tspcorr', 'temperature, salinity, and pressure corrected dissolved oxygen', ...
            'o2_sol', 'corrected oxygen solubility', ...
            'o2_sat', 'corrected oxygen saturation');

        correction_coefficients = struct('C0coef', C0coef, ...
                                         'C1coef', C1coef, ...
                                         'C2coef', C2coef, ...
                                         'C3coef', C3coef, ...
                                         'C4coef', C4coef);

        config = struct('glider_name', strGliderName,...
                        'deployment_number', strDeploymentNumber,...
                        'start_date', strStartDate,...
                        'end_date', strEndDate,...
                        'correction_coefficients', correction_coefficients,...
                        'var_descriptions', variable_description,...
                        'var_units', units);

        % set Level 1 data mat file name...
        strMatFileName = strcat(strGliderName, '_Deployment', strDeploymentNumber, '_DO_L1.mat');

        % save flight data to mat file...
        save(strMatFileName,...
             'config', ...
             'oxyw_oxygen', ...
             'oxyw_saturation', ...
             'oxyw_temp', ...
             'oxyw_dphase', ...
             ... % 'oxyw_bphase', ...
             ... % 'oxyw_rphase', ...
             ... % 'oxyw_bamp', ...
             ... % 'oxyw_bpot', ...
             ... % 'oxyw_ramp', ...
             ... % 'oxyw_rawtemp', ...
             ... % 'oxyw_time', ...
             ... % 'oxyw_installed', ...
             'ptime_ebd', ...
             'ptime_ebd_datenum', ...
             'tempi', ...
             'salini', ...
             'depthi', ...
             'o2_tcorr', ...
             'o2_tscorr', ...
             'o2_tspcorr', ...
             'o2_sol', ...
             'o2_sat');

    end
end