function r = load_raw(instclass, file, params,  varargin)
% function load_raw(instclass, file, params, ...)
%
% Writes the command (*.cmd) and binary (*.bin) files
% required by ldcodas, starting from VMDAS *.ENR or *.ENS.
% Each call to load_avg processes a single input file.
%                 ** Will change to process whole set.
% It also writes out files with (t x y) in a format suitable
% for refabs input.  The first fix from the ensemble goes to
% *.gps1, the last to *.gps2.
%
% Required arguments:
%     instclass:              instrument class (bb, nb, or os)
%     file: name of ENR or ENS file
%     params: a structure with fields
%         params.h_align      : Heading alignment in degrees;
%                                 RDI sign convention
%         params.interval     : Averaging interval in seconds (default 300)
%         params.ilayer       :  bin0:bin1, reference layer indices
%                                 for averaging.  Give the array of
%                                 indices, not just the first and last bin.
%         params.min_pings    : Minimum number of pings for an average;
%                                 only relevant at the end of file with
%                                 "finish" set to 1 (see below)
%         params.add_to_depth : meters to add to recorded depth,
%                                 if recorded transducer depth is incorrect
%         params.min_pg       : Within the reference layer, use only
%                                 the depths with at least this percent good;
%                                 if there are not at least 2 such levels,
%                                 don't use a reference layer, just calculate
%                                 a simple average.
%         params.max_ev       : Error velocity screening criterion, m/s
%         params.use_3_beam   : Not yet implemented; only 4-beam solutions
%                                 are used at present
%         params.bigtarget_amp       : 25
%         params.bigtarget_offset    :  5
%         params.bigtarget_avg       :  2
%                                          (See bigtarget.m)
%
%     >>The first three param fields are required;<<
%     the remainder have defaults.  Read the source to see what
%     they are, or eliminate all uncertainty by specifying all
%     param fields explicitly.
%
%     Other fields will likely be added as additional screening
%     methods are implemented.
%
% Optional named arguments: (keyword, value pairs)
%     "bin_path", bindir    : directory for *.cmd, *.bin
%     "nav_path", navdir    : directory for *.gps1 and *.gps2
%     "ens_range", [e0, e1] : start and end ensemble
%     "second_set", num     : 1 for second set, otherwise 0 (see below)
%     "leftover", b         : remaining profiles from last file
%     "finish", num         : 1 to finish, 0 to leave "leftover" for next call
%     "heading_offset", dh  : transducer heading offset in degrees
%
%        (second_set is normally 0, the default; set it to 1 only
%        if you have interleaved BB and NB pings and want to load
%        the NB part of the data set; leave it at 0 for the BB data,
%        or for any case when only BB or NB pings were collected.)
%
% 2000/09/22 EF
%

% 2001/01/25 JH adding instclass

global fid_bin fid_cmd

fid_cmd = 0;  % initialize so it can serve as a "files are open" flag

bin_path = './';
nav_path = './';
ens_range = [1 10000000];
second_set = 0;
leftover = [];
finish = 1;

% The following parameter structure will become an input, but
% for now it is hardwired here.
%  p.h_align =  43.28; %%% Darryl's value for Kaiyo, ky0006L3 later files
%  p.interval = 300;
%  p.ilayer = 2:40;
p.min_pg = 25;
p.max_ev = 1;
p.use_3_beam = 0;    % 3-beam not yet implemented
p.min_pings = 10;
p.add_to_depth = 0;
p.bigtarget_amp = 25;
p.bigtarget_offset = 5;
p.bigtarget_avg = 2;

if nargin < 3 | ~ischar(file) | ~isstruct(params),
   help load_raw
   fprintf('At least 3 arguments are required:\n')
   fprintf('an instrument class, a filename and a structure.')
   error('too few input arguments ')
end

if (~isfield(params, 'h_align') | ~isfield(params, 'interval') ...
      | ~isfield(params, 'ilayer')),
   help load_raw
   error('Parameter fields h_align, interval, and ilayer are required')
end

% Transfer input parameters, overriding defaults.
fn = fieldnames(params);
for ii = 1:length(fn)
   fname = fn{ii};
   p = setfield(p, fname, getfield(params, fname));
end


for iarg = 1:2:(nargin-3)
   arg = varargin{iarg};
   val = varargin{iarg+1};
   switch lower(arg)
      case 'bin_path',     bin_path = val;
      case 'nav_path',     nav_path = val;
      case 'ens_range',    ens_range = val;
      case 'second_set',   second_set = val;
      case 'leftover',     leftover = val;
      case 'finish',       finish = val;
      otherwise,           error('invalid argument name')
   end
end

[a_path, a_file, a_ext] = fileparts(file);
cmdfile = fullfile(bin_path, [a_file a_ext '.cmd']);
binfile = fullfile(bin_path, [a_file a_ext '.bin']);
nav1file = fullfile(nav_path, [a_file a_ext '.gps1']);
nav2file = fullfile(nav_path, [a_file a_ext '.gps2']);

fprintf('getting ens_ranges...\n')
ens_ranges = find_ens_ranges(instclass, leftover, file, ens_range, p.interval);
NP = size(ens_ranges, 2);

% Read the first profile to get the configuration
% and initialize the structure for the average.  This is very fast.

if second_set,
   a = read(instclass, file, 'ends', 1 ,...
            'vars', {'all', 'nav'},   'second_set', 1 );
else
   a = read(instclass, file, 'ends', 1, 'vars', {'all', 'nav'});
end

