/
estimate_bbp.m
128 lines (122 loc) · 4.45 KB
/
estimate_bbp.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
function [ bbp, beta_p, sigma_bbp, sigma_beta_p ] = estimate_bbp( beta, t, s, lambda, theta, X_p, sigma_X_p, delta )
%ESTIMATE_BBP The backscattering coefficient of particles bbp is commonly
% estimated from measurement of scattering at a single angle in the
% backward hemisphere beta.
%
% $$\beta_p(\theta) &= \beta(\theta) - \beta_{sw}(\theta)\\$$
% $$b_{bp} &= 2 \times \pi \times \chi_p(\theta) \times \beta_{p}(\theta)$$
%
% beta_sw scattering at a single angle from sea water is estimated with:
% Xiaodong Zhang, Lianbo Hu, and Ming-Xia He (2009), Scatteirng by pure
% seawater: Effect of salinity, Optics Express, Vol. 17, No. 7, 5698-5710
%
% chi_p (X_p) the particulate conversion coefficient comes from:
% James M. Sullivan, Michael S. Twardowski, J. Ronald, V. Zaneveld, and
% Casey C. Moore. Measuring optical backscattering in water. 2013.
% ISBN 978-3-540-25315-0. doi: 10.1007/3-540-37672-0.
% URL http: //www.springerlink.com/index/10.1007/3-540-37672-0.
%
%Syntax: [ bbp, beta_p ] = estimate_bbp( beta, t, s, lambda, theta, X_p, delta )
%
%Inputs:
% Required:
% beta NxM double corresponding to the values of the VSF at one angle in m^{-1}.sr^{-1}
% t Nx1 double corresponding to the temperature in deg C
% s Nx1 double corresponding to the salinity in psu
% Optional:
% lambda 1x1 or 1xM double corresponding to the wavelength in nm
% default: 700
% theta 1x1 or 1xM double corresponding to the scattering angle in deg
% default: 142 (ECO-FLBB and ECO-FLNTU)
% MCOM is 150
% other single and three channels are 124
% X_p 1x1 or 1xM double, chi_p(theta) conversion coefficient at theta
% default: interpolated from Sullivan et al. (2013)
% sigma_X_p 1x1 or 1xM double, sigma_chi_p(theta) relative error of
% conversion coefficient at theta. Used only for sigma_bbp
% default: interpolated from Sullivan et al. (2013)
% delta double corresponding to the delta for Zhang et al. (2009)
% default: 0.039
%
%Outputs:
% bbp NxM double corresponding to the particulate backscattering in m^{-1}
% beta_p NxM double corresponding to particulate beta at the angle theta in m^{-1}.sr^{-1}
% sigma_bbp NxM double corresponding to the relative error of bbp in (%)
% sigma_beta_p double corresponding to the relative error of beta_p (%)
%
%
%Examples:
% [bbp_BB9, beta_p_BB9] = estimate_bbp(bb9, t, s, wl_bb9, 117);
% [bbp_VSF, beta_p_VSF] = estimate_bbp(eco_vsf, t, s, 650, scat_angle_vsf);
%
%m-files required: betasw_ZHH2009
%
% Tested with: Matlab R2015a
% Matlab R2015b
% Matlab R2016a
%
% Author: Nils Haentjens, Ms, University of Maine
% Email: nils.haentjens@maine.edu
% Created: 11th June 2015
% Updated: 25th August 2016
% Check Nargin
if nargin > 7
error('Too many input arguments.')
elseif nargin < 3
error('Not enough input arguments.')
end;
% Set default param
if ~exist('lambda','var');
lambda = 700;
end;
if ~exist('theta','var');
theta = 142;
end;
if ~exist('delta','var');
delta = 0.039;
end;
if ~exist('X_p','var');
% Interpolate X_p with values from Sullivan et al. 2013
theta_ref = 90:10:170;
X_p_ref = [0.684 0.858 1.000 1.097 1.153 1.167 1.156 1.131 1.093];
X_p = interp1(theta_ref, X_p_ref, theta, 'spline');
end;
if ~exist('sigma_X_p','var');
% Interpolate sigma_X_p with values from Sullivan et al. 2013
theta_ref = 90:10:170;
sigma_X_p_ref = [0.034 0.032 0.026 0.032 0.044 0.049 0.054 0.054 0.057];
sigma_X_p = interp1(theta_ref, sigma_X_p_ref, theta, 'linear');
end;
% Set relative error of beta_sw from Zhang et al. 2009
sigma_beta_sw = 0.02;
% Check size of input
if size(beta) ~= size(t) | size(beta) ~= size(s)
error('beta, t and s should have same size');
end;
% Check size of lambda, theta and X_p
if size(lambda,2) == 1;
lambda = lambda * ones(size(beta,2),1);
end;
if size(theta,2) == 1;
theta = theta * ones(size(beta,2),1);
end;
if size(X_p,2) == 1;
X_p = X_p * ones(size(beta,2),1);
end;
if size(sigma_X_p,2) == 1;
sigma_X_p = sigma_X_p * ones(size(beta,2),1);
end;
% compute beta_sw
beta_p = zeros(size(beta));
sigma_beta_p = sigma_beta_sw;
bbp = zeros(size(beta));
sigma_bbp = zeros(size(beta));
for i=1:size(beta,1);
for j=1:size(beta,2);
beta_sw = betasw_ZHH2009(lambda(j), t(i), theta(j), s(i), delta);
beta_p(i,j) = beta(i, j) - beta_sw;
bbp(i,j) = 2 * pi * X_p(j) * beta_p(i,j);
sigma_bbp(i,j) = sigma_beta_sw + sigma_X_p(j);
end;
end;
end