function init()
%
% Usage:  init
%
% A script M-file to load U, V, W, E, AMP, PGOOD, and other data from .MAT files,
% generate the array YAXIS for the y-axis (bins),
% generate back-up copies of the original U,V,W,E arrays for use by RESTORE.M,
% and reset offset to 0.
%

%editglob
% Script: should not need global declaration; all should have been
% done in setup.m

global ABS_U ABS_V AMP AMP_OFFSET BADACC BADBIN BADBIN_BAK
global BADPRF BADPRF_BAK BMAX BMIN BOTM_EDIT BOTTOM BOTTOM_BAK
global CURRX DAYS DBNAME DEFAULT_NBINS DEFAULT_NPRFS DEFAULT_SEQ_FLAG
global DEPTH DISP E EDIT_MODE E_OFFSET IBLKPRF IKEY INITed 
global LASTCOMM LASTPLOT LAT LISTED LON LONLAT NBINS NPRFS
global PFLAG PGOOD PLOTBAD PMAX PMIN REFLAYER_ENDBIN REFLAYER_STARTBIN
global TIME U UV_OFFSET U_ROTATE V V_ROTATE W W_OFFSET YAXIS
global PLOTBAD_ORIGINAL
global PGOOD_THRESHOLD D2W_THRESHOLD D2UV_THRESHOLD ...
       WVAR_THRESHOLD AMP_THRESHOLD EV_THRESHOLD
global USE_LGB PLOTNAV start_bin verbose PLOTUV
global PGOOD_OFFSET LGB MAB U_SHIP V_SHIP U_MEAN V_MEAN HEADING
global GOOD_PG 


if EDIT_MODE & ~LISTED
  if ~isempty([find(BADBIN & ~mod(BADBIN,2))' ...
               find(abs(BOTTOM) < NBINS)  ...
               find(BADPRF & ~mod(BADPRF,2)) ])
    reply = 'x';
    while reply ~= 'y' & reply ~= 'Y' & reply ~= 'n' & reply ~= 'N'
      reply = input('Flagged profiles/bin not yet listed.  List?  <y|n> ', 's');
      if reply == 'y' | reply == 'Y'
        list      % in case user forgets to enter list command for previous set
      end
    end
  end
end;

fprintf('Loading matfiles...\n')
load u.mat
load v.mat
load w.mat
load e.mat
load pg.mat
load pf.mat
load amp.mat
load other.mat
load depth.mat

%% create other auxiliary arrays, as needed:
if isempty(PFLAG)
  PFLAG = zeros(NBINS,NPRFS);
end

% plotting arrays:
offset(0);
YAXIS = [-1:-1:-NBINS];
IKEY = [];
PMIN = 1;
PMAX = NPRFS;
BMIN = 1;
BMAX = NBINS;
CURRX = [];
S = 1.0;       % m/s cutoff for rotation: was 2.5 m/s (jules 97/02/20)

% flagging arrays:
IBAD_EV = [];
IBAD_D2 = [];
IBAD_WVAR = [];
IBAD_AMP = [];

% editing arrays:
BADACC = ones(NBINS,NPRFS);
GOOD_PG = BADACC;
   %% apply plot flagging criteria here:
   if ~isempty(PGOOD_THRESHOLD)
      GOOD_PG = PGOOD >= PGOOD_THRESHOLD;
      BADACC(find(~GOOD_PG)) = NaN;
      GOOD_PG = BADACC;
   end;

   %% Separate PFLAG into the glitch part and the remainder;
   %% use the remainder in BADACC, and leave only the glitch
   %% part in PFLAG.
   glitch = rem(PFLAG, 2);
   ibad = find( round(PFLAG - glitch) ~= 0);
   BADACC(ibad) = NaN;
   PFLAG = glitch;
   clear glitch
   NaN_pflag = ones(size(PFLAG));
   NaN_pflag(find(PFLAG)) = NaN;


BADPRF = zeros(1,NPRFS);
BADBIN = PFLAG;
BOTTOM = (NBINS+1) * ones(1,NPRFS);
BADPRF_BAK = BADPRF;
BADBIN_BAK = BADBIN;
BOTTOM_BAK = BOTTOM;

% saved commands:
LASTCOMM = '';


if USE_LGB
  ineg = find(LGB < 0);
  if ~isempty(ineg)
    BADPRF(ineg) = 1;
    LGB(ineg) = 0;
  end
  ibad = find(LGB < NBINS);
  for ii = 1:length(ibad)
    BADACC(1+LGB(ibad(ii)):NBINS,ibad(ii)) = NaN;
  end
  clear ineg ibad ii
end

if ~isempty(AMP_THRESHOLD)
  IBAD_AMP = flagamp(AMP .* BADACC, AMP_THRESHOLD);
  if verbose & ~isempty(IBAD_AMP)
    LISTED = 0;
    ibinprf = ndx2subs(IBAD_AMP, NBINS);
    BOTTOM(ibinprf(:,2)) = ibinprf(:,1);
    i_already = find(BOTTOM == MAB');
    BOTTOM(i_already) = BOTTOM_BAK(i_already);
    Blk_Prf_Bin_Z_Gap = [IBLKPRF(ibinprf(:,2),1)'; IBLKPRF(ibinprf(:,2),2)'; ibinprf(:,1)'; ...
      DEPTH(IBAD_AMP)'; diff([0 ibinprf(:,2)'])-1]
  end
