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 ue4 or EF RAW nb logging
% Each call to load_raw 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 RAW
%     params: a structure with fields
%         params.h_align      : Heading alignment in degrees;
%                                 RDI sign convention  (DAS OH)
%         params.h_bias       : Heading Bias           (DAS HB)
%         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.geometry     : 'convex' (OS) 'concave' (NB) 
%         params.yearbase     : yearbase (*** required *** (unlike os)
%         params.add_to_depth : transducer depth (required)  ***(unlike os)***
%
%         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.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 four 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)
%     "ldcodas_path", bindir: directory to write *.cmd, *.bin
%     "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
%     "use_external_gyro"   : def:0 (internally recorded) 1 (get from hdg file)
%
%   %%if any of the following files exist, their data will be added to the mix
%   %% before averaging
%
%     time_path             : path to [logger_dday best_dday] gbin file
%                                default '../gbinfiles/time'
%     time_ext              : extension for .???.gbin time file
%                                default 'tim'
%
%     position_path         : path to [ lon lat] gbin file
%                                default '../gbinfiles/pcode'
%     position_ext          : extension for .???.gbin position file
%                                default 'pos'
%
%     ashtech_path          : path to [heading pitch roll] gbin file
%                                default '../gbinfiles/adu2'
%     ashtech_ext           : extension for .???.gbin ashtech file
%                                default 'adu'
%
%     gyro_path             : path to [meas_heading] gbin file
%                                default '../gbinfiles/gyro'
%     gyro_ext              : extension for .???.gbin gyro heading file
%                                default 'hdg'
%
%     sound_path            : path to [meas_soundspeed] gbin file
%                                default '../gbinfiles/sndspd'
%     sound_ext             : extension for .???.gbin soundspeed file
%                                default 'snd'


%
% 2000/09/22 EF
% 2001/01/25 JH adding instclass, morphing for NBADCP; get nav fields 
%            from gbin files
%            

global fid_bin fid_cmd

if (nargin <= 4)
   help(mfilename)
   error('too few variables')
end

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

fid_log = fopen('load_raw.log','a');
cc = clock;
fprintf(fid_log,'running load_raw on %s at %s\n',...
        file, to_date(cc(1), to_day(cc(1), cc)));
clear cc

%set defaults for varargin structure 
cfg.ldcodas_path = './';
cfg.nav_path = './';     
cfg.att_path = './';
cfg.second_set = 0;
cfg.leftover = [];
cfg.finish = 1;
cfg.use_external_gyro = 0;
cfg.time_path =  '../gbinfiles/time';
cfg.time_ext =  'tim';
cfg.position_path =  '../gbinfiles/pcode';
cfg.position_ext =  'pos';
cfg.ashtech_path =  '../gbinfiles/adu2';
cfg.ashtech_ext =  'adu';
cfg.gyro_path =  '../gbinfiles/gyro';
cfg.gyro_ext =  'hdg';
cfg.sound_path =  '../gbinfiles/sndspd';
cfg.sound_ext =  'spd';


cfg = fillstruct(cfg, varargin, 0);

% rename parameters to match code
ldcodas_path = cfg.ldcodas_path;
nav_path = cfg.nav_path;     %don't write nav yet
leftover = cfg.leftover;
finish = cfg.finish;

% The following parameter structure will become an input, but
% for now it is hardwired here.
%  p.h_bias  = 0;
%  p.h_align =  43.28; %%% Darryl's value for Kaiyo, ky0006L3 later files
%  p.interval = 300;
%  p.ilayer = 2:40;
p.yearbase = 0;     % for string-formatting to_date correctly
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 < 4 | ~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, 'h_bias') | ...
    ~isfield(params, 'interval') | ...
    ~isfield(params, 'ilayer') | ...
    ~isfield(params, 'geometry') | ...
    ~isfield(params, 'add_to_depth') | ...
    ~isfield(params, 'yearbase'));
   help load_raw
   fprintf('required Parameter fields: \n');
   error('h_align, h_bias, interval, geometry, ilayer, yearbase, add_to_depth ')
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

[a_path, rawfilebase, a_ext] = fileparts(file);
cmdfile = fullfile(ldcodas_path, [rawfilebase  '.cmd']);
binfile = fullfile(ldcodas_path, [rawfilebase  '.bin']);
nav1file = fullfile(nav_path, [rawfilebase  '.gps1']);  
nav2file = fullfile(nav_path, [rawfilebase  '.gps2']);  
attfile =  fullfile(nav_path, [rawfilebase  '.attstats.bin']);  

