COMPUTECTDFLOWSPEED Compute flow speed along CTD cell.
Syntax:
FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH)
FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, PITCH)
FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, OPTIONS)
FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, OPT1, VAL1, ...)
FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, PITCH, OPTIONS)
FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, PITCH, OPT1, VAL1, ...)
Description:
FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH) computes the flow speed
through the cell of a CTD in a vertical profile sequence given by vectors
TIMESTAMP (timestamp) and DEPTH (depth), and returns it in vector FLOW.
All these vectors should have the same dimensions.
FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, PITCH) performs the same
correction but for a non-vertical profile, with pitch angle given by PITCH
in radians. PITCH may be either a vector with the same dimensions as
TIMESTAMP and DEPTH, or a scalar taken to be the constant pitch across the
whole profile sequence.
FLOW = COMPUTECTDFLOWSPEED(..., OPTIONS) and
FLOW = COMPUTECTDFLOWSPEED(..., 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:
FACTORPOLY: flow speed factor polynomial.
Vector with the coefficients of the polynomial that returns the flow
speed factor whem evaluated by POLYVAL at each surge speed value.
See note on flow speed computation below.
Default value: [0.00 0.03 1.15] (see note below).
MINVEL: minimum vertical velocity threshold.
Scalar with the minimum vertical velocity threshold below which the
aproximated surge speed value is supposed to be unreliable. On samples
with an absolute vertical velocity value less than given threshold
flow speed is invalid (NaN).
Default value: 0 (all samples are valid)
MINPITCH: minimum pitch threshold.
Scalar with the minimum pitch threshold in radians below which the
aproximated surge speed value is supposed to be unreliable. On samples
with an absolute pitch value lesser than given threshold flow speed is
invalid (NaN).
Default value: 0 (all samples are valid)
Notes:
This function is based on the flow speed computation code by Tomeu Garau in
function CORRECTTHERMALLAG. Main changes are:
- Changed vertical velocity approximation method (central instead of
forward differences).
- Added support for discarding regions with low velocity or low pitch.
The approximate flow speed is computed from the absolute value of the
vertical velocity. The vertical velocity is approximated using central
differences, but accounting for the possibly unevenly spaced samples.
In non-vertical profiles with given pitch, the vertical velocity is scaled
by the sine of the pitch angle. The absolute value of this intermediate
quantity is called surge speed. To obtain the flow speed, each surge speed
value is scaled by a factor resulting from the evaluation of the given
polynomial at the surge speed value. If no polynomial given, the surge
speed is not scaled and equals the flow speed.
The note on Garau's original code was:
% The relative coefficient between the flow speed inside and outside of
% the conductivity cell. This is still uncertain (ask Gordon for origin).
% Here are three choices for first three orders polynomial.
speed_factor_polynoms = [0.00, 0.00, 0.40; % 0th order degree.
0.00, 0.03, 0.45; % 1st order degree.
1.58, 1.15, 0.70]; % 2nd order degree.
% First order approximation, second row of the matrix.
Surge speed is supposed to be unreliable when the vertical velocity or the
pitch is too low, so the flow speed is left undefined. In October 2013
Gerd Krahmann noted that there might be a bug in the original
implementation, and the matrix of polynomial coefficients (above) should be
transposed:
speed_factor_polynoms = [0.00, 0.00, 1.58; % 0th order degree.
0.00, 0.03, 1.15; % 1st order degree.
0.40, 0.45, 0.70]; % 2nd order degree.
Examples:
% Vertical profile:
flow = computeCTDFlowSpeed(timestamp, depth, pitch)
% Tilted profile:
flow = computeCTDFlowSpeed(timestamp, depth, 25)
flow = computeCTDFlowSpeed(timestamp, depth, pitch)
% No flow speed factor scaling:
flow = computeCTDFlowSpeed(timestamp, depth, pitch, 'factorpoly', [])
Authors:
Joan Pau Beltran <joanpau.beltran@socib.cat>0001 function flow = computeCTDFlowSpeed(varargin) 0002 %COMPUTECTDFLOWSPEED Compute flow speed along CTD cell. 0003 % 0004 % Syntax: 0005 % FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH) 0006 % FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, PITCH) 0007 % FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, OPTIONS) 0008 % FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, OPT1, VAL1, ...) 0009 % FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, PITCH, OPTIONS) 0010 % FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, PITCH, OPT1, VAL1, ...) 0011 % 0012 % Description: 0013 % FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH) computes the flow speed 0014 % through the cell of a CTD in a vertical profile sequence given by vectors 0015 % TIMESTAMP (timestamp) and DEPTH (depth), and returns it in vector FLOW. 0016 % All these vectors should have the same dimensions. 0017 % 0018 % FLOW = COMPUTECTDFLOWSPEED(TIMESTAMP, DEPTH, PITCH) performs the same 0019 % correction but for a non-vertical profile, with pitch angle given by PITCH 0020 % in radians. PITCH may be either a vector with the same dimensions as 0021 % TIMESTAMP and DEPTH, or a scalar taken to be the constant pitch across the 0022 % whole profile sequence. 0023 % 0024 % FLOW = COMPUTECTDFLOWSPEED(..., OPTIONS) and 0025 % FLOW = COMPUTECTDFLOWSPEED(..., OPT1, VAL1, ...) allow passing extra 0026 % options given in key-value pairs OPT1, VAL1... or in a struct OPTIONS 0027 % with field names as option keys and field values as option values. 0028 % Recognized options are: 0029 % FACTORPOLY: flow speed factor polynomial. 0030 % Vector with the coefficients of the polynomial that returns the flow 0031 % speed factor whem evaluated by POLYVAL at each surge speed value. 0032 % See note on flow speed computation below. 0033 % Default value: [0.00 0.03 1.15] (see note below). 0034 % MINVEL: minimum vertical velocity threshold. 0035 % Scalar with the minimum vertical velocity threshold below which the 0036 % aproximated surge speed value is supposed to be unreliable. On samples 0037 % with an absolute vertical velocity value less than given threshold 0038 % flow speed is invalid (NaN). 0039 % Default value: 0 (all samples are valid) 0040 % MINPITCH: minimum pitch threshold. 0041 % Scalar with the minimum pitch threshold in radians below which the 0042 % aproximated surge speed value is supposed to be unreliable. On samples 0043 % with an absolute pitch value lesser than given threshold flow speed is 0044 % invalid (NaN). 0045 % Default value: 0 (all samples are valid) 0046 % 0047 % Notes: 0048 % This function is based on the flow speed computation code by Tomeu Garau in 0049 % function CORRECTTHERMALLAG. Main changes are: 0050 % - Changed vertical velocity approximation method (central instead of 0051 % forward differences). 0052 % - Added support for discarding regions with low velocity or low pitch. 0053 % 0054 % The approximate flow speed is computed from the absolute value of the 0055 % vertical velocity. The vertical velocity is approximated using central 0056 % differences, but accounting for the possibly unevenly spaced samples. 0057 % In non-vertical profiles with given pitch, the vertical velocity is scaled 0058 % by the sine of the pitch angle. The absolute value of this intermediate 0059 % quantity is called surge speed. To obtain the flow speed, each surge speed 0060 % value is scaled by a factor resulting from the evaluation of the given 0061 % polynomial at the surge speed value. If no polynomial given, the surge 0062 % speed is not scaled and equals the flow speed. 0063 % The note on Garau's original code was: 0064 % % The relative coefficient between the flow speed inside and outside of 0065 % % the conductivity cell. This is still uncertain (ask Gordon for origin). 0066 % % Here are three choices for first three orders polynomial. 0067 % speed_factor_polynoms = [0.00, 0.00, 0.40; % 0th order degree. 0068 % 0.00, 0.03, 0.45; % 1st order degree. 0069 % 1.58, 1.15, 0.70]; % 2nd order degree. 0070 % % First order approximation, second row of the matrix. 0071 % Surge speed is supposed to be unreliable when the vertical velocity or the 0072 % pitch is too low, so the flow speed is left undefined. In October 2013 0073 % Gerd Krahmann noted that there might be a bug in the original 0074 % implementation, and the matrix of polynomial coefficients (above) should be 0075 % transposed: 0076 % speed_factor_polynoms = [0.00, 0.00, 1.58; % 0th order degree. 0077 % 0.00, 0.03, 1.15; % 1st order degree. 0078 % 0.40, 0.45, 0.70]; % 2nd order degree. 0079 % 0080 % Examples: 0081 % % Vertical profile: 0082 % flow = computeCTDFlowSpeed(timestamp, depth, pitch) 0083 % % Tilted profile: 0084 % flow = computeCTDFlowSpeed(timestamp, depth, 25) 0085 % flow = computeCTDFlowSpeed(timestamp, depth, pitch) 0086 % % No flow speed factor scaling: 0087 % flow = computeCTDFlowSpeed(timestamp, depth, pitch, 'factorpoly', []) 0088 % 0089 % Authors: 0090 % Joan Pau Beltran <joanpau.beltran@socib.cat> 0091 0092 % Copyright (C) 2013-2016 0093 % ICTS SOCIB - Servei d'observacio i prediccio costaner de les Illes Balears 0094 % <http://www.socib.es> 0095 % 0096 % This program is free software: you can redistribute it and/or modify 0097 % it under the terms of the GNU General Public License as published by 0098 % the Free Software Foundation, either version 3 of the License, or 0099 % (at your option) any later version. 0100 % 0101 % This program is distributed in the hope that it will be useful, 0102 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0103 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0104 % GNU General Public License for more details. 0105 % 0106 % You should have received a copy of the GNU General Public License 0107 % along with this program. If not, see <http://www.gnu.org/licenses/>. 0108 0109 error(nargchk(2, 9, nargin, 'struct')); 0110 0111 0112 %% Parse mandatory input arguments. 0113 % Get numeric (non-option) arguments. 0114 narginnum = find(~cellfun(@isnumeric, varargin), 1, 'first') - 1; 0115 if isempty(narginnum) 0116 narginnum = nargin; 0117 end 0118 switch(narginnum) 0119 case 2 0120 % Vertical profile. 0121 [timestamp, depth] = varargin{1:2}; 0122 pitch = []; 0123 case 3 0124 % Tilted profile. 0125 [timestamp, depth, pitch] = varargin{1:3}; 0126 otherwise 0127 % Bad arguments. 0128 error('glider_toolbox:computeCTDFlowSpeed:InvalidArguments', ... 0129 'Invalid arguments (first 2 or 3 arguments should be numeric).'); 0130 end 0131 0132 0133 %% Configure default options. 0134 % This is NOT equivalent to the original version by Tomeu Garau. 0135 % The 1st degree flow speed factor polynomial is chosen, but using the proper 0136 % coefficients as noted by Gerd Krahmann. 0137 options.factorpoly = [0.00, 0.03, 1.15]; 0138 options.minvel = 0.0; 0139 options.minpitch = 0.0; 0140 0141 0142 %% Parse option arguments. 0143 % Get option key-value pairs in any accepted call signature. 0144 argopts = varargin(narginnum+1:end); 0145 if isscalar(argopts) && isstruct(argopts{1}) 0146 % Options passed as a single option struct argument: 0147 % field names are option keys and field values are option values. 0148 opt_key_list = fieldnames(argopts{1}); 0149 opt_val_list = struct2cell(argopts{1}); 0150 elseif mod(numel(argopts), 2) == 0 0151 % Options passed as key-value argument pairs. 0152 opt_key_list = argopts(1:2:end); 0153 opt_val_list = argopts(2:2:end); 0154 else 0155 error('glider_toolbox:findCTDFlowSpeed:InvalidOptions', ... 0156 'Invalid optional arguments (neither key-value pairs nor struct).'); 0157 end 0158 % Overwrite default options with values given in extra arguments. 0159 for opt_idx = 1:numel(opt_key_list) 0160 opt = lower(opt_key_list{opt_idx}); 0161 val = opt_val_list{opt_idx}; 0162 if isfield(options, opt) 0163 options.(opt) = val; 0164 else 0165 error('glider_toolbox:findCTDFlowSpeed:InvalidOption', ... 0166 'Invalid option: %s.', opt); 0167 end 0168 end 0169 0170 0171 %% Select full CTD rows. 0172 % The positive timestamp test is needed to deal with odd data from initial 0173 % lines in Slocum segment files. 0174 if numel(pitch) > 1 0175 valid = (timestamp(:) > 0) & ~isnan(depth(:)) & ~isnan(pitch(:)); 0176 timestamp_val = timestamp(valid); 0177 depth_val = depth(valid); 0178 pitch_val = pitch(valid); 0179 else 0180 valid = (timestamp(:) > 0) & ~isnan(depth(:)); 0181 timestamp_val = timestamp(valid); 0182 depth_val = depth(valid); 0183 pitch_val = pitch; 0184 end 0185 0186 0187 %% Compute glider surge speed from the vertical speed and the pitch if any. 0188 % For Slocum data, pitch is positive when nose is up (so positive z is down). 0189 % Use central differences scheme to compute 2nd order approximation, 0190 % but adapted to deal with irregular time sampling. 0191 vertical_velocity = zeros(size(timestamp_val)); 0192 dd = diff(depth_val); 0193 dt = diff(timestamp_val); 0194 ddt = timestamp_val(3:end) - timestamp_val(1:end-2); 0195 dd_dt = dd ./ dt; 0196 n = numel(timestamp_val); 0197 if n > 1 0198 vertical_velocity([1 end]) = dd([1 end]) ./ dt([1 end]); 0199 end 0200 if n > 2 0201 vertical_velocity(2:end-1) = ... 0202 (dt(2:end) .* dd_dt(1:end-1) + dt(1:end-1) .* dd_dt(2:end)) ./ ddt; 0203 % vertical_velocity(2:end-1) = (1.0 - dt(1:end-1)./ddt) .* dd_dt(1:end-1) ... 0204 % + (1.0 - dt(2:end )./ddt) .* dd_dt(2:end ); 0205 end 0206 if isempty(pitch) 0207 surge_speed = abs(vertical_velocity); 0208 else 0209 surge_speed = abs(vertical_velocity ./ sin(pitch_val)); 0210 end 0211 0212 0213 %% Discard unreliable surge speed estimates: 0214 low_vel = abs(vertical_velocity) < options.minvel; 0215 low_pitch = abs(pitch_val) < options.minpitch; 0216 surge_speed(low_vel | low_pitch) = nan; 0217 0218 0219 %% Compute flow speed inside cell from surge speed. 0220 % Do not scale surge speed if no polynomial is given. 0221 if isempty(options.factorpoly) 0222 flow_factor = 1; 0223 else 0224 flow_factor = polyval(options.factorpoly, surge_speed); 0225 end 0226 flow = nan(size(timestamp)); 0227 flow(valid) = flow_factor .* surge_speed; 0228 0229 end