function [data, cfg] =  agetmat(arg1, varargin)
% [data, cfg] = agetmat(config, varargin)
% two modes exist: 
% (1) set up configuration to call getmat and do so
%     usage:  [data, cfg] = agetmat(config, name1, val1, name2, val2)
% (2) just read the output of a previous getmat call (eg during
%                 shipboard adcp editing). 
%     usage: [data, cfg] = agetmat;
%
%
% --- mode: run getmat ----
% config has lots of stuff in it.
% varargin is name,value pairs which can override the values in config
% data is getmat data; cfg is titles, clips, fields
%
% required fields in config are (example of syntax given)
% end_dd         =       224;           %beginning of last decimal day
% ddrange        =       1              %number of ddays to get data
% yearbase       =       2000;          % year base
% comment        =       'ka0006';      % title
% dbpath         =       '../adcpdb/';  %with trailing slash!
%                                       % path to database
% dbname         =       'a0006';       % database name
%
%  usage: 
%      [data, cfg] = agetmat(config, name1, val1, name2, val2,...)
%
% --- mode: just load it ---
% suitable for loading getmat data which were generated by shipboard editing
% data comes out the same; cfg gets the following:
%
% usage:
%     [data, config] = agetmat;
%
% NOTE: 
% units of diagnostics are all directly from getmat
% abs water speeds are cm/s, ship  and relative water speeds are m/s
% e and w are X1000
%

%Jules 2000/11/12


% --------------------- agetmat: additional details --------------------
% if only loading, 
%   comment = '';
%   rungetmat = 0;         %generally want to run getmat
%   os_cutoff = 2;         %on-station is under 2m/s
%   rl_range = [5 12];     %reference layer range  
%   start_dd               %from DAYS (loaded)
%   ddrange                %from DAYS (loaded)
%   yearbase = -9999; comment = ''; dbname = ''; dbpath = '';
%
%
%
% optional fields in config default to:
%
%   nbins = 128;           %number of bins to extract using getmat
%   cleanerparams = [200,1]; %start  with [2,3] for more cleaner
%   comment = '';          % no additional comment
%   rungetmat = 1;         %generally want to run getmat
%   os_cutoff = 2;         %on-station is under 2m/s (to show fore/port )
%   rl_range = [5 12];     %reference layer range (should match smoothr)
%   lon_branch = 180;      % where to split longitude (use 0 if near date line)
%   refstd_window = 5;     %running std window length (for jitter calc)
%
%
% substructures containing plotting information defaults, accessed by
% field name in structure "data"
%
%   title                  %structure of titles
%   minclipval             %structure of minimum value to show
%   maxclipval             %structure of maximum value to show

% Jules 2000/11/27



program_name = 'agetmat';

if ((nargin == 0) & (missing_getmat_output))
help(program_name)
  error(...
     sprintf('-->>>> could not find sufficient getmat output in\n-->>>> %s\n', pwd))
end

if (nargin >= 1)
  if  (~isstruct(arg1))
    help(program_name)
    error(' must set single input argument structure')
  else
    num_unset = print_fields(arg1);
    if (num_unset > 0)
      help(program_name)
      error('at least one field was not set in input argument')
    end
  end
end

    
%initialize cfg
cfg = apdefs; %function is at the bottom: sets defaults

% use input structure to override defaults ------------------------------
%defaults for non-required fields of argument structure
cfg.nbins = 128;  
cfg.cleanerparams = [200,1]; % no cleaning.  would usually would start 
                             %                with [2,3] for more cleaner
cfg.refstd_window = 5;    %for jitter calculation - window length
cfg.comment = '';
cfg.rungetmat = 1;         %generally want to run getmat
cfg.os_cutoff = 2;         %on-station is under 2m/s (to show fore/port )
cfg.rl_range = [5 12];     %reference layer range
cfg.lon_branch = 180;      %split longitudes here
%these should get overwritten by input call
cfg.start_dd       =       -9999;        %bad value 
cfg.ddrange        =       -9999;        %bad value
cfg.yearbase      =       -9999;        %bad value
cfg.dbpath         =       '';           %bad value
cfg.dbname         =       '';           %bad value
cfg.comment        =       '';           %bad value


