function [] = process_cast(stn,begin_step,stop,eval_expr)

% function [] = process_cast(stn[,begin_step[,stop[,eval_expr]]])

%

% Process LADCP cast, including GPS, SADCP, and BT data.

%

% Before a cast is processed, [set_cast_params.m] is called --- cruise-

% and/or cast-specific parameters should be set there. After a cast

% is processed, [post_process_cast.m] is called (if it exists) --- cruise-

% and/or cast-specific post-processing should be carried out there.

%

% stop:  0 don't stop (default)

%	-1 stop before begin_step

%	 1 stop after begin_step

%	 2 stop after all steps

%

% eval_expr is evaluated after loading the checkpoint, just before

% 	    processing starts; it can be used to override variables

%	    saved in the checkpoint.

%

% LIST OF STEPS FOR begin_step:

%

%  1: LOAD LADCP DATA

%  2: FIX LADCP-DATA PROBLEMS

%  3: LOAD GPS DATA

%  4: GET BOTTOM-TRACK DATA

%  5: LOAD CTD PROFILE

%  6: LOAD CTD TIME SERIES

%  7: FIND SURFACE & SEA BED

%  8: APPLY PITCH/ROLL CORRECTIONS

%  9: EDIT DATA 

% 10: FORM SUPER ENSEMBLES

% 11: REMOVE SUPER-ENSEMBLE OUTLIERS

% 12: RE-FORM SUPER ENSEMBLES

% 13: (RE-)LOAD SADCP DATA

% 14: CALCULATE INVERSE SOLUTION

% 15: CALCULATE SHEAR SOLUTION

% 16: CALCULATE DIFFUSIVITY PROFILE

% 17: PLOT RESULTS & SHOW WARNINGS

% 18: SAVE OUTPUT

%======================================================================

%                    P R O C E S S _ C A S T . M 

%                    doc: Thu Jun 24 16:54:23 2004

%                    dlm: Wed Jan  7 16:48:20 2009

%                    (c) 2004 A.M. Thurnherr

%                    uE-Info: 609 0 NIL 0 0 72 2 2 8 NIL ofnI

%======================================================================

% NOTES:

%  - changing this function should not be required, except to fix bugs

%  - in order to preserve state variables across load/save calls in matlab,

%    they must be fields of the pcs structure

% HISTORY:

%  Jun 24, 2004: - created by combining laproc.m presolve.m resolve.m

%  Jun 25, 2004: - IMPROVEMENTS:

%			- save plots only if f.res defined

%			- add processing warnings to fig. 11

%  Jun 27, 2004: - removed the saving of processing structures

%		 - added 0 to debug_steps

%  Jun 28, 2004: - replaced debug_steps by stop parameter

%  Jul  2, 2004: - merged with [do_stn.m]

%  Jul  4, 2004: - added update_figures

%  Jul  9, 2004: - BUG: warnings were not shown if only pwarnp was set

%  Jul 16, 2004: - added support for fig. 14 (edit_data.m)

%  Jul 17, 2004: - changed getinv return params

%  Jul 18, 2004: - added eval_expr parameter

%  Jul 19, 2004: - changed [set_cast_params.m] from function to script

%  Jul 21, 2004: - made checkpoints user-selectable

%		 - made sure [set_cast_params.m] is executed after loading

%		   checkpoint

%  Jul 22, 2004: - changed [post_process_cast.m] from function to script

%		   to allow access to variables from [set_cast_params.m]

%  Jul 23, 2004: - moved all state variables to structure pcs.

%  Nov 23, 2004: - added handle wrapper to prevent printing of empty

%                  figures, which bombs on Matlab 7 R 14 (linux)

%  May 23, 2008: - added test after call to getdpthi() to make sure

