function [amps,compspec,sumspec] = spectral_decomp(wl,ap,wlunc,apunc,acs)
% AC-S "un-smoothing" and spectral decomposition method
%
% Ron Zaneveld, WET Labs, Inc., 2005
% Ali Chase, University of Maine, 2014
%
% See the following publication for details on the method:
% Chase, A., et al., Decomposition of in situ particulate absorption
% spectra. Methods in Oceanography (2014), http://dx.doi.org/10.1016/j.mio.2014.02.022
%
% INPUTS:
% wl - wavelengths associated with measured particulate absorption
% spectra
% ap - measured particulate absorption spectra
% wlunc - wavelengths associated with uncertainty of measured
% particulate absorption spectra
% apunc - uncertainty of measured particulate absorption spectra
% acs - 1 (for AC-S data) or 0 (for a different instrument). This code
% was designed for use with particulate absorption spectra measured
% using a WETLabs AC-S instrument, and thus it applies a spectral
% un-smoothing that is specific to that instrument. To use with
% other types of absorption data, input "0" and the correction
% applied for AC-S will be bypassed.
%
% The uncertainty values we use are the standard deviation of one-minute
% binned particulate absorption spectra. If the uncertainty is unknown,
% then all wavelenghts will be treated equally, i.e. no spectral weighting
% will be applied. In this case populate apunc with -999 values.
%
% OUTPUTS:
% amps - the amplitudes of 12 Gaussian function peaks and one non-algal
% particle (NAP) function
% compspec - the component spectra after spectral decomposition; the
% Gaussian and NAP functions multiplied by the 'amps'
% values
% sumspec - the sum of the component spectra. will be simliar to the
% original measured spectrum, but the fit may not be as good in
% parts of the spectrum with higher uncertainty.
%
% The amplitudes of the Gaussian functions represent the relative amounts
% of different phytoplankton pigments and can be compared to HPLC pigment
% concentrations to evaluate the method. See Chase et al. (2014) for
% detail on such evaluation (reference above).
if min(apunc)<0
apunc=ones(max(length(ap)),1);
end
apsize=size(ap);
if apsize(1)