% show_sigvar.m - show variation of sigma +-60m from 26.26

cmax = 7; % last CTD cast number
figure(cmax+2)

clf
mat_path = '/mnt/oc415-4a/oc415_4/ctd/process/';

tvec=[];
svec=[];
jj=1;

col=['b','r','k','y','c','m','g','y'];
for ii = [1:3 5:cmax]
  if ii <10
    fn = ['oc41540',int2str(ii)];
  elseif cast <100
    fn = ['oc4154',int2str(ii)];
  end

  %eval(['load ',mat_path,fn])
  % use 1m averages
  eval(['load ',mat_path,'d',fn,'avg'])

  aa=find(s0 < 30);
  s0(aa)=nan;
  
  sig = sw_pden(s0,t0,press,0)-1000;
  
  if 0
    % for downcast
    aa=find(press == max(press));
    aa = 1:aa(1);
  else
    % using 1m avg
    aa=1:length(s0);
  end

  % find closest sigma
  dd=find( min(abs(sig(aa)-26.26)) == abs(sig(aa)-26.26));
  % is it close enough?

  if sig(aa(dd))- 26.26 < .0001
    % find depth at closest sigma
    p26d(ii)= press(aa(dd));
  
    % 60m above
    cc=find(min(abs(press(aa)-(p26d(ii)-60))) == ...
	    abs(press(aa)-(p26d(ii)-60)));
    
    % 60m below
    ee=find(min(abs(press(aa)-(p26d(ii)+60))) == ...
	    abs(press(aa)- (p26d(ii)+60)));
    ee=ee(1);
    
    ind=[cc:ee];

    if ii == 6
      indx = 1:2
      t0(ind(indx))=nan;
      sig(ind(indx))=nan;
      press(ind(indx))=nan;
    end

    gg=find(~isnan(sig(aa(ind))));
    tvec(jj:jj+length(gg)-1) = t0(aa(ind(gg)));
    svec(jj:jj+length(gg)-1) = sig(aa(ind(gg)));
    jj=length(tvec)+1;
    
    plot(sig(aa(cc:ee)),t0(aa(cc:ee)),'.','color',col(ii))
    if ~ishold,hold on;end
    keyboard
  end
%  tvec(26481:26482)=[];
%  svec(26481:26482)=[];
%  tvec(26124:26125)=[];
%  svec(26124:26125)=[];
  
  if 0
    % for upcast
    bb=find(press == max(press));
    bb = bb(1):length(press);
  
    % find closest sigma
    uu=find( min(abs(sig(bb)-26.26)) == abs(sig(bb)-26.26));
    % find depth at closest sigma
    p26u(ii)= press(bb(uu));
  end
  
end
return
  %want to make an equation to map T from tpods to sigma
  [p,s]=polyfit(tvec,svec,2);
  

  xxx=min(tvec):.25:max(tvec);
  yyy=polyval(p,xxx);
  figure(3)
  plot(tvec,svec,'.',xxx,yyy,'r')
  
% p= -0.0088    0.0992   27.6111

% now need a calibration between tpods and sled (see t-p-pod
% directory)
% have these factors but revisit the T/Sigma calibration and 
% make a table from which to interpolate rather than the 
% polynomial above

% make t bins
%binsize = (25.1-18.9)/50; % 51 bins
binsize = (25.1-18.9)/10; % 11 bins
tbins = [18.9:binsize:25.1];
clear tline
aa=ones(length(tbins),length(tvec))* nan;
for ii = 1:length(tbins-1)
  bb=find(tvec >= tbins(ii) & tvec<= tbins(ii+1));
  aa(ii,1:length(bb))= bb;
  cc(ii)=length(bb);
  if cc(ii) >= 3
     [p,s] = polyfit(tvec(aa(ii)),svec(aa(ii,:)),1);
     tline(ii,1)=polyval(p,tbins(ii));
     tline(ii,2)=polyval(p,tbins(ii+1));  
  else
    tline(ii,1) = min(svec(aa,:));
    tline(ii,2) = max(svec(aa,:));

  end
  
end

plot(tvec,svec,'.')
hold on
plot(tline(:,1),'r')
plot(tline(:,2),'c')
