function varargout=despike2(cruiseid,staa,stab)
% DESPIKE2 -  remove spikes in the secondary sensors from CTD Data    
%
% CTD Calibration toolbox.
%
% USAGE
% 	cnv_ds2=despike2(cruiseid,staa,stab);
%
% INPUT
%	cruiseid	cruise identification name
%   
% OUTPUT
%	cnv_ds2		despiked secondary sensors
%          
% DESCRIPTION: This function uses the cnv files from the cruise database
% and remove the spikes related to the secondary sensors
%

% AUTHORS:Carlos Fonseca 
%         UM/CIMAS
%         Fri Dec 2 09:19:37 EDT 2005
% 
% VERSION: 
%	$Revision: 1.13 $
%	$Date: 2008/04/29 20:14:04 $
% 	$Id: despike2.m,v 1.13 2008/04/29 20:14:04 duncombe Exp duncombe $
%
% CHANGELOG: 
%   23-Jul-2004, Version 1.0
%	* Initial version.
%   26-Apr-2006 
%       * Fixed comments. Reads cnv_ds1 and incrementally saves to cnv_ds2
%  2007/12/12 CMDR
%  	* Get cruiseid from a file containing the identifier in the data
%  	directory.
%  2007/12/31 CMDR
%	* plot dots as well as lines, makes choosing the desired 
% 	pressures easier 
%  2008/02/17 CMDR
%	* comments and style
%  2008/03/11 CMDR
%	* enter station range (staa, stab) in arguments and provide 
%	defaults if empty
%	* provide default for cruiseid if empty.
%  2008/03/27 CMDR
%	* fixed the assumption that the dir command would provide a nicely
%	sorted array of station numbers. Now sort station numbers and loop 
%	through the station number index.
%  2008/04/17 CMDR
%	* incorporate beckman oxygen sensors
%  2008/04/29 CMDR
%	* make message more meaningful
%

%%%%%%%%%%% CMDR
if nargin==0, cruiseid=getcruiseid; end;
if isempty(cruiseid), cruiseid=getcruiseid; end;
%%%%%%%%%%%
% Retrieve the necessary preferences for this cruise
group = group_name(cruiseid);

% cruise_dir is the main calibration data directory 
if ispref(group,'cal_data_dir')
	cruise_dir=getpref(group,'cal_data_dir');
else
	error(['The cal_data_dir preference was not set for cruise: ' cruiseid '.  Run register_cruise.m']);
end

% Load data from the cruise database
db_file = fullfile(cruise_dir,[cruiseid '_db.mat']);
if exist(db_file)==2,
	ld_file=(['load ',db_file,' cnv_ds1']);
	eval(ld_file)
else,
	error(['The cruise data base file was not found on the path.  ' db_file]);
end;

n=who('cnv_ds1');
if isempty(n)==1,
 	error(['cnv_ds1 not found. run despike1.m ']);
 	return;
end;

beckman=false; 
seabird=false;

if isfield(cnv_ds1,'oxcsec'), 
	beckman=true; 
elseif isfield(cnv_ds1,'oxvsec'), 
	seabird=true; 
end;

% we need to be very sure which order the stations are being accessed, so 
% set up aux_stations so that we do not need to rely on any arbitrary order
% that the variable structure might have been created in.
[aux_stations,ind_stns]=sort(vertcat(cnv_ds1.station));

disp('****************************************** ')
disp('DESPIKE SECONDARY SENSORS')
disp('****************************************** ')
disp('The data found in cnv goes from') 
disp(['the station ',num2str(min(aux_stations)),' to the station ', ...
		num2str(max(aux_stations)),'.'] ) 
disp('****************************************** ')
%
if ~exist('staa'), staa=input('Input first station: '); end;
if ~exist('stab'), stab=input('Input last station: '); end;
if isempty(staa), staa=min(aux_stations); end;
if isempty(stab), stab=max(aux_stations); end;
%% 

