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)<apsize(2);
    ap=ap';
end

apuncsize=size(apunc);
if apuncsize(1)<apuncsize(2);
    apunc=apunc';
end

wlsize=size(wl);
if wlsize(1)<wlsize(2);
    wl=wl';
end

wluncsize=size(wlunc);
if wluncsize(1)<wluncsize(2);
    wlunc=wlunc';
end

if acs==1;
    
    % Set up filter factors at every 0.1 nm from 1 to 799 nm, with center
    % wavelength at centwavel (i.e. at the data wavelengths)
    wavelength=.1:.1:799; % Thus index of 1 nm = 10; 356 nm= 3560;
    clear filtfunc
    SIG1=(-9.845*10^-8.*wl.^3  +1.639*10^-4*wl.^2- 7.849*10^-2*wl + 25.24)/2.3547 ;
    for i=1:max(size(wl));
        for jkl=1:max(size(wavelength));
            filtfunc(jkl,i)=(1/(sqrt(2*pi)*SIG1(i)))*exp(-0.5*((wavelength(jkl)-wl(i))/SIG1(i)).^2); % First term normalizes area under the curve to 1.
        end
    end
    
    error=10^-2;    % stop iterating when error in calculated ap is less than this from measured ap.
    numbits=10 ;    % number of iterations
    
    % Convolve the measurement with the fiter factors add the difference to
    % the measured spectrum to get the first corrected spectrum and iterate
    % to find the spectrum that when convolved with the filters, will give the
    % measured spectrum. This is the corrected absorption spectrum "abscorr".
    
    minwavel=min(wl);
    maxwavel=max(wl);
    
    for ii=1:numbits;
        xixi=minwavel:.1:maxwavel;% The range of centwavel is 0.1 nm.
        yiyi=spline(wl,ap,xixi); % Spline the measured data to every 0.1 nm.
        % We need data from 0 to 799 nm to multiply by filtfac.
        absspec=zeros(size(wavelength));
        absspec(minwavel*10:maxwavel*10)=yiyi;
        absspec(1:minwavel*10-1)=absspec(minwavel*10);
        absspec(maxwavel*10+1:size(wavelength))=0;
        aspecprime=absspec';
        
        clear measur2 meassignal6
        for i=1:min(size(filtfunc));% for all filters
            measur2(:,i)=aspecprime(:,1).*filtfunc(:,i); % the measured signal for every filter factor.
            meassignal6(i)=0.1*sum(measur2(:,i)); % The measured spectrum at a wavelength i is the sum of what a filter measured at
            %all wavelengths times 0.1 as we did it at every nm..
        end
        
        abscorr=ap-meassignal6'+ap;
        if max((ap-meassignal6')/ap)<=error
            break
        end     
        
        number_of_iterations=ii;  % number of iterations
        
    end
    
else
    
    abscorr=ap;
    
end

% Peak center values ("peak_loc") determined using a interative
% approach that allows the location to vary (uses the matlab
% function LSQNONLIN), and are rounded to nearest integer.
% Sigma values ("lsqsig") are determined similarly. FWHM = sigma*2.355
peak_loc=[406,434,453,470,492,523,550,584,617,638,660,675];

lsqsig=[16,12,12,13,16,14,14,16,13,11,11,10];

onenm = 400:1:720;

%interpolate the un-smoothed ap spectra to one nm resolution
acorr2onenm = interp1(wl,abscorr,onenm,'spline');

%define the matrix of component Gaussian functions using the peaks
%and widths (sigma) above
clear coef2
for ii=1:max(size(peak_loc));
    for jj=1:max(size(onenm));
        coef2(jj,ii) = exp(-0.5.*((onenm(jj)-peak_loc(ii))/lsqsig(ii))^2);
    end
end

%define a function for non-algal particles and concatenate this to the Gaussian matrix
nap=exp(-0.01*(onenm-400));
coef2nap(:,1) = nap;
coef2=[coef2nap,coef2];

%normalize both the component functions and the measured ap
%spectrum by the uncertainty (std dev) in the ap spectra
apunc_int = interp1(wlunc,apunc,onenm,'linear','extrap');
for iC=1:length(coef2(1,:));
    coef2_new(:,iC) = coef2(:,iC)./apunc_int';
end
acorr2onenm_new = acorr2onenm./apunc_int;

%Inversion analysis
amps = lsqnonneg(coef2_new,acorr2onenm_new');

% Using the inverted amplitudes, build a matrix of new component
% spectra (compspec) and the sum of the Gaussian and nap functions (sumspec):
for ll=1:max(size(onenm));
    sumspec(ll)=0;
    for kk=1:max(size(amps));
        sumspec(ll) = sumspec(ll)+sum(amps(kk)*coef2(ll,kk));    %sum of components
        compspec(ll,kk) = amps(kk)*coef2(ll,kk);  %gaussian components
    end
end

%interpolate back to the original resolution
compspec=interp1(onenm,compspec,wl,'spline');
sumspec=interp1(onenm,sumspec,wl,'spline');