% Read the first and last profile to get the configuration, get start and 
% end dates, initialize the structure  for the average.  This is very fast.
% AND return the fid
a = read(instclass, file,  'yearbase', p.yearbase, 'ends', 1,...
         'vars', 'all'); %no nav in nb
finfo = a.info; %pass this in to subsequent reads; 



timbin_filename = ...
    fullfile(cfg.time_path, sprintf('%s.tim.gbin', rawfilebase));
tim_binfile = read_bin(timbin_filename);
if (isempty(tim_binfile))
   r = [];
   return
end

   
p.dday = tim_binfile(2,:)';  %best time


% reset leftover to [] if we're starting a new configuration segment
if (isempty(leftover))
   fprintf('leftover is empty\n')
end
if (a.ens_num(1) == 1)
   leftover = [];
   fprintf('starting a new configuration segment: resetting leftover to []\n')
end

if ~(isempty(leftover))
   if (p.dday(1) - leftover.dday(end) > .5*(p.interval/86400))
      leftover = [];
   end
end

fprintf('calculating iranges\n')
index_ranges = find_iranges(leftover,  p); 
NP = size(index_ranges, 2);

%% temporary check
sec_res = mstdgap(86400*diff(p.dday));
for ip = 1:NP
   secs = 86400*(p.dday(index_ranges(2,ip)) -p.dday(index_ranges(1,ip)));
   if (secs > p.interval + sec_res)
      warning(sprintf('averaging interval at ip = %d is too big: %f secs\n',...
                      ip, secs))
      keyboard
   end
end


% note on the proliferation of "depth" variables
% depth: temporary variable just to fill the structure a.depth
% a.depth: depth of bins for average
% r.depth (below): raw depth
% override a.config.tr_depth (starts as 0 in @nb/read.m)
a.config.tr_depth = p.add_to_depth;
a.depth = round(mk_depth(a, a.config.tr_depth)); 
NC = length(a.depth);  %from a.config.ncells (i.e. number of cells)
config1 = mk_config1(a);
c1_index =  config1_ind;
config1(c1_index.avg_interval)   = params.interval;
config1(c1_index.hd_offset)      = params.h_align;
config1(c1_index.heading_bias)   = params.h_bias;
access = mk_access(a);    % Doesn't change.

%get all the auxilliary data for this raw file
% ("data" can simply be data.info.rawfilenames and data.info.ilist)
%must replace ilist!!!
a.info.ilist = [1:length(p.dday)];
aux_struct = get_serfields(...
   a,'verbose',1, 'use_external_gyro',cfg.use_external_gyro,...
   'timpath', cfg.time_path,    'timgbin', [cfg.time_ext,    '.gbin'],...
   'pospath', cfg.position_path,'posgbin', [cfg.position_ext,'.gbin'],...
   'adupath', cfg.ashtech_path, 'adugbin', [cfg.ashtech_ext ,'.gbin'],...
   'hdgpath', cfg.gyro_path,    'hdggbin', [cfg.gyro_ext,    '.gbin'],...
   'spdpath', cfg.sound_path,   'spdgbin', [cfg.sound_ext,   '.gbin']);
aux_fnames = fieldnames(aux_struct);
att_fnames = {'aduheading', 'attdday', 'mrms', 'brms'};
if (sum(ismember(att_fnames, aux_fnames)) ~= length(att_fnames))
   noashtech = 1;
else
   noashtech = 0;
end



