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 (specify prefix)
%           eg data from waterfall adcp editing:  [data, cfg] = agetmat;
%           eg data from autoedit adcp editing:   [data, cfg] = agetmat('a_');
%        
%
%
% --- 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(file_prefix);      %file_prefix defaults to ''
%
% 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
%
%
% along-track bias parameter is gotten as: twiddle = get_atbiasp(data);
%

%Jules 2000/11/12
% 2001/03/18 added default autoedit prefix of "a_"

% --------------------- agetmat: additional details --------------------
% if only loading, 
%   comment = '';
%   prefix = '';           %load bare getmat output
%   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 [obsolete]
%   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 = [];       % (automatic: can override if this doesn't work)
%                          % where to split longitude (use 0 if near date line)
%   refstd_window = 5;     %running std window length (for jitter calc)
%   navfile = '';          % override navigation from getmat with this ascii 
%                          %   txy fix file (eg. the *.ags file from ubprint)
%                          % NOTE: this can only be accessed by the 
%                          %   [data, config] which are returned from the 
%                          %   function.  If you use agetmat to load the 
%                          %   matfiles from the disk, the navigation will 
%                          %   be the original navigation stored by the user 
%                          %   buffer, not the values loaded from this file.
%
% 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


rungetmat = 1;
prefix = '';  %the default is no prefix

if (nargin == 0)
   rungetmat = 0;
end

if (nargin == 1)
   if (ischar(arg1))
      prefix = arg1;
      rungetmat = 0;
   elseif  (~isstruct(arg1))
      help(mfilename)
      error(' must set single input argument structure')
   
   end
end
  
   
if (nargin > 1)
   if  (~isstruct(arg1))
      help(mfilename)
      error('more than one argument implies input structure and [name, val] pairs\n')
   else
      num_unset = print_fields(arg1);
      if (num_unset > 0)
         help(mfilename)
         error('at least one required input 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.prefix = prefix;
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 = rungetmat; %set above as 0 for prefix only; 1 for struct.
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 = [];       %split longitudes here if specified
%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
cfg.verbose        =        0;           %be verbose
cfg.navfile        =       '';           %don't use one
cfg.interpolate_fixes = 1;               %interpolate fixes

cfg.refsm_cfg.yearbase         = cfg.yearbase;
cfg.refsm_cfg.min_pg           = 50;
cfg.refsm_cfg.min_maxpg        = 70;
cfg.refsm_cfg.bl_half_width    = 3;
cfg.refsm_cfg.xy_bl_half_width = 6;
cfg.refsm_cfg.ref_top          = 1;
cfg.refsm_cfg.ref_mid          = 20;
cfg.refsm_cfg.ref_bot          = 20;
cfg.refsm_cfg.fix_t_mismatch   = ceil((10/300)*300); % seconds between ensemble time and fix time
cfg.refsm_cfg.min_xy_frac      = 0.05;
cfg.refsm_cfg.min_uv_frac      = 0.5;
cfg.refsm_cfg.max_dt           = 300 + cfg.refsm_cfg.fix_t_mismatch;


% for use with refsm

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


% 
if (nargin > 2)%   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, 0);
end


if (isempty(cfg.navfile))
   cfg.navfile = sprintf('../nav/%s.ags', cfg.dbname);
   fprintf('NO NAV FILE SPECIFIED: attempting to use %s\n', cfg.navfile)   
end

cfg.refsm_cfg.yearbase = cfg.yearbase;

%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)
   prefix_option = '';
   if ~isempty(cfg.prefix)
      prefix_option = sprintf(' -f %s ', cfg.prefix);
   end
   gm_opt = '-s';

   cmd = ['getmat -t ', gm_opt, ' ', ...
          prefix_option, cfg.dbpath, cfg.dbname, time_specs];
   cmd = [cmd, ' ', int2str(cfg.nbins)];
   fprintf(cmd)

   if (strcmp(computer, 'PCWIN'))
      dos([cmd])
      paused(1)
   else
      eval(sprintf('!%s',cmd))
   end
else
   fprintf('not running getmet: just loading data with prefix <%s>\n', cfg.prefix)
end

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

try
   load(sprintf('%sother.mat', cfg.prefix))
catch
   fprintf('warning: cannot load filename %sother.mat\n', cfg.prefix)
   fprintf('exiting agetmat with no data\n')
   % lasterr  %%% this could be used
   data = [];
   return
end
   
   
data.dday = DAYS(:)';
if (~isempty(cfg.navfile) & ~exist(cfg.navfile, 'file'))
   warning(sprintf('navfile %s does not exist: not loading\n', cfg.navfile))
   cfg.navfile = '';
   navfile_ok = 0;
else
   navfile_ok = 1;
end

if (length(data.dday) == 1)
   fprintf('exiting agetmat: only one dday\n')
   data = [];
   return
end   