% Assumes stations are in "order" in cnv data array (which they should be
% since they were loaded with a "dir" command that puts things in order.
% And it turns out they may not be in the expected order. So aux_stations 
% has been rearranged (sorted) above so that we know the order.
sta_min_idx=find(staa==aux_stations);
sta_max_idx=find(stab==aux_stations);

%initial setup
 cnv_ds2 = cnv_ds1;
 
% and loop through the known order of ind_stns
for i = ind_stns(sta_min_idx:sta_max_idx).',

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% OXYGEN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
% secondary oxv despiking: linearly interpolate over zeros; WOCE flag to 6

% find all good O2
   if seabird,
	jj=find(cnv_ds2(i).oxvsec~=0);  
   elseif beckman,
   	jj=find(cnv_ds2(i).oxcsec~=0);
   end;
   if isempty(jj) ==0
    
% where O2==0  and all pr >= first good O2 and all pr <= last good O2
    if seabird,
    	oxydat=abs(cnv_ds2(i).oxvsec==0);
    elseif beckman,
    	oxydat=abs(cnv_ds2(i).oxcsec==0);
    end;
    
    ii=find(oxydat & ...
        cnv_ds2(i).pressure >= cnv_ds2(i).pressure(jj(1)) & ...
        cnv_ds2(i).pressure <= cnv_ds2(i).pressure(jj(length(jj))));

    if length(ii) > 0
      if seabird,
        cnv_ds2(i).oxvsec(ii)= 	...
		interp1(cnv_ds2(i).pressure(jj),	...
			cnv_ds2(i).oxvsec(jj),	...
			cnv_ds2(i).pressure(ii));
      elseif beckman,
        cnv_ds2(i).oxcsec(ii)= ...
		interp1(cnv_ds2(i).pressure(jj), ...
			cnv_ds2(i).oxcsec(jj), ...
			cnv_ds2(i).pressure(ii));
        cnv_ds2(i).oxctsec(ii)= ...
		interp1(cnv_ds2(i).pressure(jj), ...
			cnv_ds2(i).oxctsec(jj), ...
			cnv_ds2(i).pressure(ii));
      end;
        for n=1:length(ii);
            dum=int2str(cnv_ds2(i).quality(ii(n)));
            if dum(7)=='2'
                cnv_ds2(i).quality(ii(n))=cnv_ds2(i).quality(ii(n))+4;
            end; % if dum(7)
        end; % for n=1:
    end; % for i=sta_min_idx:

% secondary oxv despiking: use cleaner program to despike deep water;
% WOCE flag to 6

    if length(cnv_ds2(i).pressure)>1000;
        ii=find(cnv_ds2(i).pressure>1000);
        factor = 1;
        order = 2;
        if seabird,
        	[pc,pm]=cleaner(cnv_ds2(i).oxvsec(ii),factor,order);
        	jj=find(pm==0);
	elseif beckman,
        	[pc,pm]=cleaner(cnv_ds2(i).oxcsec(ii),factor,order);
        	factor = 8;
        	[pct,pmt]=cleaner(cnv_ds2(i).oxctsec(ii),factor,order);
        	jj=find(pm==0 | pmt==0);
	end;
        for n=1:length(jj);
	  if seabird,
		 oxytest=abs(cnv_ds2(i).oxvsec(ii(jj(n)))-pc(jj(n)))>0.0004;
           elseif beckman,
		oxytest=(abs(cnv_ds2(i).oxcsec(ii(jj(n)))-pc(jj(n)))>0.0004 | ...
		abs(cnv_ds2(i).oxctsec(ii(jj(n)))-pct(jj(n)))>1);
	  end;
            if oxytest,
	if seabird,
                cnv_ds2(i).oxvsec(ii(jj(n)))=pc(jj(n));
	  elseif beckman,
                cnv_ds2(i).oxcsec(ii(jj(n)))=pc(jj(n));
                cnv_ds2(i).oxctsec(ii(jj(n)))=pct(jj(n));
	  end;
                dum=int2str(cnv_ds2(i).quality(ii(jj(n))));
                if dum(7)=='2'
                    cnv_ds2(i).quality(ii(jj(n)))= 	...
			cnv_ds2(i).quality(ii(jj(n)))+4;
                end;
            end;
        end;
    end;

% plot secondary ctd oxygen voltage vs pressure 0-1000 db
    ii=find(cnv_ds2(i).pressure<=1000);
%%%%%%%%%%% CMDR
fh=bigplot; %make larger plot window
    if seabird,
       figure(fh);
       plot(cnv_ds2(i).pressure,cnv_ds2(i).oxvsec,'r.-','linew',1.25);
       axis([0 1000 min(cnv_ds2(i).oxvsec(ii)) max(cnv_ds2(i).oxvsec(ii))]);
    elseif beckman,
       figure(fh);%makes a larger plot window  
       plot(cnv_ds2(i).pressure,cnv_ds2(i).oxcsec,'r.-','linew',1.25);
       axis([0 1000 min(cnv_ds2(i).oxcsec(ii)) max(cnv_ds2(i).oxcsec(ii))]);
    end;
    if (cnv_ds2(i).pressure<1000)
	if seabird,
        axis([0 max(cnv_ds2(i).pressure)+20 ...
		min(cnv_ds2(i).oxvsec(ii)) max(cnv_ds2(i).oxvsec(ii))]);
	elseif beckman,
        axis([  0 max(cnv_ds2(i).pressure)+20   ...
		min(cnv_ds2(i).oxcsec(ii)) max(cnv_ds2(i).oxcsec(ii))]);
	end;
    end
    title(['Station ' int2str(cnv_ds2(i).station) ...
		':  Secondary Oxygen Voltage Despiking'])

% remove artifacts <=1000 db interactively; WOCE flag to 6

    loop1=menu('Interpolate?','yes','no');
    while loop1==1
        title(['Station ' int2str(cnv_ds2(i).station) ...
		':  Zoom in on current area of interest, type ''return'' ' ...
		'when done']);
        zoom on
        keyboard
        zoom off
        title(['Station ' int2str(cnv_ds2(i).station)	...
		':  Click on end-points for interpolation']);
        [x,y]=ginput(2);
        x=round(x);
        kk=x(1)+1:x(2)-1;
        ll=[x(1),x(2)];
	if seabird,
        cnv_ds2(i).oxvsec(kk)=	...
		interp1(cnv_ds2(i).pressure(ll),	...
			cnv_ds2(i).oxvsec(ll),	...
			cnv_ds2(i).pressure(kk));
	elseif beckman,
        cnv_ds2(i).oxcsec(kk)=  ...
		interp1(cnv_ds2(i).pressure(ll),  ...
			cnv_ds2(i).oxcsec(ll),  ...
			cnv_ds2(i).pressure(kk));
        cnv_ds2(i).oxctsec(kk)=  ...
		interp1(cnv_ds2(i).pressure(ll),  ...
			cnv_ds2(i).oxctsec(ll),  ...
			cnv_ds2(i).pressure(kk));
	end;
        for n=1:length(kk);
            dum=int2str(cnv_ds2(i).quality(kk(n)));
            if dum(7)=='2'
                cnv_ds2(i).quality(kk(n))=cnv_ds2(i).quality(kk(n))+4;
            end
        end
%%%%%%%%%%% CMDR
	if seabird,
        bigplot;%makes a larger plot window
        plot(cnv_ds2(i).pressure,cnv_ds2(i).oxvsec,'r.-','linew',1.25);
        axis([0 1000 min(cnv_ds2(i).oxvsec(ii)) max(cnv_ds2(i).oxvsec(ii))]);
	elseif beckman,
        bigplot;%makes a larger plot window
        plot(cnv_ds2(i).pressure,cnv_ds2(i).oxcsec,'r.-','linew',1.25);
        axis([0 1000 min(cnv_ds2(i).oxcsec(ii)) max(cnv_ds2(i).oxcsec(ii))]);
	end;
        if (cnv_ds2(i).pressure<1000)
	 if seabird,
            axis([	0		...
			max(cnv_ds2(i).pressure)+20	...
			min(cnv_ds2(i).oxvsec(ii))	...
			max(cnv_ds2(i).oxvsec(ii))]);
         elseif beckman,
	    axis([	0		...
			 max(cnv_ds2(i).pressure)+20 	...
			min(cnv_ds2(i).oxcsec(ii))	...
			 max(cnv_ds2(i).oxcsec(ii))]);
	  end;
        end;
        title(['Station ' int2str(cnv_ds2(i).station)	...
		':  More Interpolating?'])
        loop1=menu('More Interpolating?','yes','no');     
    end 
    close
else
      disp(['Station ' num2str(cnv_ds2(i).station) ...
		' has no secondary oxygen data. Skipping to the salinity despiking.']);
     
   end;
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% SALINITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% salinity despiking: plot full scale, interactively interpolate over spikes,
% and backcalculate cond, theta, sigma-theta.  repeat for deep water scale.
% Ammend salinity flag to 6.

% plot secondary ctd salinity vs pressure 0-1000 db
    ii=find(cnv_ds2(i).pressure <= 1000);
%%%%%%%%%%% CMDR
    bigplot;%makes a larger plot window
    plot(cnv_ds2(i).pressure, cnv_ds2(i).salsec,'g.-','linew',1.25);
    axis([0 1000 min(cnv_ds2(i).salsec(ii)) max(cnv_ds2(i).salsec(ii))]);
    if (cnv_ds2(i).pressure<1000)
        axis([	0	...
		max(cnv_ds2(i).pressure)+20	...
		min(cnv_ds2(i).salsec(ii))	...
		max(cnv_ds2(i).salsec(ii))]);
    end
    title(['Station ' int2str(cnv_ds2(i).station)	...
		':  Secondary Salinity Despiking'])

% remove artifacts <=1000 db interactively; WOCE flag to 6

    loop1=menu('Interpolate?','yes','no');
    while loop1==1
        title(['Station ' int2str(cnv_ds2(i).station)	...
		':  Zoom in on current area of interest, type'   ...
		' ''return'' when done']);
        zoom on
        keyboard
        zoom off
        title(['Station ' int2str(cnv_ds2(i).station)	...
		':  Click on end-points for interpolation']);
        [x,y]=ginput(2);
        x=round(x);
        kk=x(1)+1:x(2)-1;
        ll=[x(1),x(2)];
        cnv_ds2(i).salsec(kk)=		...
		interp1(cnv_ds2(i).pressure(ll),	...
			cnv_ds2(i).salsec(ll),		...
			cnv_ds2(i).pressure(kk));
        cnv_ds2(i).tempsec(kk)=		...
		interp1(cnv_ds2(i).pressure(ll),	...
			cnv_ds2(i).tempsec(ll),		...
			cnv_ds2(i).pressure(kk));
        cnv_ds2(i).condsec(kk)=		...
		sw_c3515*sw_cndr(cnv_ds2(i).salsec(kk),		...
			cnv_ds2(i).tempsec(kk),		...
			cnv_ds2(i).pressure(kk))/10;
        cnv_ds2(i).thetasec(kk)=		...
		sw_ptmp(cnv_ds2(i).salsec(kk),		...
			cnv_ds2(i).tempsec(kk),		...
			cnv_ds2(i).pressure(kk), 0);
        cnv_ds2(i).sigmasec(kk)=		...
		sw_pden(cnv_ds2(i).salsec(kk),		...
			cnv_ds2(i).thetasec(kk),		...
			cnv_ds2(i).pressure(kk), 0)-1000;
        for n=1:length(kk);
            dum=int2str(cnv_ds2(i).quality(kk(n)));
            if dum(5)=='2'
                cnv_ds2(i).quality(kk(n))=cnv_ds2(i).quality(kk(n))+400;
            end
        end
%%%%%%%%%%% CMDR
        bigplot;%makes a larger plot window
        plot(cnv_ds2(i).pressure, cnv_ds2(i).salsec, 'g.-','linew',1.25);
        axis([0 1000 min(cnv_ds2(i).salsec(ii)) max(cnv_ds2(i).salsec(ii))]);
        title(['Station ' int2str(cnv_ds2(i).station)	...
		':  More Interpolating?']);
        loop1=menu('More Interpolating?','yes','no');     
    end 
    close

% plot primary ctd salinity vs pressure 500-6000 db

    if length(cnv_ds2(i).pressure)>1000;
        ii=find(cnv_ds2(i).pressure>1000);
%%%%%%%%%%% CMDR
        bigplot;%makes a larger plot window
        plot(cnv_ds2(i).pressure,cnv_ds2(i).salsec,'g.-','linew',1.25);
	axis([	500		...
		6000		...
		min(cnv_ds2(i).salsec(ii))		...
		max(cnv_ds2(i).salsec(ii))]);
        if (cnv_ds2(i).pressure<6000)
            axis([	500		...
			max(cnv_ds2(i).pressure)+20		...
			min(cnv_ds2(i).salsec(ii))		...
			max(cnv_ds2(i).salsec(ii))]);
        end
        title(['Station ' int2str(cnv_ds2(i).station)	...
		':  Secondary Salinity Despiking - Deep'])

% remove artifacts >500 db interactively; WOCE flag to 6

        loop1=menu('Interpolate?','yes','no');
        while loop1==1
            title(['Station ' int2str(cnv_ds2(i).station)	...
		':  Zoom in on current area of interest, type'	...
		' ''return'' when done']);
            zoom on
            keyboard
            zoom off
            title(['Station ' int2str(cnv_ds2(i).station)	...
		':  Click on end-points for interpolation']);
            [x,y]=ginput(2);
            x=round(x);
            kk=x(1)+1:x(2)-1;
            ll=[x(1),x(2)];
            cnv_ds2(i).salsec(kk)=		...
		interp1(cnv_ds2(i).pressure(ll),		...
			cnv_ds2(i).salsec(ll),		...
			cnv_ds2(i).pressure(kk));
           cnv_ds2(i).tempsec(kk)=		...
		interp1(cnv_ds2(i).pressure(ll),		...
			cnv_ds2(i).tempsec(ll),		...
			cnv_ds2(i).pressure(kk));
            cnv_ds2(i).condsec(kk)=		...
		sw_c3515*sw_cndr(cnv_ds2(i).salsec(kk),		...
			cnv_ds2(i).tempsec(kk),		...
			cnv_ds2(i).pressure(kk))/10;
            cnv_ds2(i).thetasec(kk)=		...
		sw_ptmp(cnv_ds2(i).salsec(kk),		...
			cnv_ds2(i).tempsec(kk),		...
			cnv_ds2(i).pressure(kk), 0);
            cnv_ds2(i).sigmasec(kk)=		...
		sw_pden(cnv_ds2(i).salsec(kk),		...
			cnv_ds2(i).thetasec(kk),		...
			cnv_ds2(i).pressure(kk), 0)-1000;
            for n=1:length(kk);
                dum=int2str(cnv_ds2(i).quality(kk(n)));
                if dum(5)=='2'
                    cnv_ds2(i).quality(kk(n))=cnv_ds2(i).quality(kk(n))+400;
                end
            end
%%%%%%%%%%% CMDR
            bigplot;%makes a larger plot window
            plot(cnv_ds2(i).pressure,cnv_ds2(i).salsec,'g.-','linew',1.25);
            axis([	500		...
			6000		...
			min(cnv_ds2(i).salsec(ii))		...
			max(cnv_ds2(i).salsec(ii))]);
            if (cnv_ds2(i).pressure<6000)
                axis([	0		...
			max(cnv_ds2(i).pressure)+20		...
			min(cnv_ds2(i).salsec(ii))		...
			max(cnv_ds2(i).salsec(ii))]);
            end
            title(['Station ' int2str(cnv_ds2(i).station)	...
		':  More Interpolating?'])
            loop1=menu('More Interpolating?','yes','no');     
        end 
    end
    close
end

%Save the conductivity summary structure.
save(db_file,'cnv_ds2','-append');

if nargout==0
	assignin('caller','ans',cnv_ds2);
else
    	varargout{1}=cnv_ds2;
end

return;