%		   results are complete (current getdpthi does not work

%		   with CLIVAR P18 leg1 casts 46, 48, 49 & 50

%  Sep 18, 2008: - BUG: p.navdata was non-existent field when no nav

%			data were loaded

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

% STEP 0: EXECUTE ALWAYS

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

if nargin < 4, eval_expr = ''; end	% by default, no variable overrides

if nargin < 3, stop = 0; end		% by default, no debug stops

if nargin < 2, begin_step = 1; end	% by default, start from scratch

pcs.begin_step 	= begin_step;		% parameters are state variables

pcs.stop 	= stop;

pcs.eval_expr 	= eval_expr;

clear f d dr p ps;			% blank slate

f.ladcpdo = ' ';			% required by [m/default.m]

default;				% load default parameters

p = setdefv(p,'checkpoints',[]);	% disable checkpointing by default

if exist('set_cast_params') ~= 2	% set per-cast parameters

  error('cannot find set_cast_params.m');

else

  set_cast_params;

end

openwindows;				% open all windows

if length(f.res) > 1			% open log file

 if exist([f.res,'.log'])==exist('loadrdi.m')

  eval(['delete ',f.res,'.log'])

 end

 diary([f.res,'.log'])

 diary on

end

disp(p.software);			% show version

if length(f.checkpoints) <= 1		% setup checkpointing

 error('Need to set f.checkpoints to write checkpoint files');

end

eval(sprintf('save %s_0 stn',f.checkpoints)); % sentinel

last_checkpoint = pcs.begin_step - 1;	% find last valid checkpoint

while ~exist(sprintf('%s_%d.mat',f.checkpoints,last_checkpoint),'file')

  last_checkpoint = last_checkpoint - 1;

end

pcs.target_begin_step = pcs.begin_step;	% backtrack to last valid checkpoint

if last_checkpoint >= 0 & pcs.begin_step > last_checkpoint+1

  disp(sprintf('Backtracking begin_step from %d to %d',...

			pcs.begin_step,last_checkpoint+1));

  pcs.begin_step = last_checkpoint + 1;

end

while last_checkpoint < 100		% remove now-invalid checkpoints

  last_checkpoint = last_checkpoint + 1;

  if exist(sprintf('%s_%d.mat',f.checkpoints,last_checkpoint),'file')

    delete(sprintf('%s_%d.mat',f.checkpoints,last_checkpoint));

  end

end

tic;					% start timer

%pcs.step_name = 'INITIALIZATION';		% allow debugging of initial params

pcs.cur_step = 0; pcs.update_figures = [];

%end_processing_step;

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

% STEP 1: LOAD LADCP DATA

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'LOAD LADCP DATA'; begin_processing_step;

  % LOAD RDI BB-raw data

  %  this is a rather complex set of functions

  %  1) load raw data

  %     apply magentic deviation if given

  %  2) merge down-up data

  %  3) do some fist order error checks

  [d,p]=loadrdi(f,p); 

  % get instrument serial number

  p=getserial(f,p);

  end_processing_step;

end % OF STEP 1: LOAD DATA

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

% STEP 2: FIX LADCP-DATA PROBLEMS

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'FIX LADCP-DATA PROBLEMS'; begin_processing_step;

  % fix problems with switched beams on instrument

  if existf(p,'beam_switch')==1, [d,p]=switchbeams(d,p); end

  % fix problems with a compass

  if p.fix_compass>0, [d,p]=fixcompass(d,p); end

  end_processing_step;

end % OF STEP 2: FIX LADCP-DATA PROBLEMS

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

% STEP 3: LOAD GPS DATA

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'LOAD GPS DATA'; begin_processing_step;

  p.navdata = 0;

  if length(f.nav)>1 & exist('loadnav')==exist('loadrdi')

    [d,p]=loadnav(f,d,p);

  else

    d.slon=NaN*d.time_jul; d.slat=d.slon;

  end

  end_processing_step;

end % OF STEP 3: LOAD GPS DATA

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

% STEP 4: GET BOTTOM-TRACK DATA

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'GET BOTTOM-TRACK DATA'; begin_processing_step;

  %  Check if hbot values are mostly ==0

  %

  %  some RDI instruments seem to have trouble reporting the distance of the bottom

  %  despite giving reasonable bottom track values.

  %  If this problem is diagnosed we make the distance ourselves.

  %

  ii1=sum(isfinite(d.hbot));

  ii0=sum(d.hbot==0);

  p.hbot_0=ii0/(ii1+1)*100;

  if p.hbot_0>80 

   p.bottomdist=1;

   disp([' WARNING found ',int2str(p.hbot_0),'% of  hbot=0  WARNING'])

  end

  [d,p]=getbtrack(d,p);  

  end_processing_step;

