processGliderData

PROCESSGLIDERDATA  Glider data processing applying conversions and derivations.

  Syntax:
    [DATA_PROC, META_PROC] = PROCESSGLIDERDATA(DATA_PRE, META_PRE)
    [DATA_PROC, META_PROC] = PROCESSGLIDERDATA(DATA_PRE, META_PRE, OPTIONS)
    [DATA_PROC, META_PROC] = PROCESSGLIDERDATA(DATA_PRE, META_PRE, OPT1, VAL1, ...)

  Description:
    DATA_PROC = PROCESSGLIDERDATA(DATA_PRE, META_PRE, ...) processes 
    preprocessed data from a glider deploymnet according to given options, 
    performing the following actions:
      - Interpolation of reference sensors: 
        Missing values of time, latitude and longitude sequences are filled 
        by interpolation if required.
      - Interpolation of optional reference sensors:
        Missing values of navigation depth, pitch, roll and heading sequences, 
        if present, are filled by interpolation if required.
      - Identification of transects:
        Transects are identified finding their boundaries at changes of 
        waypoint coordinates, if they are available.
      - Computation of distance over ground:
        The planar distance covered along the trajectory is computed cumulating
        the distance between consecutive positions with valid coordinates.
      - Pressure processing:
        Pressure is optionally filtered using a filter proposed by Seabird.
        Depth is optionally derived from pressure and longitude sequences.
      - Identification of casts:
        Upcasts and downcasts are identified finding local extrema of the
        chosen depth or pressure sequence, and the glider vertical direction
        is deduced.
      - CTD flow speed derivation:
        Flow speed through the CTD cell may be derived from selected depth,
        time and pitch sequences. A nominal pitch value may be given if the
        pitch sequence is not available.
      - Sensor lag correction:
        Any already selected sequence may be corrected from sensor lag effects.
        The sensor lag time constant may be provided as option or estimated
        from identified consecutive casts with opposite directions.
      - Thermal lag correction:
        Any temperature and conductivity sequence pair may be corrected from
        thermal lag effects. The thermal lag parameters may be provided as
        option or estimated from identified consecutive casts with opposite
        directions.
      - Salinity derivation:
        In situ salinity may be derived from any set of conductivity,
        temperature and pressure sequences already selected or produced.
      - Density derivation:
        In situ density may be derived from any set of conductivity,
        temperature and pressure sequences already selected or produced.

    DATA_PRE should be a struct in the format returned by PREPROCESSGLIDERDATA,
    where each field is a sequence of measurements of the variable with the 
    same name.

    DATA_PROC is a struct in the same format as DATA_PRE, with time sequences 
    resulting from the processing actions described above, performed according
    to the options described below.

    META_PROC is also a struct with one field per variable, adding processing 
    metadata to any existing metadata in META_PRE.

    Options may be given in key-value pairs OPT1, VAL1... or in a struct
    OPTIONS with field names as option keys and field values as option values.
    Recognized options are:
      TIME_FILLING: time interpolation switch.
        Boolean setting whether time missing values should be filled by 
        interpolation.
        Default value: false
      POSITION_FILLING: position interpolation switch.
        Boolean setting whether latitude and longitude missing values should be 
        filled by interpolation.
        Default value: false
      DEPTH_FILLING: depth interpolation switch.
        Boolean setting whether depth missing values should be filled by
        interpolation.
        Default value: false
      ATTITUDE_FILLING: attitude interpolation switch.
        Boolean setting whether roll and pitch missing values should be
        filled by interpolation.
        Default value: false
      HEADING_FILLING: heading interpolation switch.
        Boolean setting whether heading missing values should be filled by
        interpolation.
        Default value: false
      WAYPOINT_FILLING: waypoint interpolation switch.
        Boolean setting whether waypoint latitude and longitude missing values 
        should be filled with the previous valid value.
        Default value: true
      PRESSURE_FILTERING: Seabird pressure filtering switch.
        Boolean setting whether pressure should be filtered with low pass
        filter described in the Seabird Data Processing Manual.
        Default value: true
      PRESSURE_FILTER_CONSTANT: Seabird pressure filter parameter.
        Non-negative number, the time constant for the Seabird low-pass filter.
        Default value: 4 (recommended by Seabird Data Processing Manual)
      DEPTH_CTD_DERIVATION: depth from CTD pressure derivation.
        Boolean setting whether a depth sequence should be derived from CTD
        pressure readings.
        Default value: true
      PROFILING_LIST: cast identification settings.
        Struct array selecting input sequences and parameters for cast
        boundary identification, in order of preference. It should have 
        the following fields:
          DEPTH: depth or pressure sequence name.
        It may have the following optional fields (empty or missing):
          TIME: time sequence name or empty.
            Default value: []
          STALL: scalar with the maximum vertical displacement when stalled.
            Default value: 3
          SHAKE: scalar with the maximum duration of a vertical shake.
            Default value: 20
          INVERSION: scalar with the maximum depth inversion allowed during a
            cast.
            Default value: 3
          INTERRUPT: scalar with the maximum duration of stalled and/or shake 
            intervals during a cast.
            Default value: 180
          LENGTH: scalar with the minimum depth range a cast must span.
            Default value: 10
          PERIOD: scalar with the minimum duration range a cast must last.
            Default value: 0
        Each struct in the array specifies a choice of inputs for cast 
        boundary identification. The time sequence is optional and is
        relevant only if the SHAKE, INTERRUPT or PERIOD options are used.
        Identification will be performed with the first input choice whose 
        input sequences are available. If no input choice is available,
        profiles are not identified.
        Default value: struct('depth', {'depth' 'depth_ctd' 'depth_ctd'}, ...
                              'time',  {'time'  'time_ctd'  'time'});
      PROFILE_MIN_RANGE: minimum depth range allowed for a valid profile.
        Non-negative real number setting the minimum depth range threshold for
        cast validation. If the difference between the maximum and the minimum 
        depth of a valid reading in a cast is less than the given threshold,
        the cast will be discarded. Set it to 0 to prevent discarding any cast.
        Default value: 10
      PROFILE_MAX_GAP_RATIO: maximum gap ratio allowed for a valid profile.
        Real number (in [0..1]) setting the maximum gap ratio threshold for
        cast validation. A gap is a sequence of consecutive readings taken
        during a depth inversion or in which the value of at least one of the
        involved sensors is invalid. The gap ratio is the ratio of the depth 
        range covered during the gap over the total depth covered of the cast.
        If the ratio of the largest gap over the total depth range is greater
        than the given threshold, the cast will be discarded.
        Set it to 1 to prevent discarding any cast.
        Default value: 0.8
      FLOW_CTD_LIST: CTD flow speed derivation input set choices.
        Struct array selecting input sequences for CTD flow speed derivation,
        in order of preference. It should have the following fields:
          TIME: time sequence name.
          DEPTH: depth sequence name.
          PITCH: pitch sequence name.
        Each struct in the array specifies a choice of inputs for CTD flow 
        speed derivation. Pitch sequence is optional. If missing or empty, 
        the nominal pitch value may be used (see FLOW_CTD_PITCH VALUE below).
        Derivation will be performed only if the casts are properly identified,
        and using the first input choice whose time and depth sequences are 
        available, and either the pitch sequence is available or the pitch 
        nominal value is set. If no input choice is available, CTD flow speed
        is not derived.
        Default value: struct('time', {}, 'depth', {}, 'pitch', {})
      FLOW_CTD_PITCH_VALUE: nominal pitch value for CTD flow speed derivation.
        Number with the nominal pitch value (radians) to use for CTD flow
        speed derivation when no pitch sequence is available.
        Default value: [] (no default pitch value)
      FLOW_CTD_MIN_PITCH: low pitch threshold for CTD flow derivation.
        Number with the minimum absolute pitch value below which flow speed is 
        considered invalid during CTD flow speed derivation.
        Default value: 0 (all values are valid).
      FLOW_CTD_MIN_VELOCITY: low velocity threshold for CTD flow derivation.
        Number with the minimum absolute vertical velocity value below which 
        flow speed is considered invalid during CTD flow speed derivation.
        Default value: 0 (all values are valid).
      SENSOR_LAG_LIST: sensor lag correction settings.
        Struct array specifying the sequences to produce by correcting
        the sensor lag in the corresponding original sensor sequences.
        It should have the following fields:
          CORRECTED: string with the name for the corrected sequence (field in 
            struct DATA_PROC).
          ORIGINAL: string with the name of the original sequence (field in 
            struct DATA_PROC).
          PARAMETERS: non-negative number as predefined time constant,
            or string 'auto' for automatic estimation from casts.
        It may have the following optional fields (empty or missing):
          TIME: string cell array with the names of the time sequence to use 
            for estimation or correction, in order of preference.
            Default value: {'time'}
          DEPTH: string cell array with the names of the depth sequence to use 
            for estimation or correction, in order of preference.
            Default value: {'depth'}
          FLOW: string cell array with the names of the flow sequence to use
            for estimation or correction, in order of preference. This is only
            used if flow speed is not constant.
            Default value: {'flow_ctd'}
          CONSTANT_FLOW: boolean setting whether parameters are static or
            dynamic (varying with the flow speed).
            Default value: false
          ESTIMATOR: function handle or string with the name of the estimator
            to combine the parameter estimates computed for each cast pair.
            Default value: @nanmedian
          MINOPTS: struct to pass custom minimization options for estimation,
            in the format accepted by function FINDSENSORLAGPARAMS.
            Default value: struct()
        Each struct in the struct array specifies a sensor lag correction.
        It will be performed only if the casts are properly identified,
        all the input sequences are available, and the correction parameter is
        available too (either given as option or estimated from pair of casts).
        Default value: struct('corrected',  {}, ...
                              'original',   {}, ...
                              'parameters', {})
      THERMAL_LAG_LIST: thermal lag correction settings.
        Struct array specifying the temperature and conductivity sequences
        to produce by correcting the thermal lag in the corresponding 
        original sensor sequences.
        It should have the following fields:
          CONDUCTIVITY_CORRECTED: string with the name for the corrected
            conductivity sequence (field in DATA_PROC).
          TEMPERATURE_CORRECTED: string with the name for the corrected
            temperature sequence (field in DATA_PROC).
          CONDUCTIVITY_ORIGINAL: string with the name of the original
            conductivity sequence (field in DATA_PROC).
          TEMPERATURE_ORIGINAL: string with the name of the original 
            temperature sequence (field in DATA_PROC).
          PRESSURE_ORIGINAL: string with the name of the original 
            pressure sequence (field in DATA_PROC).
          PARAMETERS: numeric vector with predefined thermal lag parameters or
            string 'auto' for automatic estimation from casts. If a vector, it
            should be a 2 element array when flow speed is constant (error and 
            error time), and a 4 element array otherwise (error offset, error 
            slope, error time offset and error time slope).
        It may have the following optional fields (empty or missing):
          TIME: string cell array with the names of the time sequence to use 
            for estimation or correction, in order of preference.
            Default value: {'time_ctd' 'time'}
          DEPTH: string cell array with the names of the depth sequence to use 
            for estimation or correction, in order of preference.
            Depth is only used to ignore invalid profiles.
            Default value: {'depth_ctd' 'depth'}
          FLOW: string cell array with the names of the flow sequence to use
            for estimation or correction, in order of preference.
            This is only used if flow speed is not constant.
            Default value: {'flow_ctd'}
          CONSTANT_FLOW: boolean setting whether parameters are static or 
            dynamic (varying with the flow speed).
            Default value: false
          ESTIMATOR: function handle or string with the name of the estimator
             to combine the parameter estimates computed for each cast pair.
            Default value: @nanmedian
          MINOPTS: struct to pass custom minimization options for estimation,
            in the format accepted by function FINDTHERMALLAGPARAMS.
            Default value: struct()
        Each struct in the struct array specifies a thermal lag correction.
        It will be performed only if casts are properly identified, all the 
        input sequences are available, and the correction parameters are
        available too (either given as option or estimated from pair of casts).
        Default value: struct('conductivity_corrected', {'conductivity_corrected_thermal'}, ...
                              'temperature_corrected',  {'temperature_corrected_thermal'}, ...         
                              'conductivity_original', {'conductivity'}, ...
                              'temperature_original', {'temperature'}, ...
                              'pressure_original', {'pressure'}, ...
                              'parameters', {'auto'} )
      SALINITY_LIST: salinity derivation settings.
        Struct cell array specifying which salinity sequences to produce
        by derivation from corresponding conductivity, temperature 
        and pressure sequences. It should have the following fields:
          SALINITY: string with the name for the salinity sequence
            (field in DATA_PROC).
          CONDUCTIVITY: string with the name of the original conductivity 
            sequence (field in DATA_PROC).
          TEMPERATURE: string with the name of the original temperature 
            sequence (field in DATA_PROC).
          PRESSURE: string with the name of the original pressure sequence
            (field in DATA_PRE).
        Each struct in the struct array specifies a salinity derivation.
        It will be performed only if all the original sequences are available.
        Default value: struct('salinity',     {'salinity     salinity_corrected_thermal'}, ...
                              'conductivity', {'conductivity conductivity'}, ...
                              'temperature',  {'temperature  temperature_corrected_thermal'},
                              'pressure',     {'pressure' '  pressure});
      DENSITY_LIST: density derivation settings.
        Struct cell array specifying which salinity sequences to produce
        by derivation from corresponding salinity, temperature
        and pressure sequences. It should have the following fields:
          DENSITY: string with the name for the density sequence (field in 
            DATA_PROC).
          SALINITY: string with the name of the original salinity sequence 
            (field in DATA_PROC).
          TEMPERATURE: string with the name of the original temperature 
            sequence (field in DATA_PROC).
          PRESSURE: string with the name of the original pressure sequence (field
            in DATA_PROC).
        Each struct in the struct array specifies a density derivation.
        It will be performed only if all the or'time'iginal sequences are available.
        Default value: struct('density',     {'density'     'density_corrected_thermal'}, ...
                              'salinity',    {'salinity'    'salinity_corrected_thermal'}, ...
                              'temperature', {'temperature' 'temperature'}, ...
                              'pressure',    {'pressure'    'pressure'})

    The following options are deprecated and should not be used:
      PROFILING_SEQUENCE_LIST: sequence choices for cast identification.
        String cell array with the names of the pressure or depth sequence to
        use for cast identification, in order of preference.
        Deprecated in v1.1.0:
          Superseeded by PROFILING_LIST.
      PROFILE_JOIN_SAME_DIR: join consecutive profiles with the same direction.
        Boolean setting whether consecutive valid profiles with the same 
        vertical direction should be joined into one.
        Deprecated in v1.1.0:
          Superseeded by PROFILING_LIST.
      PROFILING_SEQUENCE_FILLING: profiling sequence interpolation switch.
        Boolean setting whether the missing values in the profiling sequence
        should be filled by interpolation before cast identification.
        Deprecated in v1.1.0:
          Not needed anymore since profile identification is consistent and 
          robust against missing values.

  Notes:
    This function is based on the previous work by Tomeu Garau. He is the true
    glider man.

  Examples:
    [data_proc, meta_proc] = processGliderData(data_pre, meta_pre, options)

  See also:
    PREPROCESSGLIDERDATA
    FILLINVALIDVALUES
    FINDTRANSECTS
    COMPUTECUMULATIVEDISTANCE
    FINDPROFILES
    VALIDATEPROFILE
    APPLYSEABIRDPRESSUREFILTER
    FINDSENSORLAGPARAMS
    CORRECTSENSORLAG
    FINDTHERMALLAGPARAMS
    CORRECTTHERMALLAG
    SW_DPTH
    SW_SALT
    SW_DENS

  Authors:
    Joan Pau Beltran  <joanpau.beltran@socib.cat>

0001 function [data_proc, meta_proc] = processGliderData(data_pre, meta_pre, varargin)
0002 %PROCESSGLIDERDATA  Glider data processing applying conversions and derivations.
0003 %
0004 %  Syntax:
0005 %    [DATA_PROC, META_PROC] = PROCESSGLIDERDATA(DATA_PRE, META_PRE)
0006 %    [DATA_PROC, META_PROC] = PROCESSGLIDERDATA(DATA_PRE, META_PRE, OPTIONS)
0007 %    [DATA_PROC, META_PROC] = PROCESSGLIDERDATA(DATA_PRE, META_PRE, OPT1, VAL1, ...)
0008 %
0009 %  Description:
0010 %    DATA_PROC = PROCESSGLIDERDATA(DATA_PRE, META_PRE, ...) processes
0011 %    preprocessed data from a glider deploymnet according to given options,
0012 %    performing the following actions:
0013 %      - Interpolation of reference sensors:
0014 %        Missing values of time, latitude and longitude sequences are filled
0015 %        by interpolation if required.
0016 %      - Interpolation of optional reference sensors:
0017 %        Missing values of navigation depth, pitch, roll and heading sequences,
0018 %        if present, are filled by interpolation if required.
0019 %      - Identification of transects:
0020 %        Transects are identified finding their boundaries at changes of
0021 %        waypoint coordinates, if they are available.
0022 %      - Computation of distance over ground:
0023 %        The planar distance covered along the trajectory is computed cumulating
0024 %        the distance between consecutive positions with valid coordinates.
0025 %      - Pressure processing:
0026 %        Pressure is optionally filtered using a filter proposed by Seabird.
0027 %        Depth is optionally derived from pressure and longitude sequences.
0028 %      - Identification of casts:
0029 %        Upcasts and downcasts are identified finding local extrema of the
0030 %        chosen depth or pressure sequence, and the glider vertical direction
0031 %        is deduced.
0032 %      - CTD flow speed derivation:
0033 %        Flow speed through the CTD cell may be derived from selected depth,
0034 %        time and pitch sequences. A nominal pitch value may be given if the
0035 %        pitch sequence is not available.
0036 %      - Sensor lag correction:
0037 %        Any already selected sequence may be corrected from sensor lag effects.
0038 %        The sensor lag time constant may be provided as option or estimated
0039 %        from identified consecutive casts with opposite directions.
0040 %      - Thermal lag correction:
0041 %        Any temperature and conductivity sequence pair may be corrected from
0042 %        thermal lag effects. The thermal lag parameters may be provided as
0043 %        option or estimated from identified consecutive casts with opposite
0044 %        directions.
0045 %      - Salinity derivation:
0046 %        In situ salinity may be derived from any set of conductivity,
0047 %        temperature and pressure sequences already selected or produced.
0048 %      - Density derivation:
0049 %        In situ density may be derived from any set of conductivity,
0050 %        temperature and pressure sequences already selected or produced.
0051 %
0052 %    DATA_PRE should be a struct in the format returned by PREPROCESSGLIDERDATA,
0053 %    where each field is a sequence of measurements of the variable with the
0054 %    same name.
0055 %
0056 %    DATA_PROC is a struct in the same format as DATA_PRE, with time sequences
0057 %    resulting from the processing actions described above, performed according
0058 %    to the options described below.
0059 %
0060 %    META_PROC is also a struct with one field per variable, adding processing
0061 %    metadata to any existing metadata in META_PRE.
0062 %
0063 %    Options may be given in key-value pairs OPT1, VAL1... or in a struct
0064 %    OPTIONS with field names as option keys and field values as option values.
0065 %    Recognized options are:
0066 %      TIME_FILLING: time interpolation switch.
0067 %        Boolean setting whether time missing values should be filled by
0068 %        interpolation.
0069 %        Default value: false
0070 %      POSITION_FILLING: position interpolation switch.
0071 %        Boolean setting whether latitude and longitude missing values should be
0072 %        filled by interpolation.
0073 %        Default value: false
0074 %      DEPTH_FILLING: depth interpolation switch.
0075 %        Boolean setting whether depth missing values should be filled by
0076 %        interpolation.
0077 %        Default value: false
0078 %      ATTITUDE_FILLING: attitude interpolation switch.
0079 %        Boolean setting whether roll and pitch missing values should be
0080 %        filled by interpolation.
0081 %        Default value: false
0082 %      HEADING_FILLING: heading interpolation switch.
0083 %        Boolean setting whether heading missing values should be filled by
0084 %        interpolation.
0085 %        Default value: false
0086 %      WAYPOINT_FILLING: waypoint interpolation switch.
0087 %        Boolean setting whether waypoint latitude and longitude missing values
0088 %        should be filled with the previous valid value.
0089 %        Default value: true
0090 %      PRESSURE_FILTERING: Seabird pressure filtering switch.
0091 %        Boolean setting whether pressure should be filtered with low pass
0092 %        filter described in the Seabird Data Processing Manual.
0093 %        Default value: true
0094 %      PRESSURE_FILTER_CONSTANT: Seabird pressure filter parameter.
0095 %        Non-negative number, the time constant for the Seabird low-pass filter.
0096 %        Default value: 4 (recommended by Seabird Data Processing Manual)
0097 %      DEPTH_CTD_DERIVATION: depth from CTD pressure derivation.
0098 %        Boolean setting whether a depth sequence should be derived from CTD
0099 %        pressure readings.
0100 %        Default value: true
0101 %      PROFILING_LIST: cast identification settings.
0102 %        Struct array selecting input sequences and parameters for cast
0103 %        boundary identification, in order of preference. It should have
0104 %        the following fields:
0105 %          DEPTH: depth or pressure sequence name.
0106 %        It may have the following optional fields (empty or missing):
0107 %          TIME: time sequence name or empty.
0108 %            Default value: []
0109 %          STALL: scalar with the maximum vertical displacement when stalled.
0110 %            Default value: 3
0111 %          SHAKE: scalar with the maximum duration of a vertical shake.
0112 %            Default value: 20
0113 %          INVERSION: scalar with the maximum depth inversion allowed during a
0114 %            cast.
0115 %            Default value: 3
0116 %          INTERRUPT: scalar with the maximum duration of stalled and/or shake
0117 %            intervals during a cast.
0118 %            Default value: 180
0119 %          LENGTH: scalar with the minimum depth range a cast must span.
0120 %            Default value: 10
0121 %          PERIOD: scalar with the minimum duration range a cast must last.
0122 %            Default value: 0
0123 %        Each struct in the array specifies a choice of inputs for cast
0124 %        boundary identification. The time sequence is optional and is
0125 %        relevant only if the SHAKE, INTERRUPT or PERIOD options are used.
0126 %        Identification will be performed with the first input choice whose
0127 %        input sequences are available. If no input choice is available,
0128 %        profiles are not identified.
0129 %        Default value: struct('depth', {'depth' 'depth_ctd' 'depth_ctd'}, ...
0130 %                              'time',  {'time'  'time_ctd'  'time'});
0131 %      PROFILE_MIN_RANGE: minimum depth range allowed for a valid profile.
0132 %        Non-negative real number setting the minimum depth range threshold for
0133 %        cast validation. If the difference between the maximum and the minimum
0134 %        depth of a valid reading in a cast is less than the given threshold,
0135 %        the cast will be discarded. Set it to 0 to prevent discarding any cast.
0136 %        Default value: 10
0137 %      PROFILE_MAX_GAP_RATIO: maximum gap ratio allowed for a valid profile.
0138 %        Real number (in [0..1]) setting the maximum gap ratio threshold for
0139 %        cast validation. A gap is a sequence of consecutive readings taken
0140 %        during a depth inversion or in which the value of at least one of the
0141 %        involved sensors is invalid. The gap ratio is the ratio of the depth
0142 %        range covered during the gap over the total depth covered of the cast.
0143 %        If the ratio of the largest gap over the total depth range is greater
0144 %        than the given threshold, the cast will be discarded.
0145 %        Set it to 1 to prevent discarding any cast.
0146 %        Default value: 0.8
0147 %      FLOW_CTD_LIST: CTD flow speed derivation input set choices.
0148 %        Struct array selecting input sequences for CTD flow speed derivation,
0149 %        in order of preference. It should have the following fields:
0150 %          TIME: time sequence name.
0151 %          DEPTH: depth sequence name.
0152 %          PITCH: pitch sequence name.
0153 %        Each struct in the array specifies a choice of inputs for CTD flow
0154 %        speed derivation. Pitch sequence is optional. If missing or empty,
0155 %        the nominal pitch value may be used (see FLOW_CTD_PITCH VALUE below).
0156 %        Derivation will be performed only if the casts are properly identified,
0157 %        and using the first input choice whose time and depth sequences are
0158 %        available, and either the pitch sequence is available or the pitch
0159 %        nominal value is set. If no input choice is available, CTD flow speed
0160 %        is not derived.
0161 %        Default value: struct('time', {}, 'depth', {}, 'pitch', {})
0162 %      FLOW_CTD_PITCH_VALUE: nominal pitch value for CTD flow speed derivation.
0163 %        Number with the nominal pitch value (radians) to use for CTD flow
0164 %        speed derivation when no pitch sequence is available.
0165 %        Default value: [] (no default pitch value)
0166 %      FLOW_CTD_MIN_PITCH: low pitch threshold for CTD flow derivation.
0167 %        Number with the minimum absolute pitch value below which flow speed is
0168 %        considered invalid during CTD flow speed derivation.
0169 %        Default value: 0 (all values are valid).
0170 %      FLOW_CTD_MIN_VELOCITY: low velocity threshold for CTD flow derivation.
0171 %        Number with the minimum absolute vertical velocity value below which
0172 %        flow speed is considered invalid during CTD flow speed derivation.
0173 %        Default value: 0 (all values are valid).
0174 %      SENSOR_LAG_LIST: sensor lag correction settings.
0175 %        Struct array specifying the sequences to produce by correcting
0176 %        the sensor lag in the corresponding original sensor sequences.
0177 %        It should have the following fields:
0178 %          CORRECTED: string with the name for the corrected sequence (field in
0179 %            struct DATA_PROC).
0180 %          ORIGINAL: string with the name of the original sequence (field in
0181 %            struct DATA_PROC).
0182 %          PARAMETERS: non-negative number as predefined time constant,
0183 %            or string 'auto' for automatic estimation from casts.
0184 %        It may have the following optional fields (empty or missing):
0185 %          TIME: string cell array with the names of the time sequence to use
0186 %            for estimation or correction, in order of preference.
0187 %            Default value: {'time'}
0188 %          DEPTH: string cell array with the names of the depth sequence to use
0189 %            for estimation or correction, in order of preference.
0190 %            Default value: {'depth'}
0191 %          FLOW: string cell array with the names of the flow sequence to use
0192 %            for estimation or correction, in order of preference. This is only
0193 %            used if flow speed is not constant.
0194 %            Default value: {'flow_ctd'}
0195 %          CONSTANT_FLOW: boolean setting whether parameters are static or
0196 %            dynamic (varying with the flow speed).
0197 %            Default value: false
0198 %          ESTIMATOR: function handle or string with the name of the estimator
0199 %            to combine the parameter estimates computed for each cast pair.
0200 %            Default value: @nanmedian
0201 %          MINOPTS: struct to pass custom minimization options for estimation,
0202 %            in the format accepted by function FINDSENSORLAGPARAMS.
0203 %            Default value: struct()
0204 %        Each struct in the struct array specifies a sensor lag correction.
0205 %        It will be performed only if the casts are properly identified,
0206 %        all the input sequences are available, and the correction parameter is
0207 %        available too (either given as option or estimated from pair of casts).
0208 %        Default value: struct('corrected',  {}, ...
0209 %                              'original',   {}, ...
0210 %                              'parameters', {})
0211 %      THERMAL_LAG_LIST: thermal lag correction settings.
0212 %        Struct array specifying the temperature and conductivity sequences
0213 %        to produce by correcting the thermal lag in the corresponding
0214 %        original sensor sequences.
0215 %        It should have the following fields:
0216 %          CONDUCTIVITY_CORRECTED: string with the name for the corrected
0217 %            conductivity sequence (field in DATA_PROC).
0218 %          TEMPERATURE_CORRECTED: string with the name for the corrected
0219 %            temperature sequence (field in DATA_PROC).
0220 %          CONDUCTIVITY_ORIGINAL: string with the name of the original
0221 %            conductivity sequence (field in DATA_PROC).
0222 %          TEMPERATURE_ORIGINAL: string with the name of the original
0223 %            temperature sequence (field in DATA_PROC).
0224 %          PRESSURE_ORIGINAL: string with the name of the original
0225 %            pressure sequence (field in DATA_PROC).
0226 %          PARAMETERS: numeric vector with predefined thermal lag parameters or
0227 %            string 'auto' for automatic estimation from casts. If a vector, it
0228 %            should be a 2 element array when flow speed is constant (error and
0229 %            error time), and a 4 element array otherwise (error offset, error
0230 %            slope, error time offset and error time slope).
0231 %        It may have the following optional fields (empty or missing):
0232 %          TIME: string cell array with the names of the time sequence to use
0233 %            for estimation or correction, in order of preference.
0234 %            Default value: {'time_ctd' 'time'}
0235 %          DEPTH: string cell array with the names of the depth sequence to use
0236 %            for estimation or correction, in order of preference.
0237 %            Depth is only used to ignore invalid profiles.
0238 %            Default value: {'depth_ctd' 'depth'}
0239 %          FLOW: string cell array with the names of the flow sequence to use
0240 %            for estimation or correction, in order of preference.
0241 %            This is only used if flow speed is not constant.
0242 %            Default value: {'flow_ctd'}
0243 %          CONSTANT_FLOW: boolean setting whether parameters are static or
0244 %            dynamic (varying with the flow speed).
0245 %            Default value: false
0246 %          ESTIMATOR: function handle or string with the name of the estimator
0247 %             to combine the parameter estimates computed for each cast pair.
0248 %            Default value: @nanmedian
0249 %          MINOPTS: struct to pass custom minimization options for estimation,
0250 %            in the format accepted by function FINDTHERMALLAGPARAMS.
0251 %            Default value: struct()
0252 %        Each struct in the struct array specifies a thermal lag correction.
0253 %        It will be performed only if casts are properly identified, all the
0254 %        input sequences are available, and the correction parameters are
0255 %        available too (either given as option or estimated from pair of casts).
0256 %        Default value: struct('conductivity_corrected', {'conductivity_corrected_thermal'}, ...
0257 %                              'temperature_corrected',  {'temperature_corrected_thermal'}, ...
0258 %                              'conductivity_original', {'conductivity'}, ...
0259 %                              'temperature_original', {'temperature'}, ...
0260 %                              'pressure_original', {'pressure'}, ...
0261 %                              'parameters', {'auto'} )
0262 %      SALINITY_LIST: salinity derivation settings.
0263 %        Struct cell array specifying which salinity sequences to produce
0264 %        by derivation from corresponding conductivity, temperature
0265 %        and pressure sequences. It should have the following fields:
0266 %          SALINITY: string with the name for the salinity sequence
0267 %            (field in DATA_PROC).
0268 %          CONDUCTIVITY: string with the name of the original conductivity
0269 %            sequence (field in DATA_PROC).
0270 %          TEMPERATURE: string with the name of the original temperature
0271 %            sequence (field in DATA_PROC).
0272 %          PRESSURE: string with the name of the original pressure sequence
0273 %            (field in DATA_PRE).
0274 %        Each struct in the struct array specifies a salinity derivation.
0275 %        It will be performed only if all the original sequences are available.
0276 %        Default value: struct('salinity',     {'salinity     salinity_corrected_thermal'}, ...
0277 %                              'conductivity', {'conductivity conductivity'}, ...
0278 %                              'temperature',  {'temperature  temperature_corrected_thermal'},
0279 %                              'pressure',     {'pressure' '  pressure});
0280 %      DENSITY_LIST: density derivation settings.
0281 %        Struct cell array specifying which salinity sequences to produce
0282 %        by derivation from corresponding salinity, temperature
0283 %        and pressure sequences. It should have the following fields:
0284 %          DENSITY: string with the name for the density sequence (field in
0285 %            DATA_PROC).
0286 %          SALINITY: string with the name of the original salinity sequence
0287 %            (field in DATA_PROC).
0288 %          TEMPERATURE: string with the name of the original temperature
0289 %            sequence (field in DATA_PROC).
0290 %          PRESSURE: string with the name of the original pressure sequence (field
0291 %            in DATA_PROC).
0292 %        Each struct in the struct array specifies a density derivation.
0293 %        It will be performed only if all the or'time'iginal sequences are available.
0294 %        Default value: struct('density',     {'density'     'density_corrected_thermal'}, ...
0295 %                              'salinity',    {'salinity'    'salinity_corrected_thermal'}, ...
0296 %                              'temperature', {'temperature' 'temperature'}, ...
0297 %                              'pressure',    {'pressure'    'pressure'})
0298 %
0299 %    The following options are deprecated and should not be used:
0300 %      PROFILING_SEQUENCE_LIST: sequence choices for cast identification.
0301 %        String cell array with the names of the pressure or depth sequence to
0302 %        use for cast identification, in order of preference.
0303 %        Deprecated in v1.1.0:
0304 %          Superseeded by PROFILING_LIST.
0305 %      PROFILE_JOIN_SAME_DIR: join consecutive profiles with the same direction.
0306 %        Boolean setting whether consecutive valid profiles with the same
0307 %        vertical direction should be joined into one.
0308 %        Deprecated in v1.1.0:
0309 %          Superseeded by PROFILING_LIST.
0310 %      PROFILING_SEQUENCE_FILLING: profiling sequence interpolation switch.
0311 %        Boolean setting whether the missing values in the profiling sequence
0312 %        should be filled by interpolation before cast identification.
0313 %        Deprecated in v1.1.0:
0314 %          Not needed anymore since profile identification is consistent and
0315 %          robust against missing values.
0316 %
0317 %  Notes:
0318 %    This function is based on the previous work by Tomeu Garau. He is the true
0319 %    glider man.
0320 %
0321 %  Examples:
0322 %    [data_proc, meta_proc] = processGliderData(data_pre, meta_pre, options)
0323 %
0324 %  See also:
0325 %    PREPROCESSGLIDERDATA
0326 %    FILLINVALIDVALUES
0327 %    FINDTRANSECTS
0328 %    COMPUTECUMULATIVEDISTANCE
0329 %    FINDPROFILES
0330 %    VALIDATEPROFILE
0331 %    APPLYSEABIRDPRESSUREFILTER
0332 %    FINDSENSORLAGPARAMS
0333 %    CORRECTSENSORLAG
0334 %    FINDTHERMALLAGPARAMS
0335 %    CORRECTTHERMALLAG
0336 %    SW_DPTH
0337 %    SW_SALT
0338 %    SW_DENS
0339 %
0340 %  Authors:
0341 %    Joan Pau Beltran  <joanpau.beltran@socib.cat>
0342 
0343 %  Copyright (C) 2013-2016
0344 %  ICTS SOCIB - Servei d'observacio i prediccio costaner de les Illes Balears
0345 %  <http://www.socib.es>
0346 %
0347 %  This program is free software: you can redistribute it and/or modify
0348 %  it under the terms of the GNU General Public License as published by
0349 %  the Free Software Foundation, either version 3 of the License, or
0350 %  (at your option) any later version.
0351 %
0352 %  This program is distributed in the hope that it will be useful,
0353 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
0354 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0355 %  GNU General Public License for more details.
0356 %
0357 %  You should have received a copy of the GNU General Public License
0358 %  along with this program.  If not, see <http://www.gnu.org/licenses/>.
0359   
0360   error(nargchk(2, 56, nargin, 'struct'));
0361   
0362   %% Configure default values for optional profile identification settings.
0363   default_profiling_time = [];
0364   default_profiling_stall = 3;
0365   default_profiling_shake = 20;
0366   default_profiling_inversion = 3;
0367   default_profiling_interrupt = 180;
0368   default_profiling_length = 10;
0369   default_profiling_period =  0;
0370   
0371   
0372   %% Configure default values for optional sensor lag settings.
0373   default_sensor_lag_time_list = {'time'};
0374   default_sensor_lag_depth_list = {'depth'};
0375   default_sensor_lag_flow_list = {'flow_ctd'};
0376   default_sensor_lag_flow_const = false;
0377   default_sensor_lag_estimator = @nanmedian;
0378   default_sensor_lag_minopts = struct();
0379   
0380   
0381   %% Configure default values for optional thermal lag settings.
0382   default_thermal_lag_time_list = {'time_ctd' 'time'};
0383   default_thermal_lag_depth_list = {'depth_ctd' 'depth'};
0384   default_thermal_lag_flow_list = {'flow_ctd'};
0385   default_thermal_lag_flow_const = false;
0386   default_thermal_lag_estimator = @nanmedian;
0387   default_thermal_lag_minopts = struct();
0388   
0389   
0390   %% Set processing options and default values.
0391   options = struct();
0392   options.time_filling = false;
0393   options.position_filling = false;
0394   options.depth_filling = false;
0395   options.attitude_filling = false;
0396   options.heading_filling = false;
0397   options.waypoint_filling = false;
0398   
0399   options.pressure_filtering = true;
0400   options.pressure_filter_constant = 4; % Recommended setting from Seabird Data Processing Manual.
0401   options.depth_ctd_derivation = true;
0402   
0403   options.profiling_list = ...
0404     struct('depth', {'depth' 'depth_ctd' 'depth_ctd'}, ...
0405            'time',  {'time'  'time_ctd'  'time'});
0406     
0407   options.profiling_sequence_list = [];    % Deprecated in v1.1.0.
0408   options.profiling_sequence_filling = []; % Deprecated in v1.1.0.
0409   options.profile_join_same_dir = [];      % Deprecated in v1.1.0.
0410   
0411   options.profile_min_range = 10;
0412   options.profile_max_gap_ratio = 0.8;
0413   
0414   options.flow_ctd_list = struct('time', {}, 'depth', {}, 'pitch', {});
0415   options.flow_ctd_pitch_value = []; % Before refactoring it was DEG2RAD(26).
0416   options.flow_ctd_min_pitch = 0;
0417   options.flow_ctd_min_velocity = 0;
0418   
0419   options.sensor_lag_list = ...
0420     struct('corrected', {}, 'original', {}, 'parameters', {});
0421   
0422   options.thermal_lag_list = ...
0423     struct('conductivity_corrected', {'conductivity_corrected_thermal'}, ...
0424            'temperature_corrected',  {'temperature_corrected_thermal'}, ...         
0425            'conductivity_original',  {'conductivity'}, ...
0426            'temperature_original',   {'temperature'}, ...
0427            'pressure_original',      {'pressure'}, ...
0428            'parameters',             {'auto'} );
0429   
0430   options.salinity_list = ...
0431     struct('salinity',     {'salinity'     'salinity_corrected_thermal'}, ...
0432            'conductivity', {'conductivity' 'conductivity'}, ...
0433            'temperature',  {'temperature'  'temperature_corrected_thermal'}, ...
0434            'pressure',     {'pressure'     'pressure'});
0435   
0436   options.density_list = ...
0437     struct('density',     {'density'     'density_corrected_thermal'}, ...
0438            'salinity',    {'salinity'    'salinity_corrected_thermal'}, ...
0439            'temperature', {'temperature' 'temperature'}, ...
0440            'pressure',    {'pressure'    'pressure'});
0441   
0442   
0443   %% Get options from extra arguments.
0444   % Parse option key-value pairs in any accepted call signature.
0445   if isscalar(varargin) && isstruct(varargin{1})
0446     % Options passed as a single option struct argument:
0447     % field names are option keys and field values are option values.
0448     option_key_list = fieldnames(varargin{1});
0449     option_val_list = struct2cell(varargin{1});
0450   elseif mod(numel(varargin), 2) == 0
0451     % Options passed as key-value argument pairs.
0452     option_key_list = varargin(1:2:end);
0453     option_val_list = varargin(2:2:end);
0454   else
0455     error('glider_toolbox:processGliderData:InvalidOptions', ...
0456           'Invalid optional arguments (neither key-value pairs nor struct).');
0457   end
0458   % Overwrite default options with values given in extra arguments.
0459   for opt_idx = 1:numel(option_key_list)
0460     opt = lower(option_key_list{opt_idx});
0461     val = option_val_list{opt_idx};
0462     if isfield(options, opt)
0463       options.(opt) = val;
0464     else
0465       error('glider_toolbox:processGliderData:InvalidOption', ...
0466             'Invalid option: %s.', opt);
0467     end
0468   end
0469   
0470   
0471   %% Handle deprecated options.
0472   % Handle deprecated options for profile identification.
0473   if ~isempty(options.profiling_sequence_list) ...
0474       || ~isempty(options.profiling_sequence_filling) ...
0475       || ~isempty(options.profile_join_same_dir)
0476     warning('glider_toolbox:findProfiles:DeprecatedOption', ...
0477             'Deprecated option: %s, %s and %s. See option: %s.', ...
0478             'profiling_sequence_list', 'profiling_sequence_filling', ...
0479             'profile_join_same_dir', 'profiling_list');
0480     options.profiling_list = ...
0481       struct('depth', options.profiling_sequence_list, ...
0482              'stall', options.profile_min_range, ...
0483              'interrupt', options.join * inf, ...
0484              'inversion', options.join * inf);
0485   end
0486   
0487   
0488   %% Initialize output variables.
0489   data_proc = data_pre;
0490   meta_proc = meta_pre;
0491   
0492   
0493   %% Fill missing time readings, if needed.
0494   % Regular sampling is assumed on time gaps.
0495   if options.time_filling && isfield(data_proc, 'time')
0496     fprintf('Filling missing time readings...\n');
0497     data_proc.time = ...
0498       fillInvalidValues(data_proc.time, 'linear');
0499     meta_proc.time.filling = 'linear';
0500   end
0501   
0502   
0503   %% Fill missing position readings, if needed.
0504   % Use linear interpolation of valid coordinate readings.
0505   if options.position_filling ...
0506       && all(isfield(data_proc, {'latitude' 'longitude'})) ...
0507       && any(isfield(data_proc, {'time' 'time_position'}))
0508     fprintf('Filling missing position readings...\n');
0509     if isfield(data_proc, 'time_position')
0510       data_proc.latitude = ...
0511         fillInvalidValues(data_proc.time_position, data_proc.latitude, ...
0512                           data_proc.time, nan(size(data_proc.time)), 'linear');
0513       data_proc.longitude = ...
0514         fillInvalidValues(data_proc.time_position, data_proc.longitude, ...
0515                           data_proc.time, nan(size(data_proc.time)), 'linear');
0516       meta_proc.latitude.sources = ...
0517         vertcat(cellstr(meta_proc.longitude.sources), ...
0518                 cellstr(meta_proc.time.sources));
0519       meta_proc.longitude.sources = ...
0520         vertcat(cellstr(meta_proc.longitude.sources), ...
0521                 cellstr(meta_proc.time.sources));
0522     else
0523       data_proc.latitude = ...
0524         fillInvalidValues(data_proc.time, data_proc.latitude, 'linear');
0525       data_proc.longitude = ...
0526         fillInvalidValues(data_proc.time, data_proc.longitude, 'linear');
0527     end
0528     meta_proc.latitude.filling = 'linear';
0529     meta_proc.longitude.filling = 'linear';
0530   end
0531   
0532   
0533   %% Fill missing depth readings, if needed.
0534   % Use linear interpolation of valid depth readings.
0535   if options.depth_filling && all(isfield(data_proc, {'time' 'depth'}))
0536     fprintf('Filling missing depth readings...\n');
0537     data_proc.depth = ...
0538       fillInvalidValues(data_proc.time, data_proc.depth, 'linear');
0539     meta_proc.depth.filling = 'linear';
0540   end
0541   
0542   
0543   %% Fill missing attitude readings, if needed.
0544   % Use linear interpolation of valid roll and pitch readings.
0545   if options.attitude_filling && ...
0546       all(isfield(data_proc, {'time' 'roll' 'pitch'}))
0547     fprintf('Filling missing attitude readings...\n');
0548     data_proc.roll = ...
0549       fillInvalidValues(data_proc.time, data_proc.roll, 'linear');
0550     data_proc.pitch = ...
0551       fillInvalidValues(data_proc.time, data_proc.pitch, 'linear');
0552     meta_proc.roll.filling = 'linear';
0553     meta_proc.pitch.filling = 'linear';
0554   end
0555   
0556   
0557   %% Fill missing heading readings, if needed.
0558   % Use linear interpolation of valid coordinate readings.
0559   if options.heading_filling && all(isfield(data_proc, {'time' 'heading'}))
0560     fprintf('Filling missing heading readings...\n');
0561     data_proc.heading = ...
0562       fillInvalidValues(data_proc.time, data_proc.heading, 'linear');
0563     meta_proc.heading.filling = 'linear';
0564   end
0565   
0566   
0567   %% Fill missing waypoint coordinate readings, if needed.
0568   % Waypoint coordinates are assumed constant until next valid waypoint
0569   % coordinate reading.
0570   if options.waypoint_filling ...
0571       && all(isfield(data_proc, {'waypoint_latitude' 'waypoint_longitude'}))
0572     fprintf('Filling missing commanded waypoint readings...\n');
0573     data_proc.waypoint_latitude = ...
0574       fillInvalidValues(data_proc.waypoint_latitude, 'prev');
0575     data_proc.waypoint_longitude = ...
0576       fillInvalidValues(data_proc.waypoint_longitude, 'prev');
0577     meta_proc.waypoint_latitude.filling = 'prev';
0578     meta_proc.waypoint_longitude.filling = 'prev';
0579   end
0580   
0581   
0582   %% Identify transect boundaries, if waypoint coordinates available.
0583   if all(isfield(data_proc, {'waypoint_latitude' 'waypoint_longitude'}))
0584     fprintf('Computing transect index by transect boundary identification...\n');
0585     data_proc.transect_index = ...
0586       findTransects(data_proc.waypoint_latitude, data_proc.waypoint_longitude);
0587     meta_proc.transect_index.sources = ...
0588       {'waypoint_latitude' 'waypoint_longitude'}';
0589     meta_proc.transect_index.method = 'findTransects';
0590   end
0591   
0592   
0593   %% Compute navigated distance over ground.
0594   if all(isfield(data_proc, {'latitude' 'longitude'}))
0595     fprintf('Computing covered horizontal distance...\n');
0596     data_proc.distance_over_ground = ...
0597       computeCumulativeDistance(data_proc.latitude, data_proc.longitude);
0598     meta_proc.distance_over_ground.sources = {'latitude' 'longitude'}';
0599     meta_proc.distance_over_ground.method = 'computeCumulativeDistance';
0600   end
0601   
0602   
0603   %% Convert and filter pressure, if pressure available and needed.
0604   if isfield(data_proc, 'pressure')
0605     % Apply pressure filter, if needed.
0606     if options.pressure_filtering
0607       if isfield(data_proc, 'time_ctd')
0608         fprintf('Filtering pressure sequence using CTD time stamp...\n');
0609         data_proc.pressure = ...
0610           applySeabirdPressureFilter(data_proc.time_ctd, data_proc.pressure, ...
0611                                      options.pressure_filter_constant);
0612         meta_proc.pressure.sources = ...
0613           vertcat(cellstr(meta_proc.pressure.sources), ...
0614                   cellstr(meta_proc.time_ctd.sources));
0615       elseif isfield(data_proc, 'time')
0616         fprintf('Filtering pressure sequence using global time stamp...\n');
0617         data_proc.pressure = ...
0618           applySeabirdPressureFilter(data_proc.time, data_proc.pressure, ...
0619                                      options.pressure_filter_constant);
0620         meta_proc.pressure.sources = ...
0621           vertcat(cellstr(meta_proc.pressure.sources), ...
0622                   cellstr(meta_proc.time.sources));
0623       end
0624       meta_proc.pressure.filter_method = 'applySeabirdPressureFilter';
0625       meta_proc.pressure.filter_parameters = options.pressure_filter_constant;
0626     end
0627   end
0628   
0629   
0630   %% Derive depth from pressure, if pressure available and needed.
0631   if options.depth_ctd_derivation ...
0632       && all(isfield(data_proc, {'pressure' 'latitude'}))
0633     fprintf('Deriving CTD depth from pressure and latitude readings...\n');
0634     data_proc.depth_ctd = sw_dpth(data_proc.pressure, data_proc.latitude);
0635     meta_proc.depth_ctd.sources = {'pressure' 'latitude'}';
0636     meta_proc.depth_ctd.method = 'sw_depth';
0637   end
0638   
0639   
0640   %% Identify start and end of profiles.
0641   % Find preferred vertical coordinate sequence
0642   % (e.g. navigation depth, CTD-derived depth, pressure...)
0643   % with optional duration sequence (navigation time, CTD timestamp, ...)
0644   % and identification parameters.
0645   profiling_avail = false;
0646   for profiling_option_idx = 1:numel(options.profiling_list)
0647     profiling_option = options.profiling_list(profiling_option_idx);
0648     profiling_depth = profiling_option.depth;
0649     profiling_time = default_profiling_time;
0650     profiling_stall = default_profiling_stall;
0651     profiling_shake = default_profiling_shake;
0652     profiling_inversion = default_profiling_inversion;
0653     profiling_interrupt = default_profiling_interrupt;
0654     profiling_length = default_profiling_length;
0655     profiling_period = default_profiling_period;
0656     if isfield(profiling_option, 'time') && ~isempty(profiling_option.time)
0657       profiling_time = profiling_option.time;
0658     end
0659     if isfield(profiling_option, 'stall') && ~isempty(profiling_option.stall)
0660       profiling_stall = profiling_option.stall;
0661     end
0662     if isfield(profiling_option, 'shake') && ~isempty(profiling_option.shake)
0663       profiling_shake = profiling_option.shake;
0664     end
0665     if isfield(profiling_option, 'inversion') ...
0666         && ~isempty(profiling_option.inversion)
0667       profiling_inversion = profiling_option.inversion;
0668     end
0669     if isfield(profiling_option, 'interrupt') ...
0670         && ~isempty(profiling_option.interrupt)
0671       profiling_interrupt = profiling_option.interrupt;
0672     end
0673     if isfield(profiling_option, 'length') && ~isempty(profiling_option.length)
0674       profiling_length = profiling_option.legnth;
0675     end
0676     if isfield(profiling_option, 'preiod') && ~isempty(profiling_option.period)
0677       profiling_period = profiling_option.period;
0678     end
0679     profiling_depth_avail = false;
0680     profiling_time_avail = false;
0681     if isfield(data_proc, profiling_depth) ...
0682         && ~all(isnan(data_proc.(profiling_depth)))
0683       profiling_depth_avail = true;
0684     end
0685     if isfield(data_proc, profiling_time) ...
0686         && any(data_proc.(profiling_time) > 0)
0687       profiling_time_avail = true;
0688     end
0689     if profiling_depth_avail
0690       if isempty(profiling_time)
0691         profiling_vars = {data_proc.(profiling_depth)};
0692         profiling_avail = true;
0693         break
0694       elseif profiling_time_avail
0695         profiling_vars = ...
0696           {data_proc.(profiling_time) data_proc.(profiling_depth)};
0697         profiling_avail = true;
0698         break
0699       end
0700     end
0701   end
0702   % Compute profile boundaries and direction if profiling sequence available.
0703   if profiling_avail
0704     fprintf('Computing vertical direction and profile index with settings:\n');
0705     fprintf('  vertical   sequence: %s\n', profiling_depth);
0706     if profiling_time_avail
0707       fprintf('  horizontal sequence: %s\n', profiling_time);
0708     end
0709     fprintf('  maximum profile stall     : %f\n', profiling_stall);
0710     fprintf('  maximum profile shake     : %f\n', profiling_shake);
0711     fprintf('  maximum profile inversion : %f\n', profiling_inversion);
0712     fprintf('  maximum profile interrupt : %f\n', profiling_interrupt);
0713     fprintf('  minimum profile length    : %f\n', profiling_length);
0714     fprintf('  minimum profile period    : %f\n', profiling_period);
0715     % Find profile directions and indices.
0716     [data_proc.profile_index, data_proc.profile_direction] = findProfiles( ...
0717         profiling_vars{:}, ...
0718         'stall', profiling_stall, 'shake', profiling_shake, ...
0719         'inversion', profiling_inversion, 'interrupt', profiling_interrupt, ...
0720         'length', profiling_length, 'period', profiling_period);
0721     meta_proc.profile_index.method = 'findProfiles';
0722     if profiling_time_avail
0723       meta_proc.profile_index.sources = {profiling_time profiling_depth};
0724     else
0725       meta_proc.profile_index.sources = {profiling_depth};
0726     end
0727     meta_proc.profile_index.shake = profiling_shake;
0728     meta_proc.profile_index.stall = profiling_stall;
0729     meta_proc.profile_index.inversion = profiling_inversion;
0730     meta_proc.profile_index.interrupt = profiling_interrupt;
0731     meta_proc.profile_index.length = profiling_length;
0732     meta_proc.profile_index.period = profiling_period;
0733     meta_proc.profile_direction.method = meta_proc.profile_index.method;
0734     meta_proc.profile_direction.sources = meta_proc.profile_index.sources;
0735     meta_proc.profile_direction.shake = meta_proc.profile_index.shake;
0736     meta_proc.profile_direction.stall = meta_proc.profile_index.stall;
0737     meta_proc.profile_direction.inversion = meta_proc.profile_index.inversion;
0738     meta_proc.profile_direction.interrupt = meta_proc.profile_index.interrupt;
0739     meta_proc.profile_direction.length = meta_proc.profile_index.length;
0740     meta_proc.profile_direction.period = meta_proc.profile_index.period;
0741   end
0742   
0743   
0744   %% Derive flow speed through CTD cell, if needed and data available.
0745   % Time and depth sequences must be present in already processed data.
0746   % Pitch may be also a sequence in processed data (preferred)
0747   % or a default pitch value when pitch sequence is not available.
0748   flow_ctd_avail = false;
0749   flow_ctd_pitch_value_avail = ~isempty(options.flow_ctd_pitch_value);
0750   for flow_ctd_option_idx = 1:numel(options.flow_ctd_list)
0751     flow_ctd_option = options.flow_ctd_list(flow_ctd_option_idx);
0752     flow_ctd_time_avail = false;
0753     flow_ctd_depth_avail = false;
0754     flow_ctd_pitch_avail = false;
0755     if isfield(data_proc, flow_ctd_option.time) ...
0756         && any(data_proc.(flow_ctd_option.time) > 0)
0757       flow_ctd_time = flow_ctd_option.time;
0758       flow_ctd_time_avail = true;
0759     end
0760     if isfield(data_proc, flow_ctd_option.depth) ...
0761         && ~all(isnan(data_proc.(flow_ctd_option.depth)))
0762       flow_ctd_depth = flow_ctd_option.depth;
0763       flow_ctd_depth_avail = true;
0764     end
0765     if isfield(flow_ctd_option, 'pitch') ...
0766         && isfield(data_proc, flow_ctd_option.pitch) ...
0767         && ~all(isnan(data_proc.(flow_ctd_option.pitch)))
0768       flow_ctd_pitch = flow_ctd_option.pitch;
0769       flow_ctd_pitch_avail = true;
0770     end
0771     if flow_ctd_time_avail && flow_ctd_depth_avail ...
0772         && (flow_ctd_pitch_avail || flow_ctd_pitch_value_avail)
0773       flow_ctd_avail = true;
0774       break
0775     end
0776   end
0777   flow_ctd_prof_avail = isfield(data_proc, 'profile_index');
0778   if flow_ctd_prof_avail && flow_ctd_avail
0779     fprintf('Deriving CTD flow speed with settings:\n');
0780     fprintf('  depth sequence: %s\n', flow_ctd_depth);
0781     fprintf('  time  sequence: %s\n', flow_ctd_time);
0782     if flow_ctd_pitch_avail
0783       fprintf('  pitch sequence: %s\n', flow_ctd_pitch);
0784     else
0785       fprintf('  pitch value   : %f\n', options.flow_ctd_pitch_value);
0786     end
0787     fprintf('  pitch minimum threshold            : %f\n', options.flow_ctd_min_pitch);
0788     fprintf('  vertical velocity minimum threshold: %f\n', options.flow_ctd_min_velocity);
0789     data_proc.flow_ctd = nan(size(data_proc.time));
0790     num_profiles = fix(max(data_proc.profile_index));
0791     for profile_idx = 1:num_profiles
0792       prof_select = (data_proc.profile_index == profile_idx);
0793       prof_time = data_proc.(flow_ctd_time)(prof_select);
0794       prof_depth = data_proc.(flow_ctd_depth)(prof_select);
0795       if flow_ctd_pitch_avail
0796         prof_pitch = data_proc.(flow_ctd_pitch)(prof_select);
0797         prof_vars = {prof_time(:) prof_pitch(:)};
0798       else
0799         prof_pitch = options.flow_ctd_pitch_value;
0800         prof_vars = {prof_time(:)};
0801       end
0802       [prof_valid, ~] = validateProfile(prof_depth(:), prof_vars{:}, ...
0803                                         'range', options.profile_min_range, ...
0804                                         'gap', options.profile_max_gap_ratio);
0805       if prof_valid
0806         data_proc.flow_ctd(prof_select) = ...
0807           computeCTDFlowSpeed(prof_time, prof_depth, prof_pitch, ...
0808                               'minpitch', options.flow_ctd_min_pitch, ...
0809                               'minvel', options.flow_ctd_min_velocity);
0810       end
0811     end
0812     if flow_ctd_pitch_avail
0813       meta_proc.flow_ctd.sources = ...
0814         {flow_ctd_time flow_ctd_depth flow_ctd_pitch 'profile_index'}';
0815     else
0816       meta_proc.flow_ctd.sources = ...
0817         {flow_ctd_time flow_ctd_depth 'profile_index'};
0818       meta_proc.flow_ctd.pitch_value = options.flow_ctd_pitch_value;
0819     end
0820     meta_proc.flow_ctd.method = 'computeCTDFlowSpeed';
0821     meta_proc.flow_ctd.min_pitch = options.flow_ctd_min_pitch;
0822     meta_proc.flow_ctd.min_vel = options.flow_ctd_min_velocity;
0823   end
0824  
0825   
0826   %% Perform sensor lag estimation and correction, if needed.
0827   % Sensor, time and depth sequences must be present in already processed data.
0828   for sensor_lag_option_idx = 1:numel(options.sensor_lag_list)
0829     % Get sensor lag arguments, setting options to default values if needed.
0830     % Name of corrected sequence must be specified in option.
0831     % Name of original sequence must be specified too.
0832     % Name of time and depth sequences may be specified as list of choices,
0833     % defaulted if missing.
0834     % Time constant may be given in option, or estimated from cast pairs.
0835     sensor_lag_option = options.sensor_lag_list(sensor_lag_option_idx);
0836     sensor_lag_cor = sensor_lag_option.corrected;
0837     sensor_lag_raw = sensor_lag_option.original;
0838     sensor_lag_params = sensor_lag_option.parameters;
0839     sensor_lag_time_list = default_sensor_lag_time_list;
0840     sensor_lag_depth_list = default_sensor_lag_depth_list;
0841     sensor_lag_flow_list = default_sensor_lag_flow_list;
0842     sensor_lag_flow_const = default_sensor_lag_flow_const;
0843     sensor_lag_estimator = default_sensor_lag_estimator;
0844     sensor_lag_minopts = default_sensor_lag_minopts;
0845     if isfield(sensor_lag_option, 'time') && ~isempty(sensor_lag_option.time)
0846       sensor_lag_time_list = sensor_lag_option.time;
0847     end
0848     if isfield(sensor_lag_option, 'depth') && ~isempty(sensor_lag_option.depth)
0849       sensor_lag_depth_list = sensor_lag_option.depth;
0850     end
0851     if isfield(sensor_lag_option, 'flow') ...
0852         && ~isempty(sensor_lag_option.flow)
0853       sensor_lag_flow_list = cellstr(sensor_lag_option.flow);
0854     end
0855     if isfield(sensor_lag_option, 'constant_flow') ...
0856         && ~isempty(sensor_lag_option.constant_flow)
0857       sensor_lag_flow_const = sensor_lag_option.constant_flow;
0858     end
0859     if isfield(sensor_lag_option, 'estimator') ...
0860         && ~isempty(sensor_lag_option.estimator)
0861       % Convert estimator function name string to function handle, if needed.
0862       if ischar(sensor_lag_option.estimator) 
0863         sensor_lag_estimator = str2func(sensor_lag_option.estimator);
0864       else
0865         sensor_lag_estimator = sensor_lag_option.estimator;
0866       end
0867     end
0868     if isfield(sensor_lag_option, 'minopts') ...
0869         && ~isempty(sensor_lag_option.minopts)
0870       sensor_lag_minopts = sensor_lag_option.minopts;
0871     end
0872     % Check if parameters are given or need to be estimated.
0873     sensor_lag_num_params = 2;
0874     if sensor_lag_flow_const
0875       sensor_lag_num_params = 1;
0876     end
0877     if isnumeric(sensor_lag_params) && ...
0878         (numel(sensor_lag_params) == sensor_lag_num_params)
0879       % Sensor lag parameters preset.
0880       sensor_lag_params_avail = true;
0881     elseif strcmpi(sensor_lag_params, 'auto')
0882       % Sensor lag parameter estimation requested.
0883       sensor_lag_params_avail = false;
0884     else
0885       % Invalid sensor lag parameter specification.
0886       error('glider_toolbox:processGliderData:InvalidSensorLagParam', ...
0887             'Invalid sensor lag settings %d: bad parameter specification.', ...
0888             sensor_lag_option_idx);
0889     end
0890     % Find input fields needed for sensor lag estimation or correction.
0891     sensor_lag_prof_avail = isfield(data_proc, 'profile_index');
0892     sensor_lag_raw_avail = false;
0893     sensor_lag_time_avail = false;
0894     sensor_lag_depth_avail = false;
0895     sensor_lag_flow_avail = false;
0896     if isfield(data_proc, sensor_lag_raw) ...
0897         && ~all(isnan(data_proc.(sensor_lag_raw)))
0898       sensor_lag_raw_avail = true;
0899     end
0900     for sensor_lag_time_idx = 1:numel(sensor_lag_time_list)
0901       sensor_lag_time = sensor_lag_time_list{sensor_lag_time_idx};
0902       if isfield(data_proc, sensor_lag_time) ...
0903           && any(data_proc.(sensor_lag_time) > 0)
0904         sensor_lag_time_avail = true;
0905         break
0906       end
0907     end
0908     for sensor_lag_depth_idx = 1:numel(sensor_lag_depth_list)
0909       sensor_lag_depth = sensor_lag_depth_list{sensor_lag_depth_idx};
0910       if isfield(data_proc, sensor_lag_depth)
0911         sensor_lag_depth_avail = true;
0912         break
0913       end
0914     end
0915     for sensor_lag_flow_idx = 1:numel(sensor_lag_flow_list)
0916       sensor_lag_flow = sensor_lag_flow_list{sensor_lag_flow_idx};
0917       if isfield(data_proc, sensor_lag_flow) ...
0918           && ~all(isnan(data_proc.(sensor_lag_flow)))
0919         sensor_lag_flow_avail = true;
0920         break
0921       end
0922     end
0923     % Perform sensor lag correction if needed input fields are there.
0924     if sensor_lag_prof_avail && sensor_lag_raw_avail ...
0925         && sensor_lag_time_avail && sensor_lag_depth_avail ...
0926         && (sensor_lag_flow_const || sensor_lag_flow_avail)
0927       num_profiles = fix(max(data_proc.profile_index));
0928       % Estimate sensor lag time constant, if needed.
0929       if sensor_lag_params_avail
0930         % Sensor lag time constant given (do not perform estimation).
0931         sensor_lag_constants = sensor_lag_params;
0932       else
0933         fprintf('Performing sensor lag parameter estimation %d with settings:\n', sensor_lag_option_idx);
0934         fprintf('  sensor sequence    : %s\n', sensor_lag_raw);
0935         fprintf('  time sequence      : %s\n', sensor_lag_time);
0936         fprintf('  depth sequence     : %s\n', sensor_lag_depth);
0937         if ~sensor_lag_flow_const
0938           fprintf('  flow speed sequence: %s\n', sensor_lag_flow);
0939         end
0940         fprintf('  estimator          : %s\n', func2str(sensor_lag_estimator));
0941         % Estimate sensor lag time constant for each pofile.
0942         sensor_lag_estimates = nan(num_profiles-1, sensor_lag_num_params);
0943         sensor_lag_exitflags = nan(num_profiles-1, 1);
0944         sensor_lag_residuals = nan(num_profiles-1, 1);
0945         for profile_idx = 1:(num_profiles-1)
0946           prof1_select = (data_proc.profile_index == profile_idx);
0947           [~, ~, prof1_dir] = ...
0948             find(data_proc.profile_direction(prof1_select), 1);
0949           prof1_raw = data_proc.(sensor_lag_raw)(prof1_select);
0950           prof1_time = data_proc.(sensor_lag_time)(prof1_select);
0951           prof1_depth = data_proc.(sensor_lag_depth)(prof1_select);
0952           prof2_select = (data_proc.profile_index == profile_idx + 1);
0953           [~, ~, prof2_dir] = ...
0954             find(data_proc.profile_direction(prof2_select), 1);
0955           prof2_raw = data_proc.(sensor_lag_raw)(prof2_select);
0956           prof2_time = data_proc.(sensor_lag_time)(prof2_select);
0957           prof2_depth = data_proc.(sensor_lag_depth)(prof2_select);
0958           if sensor_lag_flow_const
0959             prof1_vars = {prof1_time(:) prof1_depth(:) prof1_raw(:)};
0960             prof2_vars = {prof2_time(:) prof2_depth(:) prof2_raw(:)};
0961           else
0962             prof1_flow = data_proc.(sensor_lag_flow)(prof1_select);
0963             prof1_vars = ...
0964               {prof1_time(:) prof1_depth(:) prof1_raw(:) prof1_flow(:)};
0965             prof2_flow = data_proc.(sensor_lag_flow)(prof2_select);
0966             prof2_vars = ...
0967               {prof2_time(:) prof2_depth(:) prof2_raw(:) prof2_flow(:)};
0968           end
0969           [prof1_valid, ~, prof1_vars{:}] = ...
0970             validateProfile(prof1_depth(:), prof1_vars{:}, ...
0971                             'range', options.profile_min_range, ...
0972                             'gap', options.profile_max_gap_ratio);
0973           [prof2_valid, ~, prof2_vars{:}] = ...
0974             validateProfile(prof2_depth(:), prof2_vars{:}, ...
0975                             'range', options.profile_min_range, ...
0976                             'gap', options.profile_max_gap_ratio);
0977           prof_opposite_dir = (prof1_dir * prof2_dir < 0);
0978           if prof1_valid && prof2_valid && prof_opposite_dir
0979             try
0980               [sensor_lag_estimates(profile_idx, :), ...
0981                sensor_lag_exitflags(profile_idx), ...
0982                sensor_lag_residuals(profile_idx)] = ...
0983                 findSensorLagParams(prof1_vars{:}, prof2_vars{:}, sensor_lag_minopts);
0984               if sensor_lag_exitflags(profile_idx) <= 0
0985                  warning('glider_toolbox:processGliderData:SensorLagMinimizationError', ...
0986                          'Minimization did not converge for casts %d and %d, residual area: %f.', ...
0987                          profile_idx, profile_idx+1, sensor_lag_residuals(profile_idx));
0988               end
0989             catch exception
0990               fprintf('Sensor lag estimation failed for casts: %d and %d.\n', ...
0991                       profile_idx, profile_idx+1);
0992               disp(getReport(exception, 'extended'));
0993             end
0994           end
0995         end
0996         % Compute statistical estimate from individual profile estimates.
0997         sensor_lag_constants = sensor_lag_estimator(sensor_lag_estimates);
0998       end
0999       % Correct sensor lag, if possible.
1000       if any(isnan(sensor_lag_constants))
1001         fprintf('Omiting sensor lag correction %d (%s): %s.\n', ...
1002                 sensor_lag_option_idx, sensor_lag_cor, ...
1003                 'no valid parameters available');
1004       else
1005         fprintf('Performing sensor lag correction %d with settings:\n', sensor_lag_option_idx);
1006         fprintf('  output sensor sequence: %s\n', sensor_lag_cor);
1007         fprintf('  input sensor sequence : %s\n', sensor_lag_raw);
1008         fprintf('  input time sequence   : %s\n', sensor_lag_time);
1009         fprintf('  input depth sequence  : %s\n', sensor_lag_depth);
1010         if sensor_lag_flow_const
1011           fprintf('  parameters            : %f\n', sensor_lag_constants);
1012         else
1013           fprintf('  input flow sequence   : %s\n', sensor_lag_flow);
1014           fprintf('  parameters            : %f %f\n', sensor_lag_constants);
1015         end
1016         data_proc.(sensor_lag_cor) = nan(size(data_proc.(sensor_lag_raw)));
1017         for profile_idx = 1:num_profiles
1018           prof_select = (data_proc.profile_index == profile_idx);
1019           prof_raw = data_proc.(sensor_lag_raw)(prof_select);
1020           prof_time = data_proc.(sensor_lag_time)(prof_select);
1021           prof_depth = data_proc.(sensor_lag_depth)(prof_select);
1022           if sensor_lag_flow_const
1023             prof_vars = {prof_time(:) prof_raw(:)};
1024           else
1025             prof_flow = data_proc.(sensor_lag_flow)(prof_select);
1026             prof_vars = {prof_time(:) prof_raw(:) prof_flow(:)};
1027           end
1028           [prof_valid, ~, prof_vars{:}] = ...
1029             validateProfile(prof_depth(:), prof_vars{:}, ...
1030                             'range', options.profile_min_range, ...
1031                             'gap', options.profile_max_gap_ratio);
1032           if prof_valid
1033             prof_cor = correctSensorLag(prof_vars{:}, sensor_lag_constants);
1034             data_proc.(sensor_lag_cor)(prof_select) = prof_cor;
1035           end
1036         end
1037         if sensor_lag_flow_const
1038           meta_proc.(sensor_lag_cor).sources = ...
1039             {sensor_lag_raw sensor_lag_time sensor_lag_depth 'profile_index'}';
1040         else
1041           meta_proc.(sensor_lag_cor).sources = ...
1042             {sensor_lag_raw sensor_lag_time sensor_lag_depth sensor_lag_flow 'profile_index'}';
1043         end
1044         meta_proc.(sensor_lag_cor).method = 'correctSensorLag';
1045         meta_proc.(sensor_lag_cor).parameters = sensor_lag_constants;
1046         if sensor_lag_params_avail
1047           meta_proc.(sensor_lag_cor).parameter_method = 'preset';
1048         else
1049           meta_proc.(sensor_lag_cor).parameter_method = 'findSensorLagParams';
1050           meta_proc.(sensor_lag_cor).parameter_estimator = func2str(sensor_lag_estimator);
1051           meta_proc.(sensor_lag_cor).parameter_estimates = sensor_lag_estimates;
1052           meta_proc.(sensor_lag_cor).parameter_exitflags = sensor_lag_exitflags;
1053         end
1054         meta_proc.(sensor_lag_cor).profile_min_range = options.profile_min_range;
1055         meta_proc.(sensor_lag_cor).profile_gap_ratio = options.profile_max_gap_ratio;
1056       end
1057     end
1058   end
1059   
1060   
1061   %% Perform thermal lag estimation and correction, if needed.
1062   % Conductivity, temperature, pressure, and time sequences must be present in
1063   % already processed data. CTD flow speed sequence is also required for
1064   % non-constant flow (unpumped) CTDs.
1065   for thermal_lag_option_idx = 1:numel(options.thermal_lag_list)
1066     % Get thermal lag arguments, setting options to default values if needed.
1067     % Name of corrected conductivity and temperature sequences must be specified in option.
1068     % Name of original conductivity and temperature sequences must be specified too.
1069     % Name of flow speed sequence only needed when non-constant flow.
1070     % Name of time, and flow sequences may be specified as list of choices, defaulted if missing.
1071     % Parameters may be given in option, or estimated from cast pairs.
1072     thermal_lag_option = options.thermal_lag_list(thermal_lag_option_idx);
1073     thermal_lag_cond_cor = thermal_lag_option.conductivity_corrected;
1074     thermal_lag_temp_cor = thermal_lag_option.temperature_corrected;
1075     thermal_lag_cond_raw = thermal_lag_option.conductivity_original;
1076     thermal_lag_temp_raw = thermal_lag_option.temperature_original;
1077     thermal_lag_pres_raw = thermal_lag_option.pressure_original;
1078     thermal_lag_params = thermal_lag_option.parameters;
1079     thermal_lag_time_list = default_thermal_lag_time_list;
1080     thermal_lag_depth_list = default_thermal_lag_depth_list;
1081     thermal_lag_flow_list = default_thermal_lag_flow_list;
1082     thermal_lag_flow_const = default_thermal_lag_flow_const;
1083     thermal_lag_estimator = default_thermal_lag_estimator;
1084     thermal_lag_minopts = default_thermal_lag_minopts;
1085     if isfield(thermal_lag_option, 'time') ...
1086         && ~isempty(thermal_lag_option.time)
1087       thermal_lag_time_list = cellstr(thermal_lag_option.time);
1088     end
1089     if isfield(thermal_lag_option, 'depth') ...
1090         && ~isempty(thermal_lag_option.depth)
1091       thermal_lag_depth_list = cellstr(thermal_lag_option.depth);
1092     end
1093     if isfield(thermal_lag_option, 'flow') ...
1094         && ~isempty(thermal_lag_option.flow)
1095       thermal_lag_flow_list = cellstr(thermal_lag_option.flow);
1096     end
1097     if isfield(thermal_lag_option, 'constant_flow') ...
1098         && ~isempty(thermal_lag_option.constant_flow)
1099       thermal_lag_flow_const = thermal_lag_option.constant_flow;
1100     end
1101     if isfield(thermal_lag_option, 'estimator');
1102       % Convert estimator function name string to function handle, if needed.
1103       if ischar(thermal_lag_option.estimator) 
1104         thermal_lag_estimator = str2func(thermal_lag_option.estimator);
1105       else
1106         thermal_lag_estimator = thermal_lag_option.estimator;
1107       end
1108     end
1109     if isfield(thermal_lag_option, 'minopts');
1110       thermal_lag_minopts = thermal_lag_option.minopts;
1111     end
1112     % Check if parameters are given or need to be estimated.
1113     thermal_lag_num_params = 4;
1114     if thermal_lag_flow_const
1115       thermal_lag_num_params = 2;
1116     end
1117     if isnumeric(thermal_lag_params) ...
1118         && (numel(thermal_lag_params) == thermal_lag_num_params)
1119       % Thermal lag parameters preset.
1120       thermal_lag_params_avail = true;
1121     elseif strcmpi(thermal_lag_params, 'auto')
1122       % Thermal lag parameter estimation requested.
1123       thermal_lag_params_avail = false;
1124     else
1125       % Invalid thermal lag parameters given.
1126       error('glider_toolbox:processGliderData:InvalidThermalLagParam', ...
1127             'Invalid thermal lag settings %d: bad parameter specification.', ...
1128             thermal_lag_option_idx);
1129     end
1130     % Find input fields needed for thermal lag estimation or correction.
1131     thermal_lag_prof_avail = isfield(data_proc, 'profile_index');
1132     thermal_lag_cond_raw_avail = false;
1133     thermal_lag_temp_raw_avail = false;
1134     thermal_lag_pres_avail = false;
1135     thermal_lag_time_avail = false;
1136     thermal_lag_depth_avail = false;
1137     thermal_lag_flow_avail = false;
1138     if isfield(data_proc, thermal_lag_cond_raw) ...
1139         && ~all(isnan(data_proc.(thermal_lag_cond_raw))) 
1140       thermal_lag_cond_raw_avail = true;
1141     end
1142     if isfield(data_proc, thermal_lag_temp_raw) ...
1143         && ~all(isnan(data_proc.(thermal_lag_temp_raw)))
1144       thermal_lag_temp_raw_avail = true;
1145     end
1146     if isfield(data_proc, thermal_lag_pres_raw) ...
1147         && ~all(isnan(data_proc.(thermal_lag_pres_raw)))
1148       thermal_lag_pres_avail = true;
1149     end
1150     for thermal_lag_time_idx = 1:numel(thermal_lag_time_list)
1151       thermal_lag_time = thermal_lag_time_list{thermal_lag_time_idx};
1152       if isfield(data_proc, thermal_lag_time) ...
1153           && any(data_proc.(thermal_lag_time) > 0)
1154         thermal_lag_time_avail = true;
1155         break
1156       end
1157     end
1158     for thermal_lag_depth_idx = 1:numel(thermal_lag_depth_list)
1159       thermal_lag_depth = thermal_lag_depth_list{thermal_lag_depth_idx};
1160       if isfield(data_proc, thermal_lag_depth) ...
1161           && ~all(isnan(data_proc.(thermal_lag_depth)))
1162         thermal_lag_depth_avail = true;
1163         break
1164       end
1165     end
1166     for thermal_lag_flow_idx = 1:numel(thermal_lag_flow_list)
1167       thermal_lag_flow = thermal_lag_flow_list{thermal_lag_flow_idx};
1168       if isfield(data_proc, thermal_lag_flow) ...
1169           && ~all(isnan(data_proc.(thermal_lag_flow)))
1170         thermal_lag_flow_avail = true;
1171         break
1172       end
1173     end
1174     % Perform thermal lag correction if input fields are there.
1175     if thermal_lag_prof_avail && thermal_lag_cond_raw_avail ...
1176         && thermal_lag_temp_raw_avail && thermal_lag_pres_avail ...
1177         && thermal_lag_time_avail && thermal_lag_depth_avail ...
1178         && (thermal_lag_flow_const || thermal_lag_flow_avail)
1179       num_profiles = fix(max(data_proc.profile_index));
1180       % Estimate thermal lag constant, if needed.
1181       if thermal_lag_params_avail
1182         % Thermal lag parameters given.
1183         thermal_lag_constants = thermal_lag_params;
1184       else
1185         fprintf('Performing thermal lag parameter estimation %d with settings:\n', thermal_lag_option_idx);
1186         fprintf('  conductivity sequence: %s\n', thermal_lag_cond_raw);
1187         fprintf('  temperature sequence : %s\n', thermal_lag_temp_raw);
1188         fprintf('  pressure sequence    : %s\n', thermal_lag_pres_raw);
1189         fprintf('  time sequence        : %s\n', thermal_lag_time);
1190         fprintf('  depth sequence       : %s\n', thermal_lag_depth);
1191         if ~thermal_lag_flow_const
1192           fprintf('  flow speed sequence  : %s\n', thermal_lag_flow);
1193         end
1194         fprintf('  estimator            : %s\n', func2str(thermal_lag_estimator));
1195         % Estimate thermal lag time constant for each pofile.
1196         thermal_lag_estimates = nan(num_profiles-1, thermal_lag_num_params);
1197         thermal_lag_residuals = nan(num_profiles-1, 1);
1198         thermal_lag_exitflags = nan(num_profiles-1, 1);
1199         for profile_idx = 1:(num_profiles-1)
1200           prof1_select = (data_proc.profile_index == profile_idx);
1201           [~, ~, prof1_dir] = ...
1202             find(data_proc.profile_direction(prof1_select), 1);
1203           prof1_cond = data_proc.(thermal_lag_cond_raw)(prof1_select);
1204           prof1_temp = data_proc.(thermal_lag_temp_raw)(prof1_select);
1205           prof1_pres = data_proc.(thermal_lag_pres_raw)(prof1_select);
1206           prof1_time = data_proc.(thermal_lag_time)(prof1_select);
1207           prof1_depth = data_proc.(thermal_lag_depth)(prof1_select);
1208           prof2_select = (data_proc.profile_index == profile_idx + 1);
1209           [~, ~, prof2_dir] = ...
1210             find(data_proc.profile_direction(prof2_select), 1);
1211           prof2_cond = data_proc.(thermal_lag_cond_raw)(prof2_select);
1212           prof2_temp = data_proc.(thermal_lag_temp_raw)(prof2_select);
1213           prof2_pres = data_proc.(thermal_lag_pres_raw)(prof2_select);
1214           prof2_time = data_proc.(thermal_lag_time)(prof2_select);
1215           prof2_depth = data_proc.(thermal_lag_depth)(prof2_select);
1216           if thermal_lag_flow_const
1217             prof1_vars = ...
1218               {prof1_time(:) prof1_cond(:) prof1_temp(:) prof1_pres(:)};
1219             prof2_vars = ...
1220               {prof2_time(:) prof2_cond(:) prof2_temp(:) prof2_pres(:)};
1221           else
1222             prof1_flow = data_proc.(thermal_lag_flow)(prof1_select);
1223             prof1_vars = ...
1224               {prof1_time(:) prof1_cond(:) prof1_temp(:) prof1_pres(:) prof1_flow(:)};
1225             prof2_flow = data_proc.(thermal_lag_flow)(prof2_select);
1226             prof2_vars = ...
1227               {prof2_time(:) prof2_cond(:) prof2_temp(:) prof2_pres(:) prof2_flow(:)};
1228           end
1229           [prof1_valid, ~, prof1_vars{:}] = ...
1230             validateProfile(prof1_depth(:), prof1_vars{:}, ...
1231                             'range', options.profile_min_range, ...
1232                             'gap', options.profile_max_gap_ratio);
1233           [prof2_valid, ~, prof2_vars{:}] = ...
1234             validateProfile(prof2_depth(:), prof2_vars{:}, ...
1235                             'range', options.profile_min_range, ...
1236                             'gap', options.profile_max_gap_ratio);
1237           prof_opposite_dir = (prof1_dir * prof2_dir < 0);
1238           if prof1_valid && prof2_valid && prof_opposite_dir
1239             try
1240               [thermal_lag_estimates(profile_idx, :), ...
1241                thermal_lag_exitflags(profile_idx), ...
1242                thermal_lag_residuals(profile_idx)] = ...
1243                 findThermalLagParams(prof1_vars{:}, prof2_vars{:}, thermal_lag_minopts);
1244               if thermal_lag_exitflags(profile_idx) <= 0
1245                  warning('glider_toolbox:processGliderData:ThermalLagMinimizationError', ...
1246                          'Minimization did not converge for casts %d and %d, residual area: %f.', ...
1247                          profile_idx, profile_idx+1, thermal_lag_residuals(profile_idx));
1248               end
1249             catch exception
1250               fprintf('Thermal lag estimation failed for casts: %d and %d.\n', ...
1251                       profile_idx, profile_idx+1);
1252               disp(getReport(exception, 'extended'));
1253             end
1254           end
1255         end
1256         % Compute statistical estimate from individual profile estimates.
1257         % Use feval to allow estimator as either function handle or name string.
1258         thermal_lag_constants = thermal_lag_estimator(thermal_lag_estimates);
1259       end
1260       % Correct thermal lag, if possible.
1261       if any(isnan(thermal_lag_constants))
1262         fprintf('Omiting thermal lag correction %d (%s and %s): %s.\n', ...
1263                 thermal_lag_option_idx, thermal_lag_cond_cor, ...
1264                 thermal_lag_temp_cor, 'no valid parameters available');
1265       else
1266         fprintf('Performing thermal lag correction %d with settings:\n', thermal_lag_option_idx);
1267         fprintf('  output temperature sequence : %s\n', thermal_lag_temp_cor);
1268         fprintf('  output conductivity sequence: %s\n', thermal_lag_cond_cor);
1269         fprintf('  input conductivity sequence : %s\n', thermal_lag_cond_raw);
1270         fprintf('  input temperature sequence  : %s\n', thermal_lag_temp_raw);        
1271         fprintf('  input pressure sequence     : %s\n', thermal_lag_pres_raw);
1272         fprintf('  input time sequence         : %s\n', thermal_lag_time);
1273         fprintf('  input depth sequence        : %s\n', thermal_lag_depth);
1274         if thermal_lag_flow_const
1275           fprintf('  parameters                  : %f %f\n', thermal_lag_constants);
1276         else
1277           fprintf('  input flow sequence         : %s\n', thermal_lag_flow);
1278           fprintf('  parameters                  : %f %f %f %f\n', thermal_lag_constants);
1279         end
1280         data_proc.(thermal_lag_cond_cor) = ...
1281           nan(size(data_proc.(thermal_lag_cond_raw)));
1282         data_proc.(thermal_lag_temp_cor) = ...
1283           nan(size(data_proc.(thermal_lag_temp_raw)));
1284         for profile_idx = 1:num_profiles
1285           prof_select = (data_proc.profile_index == profile_idx);
1286           prof_cond_raw = data_proc.(thermal_lag_cond_raw)(prof_select);
1287           prof_temp_raw = data_proc.(thermal_lag_temp_raw)(prof_select);
1288           prof_time = data_proc.(thermal_lag_time)(prof_select);
1289           prof_depth = data_proc.(thermal_lag_depth)(prof_select);
1290           if thermal_lag_flow_const
1291             prof_vars = {prof_time(:) prof_cond_raw(:) prof_temp_raw(:)};
1292           else
1293             prof_flow = data_proc.(thermal_lag_flow)(prof_select);
1294             prof_vars = ...
1295               {prof_time(:) prof_cond_raw(:) prof_temp_raw(:) prof_flow(:)};
1296           end
1297           [prof_valid, ~, prof_vars{:}] = ...
1298             validateProfile(prof_depth(:), prof_vars{:}, ...
1299                             'range', options.profile_min_range, ...
1300                             'gap', options.profile_max_gap_ratio);
1301           if prof_valid
1302             [prof_temp_cor, prof_cond_cor] = ...
1303               correctThermalLag(prof_vars{:}, thermal_lag_constants);
1304             data_proc.(thermal_lag_temp_cor)(prof_select) = prof_temp_cor;
1305             data_proc.(thermal_lag_cond_cor)(prof_select) = prof_cond_cor;
1306           end
1307         end
1308         if thermal_lag_flow_const
1309           meta_proc.(thermal_lag_cond_cor).sources = ...
1310             {thermal_lag_cond_raw thermal_lag_temp_raw thermal_lag_pres_raw ...
1311              thermal_lag_time thermal_lag_depth 'profile_index'}';
1312         else
1313           meta_proc.(thermal_lag_cond_cor).sources = ...
1314             {thermal_lag_cond_raw thermal_lag_temp_raw thermal_lag_pres_raw ...
1315              thermal_lag_flow thermal_lag_time thermal_lag_depth 'profile_index'}';
1316         end
1317         meta_proc.(thermal_lag_cond_cor).method = 'correctThermalLag';
1318         meta_proc.(thermal_lag_cond_cor).parameters = thermal_lag_constants;
1319         if thermal_lag_params_avail
1320           meta_proc.(thermal_lag_cond_cor).parameter_method = 'preset';
1321         else
1322           meta_proc.(thermal_lag_cond_cor).parameter_method = 'findThermalLagParams';
1323           meta_proc.(thermal_lag_cond_cor).parameter_estimates = thermal_lag_estimates;
1324           meta_proc.(thermal_lag_cond_cor).parameter_exitflags = thermal_lag_exitflags;
1325           meta_proc.(thermal_lag_cond_cor).parameter_estimator = func2str(thermal_lag_estimator);
1326         end
1327         meta_proc.(thermal_lag_cond_cor).profile_min_range = options.profile_min_range;
1328         meta_proc.(thermal_lag_cond_cor).profile_gap_ratio = options.profile_max_gap_ratio;
1329         meta_proc.(thermal_lag_temp_cor) = meta_proc.(thermal_lag_cond_cor);
1330       end
1331     end
1332   end
1333   
1334   
1335   %% Derive salinity from pressure, conductivity and temperature, if available.
1336   for salinity_option_idx = 1:numel(options.salinity_list)
1337     salinity_option = options.salinity_list(salinity_option_idx);
1338     salinity_salt = salinity_option.salinity;
1339     salinity_cond = salinity_option.conductivity;
1340     salinity_temp = salinity_option.temperature;
1341     salinity_pres = salinity_option.pressure;
1342     if all(isfield(data_proc, {salinity_cond salinity_temp salinity_pres}))
1343       % Compute salinity from temperature, pressure and conductivity ratio.
1344       % Input conductivity is given in S/m (Siemens per metre),
1345       % but reference conductivity returned by sw_c3515 is in mS/cm.
1346       fprintf('Deriving salinity %d with settings:\n', salinity_option_idx);
1347       fprintf('  output salinity sequence   : %s\n', salinity_salt);
1348       fprintf('  input conductivity sequence: %s\n', salinity_cond);
1349       fprintf('  input temperature sequence : %s\n', salinity_temp);
1350       fprintf('  input pressure sequence    : %s\n', salinity_pres);
1351       data_proc.(salinity_salt) = ...
1352         sw_salt(data_proc.(salinity_cond) * (10 / sw_c3515()), ...
1353                 data_proc.(salinity_temp), data_proc.(salinity_pres));
1354       meta_proc.(salinity_salt).sources = ...
1355         {salinity_cond salinity_temp salinity_pres}';
1356       meta_proc.(salinity_salt).method = 'sw_salt';
1357     end
1358   end
1359   
1360   
1361   %% Derive density from pressure, salinity and temperature, if available.
1362   for density_option_idx = 1:numel(options.density_list)
1363     density_option = options.density_list(density_option_idx);
1364     density_dens = density_option.density;
1365     density_salt = density_option.salinity;
1366     density_temp = density_option.temperature;
1367     density_pres = density_option.pressure;
1368     if all(isfield(data_proc, {density_salt density_temp density_pres}))
1369       % Compute density from temperature, pressure and salinity.
1370       fprintf('Deriving density %d with settings:\n', density_option_idx);
1371       fprintf('  output density sequence   : %s\n', density_dens);
1372       fprintf('  input salinity sequence   : %s\n', density_salt);
1373       fprintf('  input temperature sequence: %s\n', density_temp);
1374       fprintf('  input pressure sequence   : %s\n', density_pres);
1375       data_proc.(density_dens) = ...
1376         sw_dens(data_proc.(density_salt), ...
1377                 data_proc.(density_temp), data_proc.(density_pres));
1378       meta_proc.(density_dens).sources = ...
1379         {density_salt density_temp density_pres}';
1380       meta_proc.(density_dens).method = 'sw_dens';
1381     end
1382   end
1383   
1384 end