function varargout=despike2_beckman(cruiseid,staa,stab)
% DESPIKE2_BECKMAN -  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.4 $
%	$Date: 2008/04/17 13:50:33 $
%	%Id:$
%
% CHANGELOG: 
%   23-Jul-2004, Version 1.0
%		* Initial version.
%   26-Apr-2006 
%       * Fixed comments. Reads cnv_ds1 and incrementally saves to cnv_ds2
%  2008/04/10 CMDR 
% 	* incorporate changes in despike2.m

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

% 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).',

% secondary oxv despiking: linearly interpolate over zeros; WOCE flag to 6
   
% All good O2
   jj=find(cnv_ds2(i).oxcsec~=0);
   if isempty(jj),
      disp(['Station ' num2str(cnv_ds2(i).station) ...
		' has no good oxygen data. Skipping.']);
      continue;
   end;

% where O2==0  and all pr >= first good O2 and all pr <= last good O2
    ii=find(abs(cnv_ds2(i).oxcsec==0) & ...
        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
        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));
        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:

% primary oxc despiking: use cleaner program to despike deep water;
% WOCE flag to 6

    if length(cnv_ds2(i).pressure)>500;
        ii=find(cnv_ds2(i).pressure>500);
        factor = 1;
        order = 2;
        [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);
        for n=1:length(jj);
            if (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)
                cnv_ds2(i).oxcsec(ii(jj(n)))=pc(jj(n));
                cnv_ds2(i).oxctsec(ii(jj(n)))=pct(jj(n));
                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-1000db
    ii=find(cnv_ds2(i).pressure<=1000);
%%%%%%%%%%% CMDR
    plot(cnv_ds2(i).pressure,cnv_ds2(i).oxcsec,'r.-');
    axis([0 1000 min(cnv_ds2(i).oxcsec(ii)) max(cnv_ds2(i).oxcsec(ii))]);
    if (cnv_ds2(i).pressure<1000)
        axis([  0 max(cnv_ds2(i).pressure)+20   ...
		min(cnv_ds2(i).oxcsec(ii)) max(cnv_ds2(i).oxcsec(ii))]);
    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)];
        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));
        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
        plot(cnv_ds2(i).pressure,cnv_ds2(i).oxcsec,'r.-');
        axis([0 1000 min(cnv_ds2(i).oxcsec(ii)) max(cnv_ds2(i).oxcsec(ii))]);
        if (cnv_ds2(i).pressure<1000)
            axis([0 max(cnv_ds2(i).pressure)+20 ...
		min(cnv_ds2(i).oxcsec(ii)) max(cnv_ds2(i).oxcsec(ii))]);
        end;
        title(['Station ' int2str(cnv_ds2(i).station)  ...
		':  More Interpolating?'])
        loop1=menu('More Interpolating?','yes','no');     
    end 
    close

% 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);
    plot(cnv_ds2(i).pressure, cnv_ds2(i).salsec,'g.-');
    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
        plot(cnv_ds2(i).pressure, cnv_ds2(i).salsec, 'g.-');
        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 secondary ctd salinity vs pressure bottom part db

    if length(cnv_ds2(i).pressure)>500;
        ii=find(cnv_ds2(i).pressure>500);
        plot(cnv_ds2(i).pressure,cnv_ds2(i).salsec,'g');
        axis([0 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) ...
		':  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
            plot(cnv_ds2(i).pressure,cnv_ds2(i).salsec,'g.-');
            axis([0 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;