end % OF STEP 4: GET BOTTOM-TRACK DATA

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

% STEP 5: LOAD CTD PROFILE

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'LOAD CTD PROFILE'; begin_processing_step;

  %  get processed ctd profile data

  %  We provide more than one version to support different file formats

  % 

  if length(f.ctdprof)>1 & exist('loadctdprof')==exist('loadrdi')

    [d,p]=loadctdprof(f,d,p);

  end

  end_processing_step;

end % OF STEP 5: LOAD CTD PROFILE

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

% STEP 6: LOAD CTD TIME SERIES

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'LOAD CTD TIME SERIES'; begin_processing_step;

  %  get ctd time series data

  %  We provide more than one version to support different file formats

  % 

  if length(f.ctd)>1 & exist('loadctd')==exist('loadrdi')

    pcs.update_figures = [pcs.update_figures 4];

    [d,p]=loadctd(f,d,p);

  end

  end_processing_step;

end % OF STEP 6: LOAD CTD TIME SERIES

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

% STEP 7: FIND SURFACE & SEA BED (fig. 2)

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'FIND SURFACE & SEA BED'; begin_processing_step;

  pcs.update_figures = [pcs.update_figures 2];

  %  Find depth and bottom and surface using ADCP data 

  if p.getdepth==2

   [d,p]=getdpthi(d,p);

   if length(find(~isfinite(d.izm(1,:))))

     error('Non-finite values in d.izm --- try processing with p.getdepth == 1');

   end

  else

   [d,p]=getdpth(d,p);

  end

  % Plot a summary plot of the raw data

  figure(2), clf

  p=plotraw(d,p);

  pause(.01)

  end_processing_step;

end % OF STEP 7: FIND SURFACE & SEA BED

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

% STEP 8: APPLY PITCH/ROLL CORRECTIONS

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'APPLY PITCH/ROLL CORRECTIONS'; begin_processing_step;

  if length(p.tiltcor)>1

    pd.dpit=p.tiltcor(1);

    pd.drol=p.tiltcor(2);

    d=uvwrot(d,pd,1);

  end

  if length(p.tiltcor)>2

    pu.dpit=p.tiltcor(3);

    pu.drol=p.tiltcor(4);

    d=uvwrot(d,pu,0);

 end

 end_processing_step;

end % OF STEP 8: APPLY PITCH/ROLL CORRECTIONS

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

% STEP 9: EDIT DATA

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

   pcs.step_name = 'EDIT DATA'; begin_processing_step;

   pcs.update_figures = [pcs.update_figures 14];

   d = edit_data(d,p);

  end_processing_step;

end % OF STEP 9: EDIT DATA

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

% STEP 10: FORM SUPER ENSEMBLES

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

   pcs.step_name = 'FORM SUPER ENSEMBLES'; begin_processing_step;

   pcs.update_figures = [pcs.update_figures 5 6 10];

   [di,p,d]=prepinv(d,p);

  end_processing_step;

end % OF STEP 10: FORM SUPER ENSEMBLES

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

% STEP 11: REMOVE SUPER-ENSEMBLE OUTLIERS

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'REMOVE SUPER-ENSEMBLE OUTLIERS'; begin_processing_step;

  % Reduce scatter by successively removing 1% of the data

  %  in oder to do that we need a first solution

  %

  if ps.outlier>0 | p.offsetup2down>0

     diary off

     if exist('loadsadcp')==exist('loadrdi') 

      [di,p]=loadsadcp(f,di,p);

     end

     dino=di;

     lanarrow

     diary on

  end

  end_processing_step;

end % OF STEP 11: REMOVE SUPER-ENSEMBLE OUTLIERS

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

% STEP 12: RE-FORM SUPER ENSEMBLES

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'RE-FORM SUPER ENSEMBLES'; begin_processing_step;

  %

  % once we have a first guess profile we recompute the super ensemble

  %

  if (p.offsetup2down>0 & length(d.izu)>0)

   pcs.update_figures = [pcs.update_figures 5 6 10];

   diary off

   [di,p,d]=prepinv(d,p,dr);

   diary on

  end

  end_processing_step;