for ip = 1:NP

   tic
   
   ilist = [index_ranges(1, ip) : index_ranges(2, ip)];

   r = read(instclass, finfo,  'yearbase',  p.yearbase, ...
            'ilist', ilist , 'vars', 'all');

   
   if (isempty(r.dday))
      disp('r.dday is empty: in load_raw now.  check it out')
      keyboard
   end

   
   r.depth = mk_depth(r); %a.config.tr_depth set to 0 in read.m
   r.depth = round(r.depth + p.add_to_depth);
   fprintf('%d out of %d (%s) \n', ip, NP, file)

   if (length(ilist) <= 1)
      disp('warning: ilist has 0 or 1 element: check it out')
      keyboard
   end
   
   
   aux = cutadcp(aux_struct, ilist);
   
   for fnamei = 1:length(aux_fnames);
      fname = aux_fnames{fnamei};
      r = setfield(r, fname, getfield(aux, fname));
   end


   fprintf(fid_log, ...
           '%d out of %d, indices [%d %d], dd: %7.3f - %7.3f (%s %s)\n',...
           ip, NP, index_ranges(1, ip), index_ranges(2, ip), ...
           r.dday([1 end]), to_date(p.yearbase, r.dday(1)), ...
           to_date(p.yearbase, r.dday(end)));

   if ip == 1 & ~isempty(leftover)
      if (noashtech)
         for fnamei = 1:length(att_fnames)
            if (isfield(leftover, att_fnames{fnamei}))
               leftover = rmfield(leftover, att_fnames{fnamei});
            end
            if (isfield(r, att_fnames{fnamei}))
               r = rmfield(r, att_fnames{fnamei});
            end
         end
      end
      r = catadcp(leftover, r);
   end
   
   if ip == NP
      if ~finish | (diff(index_ranges(:, NP)) < p.min_pings - 1),
         break,
      end
   end
   
   a.num_pings = length(r.dday);
   
   %get ashtech statistics
   if ~noashtech
      [attstats, att_rownames] = att_stats(aux);
   end
   
   a = avg(instclass, r, a, p);

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

   
   nav = ones(4,1) * 1e38;   
   %NAVIGATION structure wants [latitude longitude]
   %but nav structure has [time longitude latitude]
   
   % The velocity part is not used for anything.
   a.dday = r.dday(end);
   
   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', a.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(a.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.sw, 2)));
   bin_ubyte('RAW_SPECTRAL_WIDTH', a.sw);
   %% Maybe we should call it something else, but this will do
   %% for now; it is not correlation, either.
   bin_ubyte('PROFILE_FLAGS', zeros(1, NC));
   % .....
   


   %make sure there are values for navigation 
   posi = min(find(~isnan(sum(r.nav.txy1))));
   if (~isempty(posi))
      fprintf(fid_nav1, '%12.7f %12.7f %12.7f\n', r.nav.txy1(1, 1),...
              r.nav.txy1(2, posi), r.nav.txy1(3, posi));
   else
      fprintf(fid_nav1, '%%%12.7f \n', r.nav.txy1(1, 1));
   end
   posi = max(find(~isnan(sum(r.nav.txy2))));
   if (~isempty(posi))
      fprintf(fid_nav2, '%12.7f %12.7f %12.7f\n', r.nav.txy1(1, end),...
              r.nav.txy1(2, posi), r.nav.txy1(3, posi));
   else
      fprintf(fid_nav2, '%%%12.7f\n', r.nav.txy1(1, end))
   end
   
   if (~noashtech)
      if (ip == 1)
         attstruct = write_bin(attfile, attstats, 'att_stats', att_rownames); 
      else
         attstruct = write_bin(attstruct, attstats);
      end
   end

   toc
end

%close the att file
if (~noashtech)
   fclose(attstruct.fid);
   fprintf('closing the att bin file %s\n', attfile)
end



if finish,
   r = [];
elseif fid_cmd,
   fprintf(fid_cmd, '/* leftover: %d %d */\n', r.ens_num([1 end]));
end
fprintf(fid_log, '/* leftover: %d %d */\n', r.ens_num([1 end]));
fprintf(fid_log, 'leftover dday range is %f - %f (%s - %s)\n', ...
        r.dday([1 end]), to_date(r.config.yearbase,r.dday(1)), ...
        to_date(r.config.yearbase,r.dday(end)));


if fid_cmd
   fclose(fid_nav1);      %don't write nav now
   fclose(fid_nav2);      %don't write nav now
   fclose(fid_cmd);
   fclose(fid_bin);
end

fclose(fid_log);
%====================================================

function ir =  find_iranges( leftover,  p);
gapfrac = .5;  %fraction of delta to allow between leftover(end) and dday(1)



yearbase = p.yearbase;
interval = p.interval;


dday = p.dday;

delta = interval/86400;

if (isempty(leftover) | (  (dday(1) - leftover.dday(end)) > gapfrac)  )
   d0 =dday(1) : delta : dday(end);
else
   d0 = leftover.dday(1) : delta : dday(end);
end

d0(find(d0 < dday(2))) = [];   %throw out times where no dday exist
   
n = length(d0);


ir = zeros(2,n+1);

ir(1,1) = 1; 
for ii = 1:n
   ir(1, ii+1) = max(find(dday <= d0(ii)));;
end

ir(2,:) = [ir(1,2:end)-1, length(dday)];






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


