function [K, Z] = oms_k2(infile, lambda)
% OMS_K2    - Calculate diffuse attenuation profiles from Oasis Mooring Spectroradiometer data
%
% The regime used for calculating the coefficients uses the Ed[3m+] and the 
% subsurface Ed to calculate the attenuation coefficient
%
% Use as: [k Z] = oms_k2(infile, lambda)
%
% Inputs: infile = A netcdf file object 
%         lambda = Wavelength of desired attenuation coeffiecient (nanometers)(a scalar)
%
% Output: K = attenuation coefficients
%         Z = depth of Ed sensor used to calculated k
%
% Note K and K may be MxN matrices if more than one subsurface Ed sensor using
% wavelength lambda is found

% Brian Schlining
% 09 Dec 1999

OK = oms_checkfile(infile);

if ~OK
   K = [];
   Z = [];
   return
end


v = var(infile);									% Get all the variables

for i = 1:length(v)
   s = lower(v{i}.symbol(:));
   if (isempty(s))
      symbol{i} = 'null';
   else
      symbol{i} = s;
   end
   
   wl = v{i}.wavelength(:);
   if (isempty(wl) | wl == v{i}.FillValue_(:))
      wavelength(i) = NaN;
   else
      wavelength(i) = wl;
   end
   
   b = v{i}.oasis_bank(:);
   if (isempty(b) | b == v{i}.FillValue_(:))
      bank(i) = NaN;
   else
      bank(i) = b;
   end
end
S = oms_getinfo(infile);

%===================
% Find Es at lambda
%===================
iEs = strmatch('es',S.symbol);					% Find Es
iWp = find(S.wavelength(iEs) == lambda);		% Find Es at the desired wavelength
goodEs = iEs(iWp);
if isempty(goodEs)
   K = msg(['Es[' num2str(lambda) 'nm] (Surface irradiance) sensors were not found']);
   Z = [];
   return
end
vEs = v{goodEs};									% Get the correct Es variable object
Es = vEs(:);
bad = find(Es == vEs.FillValue_(:));
Es(bad) = NaN;

%===================
% Find Ed at lambda
%===================
iEd = strmatch('ed',S.symbol);					% Find Ed
iWm = find(S.wavelength(iEd) == lambda);		% Find Ed at the desired wavelength
goodEd = iEd(iWm);
if isempty(goodEd)
   K = msg(['Ed[' num2str(lambda) 'nm] (Sub-Surface irradiance) sensors were not found']);
   Z = [];
   return
end

%======================================
% Get the depths and sort Eds by depth
%======================================
for i = 1:length(goodEd);
   Z(:,i) = oms_getdepth(infile, v{goodEd(i)});
   
end
[buf order] = sort(Z(1,:));
order  = fliplr(order);
goodEd = goodEd(order);
Z      = Z(:,order);


% There may be more than one Ed sensor at the desired wavelength
for i = 1:length(goodEd)
   
   
   if i == 1
      dz = Z(1,1);
   elseif i > 1
      Es = Edd; 
      dz = Z(1,i) - dz;
   end
   
   n = goodEd(i);
   vEd = v{n};									% Get the correct Es variable object
    
   Edd = vEd(:);
   bad = find(Edd == vEd.FillValue_(:));
   Edd(bad) = NaN;
   
   %===============================================
   % Calculate the diffuse attenuation Coefficient
   %===============================================
   K(:,i) = calc_k(Es, Edd, dz);
end


function k = msg(s)
warning(s)
k = [];