depth = round(a.depth + p.add_to_depth);
ND = length(depth);
config1 = mk_config1(instclass, a);
ii = config1_ind;
config1(ii.hd_offset) = p.h_align;
access = mk_access(instclass, a);    % Doesn't change.


for ip = 1:NP
   fprintf('%d out of %d (%s)\n', ip, NP, file)
   tic

   
   if second_set,
      r = read(instclass, file, 'range', ens_ranges(:, ip),...
               'vars', {'all', 'nav'},   'second_set', 1 );
   else
      r = read(instclass, file, 'range', ens_ranges(:, ip),...
               'vars', {'all', 'nav'});
   end


   if ip == 1 & ~isempty(leftover)
      r = catadcp(instclass, leftover, r);
   end

   if ip == NP
      if ~finish | (diff(ens_ranges(:, NP)) < p.min_pings - 1),
         break,
      end
   end

   a.num_pings = length(r.dday);
   
   a = avg(instclass, r, a, p);
   toc

   % all the structures:
   ancil1 = mk_ancil1(instclass, r, 1);
   ancil2 = mk_ancil2(instclass, r, 1);
   bt = mk_bt(instclass, r, 1);


   %%%% nav and dday are now assuming ens files with good fixes!
   nav = ones(4,1) * 1e38;
   % NAVIGATION structure from codas3/adcp/include/adcp.h is [lat lon]
   % Positions shown on the editing stagger plots.
   nav(1:2) = r.nav.txy2([3 2],end);  
   % The velocity part is not used for anything.

   [mm, ss, nn, mdt] = mstdgap(r.nav.sec_pc_minus_utc);
   a.dday = r.dday(end) - mdt / 86400;

   if ip == 1,
      fid_nav1 = fopen(nav1file, 'w');
      fid_nav2 = fopen(nav2file, 'w');
      fid_cmd = fopen(cmdfile, 'w');
      fid_bin = fopen(binfile, 'w');
      if fid_nav1 == -1, error(['Can''t open ', nav1file]); end
      if fid_nav2 == -1, error(['Can''t open ', nav2file]); end
      if fid_cmd == -1, error(['Can''t open ', cmdfile]); end
      if fid_bin == -1, error(['Can''t open ', binfile]); end
      fprintf(fid_cmd, 'binary_file: %s\n', binfile);
   end

   if ip == 1 | (rem(ip, 400) == 1 &  (NP - ip) > 10),
      fprintf(fid_cmd, 'new_block\n');
      fprintf(fid_cmd, 'dp_mask: 1\n');

      bin_float('DEPTH', depth);
      bin_double('CONFIGURATION_1', config1);
   end
   fprintf(fid_cmd, ['new_profile: ', to_date(a.config.yearbase, a.dday),...
              '  /* %d %d */\n'], r.ens_num([1 end]));
   fprintf(fid_cmd, 'depth_range: %d %d\n', round(depth([1 end])));
   bin_double('ACCESS_VARIABLES', access);
   bin_double('ANCILLARY_1', ancil1);
   bin_double('ANCILLARY_2', ancil2);
   bin_double('BOTTOM_TRACK', bt);
   bin_double('NAVIGATION', nav);
   bin_float('U', a.vel(:, 1));
   bin_float('V', a.vel(:, 2));
   bin_float('W', a.vel(:, 3));
   bin_float('ERROR_VEL', a.vel(:, 4));
   bin_ubyte('AMP_SOUND_SCAT', round(mean(a.amp, 2)));
   bin_ubyte('RAW_AMP', a.amp);
   bin_ubyte('PERCENT_GOOD', a.pg(:, 4));
   bin_ubyte('PERCENT_3_BEAM', a.pg(:, 3));
   bin_ubyte('SPECTRAL_WIDTH', round(mean(a.cor, 2)));
   bin_ubyte('RAW_SPECTRAL_WIDTH', a.cor);
   %% Maybe we should call it something else, but this will do
   %% for now; it is not correlation, either.
   bin_ubyte('PROFILE_FLAGS', zeros(1, ND));
   % ......

   fprintf(fid_nav1, '%12.7f %12.7f %12.7f\n', r.nav.txy1(:, 1));
   fprintf(fid_nav2, '%12.7f %12.7f %12.7f\n', r.nav.txy2(:, end));


end

if finish,
   r = [];
elseif fid_cmd,
   fprintf(fid_cmd, '/* leftover: %d %d */\n', r.ens_num([1 end]));
end

if fid_cmd
   fclose(fid_nav1);
   fclose(fid_nav2);
   fclose(fid_cmd);
   fclose(fid_bin);
end

%====================================================

function er = find_ens_ranges(instclass, leftover, file, ens_range, interval)

r = read(instclass, file, 'range', ens_range);
if ~isempty(leftover)
   dday = [leftover.dday; r.dday];
   ens = [leftover.ens_num; r.ens_num];
else
   dday = r.dday;
   ens = r.ens_num;
end

delta = interval/86400;
d0 = dday(1):delta:dday(end);

n = length(d0);

if n == 1,
   er = [ens(1); ens(end)];
else
   er = ones(2, n);

   for ii = 2:n
      er(1, ii) = ens(max(find(dday <= d0(ii))));
   end

   er(2, 1:(n-1)) = er(1, 2:n) - 1;
   er(2, end) = ens(end);
end


%-----------------------------------------------------------

function mask = mask_below(depth, bottom)
np = length(bottom);
nd = length(depth);
depth = depth(:, ones(np, 1));
bottom = bottom(ones(nd, 1), :);

mask = depth > bottom;


