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

dep=[-1 01];
PT=[10 11];
figure;

plot(PT,dep);
axis([trange deprange]);

title('Average Potential Temperature vs depth');
xlabel('Potential Temperature (C)');
ylabel('Depth (m)');
text(trange(1) + .1*(trange(2)-trange(1)), .75*deprange(1),InfoString);
text(trange(1) + .1*(trange(2)-trange(1)),.81*deprange(1),CruiseName);
hold on;

%              set interpolation depths at 10 m interval
zi=deprange(1):10:deprange(2)-10;
SumPT=zeros(size(zi));
SumSqPT=zeros(size(zi));
NumberPT=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;
  ti=interp1(dep,pt,zi,'linear');
  i=find(~isnan(ti));
  SumPT(i)=SumPT(i)+ti(i);
  SumSqPT(i)=SumSqPT(i)+ti(i).^2;
  NumberPT(i)=NumberPT(i)+1;
end
AvgPT=SumPT./NumberPT;
VarPT=SumSqPT./NumberPT - (SumPT./NumberPT).^2;
SdPT=sqrt(VarPT);
lowerPT=AvgPT - SdPT;
upperPT=AvgPT + SdPT;

plot(AvgPT,zi,'b',lowerPT,zi,'r--',upperPT,zi,'r--'); 

hold off;