if (nargin == 0)
  cfg.rungetmat      =        0;           %don't run it
end


%update (replace or apppend) as necessary, fields from arg1 to cgf
if (nargin >= 1)
  cfg = fillstruct(cfg, arg1); 
end

% 
if (nargin > 1)%   we'll assume the rest is a set of name,value pairs
   %override input structure with name,value pairs
   fprintf('overriding some configuration settings with command line values\n')
   cfg = fillstruct(cfg, varargin )
end


%set up and run getmat ------------------------------
% extract data with getmat; clean up workspace
if (cfg.rungetmat)
   time_specs = [' ',num2str(cfg.yearbase),' ',num2str(cfg.start_dd),...
                 ' ',num2str(cfg.ddrange)];
   
   %set up getmat; (get everything)
   gm_opt = '-s';
   %cmd = ['!getmat -t ', gm_opt, ' ', cfg.dbpath, cfg.dbname, time_specs];
   cmd = ['!getmat -t ', gm_opt, ' ', cfg.dbname, time_specs];
   cmd = [cmd, ' ', int2str(cfg.nbins)]
   eval(cmd)
else
   fprintf('not running getmet: just loading data\n')
end

% load everything now ------------------------------
%load navigation information and trim bins as necessary before loading the rest

try
   load other
catch
   fprintf('warning: cannot load filename other.mat\n')
   fprintf('exiting agetmat with no data\n')
   % lasterr  %%% this could be used
   data = [];
   return
end
   
   
data.dday = DAYS(:)';
  %split lonitude here 
