function [suspect, ends] = inversions(z, rhoz, pflag, cst)
% function [suspect, ends] = inversions(z, rhoz, pflag, cst)
%
% Overturn Toolbox
% 2005 Mar 08 - last update 
% Original Code by David S. Ciochetto
% Department of Oceanography at Dalhousie University
%
% A function to take a marine density profile (rhoz) and pick out suspected
% overturning regions. The vertical axis (z) is assumed to be pointed
% positive upward with the origin at the water surface. It will pick overturns
% for an entire profile as described in Galbraith and Kelley (1996). Starting
% at the surface, it finds the first point where all the points above weigh
% less then it finds the next point where all the data below weigh more. This
% is repeated for the entire profile. The profiles that are sent should not have
% depth loops, i.e. z should be monotonically increasing. The suspected overturns
% are indicated as true values in a logical array returned with the same size
% as the input rhoz and z arrays.
%
% If you are using this function independent of the Overturn Toolbox, then you
% should feed the raw CTD or density profile data to this function rather than
% smoothed data. It still should not contain vertical loops as this function does
% no sorting of the data sent to it.
%
% INPUT:
%		z			= depth variable, should be increasingly negative and monotonic
%		rhoz	= density profile from surface to max depth matching z
%		pflag	= plot flag, send 'y' to make a plot highlighting the suspect overturns
%		cst		= cast number to append to the check plot title, can also add any string
%
% OUTPUT:
%		suspect	= logical array indicating true (1) for the indices of suspect overturns
%							in an array the same size as the z and rhoz input data
%		ends = the indicies for the ends of the suspect overturns, [start stop]
%							in a two column array the same length as the z and rhoz input data
%
% NOTES:
%		0. Please send the data vectors in columns, with z=0 at index 1 (the top)
%		1.  The function works with z becoming more negative with depth. If the data sent
%		does not fit this, it is converted internally to match the routine. The output is
%   then flipped back to match the input data.
%		2. Default values for pflag will be set to zero to avoid unwanted plots.
%		3. This function is part of the package to calculate overturns based on the method
%			of Galbraith and Kelley (1996).
%		4. The input should be a vector, do not send z and rhoz as matricies.
%   5. THIS IS DONE AS SPECIFIED IN THE PAPER. IN THE reorange code it is done differently.
%      THUS IT WILL NOT BE DEVELOPED FURTHER OR INCLUDED IN THE PACKAGE.
%
% REFERENCES:
%		Galbraith, P. S. and Kelley, D. E., 1996 "Identifying overturns in CTD profiles,"
%			 Journal of Atmospheric and Oceanic Technology, 13, 688-702.
%
% USES: ooz.m, linkedzoom.m
%
% SEE ALSO: 
%

% Did you forget to send data?
if ~exist('z')
	disp('inversions.m:  Missing z data')
	return
end

if ~exist('rhoz')
	disp('inversions.m:  Missing rho(z) data')
	return
end

z=z(:);															% Convert to a column if not.
rhoz = rhoz(:);

if ~exist('pflag'), pflag='n', end	% Default values
if ~exist('cst'), cst=[], end

% Make sure that z is increasingly negative downward (i.e. as index increases)
zpos = 'n';													% Keep track if z data is (+) or (-)
if sum(z>0)/length(z) > 0.95				% If > 95% of the data are (+) tag the profile as positive.
	zpos = 'y';
	z = -z;														% Want z to be increasingly negative
end
zflip = 'n';												% What is the orientation of the z data with index? Keep track.
if sum(diff(z)>0)/(length(z)-1) > 0.95	% If 95% of the values are decreasing, then flip it.
	zflip = 'y';
	z = flipud(z);
	rhoz = flipud(rhoz);
end

% %%% Check the procedure: (COMMENT THESE LINES OUT WHEN RUNNING WITH REAL DATA)
% %% Run Zoverturn first to generate a fake profile if needed.
% clear all;close all;Zoverturn; clear H L Lz N SR SS T a1 a2 dr drdz drho g idx ot r res rho0 t1 t2 y; figure;
% rhoz=rho; pflag='y'; cst='Linear Density Gradient with 1 Z-overturn and uniform random noise';
% %%% End (generate fake signal)

suspect = zeros(size(z));						% Initialize output variable(s)
ends = [];
Q = 0; T = 0;												% Initialize internal variables
QL=1; TL=1;
slopes = suspect;
sample = suspect;
chkrho=rhoz; chkz=z;

% Calculate which way the changes in density are, slopes, to indicate d(rhoz)/dz for point pairs
slopes = [diff(rhoz); 0];
sample = slopes;										% Will use this varible to loop through the profile in suspect chunks
while sum(sample<0) > 0							% If negative slopes still exist, there may be an inversion
	% Find the index of the shallowest (-) slope
	Q = min(find(sample<0));
	% Find the index for the first point heavier than our flagged point
	T = min(find(chkrho>chkrho(Q)));
	% Get the indicies with respect to the non-shrinking reference z
	% BEWARE multiple z - values, find the exact match with respect to rho and z
	QL=find( z==chkz(Q) & rhoz==chkrho(Q) );
	% Watch out for points with duplicate values; take the first occurance
	if length(QL) > 1, QL=min(QL); end
	TL=find( z==chkz(T) & rhoz==chkrho(T) );
	if length(TL) > 1, TL=min(TL); end
