%  Sigsection.m  - contour potential density from ctd on a vertical section.
%
%   need variables set: RefPress xrange deprange
%                       DataDir CruiseName DepthResolution Line
%                       SectionNx SectionNz ctdcastlist
%                       InfoX InfoZ InfoString
%   sets variables: none
%   calls: readctd swstate
%

figtitle=['Potential Density for ' Line ' Section'];

%        open an new graphics window.
figure;

%              set interpolation depths
zi=0:SectionNz-1;
zi=zi./(SectionNz-1);
zi=deprange(1) + (deprange(2)-deprange(1))*zi;

Darray=ones(SectionNz,SectionNx);
Zarray=ones(SectionNz,SectionNx);
Xarray=ones(SectionNz,SectionNx);
bottom=ones(1,SectionNx);
X=ones(1,SectionNx);
Y=ones(1,SectionNx);

[n,stations]=size(ctdcastlist);
Casts=0;
for station=1:stations
%           count number of casts
   Casts=Casts+1;
   cast=ctdcastlist(station);
   file= [DataDir '/' CruiseName '/Ctd/' DepthResolution '/ctd' castnum(cast) '.dat'];
   [nlines,ndepths,desc,p,t,s,pt,r0,r1,bv,t0,t1,s0,s1,o2,par,tr,fl] = readctd(file);
   x=ctdxgrid(cast);
   y=ctdygrid(cast);
   bdepth=ctdtotdepth(cast);
%          convert pressure to depth
   Z=-p;
   [SVAN,Pden]=swstate(s,pt,RefPress);
%              fill array with observations
   bottom(1,station)=-bdepth;
   X(1,station)=x;
   Y(1,station)=y;
   [n,nz]=size(Z);
   ri=interp1(Z,Pden,zi,'linear');
   Darray(1:SectionNz,station)=ri';
   Zarray(1:SectionNz,station)=zi';
end
dx=X(1,1:station)-X(1,1);
dy=Y(1,1:station)-Y(1,1);
dx2=dx.*dx;
dy2=dy.*dy;
if Line(1) == 'x'
%             cross shelf arrays
   Xarray(1,1:station)=Y(1:station);
elseif Line(1) == 'y'
%             along shelf arrays
   Xarray(1,1:station)=X(1:station);
else
   Xarray(1,1:station)=sqrt( (X(1,1:station)-X(1,1)).^2 ...
                            +(Y(1,1:station)-Y(1,1)).^2);
end
for i=2:SectionNz
   Xarray(i,1:station)=Xarray(1,1:station);
end

[C,h]=contour(Xarray(:,1:station),Zarray(:,1:station),Darray(:,1:station),cl);
clabel(C,h);
title(figtitle);
xlabel('distance (km)');
ylabel('depth(m)');
text(InfoX,InfoZ,InfoString);
text(InfoX,1.1*InfoZ,CruiseName);
axis([xrange deprange]);
hold on;
plot(Xarray(1,1:station),bottom(1:station));
for station=1:stations
  text(Xarray(1,station),-20,'T');
end
hold off;

%disp('Sigsection is done')
%    end of Sigsection.m