data.lon = rem( LON(:)' - cfg.lon_branch + 360, 360) + cfg.lon_branch;
data.lat = LAT(:)';
data.iblkprf = IBLKPRF;
data.heading = HEADING(:)';
data.time = TIME;
load depth;

if (cfg.rungetmat == 0)
  cfg.nbins = NBINS;
  cfg.start_dd = min(data.dday);
  cfg.ddrange = max(data.dday) - min(data.dday);
end
 

disp(sprintf('Bins requested: %d, available: %d.\n', cfg.nbins, NBINS))
cfg.nbins=min(cfg.nbins, NBINS); %nbins was set earlier, NBINS is from the data


bins=[1:cfg.nbins]';
clear DAYS      IBLKPRF   LGB       MAB       NPRFS     
clear HEADING   LAT       LON       NBINS     TIME     



% get the loadable ones --------------------

data.pg           =      load_file('pg', bins);
data.depth        =      load_file('depth',bins);             
data.pflag        =      load_file('pflag', bins);
data.amp          =      load_file('amp',bins);           
data.sw           =      load_file('sw',bins);            
data.umeas        =      load_file('umeas',bins);         
data.uship_sm     =      load_file('uship_sm',bins);      
data.vmeas        =      load_file('vmeas',bins);         
data.vship_sm     =      load_file('vship_sm',bins);      
data.w            = 1000*load_file('w',bins);             
data.wmean        = 1000*load_file('wmean',bins);         
data.e            = 1000*load_file('e',bins);             


badi = find(isnan(data.lon + data.lat));  %JUNK it if it's bad navigation
data.dday(badi)=[];            % 1 x NPRFS
data.lon(badi)=[];    
data.lat(badi)=[];    
data.heading(badi)=[]; 
data.uship_sm(badi)=[];
data.vship_sm(badi)=[];
data.wmean(badi)=[];
%
data.iblkprf(badi,:)=[];         %NPRFS x 2
%
data.pg(:,badi)=[];             %NBINS x NPRFS
data.depth(:,badi)=[]; 
data.pflag(:,badi)=[]; 
data.amp(:,badi)=[]; 
data.sw(:,badi)=[]; 
data.umeas(:,badi)=[]; 
data.vmeas(:,badi)=[]; 
data.w(:,badi)=[]; 
data.e(:,badi)=[]; 


%ship speeds: --------------------
%from smoothr navigation
data.smmps        = abs(data.uship_sm + i*data.vship_sm);
% from gps navigation
[mps, x, y] = ll2mps(data.dday, data.lon, data.lat);
data.clmps       = cleaner(mps, cfg.cleanerparams(1), cfg.cleanerparams(2) );
data.x = x;
data.y = y;

%get ship's speed in u and v (on dday times)
data.uship_g   = [0 diff(x)./diff(data.dday * 86400)];
data.vship_g   = [0 diff(y)./diff(data.dday * 86400)];

% referencing velocities --------------------
dummy1 = ones(size(data.depth(:,1)));
rads = (90-data.heading)*pi/180;
reflbins = cfg.rl_range(1):cfg.rl_range(2);
% in East (u) and North (v)  directions ...
% from smoothr nav
data.uabs_sm        = data.umeas + dummy1*data.uship_sm;
data.vabs_sm        = data.vmeas + dummy1*data.vship_sm;       
data.uabs_refl_sm   = mean(data.uabs_sm(reflbins,:));
data.vabs_refl_sm   = mean(data.vabs_sm(reflbins,:));
% from gps (i.e. other) navigation
data.uabs_g         = data.umeas + dummy1*data.uship_g;
data.vabs_g         = data.vmeas + dummy1*data.vship_g;
data.uabs_refl_g    = mean(data.uabs_g(reflbins,:));
data.vabs_refl_g    = mean(data.vabs_g(reflbins,:));

% and fore, port directions ...
% from smoothr nav

% mask to dot-multipy to turn fp into NaN when on station
osmask_sm = ones(size(data.dday));    %go sideways
osmask_sm(find(abs(data.uship_sm + i*data.vship_sm) <= cfg.os_cutoff)) =NaN;
osmask_g = ones(size(data.dday));    %go sideways
osmask_g(find(abs(data.uship_g + i*data.vship_g) <= cfg.os_cutoff)) =NaN;


%foreport of measured vel
fp = (data.umeas + i*data.vmeas).*exp(-i*dummy1*rads);
data.fmeas = (dummy1*osmask_sm).*real(fp);
data.pmeas = (dummy1*osmask_sm).*imag(fp);

% abs: from smoothr nav
fp =    (data.uabs_sm + i*data.vabs_sm).*exp(-i*dummy1*rads);
data.fabs_sm      = (dummy1*osmask_sm).*real(fp);
data.pabs_sm      = (dummy1*osmask_sm).*imag(fp);

fp =    (data.uabs_refl_sm + i*data.vabs_refl_sm).*exp(-i*rads);
data.fabs_refl_sm = osmask_sm.*real(fp);
data.pabs_refl_sm = osmask_sm.*imag(fp);

% abs: from gps (i.e. other) navigation
fp =      (data.uabs_g + i*data.vabs_g).*exp(-i*dummy1*rads);
data.fabs_g       = (dummy1*osmask_g).*real(fp);
data.pabs_g       = (dummy1*osmask_g).*imag(fp);

fp =     (data.uabs_refl_g + i*data.vabs_refl_g).*exp(-i*rads);
data.fabs_refl_g  = osmask_g.*real(fp);
data.pabs_refl_g  = osmask_g.*imag(fp);


% (out = back to meters) %%% rescale to cm/s
% data.uabs_sm       =100*data.uabs_sm       ;
% data.uabs_g        =100*data.uabs_g        ;
% data.uabs_refl_sm  =100*data.uabs_refl_sm  ;
% data.uabs_refl_g   =100*data.uabs_refl_g   ;
% 
% data.vabs_sm       =100*data.vabs_sm       ;
% data.vabs_g        =100*data.vabs_g        ;
% data.vabs_refl_sm  =100*data.vabs_refl_sm  ;
% data.vabs_refl_g   =100*data.vabs_refl_g   ;
% 
% data.fabs_sm     =100*data.fabs_sm     ;
% data.fabs_g      =100*data.fabs_g      ;
% data.fabs_refl_sm=100*data.fabs_refl_sm;
% data.fabs_refl_g =100*data.fabs_refl_g ;
% 
% data.pabs_sm     =100*data.pabs_sm     ;
% data.pabs_g      =100*data.pabs_g      ;
% data.pabs_refl_sm=100*data.pabs_refl_sm;
% data.pabs_refl_g =100*data.pabs_refl_g ;

data.bins = [1:length(data.depth(:,1))];



%jitter calculation
ww = cfg.refstd_window;  %from aflagit_setup.m

uvdiff = abs( (data.uabs_refl_g - data.uabs_refl_sm)  + ...
              i*( data.vabs_refl_g - data.vabs_refl_sm));


%initialize:
ref_std = zeros(1,length(data.dday) - ww + 1 );
%floating std
for ii=1:length(ref_std)
   ref_std(ii) = std(uvdiff(ii+[0:ww-1]));
end
data.ref_std = [zeros(1,floor(ww/2))  ref_std  zeros(1,floor(ww/2))];
clear ww ddu ddv uvdiff;


%-----------------------------------------------------------------------
function x = load_file(var, bins)

% matlab file to load         %variable name in file
fname.pg         =   'pg';    vname.pg          = 'PGOOD';
fname.pflag      =   'pf';    vname.pflag       = 'PFLAG';
fname.amp        =  'amp';    vname.amp         = 'AMP';
fname.sw         =   'sw';    vname.sw          = 'SPECTRAL_WIDTH';
fname.umeas      =    'u';    vname.umeas       = 'U';
fname.umeas_refl =    'u';    vname.umeas_refl  = 'U_MEAN';
fname.uship_sm    =   'u';    vname.uship_sm    = 'U_SHIP';
fname.vmeas      =    'v';    vname.vmeas       = 'V';
fname.vmeas_refl =    'v';    vname.vmeas_refl  = 'V_MEAN';
fname.vship_sm   =    'v';    vname.vship_sm    = 'V_SHIP';
fname.w          =    'w';    vname.w           = 'W';
fname.wmean      =    'w';    vname.wmean       = 'W_MEAN';
fname.e          =    'e';    vname.e           = 'E';
fname.depth      ='depth';    vname.depth       = 'DEPTH';

% We can add more later...
if (~isfield(fname, var))
   error(sprintf(...
      'cannot find field %s for file (check inside agetmat)\n', var));
end
filename = sprintf('%s.mat', getfield(fname, var));

if (~isfield(vname, var))
   error(sprintf(...
      'cannot find field %s for variable (check inside agetmat)\n', var))
end
varname= getfield(vname, var);

x = []; % might help the parser by initializing x
if (exist(filename))
   load(filename)
   eval(['x = ', varname, ';']);
else
   fprintf('cannot find filename %s: returning [] for %s\n',...
           filename, varname)
end

   

%----------------------------------------------------------------------
function num_missing_files = missing_getmat_output;
% return number of missing fields

reqfilenames = {'amp', 'e', 'pf', 'v', 'depth', 'other', 'pg', 'u', 'w'};

num_missing_files = 0;
for fni = 1:length(reqfilenames)
  if ~exist(sprintf('./%s.mat', reqfilenames{fni}), 'file')
     num_missing_files = num_missing_files + 1;
     fprintf('missing getmat output: %s.mat\n', reqfilenames{fni})
  end
end



%------------------------------------------------------------
function num_unset = print_fields(arg1)
% run through required fields and check to see if all are initialized

required_fields = ...
{'start_dd',
'ddrange',
'yearbase',
'comment',
'dbpath',
'dbname'}       ;
       
num_unset = 0;
for fnamei = 1:length(required_fields)
   if ~isfield(arg1, required_fields{fnamei})
      warning(sprintf('field <%s> not set\n', required_fields{fnamei}))
      num_unset = num_unset + 1;
   end
end

%----------------------------------------------------------------
function cfg = apdefs();

%%NOTE: this first section (titles) has to contail all the variables in 'data'
%%% structure


%  NBINSxNPRFS loaded variables and their titles  (could be fewer bins)
cfg.title.pg            = 'PERCENT GOOD';
cfg.title.depth         = 'DEPTH';
cfg.title.pflag         = 'PROFILE_FLAGS';
cfg.title.amp           = 'AMPLITUDE';
cfg.title.sw            = 'SPECTRAL WIDTH';
cfg.title.umeas         = 'MEASURED U';         
cfg.title.umeas_refl    = 'MEASURED REFL (U)';  
cfg.title.uship_sm      = 'SHIP U (smoothr)';   
cfg.title.vmeas         = 'MEASURED V';         
cfg.title.vmeas_refl    = 'MEASURED REFL (V)';  
cfg.title.vship_sm      = 'SHIP V (smoothr)';   
cfg.title.w             = 'W';
cfg.title.wmean         = 'REFL W';
cfg.title.e             = 'ERROR VELOCITY';
% NBINSxNPRFS calculated variables and their titles  (could be fewer bins)
cfg.title.uabs_sm           = 'UABS  (smoothr)';
cfg.title.uabs_g            = 'UABS  (gps nav)';
cfg.title.uabs_refl_sm      = 'UABS  (REFL, smoothr)';
cfg.title.uabs_refl_g       = 'UABS  (REFL, gps nav)';

cfg.title.vabs_sm           = 'VABS  (smoothr)';
cfg.title.vabs_g            = 'VABS  (gps nav)';
cfg.title.vabs_refl_sm      = 'VABS  (REFL, smoothr)';
cfg.title.vabs_refl_g       = 'VABS  (REFL, gps nav)';

cfg.title.fabs_sm         = 'FORE ABS (smoothr)';
cfg.title.fabs_g          = 'FORE ABS (gps nav)';
cfg.title.fabs_refl_sm    = 'FORE ABS (REFL, smoothr)';
cfg.title.fabs_refl_g     = 'FORE ABS (REFL, gps nav)';

cfg.title.pabs_sm         = 'PORT ABS (smoothr)';
cfg.title.pabs_g          = 'PORT ABS (gps nav)';
cfg.title.pabs_refl_sm    = 'PORT ABS (REFL, smoothr)';
cfg.title.pabs_refl_g     = 'PORT ABS (REFL, gps nav)';

% NPRFSx1 loaded or calculated variables and their titles
cfg.title.smmps         = 'SHIP SPEED (smoothr)';
cfg.title.uship_g       = 'SHIP SPEED (U: gps)';
cfg.title.vship_g       = 'SHIP SPEED (U: gps)';
cfg.title.lon           = 'LONGITUDE';
cfg.title.iblkprf       = 'BLOCK AND PROFILE';
cfg.title.lat           = 'LATITUDE';
cfg.title.dday          = 'DECIMAL DAY';
cfg.title.heading       = 'SHIP HEADING';
cfg.title.clmps         = 'SHIP SPEED (gps)';

cfg.title.bins          = 'bins';

%============================================================
% the rest of these defaults are just set to use if they exist

%----------------------------------------------------------------------
% default minimum and maximum clip values

cfg.minclipval.pg    = 30;         cfg.maxclipval.pg    = 100; 
cfg.minclipval.amp   = 10;         cfg.maxclipval.amp   = 230; 
cfg.minclipval.sw    = 30;         cfg.maxclipval.sw    = 260; 
cfg.minclipval.e     = -70;        cfg.maxclipval.e     = 70;
cfg.minclipval.uv    = -1;         cfg.maxclipval.uv    = 1; 
cfg.minclipval.mps   = 0;          cfg.maxclipval.mps   = 7;  