if (rungetmat) %try to get info from nav file
   if (navfile_ok)
      [pp_nav, nn_nav, ee_nav] = fileparts(cfg.navfile);
      if (strcmp(ee_nav, '.mat'))
         fprintf('loading nav data from (adcpsect output) matlab file %s\n', ...
                 cfg.navfile)
         load(cfg.navfile);  
         txy = xyt([3 1 2],:)';
      else
         fprintf('loading nav data from (presumed asciii) file %s\n', cfg.navfile)
         txy = read_asc(cfg.navfile);
      end
      data.lon = table1([txy(:,1) txy(:,2)], data.dday);data.lon = data.lon(:)';
      data.lat = table1([txy(:,1) txy(:,3)], data.dday);data.lat = data.lat(:)';
   else
      %split longitude here 
      if (~isempty(cfg.lon_branch))
         data.lon = ...
             rem(rem( LON(:)' - cfg.lon_branch + 360, 360) + cfg.lon_branch, 360);
      else
         data.lon = unwrap_lon(LON(:)');
      end
      data.lat = LAT(:)';
   end
end 


if (rungetmat & cfg.interpolate_fixes)
   fixes.dd = data.dday;  %from getmat
   fixes.x =  data.lon;   %from getmat
   fixes.y =  data.lat;   %from getmat
   [data.lon, data.lat] = interpolate_fixes(cfg, data.dday, fixes);
end


data.iblkprf = IBLKPRF;
data.heading = HEADING(:)';
data.time = TIME;
load(sprintf('%sdepth.mat', cfg.prefix))

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',      cfg);
data.depth        =      load_file('depth',   cfg);             
data.pflag        =      load_file('pflag',   cfg);
data.amp          =      load_file('amp',     cfg);           
data.sw           =      load_file('sw',      cfg);            
data.umeas        =      load_file('umeas',   cfg);         
data.vmeas        =      load_file('vmeas',   cfg);         
data.w            = 1000*load_file('w',       cfg);             
data.wmean        = 1000*load_file('wmean',   cfg);         
data.e            = 1000*load_file('e',       cfg);             


%%%%%%%%%% impliment adcpnav
if (rungetmat)
   profs.dd = data.dday(:)';
   profs.u =  data.umeas;
   profs.v =  data.vmeas;
   profs.pg = data.pg;
   profs.pf = data.pflag;
   
   if (navfile_ok)
      fixes.dd = txy(:,1)';
      fixes.x =  txy(:,2)';
      fixes.y =  txy(:,3)';
   else %use the data from the nav structure
      fixes.dd = data.dday;  %from getmat
      fixes.x =  data.lon;   %from getmat
      fixes.y =  data.lat;   %from getmat
   end   

   [nav, uv] = adcp_nav(fixes, profs, cfg.refsm_cfg);
else
   load(sprintf('%srefsmnav', cfg.prefix))  %load variables 'nav' and 'uv'
end


%%% plug in the holes with the original values
%%% maybe this should be in adcp_nav?

badi = find(isnan(mstdgap(uv.adj)));
numbins = length(data.amp(:,1));
uv.raw(:,badi) = ...
        (data.umeas(:,badi) + repmat(nav.u(badi), numbins, 1) ) + ...
    i * (data.vmeas(:,badi) + repmat(nav.v(badi), numbins, 1) ) ;
uv.adj(:,badi) = uv.raw(:,badi);


data.uabs_sm = real(uv.adj);
data.vabs_sm = imag(uv.adj);


data.uship = nav.u;
data.vship = nav.v;
data.uship_sm = nav.u_adj;
data.vship_sm = nav.v_adj;


data.lon = nav.x_adj;
data.lat = nav.y_adj;

data.uv = uv;
data.nav = nav;

%jitter calculation
%initialize smoothr-gps deviation
jitter = get_jitter(data.dday, data.uv);
data.ref_std = jitter(1,:); %mean


%%%%%%%%%%  back to getmat
%need cleaned ship speed over ground

%shouldn't need to clean this if adcp_nav is working.
[data.smmps, data.x, data.y] = ll2mps(nav.dd, nav.x_adj, nav.y_adj);

dummy1 = ones(size(data.depth(:,1)));
rads = (90-data.heading)*pi/180;

% on-station 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;

%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);


%%% rescale to cm/s
data.uabs_sm       =100*data.uabs_sm       ;
data.vabs_sm       =100*data.vabs_sm       ;
data.fabs_sm     =100*data.fabs_sm     ;
data.pabs_sm     =100*data.pabs_sm     ;


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

if (rungetmat)
   refsmnavfile = sprintf('%srefsmnav.mat', cfg.prefix);
   fprintf('saving nav info to %s\n', refsmnavfile)
   eval(['save ' refsmnavfile '  nav uv']);
end



%----------------------------------------------------------------
function [jitter, jtitles] = get_jitter(dday, uv)


   nbins = length(uv.raw(:,1));
   dummy= NaN*uv.raw(:,1);
   dt = mstdgap(86400*diff(dday));
   d2uvraw = [dummy  (diff(diff(uv.raw'))/dt)'  dummy];


   [d2uvraw_mean, d2uvraw_std, d2uvraw_n] = mstdgap(d2uvraw);
   d2uvraw_n(d2uvraw_n == 0) = NaN;
   d2uvraw_sem = d2uvraw_std./sqrt(d2uvraw_n);
   
   jitter = 200*[ 50*abs(d2uvraw_mean);              
                  100*abs(d2uvraw_std);
                  500*abs(d2uvraw_sem);
                  30*abs(d2uvraw_mean)];

   jtitles = {'vertical mean (d2h(u))',
                'vertical std (d2h(u))',         
                'vertical sem (d2h(u))',
                'vertical mean/pg (d2h(u))'};
   
   
%------------------------------------------------------------------

function [newx, newy] = interpolate_fixes(cfg, dday, fixes)
%dday is ensemble times
%fixes is fixes.dd, fixes.x, fixes.y from navigation

nt = length(dday);
% For each profile time, find the fix with the closest time.
ifix = zeros(nt, 1);      % indices into fixes.*
abs_min = ifix;           % corresponding magnitude of time mismatch
for ii = 1:nt,
   absdt = abs(fixes.dd - dday(ii));
   i_closest = find(absdt == min(absdt));
   ifix(ii) = i_closest(1);
   abs_min(ii) = absdt(i_closest);
end

% All vectors in fix structure match the columns in profs.
nav.dd = fixes.dd(ifix);
nav.x = fixes.x(ifix);
nav.y = fixes.y(ifix);


% A fix is "missing" if the time mismatch exceeds a threshold.
i_missing_fix = abs_min > (cfg.refsm_cfg.fix_t_mismatch / 86400);
ibad = i_missing_fix;  % Shorter name for temp use.
nav.dd(ibad) = dday(ibad);   % Assign profile time to time of
nav.x(ibad) = NaN;                %   missing nav.
nav.y(ibad) = NaN;

if 1 & (sum(ibad) > 0),  % interpolate missing fixes.    Dangerous!
   nav.x(ibad) = table1([nav.dd(~ibad).', nav.x(~ibad).'], ...
                          nav.dd(ibad)).';
   nav.y(ibad) = table1([nav.dd(~ibad).', nav.y(~ibad).'], ...
                          nav.dd(ibad)).';
   disp(['Interpolated ', num2str(sum(ibad)), ' missing fixes.']);
end

newx = nav.x;
newy = nav.y;

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

% matlab file to load         %variable name in file
fname.pg         = sprintf('%spg',cfg.prefix);   vname.pg          = 'PGOOD';
fname.pflag      = sprintf('%spf',cfg.prefix);   vname.pflag       = 'PFLAG';
fname.amp        = sprintf('%samp',cfg.prefix);  vname.amp         = 'AMP';
fname.sw         = sprintf('%ssw',cfg.prefix);   vname.sw          = 'SPECTRAL_WIDTH';
fname.umeas      = sprintf('%su',cfg.prefix);    vname.umeas       = 'U';
fname.uship_sm   = sprintf('%su',cfg.prefix);    vname.uship_sm    = 'U_SHIP';
fname.vmeas      = sprintf('%sv',cfg.prefix);    vname.vmeas       = 'V';
fname.vship_sm   = sprintf('%sv',cfg.prefix);    vname.vship_sm    = 'V_SHIP';
fname.w          = sprintf('%sw',cfg.prefix);    vname.w           = 'W';
fname.wmean      = sprintf('%sw',cfg.prefix);    vname.wmean       = 'W_MEAN';
fname.e          = sprintf('%se',cfg.prefix);    vname.e           = 'E';
fname.depth      = sprintf('%sdepth',cfg.prefix);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 
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.uship_sm      = 'SHIP U (smoothr)';   
cfg.title.vmeas         = 'MEASURED 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 
cfg.title.uabs_sm           = 'UABS  (smoothr)';
cfg.title.vabs_sm           = 'VABS  (smoothr)';
cfg.title.fabs_sm         = 'FORE ABS (smoothr)';
cfg.title.pabs_sm         = 'PORT ABS (smoothr)';

% NPRFSx1 loaded or calculated variables and their titles
cfg.title.smmps         = 'SHIP SPEED (smoothed)';
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.bins          = 'bins';

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

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

cfg.minclipval.pg    = 0;          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;  


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


%  NBINSxNPRFS loaded variables and their cmapss 
cfg.cmaps.pg            = 'pg3080';
cfg.cmaps.pflag         = 'pg3080';
cfg.cmaps.amp           = 'jet';
cfg.cmaps.sw            = 'jet';
cfg.cmaps.umeas         = 'rg_vel';
cfg.cmaps.vmeas         = 'rg_vel'; 
cfg.cmaps.w             = 'rg_vel'; 
cfg.cmaps.e             = 'rg_vel'; 
cfg.cmaps.uabs_sm       = 'gautoedit'; 
cfg.cmaps.vabs_sm       = 'gautoedit'; 
cfg.cmaps.fabs_sm       = 'gautoedit'; 
cfg.cmaps.pabs_sm       = 'gautoedit'; 
