function statrtn=cvt_nep2sed(infile,outfile, extmeta)
% cvt_nep2sed: computes sediment concentration from nep (raw volts)
%  uses do_hydraoptic equation to correct or compute sed.conc. from OBS.
% works on stats file at the moment, not bursts.
%
%  program in progress- use with caution
% emontgomery 3/3/09

%%% START USGS BOILERPLATE -------------
% Use of this program is self-described.
% Programs written in Matlab- as early as 7.2.0.232 (R2006a) and 
% most recent modifications in  v7.8.0 (R2009a)
% Programs ran on PC with Windows XP Professional OS, and RHEL4 linux.
%
% "Although this program has been used by the USGS, no warranty, 
% expressed or implied, is made by the USGS or the United States 
% Government as to the accuracy and functioning of the program 
% and related program material nor shall the fact of distribution 
% constitute any such warranty, and no responsibility is assumed 
% by the USGS in connection therewith."
%%% END USGS BOILERPLATE --------------

if exist(infile,'file')
    nc = netcdf(infile);
else
    [infile,inpath] = uigetfile('*.cdf','Choose input file');
    nc = netcdf(infile);
    if isempty(nc)
        disp (['Cannot open input file ' infile]);
        statrtn=0;
        disp (['Execution ends.']);
        return;
    end
end

[nvarin,chosen] = em_dollyvar(nc,'y');
nvarout=nvarin;

% Get the global attributes
keepatts = att(nc);
ngatts=length(keepatts);

% get the variablenames
for ivar = 1:nvarin;
    invar = nc{chosen{ivar}};
    invar = autonan(invar,1);
    inname{ivar} = name(invar);
end

% get the length of the data- this won't change
ntimes=length(nc{'time'}(:));

%get existing neph data 
% so far we only have one OBS/logger, so we don't need to handle more than
% one thing in loc (yet)
loc=strmatch('NEP',chosen);
loc2=strmatch('Sed',chosen);

extvolts=nc{chosen{loc}}(:);
  % find the number of the ext connection
  extno=chosen{loc}(4);

% this is the part from sohydraoptic
if ~isempty(findstr('OBS',extmeta.serial)),
      adj_evolts=extvolts-min(extvolts);
    if isfield(extmeta.cals,'coef') && ~isempty(extmeta.cals.coef),
        % user has supplied coefficients
         optic = min(extvolts).*ones(size(extvolts));      
        for n = 2:length(extmeta.cals.coef),
            optic = optic + extmeta.cals.coef(n).*(extvolts.^(n-1));
        end
    else
        % No coefficients, no conversion
        %disp('dohydraoptic: No coefficients supplied for the OBS')
        optic = [];
    end
else
    optic = adj_evolts;
end
% end part from do_hydraoptic

% Open output cdf
outc = netcdf(outfile,'noclobber');
if isempty(outc)
    [outfile,outpath] = uiputfile('*.nc','Choose output file');
    outc = netcdf(fullfile(outpath,outfile),'noclobber'); % MM 11/18/04 add path
    if isempty(outc)
        disp (['Cannot open output file ' outfile]);
        disp (['Data will be parked in temporary.nc']);
        outc = netcdf('temporary.nc','clobber');
    end
end

%Update global attributes and variables in output file.
history =['Sed_981 created or modified by cvt_nep2sed :' nc.history(:)];
outc.history = history;

vardesc=nc.VAR_DESC(:);
if isempty(strfind(vardesc,'sed'))  % add sed to the list, if missing
    vardesc = [vardesc ':sed' extno];
    nvarout=nvarout+1;
    loc2=nvarout;        % new cases of sed will be added to the end
end
outc.VAR_DESC = ncchar(vardesc);
outc.CREATION_DATE = ncchar(datestr(now,0));

% Dimensions:
outc('time') = 0;
outc('depth') = 1;
outc('lon') = 1;
outc('lat') = 1;

% Variables and attributes:
copy(nc{'time'},outc,0,1);
copy(nc{'time2'},outc,0,1);
copy(nc{'depth'},outc,0,1);
copy(nc{'lon'},outc,0,1);
copy(nc{'lat'},outc,0,1);

% Copy global attributes from input file to output file.
for i = 1 : length(keepatts);
    copy(keepatts{i},outc);
end

for i = 1: nvarout;
    disp (['creating output variable ' inname{i}]);
    ivar = nc{inname{i}};
    if loc2 == nvarout     % instantiate, if no existing sed variable
        sedname=['Sed' extno '_981'];
        nc{sedname} = ncfloat('time', 'depth', 'lat', 'lon')
        nc{sedname}.name = ncchar('Sed');
        nc{sedname}.long_name = ncchar('Sediment concentration ');
        nc{sedname}.generic_name = ncchar('sed');
        nc{sedname}.FORTRAN_format = ncchar('f10.2');
        nc{sedname}.units = ncchar('g/l');
        nc{sedname}.epic_code = nclong(981);
        nc{sedname}.code = ncchar('R');
        nc{sedname}.sensor_depth = ncfloat(0);
        nc{sedname}.serial_number = ncchar(' ');
        nc{sedname}.minimum = ncfloat(0);
        nc{sedname}.maximum = ncfloat(0);
        nc{sedname}.valid_range = ncfloat([0 1500]);
        nc{sedname}.FillValue_ = 1.0e35;
    else
        if ~isempty(ivar)
            copy(ivar,outc,0,1);
        end
    end
end

%disable define mode
    endef(outc)
    
%now put stuff in the variables
    outc{'lat'}(1) = nc{'lat'}(1);
    outc{'lon'}(1) = nc{'lon'}(1);
    outc{'depth'}(1) = nc{'depth'}(1);
    outc{'time'}(1:ntimes) = nc{'time'}(1:ntimes);
    outc{'time2'}(1:ntimes) = nc{'time2'}(1:ntimes);

for ik = 1:nvarout;
    if ik==loc2            %loc is index sed
        outc{chosen{ik}}(1:ntimes) = optic;
        outc{chosen{ik}}.minimum = ncfloat(min(optic));
        outc{chosen{ik}}.maximum = ncfloat(max(optic));
    else
        eval('outc{chosen{ik}}(1:ntimes) = nc{chosen{ik}}(1:ntimes);');
    end
end
% now close the ncfiles
close (outc);
close (nc);
statrtn=1;