function varargout=despike1_beckman(cruiseid,staa,stab)
% DESPIKE1_BECKMAN -  Remove spikes in the primary sensors
%             from CTD Data    
%
% CTD Calibration toolbox.
%
% USAGE
%	cnv_ds1=despike1(cruiseid,staa,stab);
%
% INPUT
%	cruiseid	cruise identification name
%   
% OUTPUT
%	cnv_ds1  
%          
% DESCRIPTION: This function uses the cnv files from the cruise database
% and remove the spikes related to the primary sensors
%

% AUTHORS:Carlos Fonseca 
%         UM/CIMAS
%         Fri Dec 2 09:19:37 EDT 2005
%
% VERSION: 
%	$Revision: 1.4 $
%	$Date: 2008/04/16 21:36:13 $
%	$Id: despike1_beckman.m,v 1.4 2008/04/16 21:36:13 duncombe Exp duncombe $
%
% CHANGELOG: 
%   23-Jul-2004, Version 1.0
%		* Initial version.
%   24-Mar-2006 Molly Baringer
%       * Fixed comments.  Reads cnv_fill and incrementally saves to
%       cnv_ds1
%  2008/04/10 CMDR
%	* incorporate changes in despike1.m
%  2008/04/16 CMDR
%	* beckman processing now in despike1.m This function is now redundant.
% 
%

%%%%%%%%%%% 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]);
    eval(ld_file)
else
	error(['The cruise data base file was not found on the path.  ' db_file]);
end

