%  plotAvgSigz.m - plot horizontal average of sigma zero versus
%                  depth for a set of stations.
%                  Includes one standard deviation envelopes.
%    needs variables set:  DataDir CruiseName DepthResolution ctdcastlist 
%                          srange deprange InfoString
%    sets variables: none
%    calls: readctd

Pden=[27 28];
dep=[0 1];
figure;

plot(Pden,dep);
axis([sigrange deprange]);
title('Average Potential Density vs depth');
xlabel('Sigma 0');
ylabel('Depth (m)');
text(sigrange(1) + .1*(sigrange(2)-sigrange(1)),.75*deprange(1),InfoString);
text(sigrange(1) + .1*(sigrange(2)-sigrange(1)),.81*deprange(1),CruiseName);
hold on;

%              set interpolation depths at 10 m interval
zi=deprange(1):10:deprange(2)-10;
SumS=zeros(size(zi));
SumSqS=zeros(size(zi));
NumberS=zeros(size(zi));

for cast = ctdcastlist
  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);
  dep=-p;
  si=interp1(dep,r0,zi,'linear');
  i=find(~isnan(si));
  SumS(i)=SumS(i)+si(i);
  SumSqS(i)=SumSqS(i)+si(i).^2;
  NumberS(i)=NumberS(i)+1;
end
AvgS=SumS./NumberS;
VarS=SumSqS./NumberS - (SumS./NumberS).^2;
SdS=sqrt(VarS);
lowerS=AvgS - SdS;
upperS=AvgS + SdS;

plot(AvgS,zi,'b',lowerS,zi,'r--',upperS,zi,'r--'); 

hold off;

[m,nz]=size(zi);
%               print values to a file
fid=fopen(['AvgSigz.' CruiseName '.asc'],'w');
for i=nz:-1:1
   fprintf(fid,'%9.4f  %9.4f  \n',zi(i),AvgS(i));
end 
fclose(fid);

