function load_avg(instclass, file,  varargin)
% function load_avg(instclass, file, ...)
%
% Writes the command (*.cmd) and binary (*.bin) files
% required by ldcodas, starting from VMDAS *.STA or *.LTA.
% Each call to load_avg processes a single input file.
%
% If data are from an Ocean Surveyoror a Broad Band adcp, it also writes 
% out files with (t x y) in a formula suitable for refabs input.  
% The first fix from the  ensemble goes to *.gps1, the last to *.gps2.
%
% Required arguments:
%     instclass: (bb or os)
%     file: name of STA or LTA file
%
% Optional named arguments:
%     "bin_path", bindir    : directory for *.cmd, *.bin
%     "nav_path", navdir    : directory for *.gps1 and *.gps2
%     "irange", [starti  endi] : start and end index
%     "add_to_depth", dplus : meters to add to recorded depths
%     "second_set", num     : 1 for second set, otherwise 0  (OS only)
%
%
%
% 2000/09/22 EF, written for os
%
% 2001/01/25 JH adding instclass
% 2001/06/05 JH adding bb


global fid_bin fid_cmd

bin_path = './';
nav_path = './';
irange = [1 10000000];
dplus = 0;
second_set = 0;

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 'irange',    irange = val;
      case 'add_to_depth', dplus = val;
      case 'second_set',   second_set = 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']);


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

%if length(a.dday) < 2,
%   disp(['Warning: ', file, ' has less than 2 ensembles and will be ignored'])
%   return
%end

%??? see a.dday below
if length(a.dday) < 3,
   disp(['Warning: ', file, ' has less than 3 ensembles and will be ignored'])
   return
end

a.dday = a.dday - a.nav.sec_pc_minus_utc / 86400;
% This is the ensemble start time, but we need the end time:
a.dday = [a.dday(2:end); 1];
a.dday(end) = 2 * a.dday(end-1) - a.dday(end-2); %% ??? WHAT IA THIS???

NP = length(a.heading);

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

fprintf(fid_nav1, '%12.7f %12.7f %12.7f\n', a.nav.txy1);
fprintf(fid_nav2, '%12.7f %12.7f %12.7f\n', a.nav.txy2);

fclose(fid_nav1);
fclose(fid_nav2);

% all the structures:
access = mk_access(a);
config1 = mk_config1(a);
ancil1 = mk_ancil1(a);
ancil2 = mk_ancil2(a);
depth = round(a.depth + dplus);
ND = length(depth);
bt = mk_bt(a);

nav = zeros(4,NP) * NaN;
% NAVIGATION structure from codas3/adcp/include/adcp.h is [latitude longitude]
% Positions shown on the editing stagger plots.
nav(1:2,:) = a.nav.txy2([3 2],:);  
% The velocity part is not used for anything.

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

   fprintf(fid_cmd, ['new_profile: ', to_date(a.config.yearbase, a.dday(ip)), '\n']);
   fprintf(fid_cmd, 'depth_range: %d %d\n', round(depth([1 end])));
   bin_double('ACCESS_VARIABLES', access);
   bin_double('ANCILLARY_1', ancil1(:,ip));
   bin_double('ANCILLARY_2', ancil2(:,ip));
   bin_double('BOTTOM_TRACK', bt(:,ip));
   bin_double('NAVIGATION', nav(:,ip));
   bin_float('U', a.vel(:, 1, ip));
   bin_float('V', a.vel(:, 2, ip));
   bin_float('W', a.vel(:, 3, ip));
   bin_float('ERROR_VEL', a.vel(:, 4, ip));
   bin_ubyte('AMP_SOUND_SCAT', round(mean(a.amp(:, :, ip), 2)));
   bin_ubyte('RAW_AMP', squeeze(a.amp(:, :, ip)));
   bin_ubyte('PERCENT_GOOD', a.pg(:, 4, ip));
   bin_ubyte('PERCENT_3_BEAM', a.pg(:, 3, ip));
   bin_ubyte('SPECTRAL_WIDTH', round(mean(a.cor(:, :, ip), 2)));
   bin_ubyte('RAW_SPECTRAL_WIDTH', squeeze(a.cor(:, :, ip)));
   %% 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));
   % ......

end

fclose(fid_cmd);
fclose(fid_bin);


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