function out = isopycs(handle,P,isopycnals)
% ISOPYCS  This function plots a temp vs. salinity diagram,
%            contouring the isopycnals at at selected pressure.
%
%            ISOPYCS(HANDLE,PR,SIGMA) draws contours lines of density 
%            SIGMA (kg/m^3) at pressure PR (dbars) anomaly on the current
%            axis of figure(HANDLE).  The freezing point (if visible) will 
%            be indicated.
%
%            ISOPYCS(HANDLE,PR) draws several randomly chosen contours.

% Part of the staplot package. For instructions, see sta_info.m.
%
% from: Oceanography toolbox  v1.4 (WHOI)
%        - R. Pawlowicz 14/Mar/94 (rich@boreas.whoi.edu)
% Notes: RP (WHOI) 9/Dec/91
%                 7/Nov/92 Changed for Matlab 4.0
%                 14/Mar/94 Made P optional.
%
% slightly revised by D Byrne, L-DEO, 6/Nov/95: 
%              eliminated specification of T,S.

% get current axis limits
figure(handle)
S = get(gca,'XLim'); T = get(gca,'YLim');

if (nargin == 2),
  isopycnals=6;
else
  % Find levels and convert to columns to be on the safe side
  isopycnals=isopycnals(find(~isnan(isopycnals)));
  isopycnals=isopycnals(:);
end;

% grid points for contouring
Sg=S(1)+[0:30]/30*(S(2)-S(1));
Tg=T(1)+[0:30]'/30*(T(2)-T(1));

[SV,SG] = swstate(ones(size(Tg))*Sg,Tg*ones(size(Sg)),P(1));

if ~isstr(get(gca,'color'))
  if all(get(gca,'color') == [1 1 1])
    [CS,h] = contour(Sg, Tg, SG, isopycnals, 'k:');
  else
    [CS,h] = contour(Sg, Tg, SG, isopycnals, 'w:');
  end
else
  if strcmp(get(gca,'color'),'none')
    if all(get(gcf,'color') == [1 1 1])
      [CS,h] = contour(Sg, Tg, SG, isopycnals, 'k:');
    else
      [CS,h] = contour(Sg, Tg, SG, isopycnals, 'w:');
    end
  else
    [CS,h] = contour(Sg, Tg, SG, isopycnals, 'w:');
  end
end
%[CS,h] = contour(Sg, Tg, SG, isopycnals, 'w:');
out = h; out = out(:);
axis([S(1) S(2) T(1) T(2)]);

%plot freezing temp.
freezeT = swfreeze(S, P(1));
if ~isstr(get(gca,'color'))
  if all(get(gca,'color') == [1 1 1])
    h = line(S,freezeT,'LineStyle','--','color','k');
  else
    h = line(S,freezeT,'LineStyle','--');
  end
else
  if strcmp(get(gca,'color'),'none')
    if all(get(gcf,'color') == [1 1 1])
      h = line(S,freezeT,'LineStyle','--','color',[0 0 0]);
    else
      h = line(S,freezeT,'LineStyle','--');
    end
  else
    h = line(S,freezeT,'LineStyle','--');
  end
end
out = [out; h];

if (max(size(isopycnals))==1), 
  h = isolabel(CS);
else 
  h = isolabel(CS,isopycnals);
end;
out = [out; h];

% Label isopycnals

h = text(1,1.1,[' Pressure = ' int2str(P(1)) ' dbars'],...
    'units','normalized',...
    'HorizontalAlignment','right','VerticalAlignment','middle');
out = [out; h];
return