n=who('cnv_fill');
if isempty(n)==1
    error(['cnv_fill not found. run fill_surface.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_fill.station));

disp('****************************************** ')
disp('DESPIKE PRIMARY 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([cnv_fill.station]); end;
if isempty(stab), stab=max([cnv_fill.station]); 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
% loop through the known order of ind_stns
for i=ind_stns(sta_min_idx:sta_max_idx).',
    cnv_ds1(i) = cnv_fill(i);

% primary oxv despiking: linearly interpolate over zeros; WOCE flag to 6
    % All good O2
    jj=find(cnv_ds1(i).oxcpri~=0);                                      
	% O2==0 ;  all pr >= first good O2 ; all pr <= last good O2
    % but only if there are good oxygens
    if ~isempty(jj),
    ii=find(abs(cnv_ds1(i).oxcpri==0) & ...
        cnv_ds1(i).pressure >= cnv_ds1(i).pressure(jj(1)) & ...
        cnv_ds1(i).pressure <= cnv_ds1(i).pressure(jj(length(jj))));
    if length(ii) > 0,
        cnv_ds1(i).oxcpri(ii)=interp1(cnv_ds1(i).pressure(jj), ...
		cnv_ds1(i).oxcpri(jj),cnv_ds1(i).pressure(ii));
        cnv_ds1(i).oxctpri(ii)=interp1(cnv_ds1(i).pressure(jj), ...
		cnv_ds1(i).oxctpri(jj),cnv_ds1(i).pressure(ii));
        for n=1:length(ii);
            dum=int2str(cnv_ds1(i).quality(ii(n)));
            if dum(6)=='2'
                cnv_ds1(i).quality(ii(n))=cnv_ds1(i).quality(ii(n))+40;
              end; % if dum(6)=='2'
          end; % for n=1:length(ii);
      end; % if length(ii) > 0,
    end; % if ~isempty(jj),

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

    if length(cnv_ds1(i).pressure)>1000;
        ii=find(cnv_ds1(i).pressure>1000);
        factor = 1;
        order = 2;
        [pc,pm]=cleaner(cnv_ds1(i).oxcpri(ii),factor,order);
        factor=8;
         [pct,pmt]=cleaner(cnv_ds1(i).oxctpri(ii),factor,order);
        jj=find(pm==0);
        for n=1:length(jj),
            if (abs(cnv_ds1(i).oxcsec(ii(jj(n)))-pc(jj(n)))>0.0004 | ...
		abs(cnv_ds1(i).oxctsec(ii(jj(n)))-pct(jj(n)))>1),
                cnv_ds1(i).oxcpri(ii(jj(n)))=pc(jj(n));
                cnv_ds1(i).oxctpri(ii(jj(n)))=pct(jj(n));
                dum=int2str(cnv_ds1(i).quality(ii(jj(n))));
                if dum(6)=='2',
                    cnv_ds1(i).quality(ii(jj(n)))= ...
			cnv_ds1(i).quality(ii(jj(n)))+40;
                end;
            end;
        end;
    end;

% plot primary ctd oxygen voltage vs pressure 0-1000 db
    ii=find(cnv_ds1(i).pressure<=1000);
%%%%%%%%%%% CMDR
fh=figure; 
    plot(cnv_ds1(i).pressure,cnv_ds1(i).oxcpri,'r.-');
% we test if there are valid data for oxcpri
    if max(cnv_ds1(i).oxcpri(ii))-min(cnv_ds1(i).oxcpri(ii)) <= eps,
    axis([0 1000 min(cnv_ds1(i).oxcpri(ii))-1 max(cnv_ds1(i).oxcpri(ii))+1]);
    else
    axis([0 1000 min(cnv_ds1(i).oxcpri(ii)) max(cnv_ds1(i).oxcpri(ii))]);
    end;
    if (cnv_ds1(i).pressure<1000)
	if max(cnv_ds1(i).oxcpri(ii))-min(cnv_ds1(i).oxcpri(ii))<=eps,
        axis([0 max(cnv_ds1(i).pressure)+20  ...
		min(cnv_ds1(i).oxcpri(ii))-1 max(cnv_ds1(i).oxcpri(ii))+1]);
	else
        axis([0 max(cnv_ds1(i).pressure)+20  ...
		min(cnv_ds1(i).oxcpri(ii)) max(cnv_ds1(i).oxcpri(ii))]);
	end;
    end
    title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Primary Oxygen Voltage Despiking'])

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

    loop1=menu('Interpolate?','yes','no');
    while loop1==1
        title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Zoom in on current area of interest, type ' ...
		'''return'' when done']);
        zoom on
        keyboard
        zoom off
        title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Click on end-points for interpolation']);
        [x,y]=ginput(2);
        x=round(x);
%%%%%% CMDR interp complains if the points are the same;
%%%%%% and dies if the order is reversed.
	if diff(x)<0, x=flipud(x); end;
	if diff(x)>1,
        kk=x(1)+1:x(2)-1;
        ll=[x(1),x(2)];
        cnv_ds1(i).oxcpri(kk)=interp1(cnv_ds1(i).pressure(ll), ...
		cnv_ds1(i).oxcpri(ll),cnv_ds1(i).pressure(kk));
        cnv_ds1(i).oxctpri(kk)=interp1(cnv_ds1(i).pressure(ll), ...
		cnv_ds1(i).oxctpri(ll),cnv_ds1(i).pressure(kk));
        for n=1:length(kk);
            dum=int2str(cnv_ds1(i).quality(kk(n)));
            if dum(6)=='2',
                cnv_ds1(i).quality(kk(n))=cnv_ds1(i).quality(kk(n))+40;
            end;
        end;
%%%%%%%%%%% CMDR
figure(fh);
        plot(cnv_ds1(i).pressure,cnv_ds1(i).oxcpri,'r.-');
        axis([0 1000 min(cnv_ds1(i).oxcpri(ii)) max(cnv_ds1(i).oxcpri(ii))]);
        if (cnv_ds1(i).pressure<1000)
            axis([0 max(cnv_ds1(i).pressure)+20 ...
		min(cnv_ds1(i).oxcpri(ii)) max(cnv_ds1(i).oxcpri(ii))]);
        end;
else
	disp('Choose distinct x values');
end;
        title(['Station ' int2str(cnv_fill(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 primary ctd salinity vs pressure 0-1000 db
    ii=find(cnv_ds1(i).pressure <= 1000);
%%%%%%%%%%% CMDR
fh=figure; 
    plot(cnv_ds1(i).pressure, cnv_ds1(i).salpri,'g.-');
    axis([0 1000 min(cnv_ds1(i).salpri(ii)) max(cnv_ds1(i).salpri(ii))]);
    if (cnv_ds1(i).pressure<1000)
        axis([ 0 max(cnv_ds1(i).pressure)+20 ...
		min(cnv_ds1(i).salpri(ii)) max(cnv_ds1(i).salpri(ii))]);
    end
    title(['Station ' int2str(cnv_fill(i).station) ...
		':  Primary Salinity Despiking'])

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

    loop1=menu('Interpolate?','yes','no');
    while loop1==1
        title(['Station ' int2str(cnv_fill(i).station) ...
		':  Zoom in on current area of interest, type ' ...
		'''return'' when done']);
        zoom on
        keyboard
        zoom off
        title(['Station ' int2str(cnv_fill(i).station) ...
		':  Click on end-points for interpolation']);
        [x,y]=ginput(2);
        x=round(x);
%%%%%% CMDR interp complains if the points are the same; 
%%%%%% and dies if the order is reversed.
	if diff(x)<0, x=flipud(x); end;
	if diff(x)>1,
        kk=x(1)+1:x(2)-1;
        ll=[x(1),x(2)];
          cnv_ds1(i).salpri(kk)=interp1(cnv_ds1(i).pressure(ll), ...
		cnv_ds1(i).salpri(ll),cnv_ds1(i).pressure(kk));
          cnv_ds1(i).temppri(kk)=interp1(cnv_ds1(i).pressure(ll), ...
		cnv_ds1(i).temppri(ll),cnv_ds1(i).pressure(kk));
          cnv_ds1(i).condpri(kk)=sw_c3515*sw_cndr(cnv_ds1(i).salpri(kk), ...
		cnv_ds1(i).temppri(kk),cnv_ds1(i).pressure(kk))/10;
          cnv_ds1(i).thetapri(kk)=sw_ptmp(cnv_ds1(i).salpri(kk), ...
		cnv_ds1(i).temppri(kk),cnv_ds1(i).pressure(kk),0);
          cnv_ds1(i).sigmapri(kk)=sw_pden(cnv_ds1(i).salpri(kk), ...
		cnv_ds1(i).thetapri(kk),cnv_ds1(i).pressure(kk),0)-1000;
          for n=1:length(kk),
            dum=int2str(cnv_ds1(i).quality(kk(n)));
            if dum(4)=='2',
                cnv_ds1(i).quality(kk(n))=cnv_ds1(i).quality(kk(n))+4000;
            end;
          end;
%%%%%%%%%%% CMDR
figure(fh);
        plot(cnv_ds1(i).pressure, cnv_ds1(i).salpri, 'g.-');
        axis([0 1000 min(cnv_ds1(i).salpri(ii)) max(cnv_ds1(i).salpri(ii))]);
else
	disp('Choose distinct x values');
end;
        title(['Station ' int2str(cnv_fill(i).station)  ...
		':  More Interpolating?'])
        loop1=menu('More Interpolating?','yes','no');     
    end 
    close

% plot primary ctd salinity vs pressure bottom db

    if length(cnv_ds1(i).pressure)>500;
        ii=find(cnv_ds1(i).pressure>500);
%%%%%%%%%%% CMDR
fh=figure; 
        plot(cnv_ds1(i).pressure,cnv_ds1(i).salpri,'g.-');
        axis([500 6000 min(cnv_ds1(i).salpri(ii)) max(cnv_ds1(i).salpri(ii))]);
        if (cnv_ds1(i).pressure<6000)
            axis([0 max(cnv_ds1(i).pressure)+20  ...
		min(cnv_ds1(i).salpri(ii)) max(cnv_ds1(i).salpri(ii))]);
        end
        title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Primary Salinity Despiking - Deep'])

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

        loop1=menu('Interpolate?','yes','no');
        while loop1==1
            title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Zoom in on current area of interest, type ' ...
		'''return'' when done']);
            zoom on
            keyboard
            zoom off
            title(['Station ' int2str(cnv_fill(i).station)  ...
		':  Click on end-points for interpolation']);
            [x,y]=ginput(2);
            x=round(x);
%%%%%% CMDR interp complains if the points are the same;
%%%%%% and dies if the order is reversed.
	if diff(x)<0, x=flipud(x); end;
	if diff(x)>1,
            kk=x(1)+1:x(2)-1;
            ll=[x(1),x(2)];
            cnv_ds1(i).salpri(kk)=	...
		interp1(cnv_ds1(i).pressure(ll),	...
			cnv_ds1(i).salpri(ll),	...
			cnv_ds1(i).pressure(kk));
            cnv_ds1(i).temppri(kk)=	...
		interp1(cnv_ds1(i).pressure(ll),	...
			cnv_ds1(i).temppri(ll),	...
			cnv_ds1(i).pressure(kk));
            cnv_ds1(i).condpri(kk)=	...
		sw_c3515*sw_cndr(	...
			cnv_ds1(i).salpri(kk),	...
			cnv_ds1(i).temppri(kk),	...
			cnv_ds1(i).pressure(kk))/10;
            cnv_ds1(i).thetapri(kk)=	...
		sw_ptmp(cnv_ds1(i).salpri(kk),	...
			cnv_ds1(i).temppri(kk),	...
			cnv_ds1(i).pressure(kk),0);
            cnv_ds1(i).sigmapri(kk)=	...
		sw_pden(cnv_ds1(i).salpri(kk),	...
			cnv_ds1(i).thetapri(kk),	...
			cnv_ds1(i).pressure(kk),0)-1000;
            for n=1:length(kk);
                dum=int2str(cnv_ds1(i).quality(kk(n)));
                if dum(4)=='2'
                    cnv_ds1(i).quality(kk(n))=cnv_ds1(i).quality(kk(n))+4000;
                end
            end
%%%%%%%%%%% CMDR
figure(fh); 
            plot(cnv_ds1(i).pressure,cnv_ds1(i).salpri,'g.-');
            axis([	500 ...
			6000 ...
			min(cnv_ds1(i).salpri(ii)) ...
			max(cnv_ds1(i).salpri(ii))]);
            if (cnv_ds1(i).pressure<6000)
                axis([  0 ...
			max(cnv_ds1(i).pressure)+20 ...
			min(cnv_ds1(i).salpri(ii)) ...
			max(cnv_ds1(i).salpri(ii))]);
            end
else
	disp('Choose distinct x values');
end;
            title(['Station ' int2str(cnv_fill(i).station) ...
		':  More Interpolating?'])
            loop1=menu('More Interpolating?','yes','no');     
        end 
    end
    close
end

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

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

  
