function plotTSRho(data,depthLim,tempLim,salLim,densLim)
%  Data format:  depth, temp, salinity, density

%  Set up the physical limits for the plot axis.
left = 0.15;
bot = 0.1;
width = 1 - (2*left);
height = 0.65;
% add code to handle that fact that there may be NaN's 
% in the data.
vt=isnan(data(:,2))==0;
vs=isnan(data(:,3))==0;
vd=isnan(data(:,4))==0;
sst=size(data(vt,2),1);
sss=size(data(vs,3),1);
ssd=size(data(vd,4),1);
temp=data(vt,2);
salt=data(vs,3);
dens=data(vd,4);

%  Set up the data limits.
if nargin < 2
   depthLim = [min(data(:,1)), max(data(:,1))];
end
depthLim=[0 300];
if max(data(:,1)) > 320,
	depthLim=[0 360];
end
if (size(data,1) > 10),
	depthLim=[0 360];
end
if depthLim(2) > 0
   depthLim = [-depthLim(2), -depthLim(1)];
   data(:,1) = -1 * data(:,1);
end
if length(temp)==0,
	xfil=[0 1 1 0];
	yfil=[0 0 1 1];
	plot(xfil,yfil);
	axis([0 1 0 1]);
	text(0.4,0.5,'NO DATA','FontWeight','bold','FontSize',20);
	return;
end
if nargin < 3,
	if length(temp)==1,
		tempLim=[temp temp+0.1];
	else,
		tempLim = [min(temp), max(temp)];
	end
end
if nargin < 4
   if sss > 0,
	  if sss==1,
		salLim=[salt salt+0.1];
	  else,
	   	salLim = [min(salt), max(salt)];
	  end
   else,
	   salLim = [NaN NaN];
   end
end
if nargin < 5
   if ssd > 0,
	if ssd==1,
		densLim=[dens-0.1 dens];
	else,
   		densLim = [min(dens), max(dens)];
	end
   else,
	densLim = [NaN NaN];
   end
end
   
axisLim1 = [tempLim(1), tempLim(2), depthLim(1), depthLim(2)];
axisLim2 = [salLim(1), salLim(2), depthLim(1), depthLim(2)];
axisLim3 = [densLim(1), densLim(2), depthLim(1), depthLim(2)];

%%axisLim1 = [1,14.5,-3000,0];
%%axisLim2 = [33.3,34.7,axisLim1(3),axisLim1(4)];
%%axisLim3 = [24.75,27.75,axisLim1(3),axisLim1(4)];
if ssd > 0,
	ax3 = axes('position',[left,bot,width,height+0.1]);
	set(ax3,'XAxisLocation','top')
	cl=['\fontsize{20}\sigma_\theta'];
	xlabel(cl,'color','r')
	set(ax3,'ytick',[],'ycolor',[0.999,0.999,0.999],'tickdir','out');
	set(ax3,'xlim',[axisLim3(1),axisLim3(2)],'ylim',[axisLim3(3),axisLim3(4)]);
%%set(ax3,'xtick',axisLim3(1):0.25:axisLim3(2));
	set(ax3,'FontWeight','bold','FontSize',20);
	set(ax3,'xcolor','r');
end
if sss > 0,
	ax2=axes('position',[left,bot,width,height+0.02]);
	set(ax2,'XAxisLocation','top')
	set(ax2,'FontWeight','bold','FontSize',20);
	xlabel('Salinity','color','g','FontWeight','bold','FontSize',20);
	set(ax2,'xlim',[axisLim2(1),axisLim2(2)],'ylim',[axisLim2(3),axisLim2(4)]);
	set(ax2,'ytick',[],'ycolor',[0.999,0.999,0.999],'tickdir','out');
%%set(ax2,'xtick',axisLim2(1):0.1:axisLim2(2));
	set(ax2,'xcolor','g');
end
if sst > 0,
	ax1=axes('position',[left,bot,width,height],'tickdir','out');
	set(ax1,'xlim',[axisLim1(1),axisLim1(2)],'ylim',[axisLim1(3),axisLim1(4)]);
%%set(ax1,'xtick',axisLim1(1):1:axisLim1(2), ...
%%        'ytick',axisLim1(3):250:axisLim1(4));
	xlabel('Temperature ( ^oC)','color','b','FontWeight','bold','FontSize',20)
	ylabel('Depth (m)','FontWeight','bold','FontSize',20)
	set(ax1,'xcolor','b');
end


% Do the data
% Line for axis 1
L1=line(data(:,2),data(:,1),'linestyle','-','color','b');
line(data(:,2),data(:,1),'marker','o','color','b');
set(gca,'Box','on');
if sss > 0,
% Line for axis 2 - remap the salinity data into temperature space
	Map = tempLim(1):0.001:tempLim(2);
	Map = Map(:);
	mapSal = remap(data(vs,3),Map,[salLim(1),salLim(2)]);
	mps(vs)=mapSal;
	mps(~vs)=NaN;
	line(mps,data(:,1),'linestyle','-','color','g') ;
	line(mps,data(:,1),'marker','o','color','g');
end
if ssd > 0,
% Line for axis 3 - remap the density data into temperature space
	mapDen = remap(data(vd,4),Map,[densLim(1),densLim(2)]);
	mpd(vd)=mapDen;
	mpd(~vd)=NaN;
	line(mpd,data(:,1),'linestyle','-','color','r');
	line(mpd,data(:,1),'marker','*','color','r');
end
return;


% ------------------------------------------------------------------------
function mappedData = remap(data,map,dataLimits);

data = data(:);

if nargin < 2
   colormap('jet');
   map = colormap;
end

if nargin < 3
   dataLimits(1) = min(data);
   dataLimits(2) = max(data);
end

numDataPts = size(data,1);
maxMapIndex = size(map,1);
minMapIndex = 1;
map(maxMapIndex+1,1:end) = NaN;

% equation of line: y = m*x + b.  Here b1 = m, and b0 = b.
% max - min index number of map values / max - min values of data for 
% for slope calculation.
b1 = (maxMapIndex - minMapIndex) / (dataLimits(2) - dataLimits(1));
b0 = (-b1 * dataLimits(1)) + 1;
mapIndex = round(b1 .* data + b0);

index = find(mapIndex < 1 | mapIndex > maxMapIndex);
mapIndex(index) = maxMapIndex + 1;
mappedData = map(mapIndex,1:end);

return;
