function [K, Z] = oms_k1(infile, lambda,flag)
% OMS_K1    - Calculate diffuse attenuation coefficients 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. Reflectance by
% The waters surface is accouted for when propagatin Ed[3m+] to Ed[0-]
%
% Use as: [k Z] = oms_k1(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

S = oms_getinfo(infile);
v = var(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

% There may be more than one Ed sensor at the desired wavelength
for i = 1:length(goodEd)
   n = goodEd(i);
   ntest=char(ncnames(infile{n}));
   if strcmp(ntest(1:2),'Ed')
   vEd = v{n};									% Get the correct Es variable object
   
   %================
   % Get the depths
   %================
   if exist('flag','var')
      Z(:,i) = oms_getdepth(infile, vEd,0);
   else
      Z(:,i) = oms_getdepth(infile, vEd);
      %Z(:,i) = vEd.instrument_depth(:);
   end;
   
   Ed = vEd(:);
   bad = find(Ed == vEd.FillValue_(:));
   Ed(bad) = NaN;
   
   %===============================================
   % Calculate the diffuse attenuation Coefficient
   %===============================================
   % R is the % transmittance through the air-water interface. It is based on solar zenith 
   % angle only as does not take into account winds stress.
   solZen = solarzenith(nctime(infile), infile{'latitude'}(:), infile{'longitude'}(:));
   bad = find(solZen > 85);
   R = 1 - fresnel_(solZen);
   R(bad) = NaN;      
   R(find(R < 0)) = NaN;
   %R = 0.98
   K(:,i) = calc_k(Es.*R, Ed, Z(:,i)); 
   end;
end


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