% 	%%% Check what is going on.
% 	clf; subplot 121, plot(chkrho(Q:T),chkz(Q:T),'-y','LineWidth',6);set(gca,'FontSize',12,'FontWeight','bold');
% 	ylabel('z [m]');xlabel('\rho(z) [kg m^{-3}]');hold on;
% 	plot(rhoz,z,'.-k',rhoz(QL),z(QL),'co',rhoz(TL),z(TL),'bo','MarkerFaceColor','r');
% 	plot(chkrho(Q),chkz(Q),'gv','MarkerFaceColor','g');plot(chkrho(T),chkz(T),'r^','MarkerFaceColor','r');hold off;
% 	subplot 122, plot(sample(Q:T),chkz(Q:T),'-y','LineWidth',6);hold on;
% 	plot(slopes,z,'.-k',slopes(QL),z(QL),'co',slopes(TL),z(TL),'bo','MarkerFaceColor','r');
% 	plot(sample(Q),chkz(Q),'gv','MarkerFaceColor','g');plot(sample(T),chkz(T),'r^','MarkerFaceColor','r');hold off;
% 	set(gca,'FontSize',12,'FontWeight','bold');xlabel('\Delta\rho(z) [kg m^{-3}]');yl=get(gca,'YLim');
% 	line([0 0],[yl]);linkedzoom ony; pause(0.01);clear yl; pause
% 	%% Don't Panic, it seems that all noise fluctuations will be labeled as an overturn and my returned index will
% 	%% be all ones. I think that Galbraith and Kelley (1996) gave recommendations on the smallest run length that
% 	%% one should expect. I can filter those small than that out here or wait til the main program. It may be
% 	%% better to impliment it here such that a vector of all true for suspect overturn is not sent back. This
% 	%% function really depends on a zero or false delineating between the suspect overturns. I may need to return
% 	%% a vector of start and stop indicies rather than true and false. Will consider this.
% 	%%% End check and comment section.
	suspect(QL:TL) = 1;								% Make the indicies in suspect to match the entire range true (1)
	sample = sample(T:end);						% Take a new sample from the profile slopes to seek remaining overturns
	chkrho = chkrho(T:end);
	chkz = chkz(T:end);
	Q = 0; T = 0;											% Reset Q and T indicies for shrinking sample vector
	ends = [ends; QL TL];
	% Note: no need for a stop condition due to length, if the final sample is negative, sample is empty and it stops
	% Therefore it should not loop indefinitely. I tested this several times and this is the conclusion.
end
clear chkz chkrho QL TL Q T sample slopes ans	% Clean up

% Make the result a logical array
suspect = logical(suspect);
% %% Check results
% figure; bar(suspect);

% If the data sent were flipped in a differnt direction from what was sent, then flip them back
% both so the output matches what was sent but also so the plot will look familiar to the user.
if zflip == 'y'
	suspect = flipud(suspect);
	ends = fliplr(flipud(length(z)-ends + 1));
	rhoz = flipud(rhoz);
	z = flipud(z);
	% zflip = 'n';
end
if zpos == 'y'
	z = -z; 
	% zpos = 'n';
end

% Make a summary plot if it is desired
if lower(pflag) == 'y'
	GS = 'n';		% Greyscale plot?
	ME = 'y';		% Mark the ends? Really only good if you zoom in.
	% Default end marks
	mrk1 = 'kv';
	mrk2 = 'k^';
	% Add NaN values to the logical indicies to allow the OT highlights to not be continuous
	II = suspect;
	II(II<1)=NaN;
	% Plot the highlight
	OT = plot(rhoz.*II,z.*II,'-r','linewidth',9);
	set(gca,'FontSize',14,'FontWeight','bold')
	hold on;
	LN = plot(rhoz,z,'.-b');
	if lower(ME) == 'y'
		if zpos == 'y'
			mrk1 = 'k^';
			mrk2 = 'kv';
		end
		SA = plot(rhoz(ends(:,1)),z(ends(:,1)),mrk1,'MarkerFaceColor','g');
	  FA = plot(rhoz(ends(:,2)),z(ends(:,2)),mrk2,'MarkerFaceColor','c');
	end
	xlabel('\rho(z) [kg m^{-3}]')
	ylabel('z [m]')
	ttl='Suspect Overturns: ';
	if isstr(cst)
		title([ttl cst])
	else
		title([ttl 'Cast' ooz(cst,3)])
	end
	if lower(GS) == 'y'
		gs = 0.50;		% Select level of grey for Overturn highlights
		g1 = 0.70;		% Select level of grey for Overturn start point
		g2 = 0.90;		% Select level of grey for Overturn ending point
		set(OT,'color',[gs gs gs]);
		set(LN,'color','k');
		if lower(ME) == 'y'
			set(SA,'color','k','MarkerFaceColor',[g1 g1 g1]);
			set(FA,'color','k','MarkerFaceColor',[g2 g2 g2]);
		end
	end
	zoom(2); pause(0.01); zoom out  % Hack to smooth out the lines
	hold off; zoom on;
	clear GS II OT LN SA FA ttl ME mrk1 mrk2 g1 g2 gs
end
	