end

NaN_vel = BADACC .* NaN_pflag;

if ~isempty(EV_THRESHOLD)
  IBAD_EV = find(abs(E .* NaN_vel) > EV_THRESHOLD/1000);
  BADBIN(IBAD_EV) = BADBIN(IBAD_EV) + 4;  % 4 = BAD EV FLAG
  if verbose & ~isempty(IBAD_EV)
    LISTED = 0;
    ibinprf = ndx2subs(IBAD_EV, NBINS);
    Blk_Prf_Bin_EV = [IBLKPRF(ibinprf(:,2),1)'; IBLKPRF(ibinprf(:,2),2)'; ...
      ibinprf(:,1)'; round(E(IBAD_EV)*1000)']
  end
  clear ibinprf Blk_Prf_Bin_EV
end

if ~isempty(D2UV_THRESHOLD) & ~isempty(D2W_THRESHOLD)
  d2u = 1000 * diff(diff(U .* NaN_vel));
  d2v = 1000 * diff(diff(V .* NaN_vel));
  d2w = 1000 * diff(diff(W .* NaN_vel));
  imask = (~isnan(d2u) & abs(d2w) > D2W_THRESHOLD & ...
            (abs(d2u) > D2UV_THRESHOLD | abs(d2v) > D2UV_THRESHOLD));
  ibad_d2  = find(imask > 0);
  if ~isempty(ibad_d2)
    LISTED = 0;
    imask = [zeros(2,NPRFS); imask];   % need to shift subscripts back to 0 diff
    ibinprf = ndx2subs(find(imask > 0), NBINS);
    ibinprf(:,1) = ibinprf(:,1) - 1;
    IBAD_D2 = ibinprf(:,1) + (ibinprf(:,2)-1) * NBINS;  % vectorize the bin-profile subscripts
    BADBIN(IBAD_D2) = BADBIN(IBAD_D2) + 2;  % 2 = BAD D2 FLAG
    if verbose & ~isempty(IBAD_D2)
      Blk_Prf_Bin_D2U_V_W = round([IBLKPRF(ibinprf(:,2),1)'; IBLKPRF(ibinprf(:,2),2)'; ibinprf(:,1)';...
        d2u(ibad_d2)'; d2v(ibad_d2)'; d2w(ibad_d2)'])
    end
  end
  clear d2u d2v d2w Blk_Prf_Bin_D2U_V_W ibinprf ibad_d2 imask
end

%keyboard

if ~isempty(WVAR_THRESHOLD)
  imask = ~isnan(W .* NaN_vel);
  if start_bin > 1  % treat the top bins as NaN also
    imask(1:start_bin-1,:) = logical(0);
  end
  w = imask+0.0;   % JH convert to double so matlab isn't confused about types
  w(imask) = W(imask) * 1000.0;  % m/s to mm/s
  n = sum(imask);
  wvar = zeros(size(n));
  ii = find(n > 1);

  if ~isempty(ii),
     wvar(ii) = (sum(w(:,ii) .^ 2) - (sum(w(:,ii)) .^ 2)./n(ii)) ./ (n(ii) - 1);
  end
  IBAD_WVAR = find(wvar > WVAR_THRESHOLD);
  BADPRF(IBAD_WVAR) = 2 + BADPRF(IBAD_WVAR);          % 2 = high W Variance flag
  if verbose & ~isempty(IBAD_WVAR)
    LISTED = 0;
    Blk_Prf_WVar = [IBLKPRF(IBAD_WVAR,1)'; IBLKPRF(IBAD_WVAR,2)'; round(wvar(IBAD_WVAR))]
  end
  clear w imask w n wvar Blk_Prf_WVar
end

% Now subtract the reference layer means:
if REFLAYER_STARTBIN > 0

  U = U - U_MEAN(ones(NBINS, 1), :);
  V = V - V_MEAN(ones(NBINS, 1), :);
  W = W - W_MEAN(ones(NBINS, 1), :);
end

% Calculate velocity profiles rotated into foreward and port components.
if REFLAYER_STARTBIN > 0
   uvc = U_MEAN + i * V_MEAN;
else
   uvc = ones(1, NPRFS) * S/2;
end
heading_radians = angle(-uvc); % approx heading based on vel rel to water
speed = abs(uvc);
fast_enough = speed >= S;
xx = exp(-i*heading_radians) .* fast_enough;
uvcprof_rotated = (U + i * V) .* xx(ones(NBINS,1),:);
% uvcprof_rotated should now be the u + iv in coordinates with x
% in the direction of the ship's heading, unless the speed is less
% than the threshold S, in which case u + iv is zero.
U_ROTATE = real(uvcprof_rotated);
V_ROTATE = imag(uvcprof_rotated);


editnav
INITed = 1;
% plot navigation & U & V if available and requested
if PLOTNAV & ~isempty(LON)
  figure(1);
  clf
  plotnav
end
if PLOTUV
  figure(2);
  clf
  savelast = LASTPLOT;
  uv
  LASTPLOT = savelast;
end

figure(3); % all user-requested plots go here
clf

