function load_raw(file,  varargin)
% function load_raw(file, ...)
%
% 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 argument:
%     file: name of ENR or ENS file
%
% Optional named arguments:
%     "bin_path", bindir    : directory for *.cmd, *.bin
%     "nav_path", navdir    : directory for *.gps1 and *.gps2
%     "ens_range", [e0, e1] : start and end ensemble
%     "add_to_depth", dplus : meters to add to recorded depths
%     "second_set", num     : 1 for second set, otherwise 0
%     "interval", interval  : averaging interval in seconds
%
%
%
% 2000/09/22 EF
%


global fid_bin fid_cmd

bin_path = './';
nav_path = './';
ens_range = [1 10000000];
dplus = 0;
second_set = 0;
interval = 300;

p.max_ev = 1;
p.max_dev = 3;
p.use_3_beam = 0;    % 3-beam not yet implemented
p.ilayer = 5:40;

for iarg = 1:2:(nargin-2)
   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 'add_to_depth', dplus = val;
      case 'second_set',   second_set = val;
      case 'interval',     interval = 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']);

ens_ranges = find_ens_ranges(file, interval);
if isempty(ens_ranges), return; end
NP = size(ens_ranges, 2);

% Read the first profile to get the configuration
% and initialize the structure.  This is very fast.
if second_set,
   a = read_os(file, ens_range, [1 1], 'all', 'nav', 'second_set');
else
   a = read_os(file, ens_range, [1 1], 'all', 'nav');
end

depth = round(a.depth + dplus);
ND = length(depth);
config1 = mk_config1(a);
access = mk_access(a);    % Doesn't change.



fid_nav1 = fopen(nav1file, 'w');
fid_nav2 = fopen(nav2file, 'w');
if fid_nav1 == -1, error(['Can''t open ', nav1file]); end
if fid_nav2 == -1, error(['Can''t open ', nav2file]); end

fid_cmd = fopen(cmdfile, 'w');
fid_bin = fopen(binfile, 'w');
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);

for ip = 1:NP
   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
   tic
   if second_set,
      r = read_os(file, ens_ranges(:, ip), 'all', 'nav', 'second_set');
   else
      r = read_os(file, ens_ranges(:, ip), 'all', 'nav');
   end
   a.num_pings = length(r.dday);
   toc
   tic
   a = avg(r, a, p);
   toc

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


   %%%% nav and dday are now assuming ens files with good fixes!
   nav = ones(4,1) * 1e38;
   nav(1:2) = r.nav.txy2(2:3,end);  % Positions shown on the editing stagger plots.
   % 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;


   fprintf(fid_cmd, ['new_profile: ', to_date(a.config.yearbase, a.dday), '\n']);
   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

fclose(fid_nav1);
fclose(fid_nav2);


fclose(fid_cmd);
fclose(fid_bin);


%----------------------------------------------------
function bin_float(name, array)

global fid_bin fid_cmd

offset = ftell(fid_bin);
fwrite(fid_bin, array, 'float32');
fprintf(fid_cmd, 'binary_data: %s FLOAT %d %d\n', ...
                 name, length(array(:)), offset);

%---------------------------------------------------
function bin_double(name, array)

global fid_bin fid_cmd
offset = ftell(fid_bin);
fwrite(fid_bin, array, 'float64');
fprintf(fid_cmd, 'binary_data: %s DOUBLE %d %d\n', ...
                 name, length(array(:)), offset);

%---------------------------------------------------
function bin_ubyte(name, array)

global fid_bin fid_cmd
offset = ftell(fid_bin);
fwrite(fid_bin, array, 'uint8');
fprintf(fid_cmd, 'binary_data: %s UBYTE %d %d\n', ...
                 name, length(array(:)), offset);

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

function er = find_ens_ranges(file, interval)
          %** Modify for multiple files.

r = read_os(file);

np = length(r.dday);
if np < 11,   er = [];   return; end

delta = interval/86400;
d0 = r.dday(1):delta:r.dday(end-10); % Require at least 10.
n = length(d0);
if n == 1, er = [r.ens_num(1); r.ens_num(end)]; return; end

er = ones(2, n);

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

er(2, 1:(n-1)) = er(1, 2:n) - 1;
er(2, end) = r.ens_num(end);

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

function a = avg(r, a, p)
% r is a raw data structure; a is initialized for one profile
% p is a structure of parameters that control the averaging

enu = beam_enu(r.vel, r.heading, 0, 0, 'error', 1);

no_vel_mask = isnan(squeeze(enu(:, 1, :)));

ev_mask = abs(enu(:, 4, :)) > p.max_ev;
% Calculate other masks and combine with logical or.
bad_mask = squeeze(ev_mask);
% The following could probably be made more efficient by
% simply passing the bad mask into refavg, and optimizing
% refavg to use a single set of masks and counts for all
% iterations.

% maybe make refavg also operate
% on all 4 components at once.

for ii = 1:4
   uu = squeeze(enu(:, ii, :));
   uu(bad_mask) = NaN;
   a.vel(:,ii) = refavg(uu, p.ilayer, 5);
end

a.pg(:,:) = 0;
a.pg(:, 4) = round(100 *  sum(~(bad_mask | no_vel_mask), 2)  / a.num_pings);

a.cor = round(mean(r.cor, 3)); % Or, could use the masks to
                               % get the average of the values actually used.

a.amp = round(mean(r.amp, 3));