end % OF STEP 12: RE-FORM SUPER ENSEMBLES

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

% STEP 13: (RE-)LOAD SADCP DATA

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = '(RE-)LOAD SADCP DATA'; begin_processing_step;

  if exist('loadsadcp')==exist('loadrdi') 

    pcs.update_figures = [pcs.update_figures 9];

    di=loadsadcp(f,di,p);

  end

  end_processing_step;

end % OF STEP 13: (RE-)LOAD SADCP DATA

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

% STEP 14: CALCULATE INVERSE SOLUTION

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'CALCULATE INVERSE SOLUTION'; begin_processing_step;

  pcs.update_figures = [pcs.update_figures 3 7 12 13 14];

  % 

  %  take advantage of presolve if it existed

  %  call the main inversion routine

  %

  [p,dr,ps,de]=getinv(di,p,ps,dr,1);

  %

  % check inversion constraints

  % 

  p=checkinv(dr,di,de,der,p,ps);

  if existf(de,'bvel'), p=checkbtrk(d,di,de,dr,p); end

  end_processing_step;

end % OF STEP 14: CALCULATE INVERSE SOLUTION

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

% STEP 15: CALCULATE SHEAR SOLUTION

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'CALCULATE SHEAR SOLUTION'; begin_processing_step;

  % Compute 'old fashioned' shear based solution 

  %  two choices, fisrt us all data

  %  second use super ensemble data

  ps=setdefv(ps,'shear',1);

  if ps.shear>0

   if ps.shear==1

    [ds,dr,ps,p]=getshear2(d,p,ps,dr);

   else

    [ds,dr,ps,p]=getshear2(di,p,ps,dr);

   end

  end

  end_processing_step;

end % OF STEP 15: CALCULATE SHEAR SOLUTION

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

% STEP 16: CALCULATE DIFFUSIVITY PROFILE

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'CALCULATE DIFFUSIVITY PROFILE'; begin_processing_step;

  % 

  %  prepare ctd data for output

  %

  if existf(d,'ctdprof_p')

   % first save CTD data with profile

   dr.ctd_t=interp1q(d.ctdprof_z,d.ctdprof_t,dr.z);

   dr.ctd_s=interp1q(d.ctdprof_z,d.ctdprof_s,dr.z);

  end

  if existf(d,'ctdprof_ss')

   dr.ctd_ss=interp1q(d.ctdprof_z,d.ctdprof_ss,dr.z);

  end

  %

  % estimate vertical diffusivity if CTD N2 profile exist

  %  this is totally experimental and might be WRONG!

  %

  if existf(d,'ctdprof_N2') & exist('pk','var')

   pcs.update_figures = [pcs.update_figures 8];   

   disp(' compute vertical diffusivity')

   dr.ctd_N2=interp1q(d.ctdprof_z,d.ctdprof_N2,dr.z);

   [dk,pk]=getkzprof(dr,pk,dr.ctd_N2);

   if existf(dk,'z')

    if length(dk.z)>1

     zz=[dk.z(1)-dk.zrange/4,dk.z,dk.z(end)+dk.zrange/4]';

     Kz=dk.Kz([1 1:end end])';

     dr.Kz=interp1q(zz,Kz,dr.z);

    end

   end

  end

  end_processing_step;

end % OF STEP 16: CALCULATE DIFFUSIVITY PROFILE

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

% STEP 17: PLOT RESULTS & SHOW WARNINGS

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

pcs.cur_step = pcs.cur_step + 1;

if pcs.begin_step <= pcs.cur_step

  pcs.step_name = 'PLOT RESULTS & SHOW WARNINGS'; begin_processing_step;

  pcs.update_figures = [pcs.update_figures 1 11];

  % Plot final results

  figure(1), clf

  plotinv(dr,d,p,ps)

  % Convert p.warn to one line of text with newline characters

  p.warnings = [];

  for i = 1:size(p.warnp,1)

     p.warnings = [p.warnings deblank(p.warnp(i,:)) char(10)];

  end

  figure(11)

  clf

  % experimental diagnostic of battery voltage

  %

  p=battery(p);

  %

  % complete task by repeating the most important warnings

  %

  if size(p.warn,1) + size(p.warnp,1) > 2