FINDTHERMALLAGPARAMS Thermal lag parameter estimation for profiles.
Syntax:
PARAMS = FINDTHERMALLAGPARAMS(TIME1, COND1, TEMP1, PRES1, TIME2, COND2, TEMP2, PRES2)
PARAMS = FINDTHERMALLAGPARAMS(TIME1, COND1, TEMP1, PRES1, FLOW1, TIME2, COND2, TEMP2, PRES2, FLOW2)
PARAMS = FINDTHERMALLAGPARAMS(..., OPTIONS)
PARAMS = FINDTHERMALLAGPARAMS(..., OPT1, VAL1, ...)
[PARAMS, EXITFLAG, RESIDUAL] = FINDTHERMALLAGPARAMS(...)
PARAMS = FINDTHERMALLAGPARAMS(TIME1, COND1, TEMP1, PRES1, TIME2, COND2, TEMP2, PRES2)
finds the thermal lag parameters for two CTD profiles with constant flow
speed (pumped CTD) given by sequences of time (s), conductivity (S/m),
temperature (deg C) and pressure (dbar) in respective vectors TIME1, COND1,
TEMP1 and PRES1, and TIME2, COND2, TEMP2 and PRES2. The profiles are
supposed to measure the same column of water in opposite directions.
The computed parameters are returned in a two element vector PARAMS,
with the error magnitude (alpha), and the error time constant (tau).
A detailed description of these parameters may be found in the references
listed below (Lueck 1990).
Under the assumption that both profiles should be as similar as possible,
the function solves the minimization problem of finding the thermal lag
parameters such that the area between profiles of temperature and salinity
is minimal; where salinity is derived from temperature, conductivity and
pressure sequences using SW_SALT with the corrected temperature sequence
returned by CORRECTTHERMALLAG. This problem is solved by function FMINCON,
using default values for the initial guess and the parameter bounds.
See OPTIONS description below.
PARAMS = FINDTHERMALLAGPARAMS(TIME1, COND1, TEMP1, PRES1, FLOW1, TIME2, DEPTH2, TEMP2, PRES2, FLOW2, ...)
performs the same estimation but for a pair of CTD profiles with variable
flow speed (unpumped CTD), given by respective vectors FLOW1 and FLOW2.
The estimated parameters are returned in a four element vector PARAMS,
with the offset and the slope of the error magnitude (alpha_o and alpha_s)
and the offset and the slope of the error time (tau_o and tau_s). Details
on these parameters may be found in references below (Morison 1994).
PARAMS = FINDTHERMALLAGPARAMS(..., OPTIONS) and
PARAMS = FINDTHERMALLAGPARAMS(..., OPT1, VAL1, ...) allow passing extra
options 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:
GRAPHICS: whether graphic output should be produced.
A boolean. If true a nice figure showing the minimization process will
be displayed. It includes the parameter values, the objective function
value, a temperature-salinity diagram, a pressure-salinity diagram,
a temperature-time plot, and a conductivity-time plot.
Default value: false.
GUESS: initial guess for minimization function FMINCON.
A two or four element vector with the initial guess for each parameter.
Default value:
For constant flow speed: [0.0677 11.1431] (see note below)
For variable flow speed: [0.0135 0.0264 7.1499 2.7858] (Morison 1994)
LOWER: lower bounds of parameters for minimization function FMINCON.
A two or four element vector with the lower bound for each parameter.
Default value:
For constant flow speed: [0 0] (no correction)
For variable flow speed: [0 0 0 0] (no correction)
UPPER: upper bounds of parameters for minimization function FMINCON.
A two or four element vector with the upper bound for each parameter.
Default value:
For constant flow speed: [4 2.5*RANGE(TIME1)]
For variable flow speed: [2 1 RANGE(TIME1) RANGE(TIME1)/2]
OPTIMOPTS: extra options for the minimization function FMINCON.
An option struct as needed by the function FMINCON.
Default value: default options for FMINCON, except for:
'Algorithm': 'interior-point'
'FinDiffType': 'central'
'TolFun': 1e-4
'TolCon': 1e-5
'TolX': 1e-5
'Display': 'off'
[PARAMS, EXITFLAG, RESIDUAL] = FINDTHERMALLAGPARAMS(...) also returns the
exit code of the minimization function FMINCON in EXITFLAG, and the
residual area in RESIDUAL. EXITFLAG is positive when minimization succeeds.
Notes:
This function is an improved version of a previous function by Tomeu Garau,
called ADJUSTTHERMALLAGPARAMS. He is the true glider man. Main changes are:
- Support for CTD profiles with constant flow speed (pumped CTD).
- Different minimization algorithm (interior-point instead of active-set).
- Support for custom minimization options.
- Optional predefined graphical output.
It remains to be assessed that the parameters found minimize the area
between the corrected profiles globally, because the solver might look for
local minimizers only.
Parameters' initial guess values and bounds are taken from original
implementation. Since it did not handle constant flow speed profiles,
the values for this case has been set according to the same reference,
aproximating the constant CTD flow speed by v = 0.4867 m/s:
alpha = 0.0135 + 0.0264 / 0.4867 = 0.0677
tau = 7.1499 + 2.7858 / sqrt(11.1431) = 11.1431
The flow speed value is estimated from reported Seabird Glider Payload CTD
(GPCTD) specifications as:
V = Q / U * L = 0.010 / 0.003 * 0.146 = 0.4867 m/s
where:
Q = 0.010 l/s is the flow rate.
U = 0.003 l is the cell volume.
L = 0.146 m is the cell length.
References:
Garau, B.; Ruiz, S.; G. Zhang, W.; Pascual, A.; Heslop, E.;
Kerfoot, J.; and Tintoré, J.; 2011:
Thermal Lag Correction on Slocum CTD Glider Data.
Journal of Atmospheric and Oceanic Technology, vol. 28, pages 1065-1071.
Morison, J.; Andersen, R.; Larson, N.; D'Asaro, E.; and Boyd, T.; 1994:
The Correction for Thermal-Lag Effects in Sea-Bird CTD Data.
Journal of Atmospheric and Oceanic Technology, vol. 11, pages 1151-1164.
Lueck, R. G.; and Picklo, J. J.; 1990:
Thermal Inertia of Conductivity Cells: Observations with a Sea-Bird cell.
Journal of Atmospheric and Oceanic Technology, vol. 7, pages 756–768.
Lueck, R. G.; 1990:
Thermal Inertia of Conductivity Cells: Theory.
Journal of Atmospheric and Oceanic Technology, vol. 7, pages 741–755.
Examples:
% Constant flow speed profiles (pumped CTD):
data = load('private/test/ctd_pumped.dat');
profile = data(:, 10);
time1 = data(profile==1, 1);
cond1 = data(profile==1, 2);
temp1 = data(profile==1, 3);
pres1 = data(profile==1, 4);
time2 = data(profile==2, 1);
cond2 = data(profile==2, 2);
temp2 = data(profile==2, 3);
pres2 = data(profile==2, 4);
params = findThermalLagParams(time1, cond1, temp1, pres1, ...
time2, cond2, temp2, pres2)
% Variable flow speed profiles (unpumped CTD):
data = load('private/test/ctd_unpumped.dat');
profile = data(:, 10);
time1 = data(profile==1, 1);
cond1 = data(profile==1, 2);
temp1 = data(profile==1, 3);
pres1 = data(profile==1, 4);
ptch1 = data(profile==1, 5);
lati1 = data(profile==1, 8);
dpth1 = sw_dpth(pres1, lati1);
flow1 = computeCTDFlowSpeed(time1, dpth1, ptch1, 'minpitch', deg2rad(11));
time2 = data(profile==2, 1);
cond2 = data(profile==2, 2);
temp2 = data(profile==2, 3);
pres2 = data(profile==2, 4);
ptch2 = data(profile==2, 5);
lati2 = data(profile==2, 8);
dpth2 = sw_dpth(pres2, lati2);
flow2 = computeCTDFlowSpeed(time2, dpth2, ptch2, 'minpitch', deg2rad(11));
params = findThermalLagParams(time1, cond1, temp1, pres1, flow1, ...
time2, cond2, temp2, pres2, flow2)
% Variable flow speed profiles with exit code, residual and extra options.
[params, exitflag, residual] = ...
findThermalLagParams(time1, cond1, temp1, pres1, flow1, ...
time2, cond2, temp2, pres2, flow2, ...
'graphics', 'true')
See also:
FMINCON
OPTIMSET
CORRECTTHERMALLAG
SW_SALT
SW_DENS
COMPUTECTDFLOWSPEED
Authors:
Joan Pau Beltran <joanpau.beltran@socib.cat>0001 function [params, exitflag, residual] = findThermalLagParams(varargin) 0002 %FINDTHERMALLAGPARAMS Thermal lag parameter estimation for profiles. 0003 % 0004 % Syntax: 0005 % PARAMS = FINDTHERMALLAGPARAMS(TIME1, COND1, TEMP1, PRES1, TIME2, COND2, TEMP2, PRES2) 0006 % PARAMS = FINDTHERMALLAGPARAMS(TIME1, COND1, TEMP1, PRES1, FLOW1, TIME2, COND2, TEMP2, PRES2, FLOW2) 0007 % PARAMS = FINDTHERMALLAGPARAMS(..., OPTIONS) 0008 % PARAMS = FINDTHERMALLAGPARAMS(..., OPT1, VAL1, ...) 0009 % [PARAMS, EXITFLAG, RESIDUAL] = FINDTHERMALLAGPARAMS(...) 0010 % 0011 % PARAMS = FINDTHERMALLAGPARAMS(TIME1, COND1, TEMP1, PRES1, TIME2, COND2, TEMP2, PRES2) 0012 % finds the thermal lag parameters for two CTD profiles with constant flow 0013 % speed (pumped CTD) given by sequences of time (s), conductivity (S/m), 0014 % temperature (deg C) and pressure (dbar) in respective vectors TIME1, COND1, 0015 % TEMP1 and PRES1, and TIME2, COND2, TEMP2 and PRES2. The profiles are 0016 % supposed to measure the same column of water in opposite directions. 0017 % The computed parameters are returned in a two element vector PARAMS, 0018 % with the error magnitude (alpha), and the error time constant (tau). 0019 % A detailed description of these parameters may be found in the references 0020 % listed below (Lueck 1990). 0021 % 0022 % Under the assumption that both profiles should be as similar as possible, 0023 % the function solves the minimization problem of finding the thermal lag 0024 % parameters such that the area between profiles of temperature and salinity 0025 % is minimal; where salinity is derived from temperature, conductivity and 0026 % pressure sequences using SW_SALT with the corrected temperature sequence 0027 % returned by CORRECTTHERMALLAG. This problem is solved by function FMINCON, 0028 % using default values for the initial guess and the parameter bounds. 0029 % See OPTIONS description below. 0030 % 0031 % PARAMS = FINDTHERMALLAGPARAMS(TIME1, COND1, TEMP1, PRES1, FLOW1, TIME2, DEPTH2, TEMP2, PRES2, FLOW2, ...) 0032 % performs the same estimation but for a pair of CTD profiles with variable 0033 % flow speed (unpumped CTD), given by respective vectors FLOW1 and FLOW2. 0034 % The estimated parameters are returned in a four element vector PARAMS, 0035 % with the offset and the slope of the error magnitude (alpha_o and alpha_s) 0036 % and the offset and the slope of the error time (tau_o and tau_s). Details 0037 % on these parameters may be found in references below (Morison 1994). 0038 % 0039 % PARAMS = FINDTHERMALLAGPARAMS(..., OPTIONS) and 0040 % PARAMS = FINDTHERMALLAGPARAMS(..., OPT1, VAL1, ...) allow passing extra 0041 % options given in key-value pairs OPT1, VAL1... or in a struct OPTIONS with 0042 % field names as option keys and field values as option values. 0043 % Recognized options are: 0044 % GRAPHICS: whether graphic output should be produced. 0045 % A boolean. If true a nice figure showing the minimization process will 0046 % be displayed. It includes the parameter values, the objective function 0047 % value, a temperature-salinity diagram, a pressure-salinity diagram, 0048 % a temperature-time plot, and a conductivity-time plot. 0049 % Default value: false. 0050 % GUESS: initial guess for minimization function FMINCON. 0051 % A two or four element vector with the initial guess for each parameter. 0052 % Default value: 0053 % For constant flow speed: [0.0677 11.1431] (see note below) 0054 % For variable flow speed: [0.0135 0.0264 7.1499 2.7858] (Morison 1994) 0055 % LOWER: lower bounds of parameters for minimization function FMINCON. 0056 % A two or four element vector with the lower bound for each parameter. 0057 % Default value: 0058 % For constant flow speed: [0 0] (no correction) 0059 % For variable flow speed: [0 0 0 0] (no correction) 0060 % UPPER: upper bounds of parameters for minimization function FMINCON. 0061 % A two or four element vector with the upper bound for each parameter. 0062 % Default value: 0063 % For constant flow speed: [4 2.5*RANGE(TIME1)] 0064 % For variable flow speed: [2 1 RANGE(TIME1) RANGE(TIME1)/2] 0065 % OPTIMOPTS: extra options for the minimization function FMINCON. 0066 % An option struct as needed by the function FMINCON. 0067 % Default value: default options for FMINCON, except for: 0068 % 'Algorithm': 'interior-point' 0069 % 'FinDiffType': 'central' 0070 % 'TolFun': 1e-4 0071 % 'TolCon': 1e-5 0072 % 'TolX': 1e-5 0073 % 'Display': 'off' 0074 % 0075 % [PARAMS, EXITFLAG, RESIDUAL] = FINDTHERMALLAGPARAMS(...) also returns the 0076 % exit code of the minimization function FMINCON in EXITFLAG, and the 0077 % residual area in RESIDUAL. EXITFLAG is positive when minimization succeeds. 0078 % 0079 % Notes: 0080 % This function is an improved version of a previous function by Tomeu Garau, 0081 % called ADJUSTTHERMALLAGPARAMS. He is the true glider man. Main changes are: 0082 % - Support for CTD profiles with constant flow speed (pumped CTD). 0083 % - Different minimization algorithm (interior-point instead of active-set). 0084 % - Support for custom minimization options. 0085 % - Optional predefined graphical output. 0086 % 0087 % It remains to be assessed that the parameters found minimize the area 0088 % between the corrected profiles globally, because the solver might look for 0089 % local minimizers only. 0090 % 0091 % Parameters' initial guess values and bounds are taken from original 0092 % implementation. Since it did not handle constant flow speed profiles, 0093 % the values for this case has been set according to the same reference, 0094 % aproximating the constant CTD flow speed by v = 0.4867 m/s: 0095 % alpha = 0.0135 + 0.0264 / 0.4867 = 0.0677 0096 % tau = 7.1499 + 2.7858 / sqrt(11.1431) = 11.1431 0097 % The flow speed value is estimated from reported Seabird Glider Payload CTD 0098 % (GPCTD) specifications as: 0099 % V = Q / U * L = 0.010 / 0.003 * 0.146 = 0.4867 m/s 0100 % where: 0101 % Q = 0.010 l/s is the flow rate. 0102 % U = 0.003 l is the cell volume. 0103 % L = 0.146 m is the cell length. 0104 % 0105 % References: 0106 % Garau, B.; Ruiz, S.; G. Zhang, W.; Pascual, A.; Heslop, E.; 0107 % Kerfoot, J.; and Tintoré, J.; 2011: 0108 % Thermal Lag Correction on Slocum CTD Glider Data. 0109 % Journal of Atmospheric and Oceanic Technology, vol. 28, pages 1065-1071. 0110 % 0111 % Morison, J.; Andersen, R.; Larson, N.; D'Asaro, E.; and Boyd, T.; 1994: 0112 % The Correction for Thermal-Lag Effects in Sea-Bird CTD Data. 0113 % Journal of Atmospheric and Oceanic Technology, vol. 11, pages 1151-1164. 0114 % 0115 % Lueck, R. G.; and Picklo, J. J.; 1990: 0116 % Thermal Inertia of Conductivity Cells: Observations with a Sea-Bird cell. 0117 % Journal of Atmospheric and Oceanic Technology, vol. 7, pages 756–768. 0118 % 0119 % Lueck, R. G.; 1990: 0120 % Thermal Inertia of Conductivity Cells: Theory. 0121 % Journal of Atmospheric and Oceanic Technology, vol. 7, pages 741–755. 0122 % 0123 % Examples: 0124 % % Constant flow speed profiles (pumped CTD): 0125 % data = load('private/test/ctd_pumped.dat'); 0126 % profile = data(:, 10); 0127 % time1 = data(profile==1, 1); 0128 % cond1 = data(profile==1, 2); 0129 % temp1 = data(profile==1, 3); 0130 % pres1 = data(profile==1, 4); 0131 % time2 = data(profile==2, 1); 0132 % cond2 = data(profile==2, 2); 0133 % temp2 = data(profile==2, 3); 0134 % pres2 = data(profile==2, 4); 0135 % params = findThermalLagParams(time1, cond1, temp1, pres1, ... 0136 % time2, cond2, temp2, pres2) 0137 % 0138 % % Variable flow speed profiles (unpumped CTD): 0139 % data = load('private/test/ctd_unpumped.dat'); 0140 % profile = data(:, 10); 0141 % time1 = data(profile==1, 1); 0142 % cond1 = data(profile==1, 2); 0143 % temp1 = data(profile==1, 3); 0144 % pres1 = data(profile==1, 4); 0145 % ptch1 = data(profile==1, 5); 0146 % lati1 = data(profile==1, 8); 0147 % dpth1 = sw_dpth(pres1, lati1); 0148 % flow1 = computeCTDFlowSpeed(time1, dpth1, ptch1, 'minpitch', deg2rad(11)); 0149 % time2 = data(profile==2, 1); 0150 % cond2 = data(profile==2, 2); 0151 % temp2 = data(profile==2, 3); 0152 % pres2 = data(profile==2, 4); 0153 % ptch2 = data(profile==2, 5); 0154 % lati2 = data(profile==2, 8); 0155 % dpth2 = sw_dpth(pres2, lati2); 0156 % flow2 = computeCTDFlowSpeed(time2, dpth2, ptch2, 'minpitch', deg2rad(11)); 0157 % params = findThermalLagParams(time1, cond1, temp1, pres1, flow1, ... 0158 % time2, cond2, temp2, pres2, flow2) 0159 % 0160 % % Variable flow speed profiles with exit code, residual and extra options. 0161 % [params, exitflag, residual] = ... 0162 % findThermalLagParams(time1, cond1, temp1, pres1, flow1, ... 0163 % time2, cond2, temp2, pres2, flow2, ... 0164 % 'graphics', 'true') 0165 % 0166 % See also: 0167 % FMINCON 0168 % OPTIMSET 0169 % CORRECTTHERMALLAG 0170 % SW_SALT 0171 % SW_DENS 0172 % COMPUTECTDFLOWSPEED 0173 % 0174 % Authors: 0175 % Joan Pau Beltran <joanpau.beltran@socib.cat> 0176 0177 % Copyright (C) 2013-2016 0178 % ICTS SOCIB - Servei d'observacio i prediccio costaner de les Illes Balears 0179 % <http://www.socib.es> 0180 % 0181 % This program is free software: you can redistribute it and/or modify 0182 % it under the terms of the GNU General Public License as published by 0183 % the Free Software Foundation, either version 3 of the License, or 0184 % (at your option) any later version. 0185 % 0186 % This program is distributed in the hope that it will be useful, 0187 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0188 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0189 % GNU General Public License for more details. 0190 % 0191 % You should have received a copy of the GNU General Public License 0192 % along with this program. If not, see <http://www.gnu.org/licenses/>. 0193 0194 error(nargchk(8, 20, nargin, 'struct')); 0195 0196 0197 %% Parse basic input arguments. 0198 % Get numeric (non-option) arguments. 0199 nargnum = find(~cellfun(@isnumeric, varargin), 1, 'first') - 1; 0200 if isempty(nargnum) 0201 nargnum = nargin; 0202 end 0203 switch(nargnum) 0204 case 8 0205 % Constant flow speed (pumped CTD). 0206 [time1, cond1, temp1, pres1] = varargin{1:4}; 0207 [time2, cond2, temp2, pres2] = varargin{5:8}; 0208 constant_flow = true; 0209 case 10 0210 % Variable flow speed (unpumped CTD). 0211 [time1, cond1, temp1, pres1, flow1] = varargin{1:5}; 0212 [time2, cond2, temp2, pres2, flow2] = varargin{6:10}; 0213 constant_flow = false; 0214 end 0215 0216 0217 %% Configure default options. 0218 % For the case of variable flow speed (unpumped CTD) 0219 % this is equivalent to the original version by Tomeu Garau, 0220 % except for the method (it was active-set). 0221 time_range = min(max(time1) - min(time1), max(time2) - min(time2)); 0222 options.graphics = false; 0223 if constant_flow 0224 options.guess = [0.0677 11.1431]; % above parameters applied to GPCTD flow speed. 0225 options.lower = [0 0]; % no correction. 0226 options.upper = [4 2.5*time_range]; % above parameters applied to GPCTD flow speed. 0227 else 0228 options.guess = [0.0135 0.0264 7.1499 2.7858]; % from Morrison (1994). 0229 options.lower = [0 0 0 0]; % no correction. 0230 options.upper = [2 1 time_range time_range/2]; % from old version. 0231 end 0232 options.optimopts = optimset(optimset('fmincon'), ... 0233 'Algorithm', 'interior-point', ... 0234 'TolFun', 1e-4, 'TolCon', 1e-5, 'TolX', 1e-5, ... 0235 'FinDiffType', 'central', ... 0236 'Display', 'off'); 0237 0238 0239 %% Parse option arguments. 0240 % Get option key-value pairs in any accepted call signature. 0241 argopts = varargin(nargnum+1:end); 0242 if isscalar(argopts) && isstruct(argopts{1}) 0243 % Options passed as a single option struct argument: 0244 % field names are option keys and field values are option values. 0245 opt_key_list = fieldnames(argopts{1}); 0246 opt_val_list = struct2cell(argopts{1}); 0247 elseif mod(numel(argopts), 2) == 0 0248 % Options passed as key-value argument pairs. 0249 opt_key_list = argopts(1:2:end); 0250 opt_val_list = argopts(2:2:end); 0251 else 0252 error('glider_toolbox:findThermalLagParams:InvalidOptions', ... 0253 'Invalid optional arguments (neither key-value pairs nor struct).'); 0254 end 0255 % Overwrite default options with values given in extra arguments. 0256 for opt_idx = 1:numel(opt_key_list) 0257 opt = lower(opt_key_list{opt_idx}); 0258 val = opt_val_list{opt_idx}; 0259 if isfield(options, opt) 0260 options.(opt) = val; 0261 else 0262 error('glider_toolbox:findThermalLagParams:InvalidOption', ... 0263 'Invalid option: %s.', opt); 0264 end 0265 end 0266 0267 0268 %% Enable graphical output, if needed. 0269 % The functions are defined below. 0270 if options.graphics 0271 if constant_flow 0272 [temp_cor1, cond_cor1] = ... 0273 correctThermalLag(time1, cond1, temp1, options.guess); 0274 [temp_cor2, cond_cor2] = ... 0275 correctThermalLag(time2, cond2, temp2, options.guess); 0276 else 0277 [temp_cor1, cond_cor1] = ... 0278 correctThermalLag(time1, cond1, temp1, flow1, options.guess); 0279 [temp_cor2, cond_cor2] = ... 0280 correctThermalLag(time2, cond2, temp2, flow2, options.guess); 0281 end 0282 salt_cor1 = sw_salt(cond1 * (10 / sw_c3515()), temp_cor1, pres1); 0283 salt_cor2 = sw_salt(cond2 * (10 / sw_c3515()), temp_cor2, pres2); 0284 salt1 = sw_salt(cond1 * (10 / sw_c3515()), temp1, pres1); 0285 salt2 = sw_salt(cond2 * (10 / sw_c3515()), temp2, pres2); 0286 options.optimopts.OutputFcn = @optimoutUpdateCorrectedData; 0287 options.optimopts.PlotFcns = ... 0288 {@optimplotx @optimplotfval ... 0289 @optimplotTempSalt @optimplotPresSalt ... 0290 @optimplotTempTime @optimplotCondTime}; 0291 end 0292 0293 0294 %% Perform estimation through minimization. 0295 % Definition of minimization objective function. 0296 objective_function = @optimobjTSArea; 0297 0298 % Run minimization procedure. 0299 [params, residual, exitflag] = ... 0300 fmincon(objective_function, options.guess, ... 0301 [], [], [], [], options.lower, options.upper, [], ... 0302 options.optimopts); 0303 0304 0305 %% Definition of auxiliary objective and plotting functions. 0306 % They should be nested to access cast data. 0307 function area = optimobjTSArea(params) 0308 %OPTIMOBJTSAREA Compute area enclosed by profiles in TS diagram. 0309 if constant_flow 0310 temp_cor1 = correctThermalLag(time1, cond1, temp1, params); 0311 temp_cor2 = correctThermalLag(time2, cond2, temp2, params); 0312 else 0313 temp_cor1 = correctThermalLag(time1, cond1, temp1, flow1, params); 0314 temp_cor2 = correctThermalLag(time2, cond2, temp2, flow2, params); 0315 end 0316 salt_cor1 = sw_salt(cond1 * (10 / sw_c3515()), temp_cor1, pres1); 0317 salt_cor2 = sw_salt(cond2 * (10 / sw_c3515()), temp_cor2, pres2); 0318 dens_cor1 = sw_dens(salt_cor1, temp1, pres1); 0319 dens_cor2 = sw_dens(salt_cor2, temp2, pres2); 0320 dens_min = max(min(dens_cor1), min(dens_cor2)); 0321 dens_max = min(max(dens_cor1), max(dens_cor2)); 0322 dens_mask1 = (dens_min <= dens_cor1) & (dens_cor1 <= dens_max); 0323 dens_mask2 = (dens_min <= dens_cor2) & (dens_cor2 <= dens_max); 0324 min_idx1 = find(dens_mask1, 1, 'first'); 0325 min_idx2 = find(dens_mask2, 1, 'first'); 0326 max_idx1 = find(dens_mask1, 1, 'last'); 0327 max_idx2 = find(dens_mask2, 1, 'last'); 0328 area = ... 0329 profileArea(salt_cor1(min_idx1:max_idx1), temp1(min_idx1:max_idx1), ... 0330 salt_cor2(min_idx2:max_idx2), temp2(min_idx2:max_idx2)); 0331 end 0332 0333 function stop = optimoutUpdateCorrectedData(params, ~, state) 0334 %OPTIMOUTUPDATEPLOTDATA Update corrected data sequences. 0335 switch state 0336 case 'init' 0337 case 'iter' 0338 if constant_flow 0339 [temp_cor1, cond_cor1] = ... 0340 correctThermalLag(time1, cond1, temp1, params); 0341 [temp_cor2, cond_cor2] = ... 0342 correctThermalLag(time2, cond2, temp2, params); 0343 else 0344 [temp_cor1, cond_cor1] = ... 0345 correctThermalLag(time1, cond1, temp1, flow1, params); 0346 [temp_cor2, cond_cor2] = ... 0347 correctThermalLag(time2, cond2, temp2, flow2, params); 0348 end 0349 salt_cor1 = sw_salt(cond1 * (10 / sw_c3515()), temp_cor1, pres1); 0350 salt_cor2 = sw_salt(cond2 * (10 / sw_c3515()), temp_cor2, pres2); 0351 case 'interrupt' 0352 case 'done' 0353 end 0354 stop = false; 0355 end 0356 0357 % Temperature-salinity diagram. 0358 function stop = optimplotTempSalt(~, ~, state) 0359 %OPTIMPLOTTEMPSALT Temperature-salinity diagram for pair of casts. 0360 switch state 0361 case 'init' 0362 case 'iter' 0363 valid_cor1 = ~(isnan(salt_cor1) | isnan(temp1)); 0364 valid_cor2 = ~(isnan(salt_cor2) | isnan(temp2)); 0365 valid1 = ~(isnan(salt1) | isnan(temp1)); 0366 valid2 = ~(isnan(salt2) | isnan(temp2)); 0367 plot(salt1(valid1), temp1(valid1), ':r', ... 0368 salt2(valid2), temp2(valid2), ':b', ... 0369 salt_cor1(valid_cor1), temp1(valid_cor1), '-r', ... 0370 salt_cor2(valid_cor2), temp2(valid_cor2), '-b'); 0371 title('Temperature-Salinity diagram'); 0372 xlabel('Salinity'); 0373 ylabel('Temperature'); 0374 case 'interrupt' 0375 case 'done' 0376 end 0377 stop = false; 0378 end 0379 0380 % Pressure-salinity diagram. 0381 function stop = optimplotPresSalt(~, ~, state) 0382 %OPTIMPLOTPRESSALT Pressure-salinity diagram for pair of casts. 0383 switch state 0384 case 'init' 0385 case 'iter' 0386 valid_cor1 = ~(isnan(salt_cor1) | isnan(pres1)); 0387 valid_cor2 = ~(isnan(salt_cor2) | isnan(pres2)); 0388 valid1 = ~(isnan(salt1) | isnan(pres1)); 0389 valid2 = ~(isnan(salt2) | isnan(pres2)); 0390 plot(salt1(valid1), pres1(valid1), ':r', ... 0391 salt2(valid2), pres2(valid2), ':b', ... 0392 salt_cor1(valid_cor1), pres1(valid_cor1), '-r', ... 0393 salt_cor2(valid_cor2), pres2(valid_cor2), '-b'); 0394 title('Pressure-Salinity diagram'); 0395 xlabel('Salinity'); 0396 ylabel('Pressure'); 0397 set(gca, 'YDir', 'reverse') 0398 case 'interrupt' 0399 case 'done' 0400 end 0401 stop = false; 0402 end 0403 0404 % Temperature-time diagram. 0405 function stop = optimplotTempTime(~, ~, state) 0406 %OPTIMPLOTTEMPTIME Temperature-time diagram for pair of casts. 0407 switch state 0408 case 'init' 0409 case 'iter' 0410 valid_cor1 = (time1 > 0) & ~isnan(temp_cor1); 0411 valid_cor2 = (time2 > 0) & ~isnan(temp_cor2); 0412 valid1 = (time1 > 0) & ~isnan(temp1); 0413 valid2 = (time2 > 0) & ~isnan(temp2); 0414 time_offset = min(min(time1(valid1)), min(time2(valid2))); 0415 plot(time1(valid1) - time_offset, temp1(valid1), ':r', ... 0416 time2(valid2) - time_offset, temp2(valid2), ':b', ... 0417 time1(valid_cor1) - time_offset, temp_cor1(valid_cor1), '-r', ... 0418 time2(valid_cor2) - time_offset, temp_cor2(valid_cor2), '-b'); 0419 title('Temperature-Time diagram'); 0420 xlabel('Time'); 0421 ylabel('Temperature'); 0422 case 'interrupt' 0423 case 'done' 0424 end 0425 stop = false; 0426 end 0427 0428 % Conductivity-time diagram. 0429 function stop = optimplotCondTime(~, ~, state) 0430 %OPTIMPLOTCONDTIME Conductivity-time diagram for pair of casts. 0431 switch state 0432 case 'init' 0433 case 'iter' 0434 valid_cor1 = (time1 > 0) & ~isnan(cond_cor1); 0435 valid_cor2 = (time2 > 0) & ~isnan(cond_cor2); 0436 valid1 = (time1 > 0) & ~isnan(cond1); 0437 valid2 = (time2 > 0) & ~isnan(cond2); 0438 time_offset = min(min(time1(valid1)), min(time2(valid2))); 0439 plot(time1(valid1) - time_offset, cond1(valid1), ':r', ... 0440 time2(valid2) - time_offset, cond2(valid2), ':b', ... 0441 time1(valid_cor1) - time_offset, cond_cor1(valid_cor1), '-r', ... 0442 time2(valid_cor2) - time_offset, cond_cor2(valid_cor2), '-b'); 0443 title('Conductivity-Time diagram'); 0444 xlabel('Time'); 0445 ylabel('Conductivity'); 0446 case 'interrupt' 0447 case 'done' 0448 end 0449 stop = false; 0450 end 0451 0452 end