% plotdat.m - to plot cast data
% this one for EDDIES 07/05 OC415-2
% for cast 15 with 4 parts [],b,c,d

printit=1;

set(0,'defaultaxesfontname','timesroman')
set(0,'defaultaxesfontsize',12)

% keep list of starting points for CTD downcast
% cast   1    2     3    4    5   6    7  8    9    10   11    12  13
%st_ind=[1060 1683 4483 3062 2589 817 1238 862 1184 1205 2415 2758 ...
%	1702 2459];

% casts are  deep to shallow
if cast == 15
  % for cast 15 ,15b, 15c, 15d
  st_ind=[1941 2157 2039 2392] ;
  letter=[' ','b','c','d'];
elseif cast == 17
  st_ind=[991 942 2316 2086];
  letter=['d','c','b','a'];
elseif cast == 20
  st_ind=[334 3389 2192 2382 2140];
  letter=['e','d','c','b','a'];
elseif cast == 22
  st_ind=[384 184 247 328];
  letter=['d','c','b','a'];
elseif cast == 24
  st_ind=[243 210 415 459];
  letter=['d','c','b','a'];
end
mat_path = '/mnt/oc415-2a/oc415_2/ctd/process/';

for ii = 1:length(letter)
  %cast = input('Enter cast number: ');

  minth=18;maxth=30;
  minsal=36.6;maxsal=37;
  minsig=23;maxsig=26.5;

  fn = ['oc4152',int2str(cast),letter(ii)];

  eval(['load ',mat_path,fn])

  [Lat,Lon,time,date]= readhdr_bcd(cast,letter(ii));

  % compute quantities
  theta0=sw_ptmp(s0,t0,press,0);
  sigma0=sw_pden(s0,t0,press,0);
    
  % initial plot of pressure vs time
  if 0
    figure
    plot(scan,press)                    
    set(gca,'ydir','reverse')
    orient landscape
    wysiwyg
    xlabel('Scan')
    ylabel('Pressure (db)')
    %title(['Cast: ',fn,' - ',dirstr(ii,:),' cast'])
    title(['Cast: ',fn])
    ax = axis;
    axis([scan(1) scan(length(scan)) -1 1010])
  end
  
  
  % this plots all 3 quantites on one plot
  figure(ii)
  clf
  %hold off
  bb = find(max(press) == press);
  bb = st_ind(ii):bb(1);

  cc=find(press(bb)>0 & s0(bb)>20 & theta0(bb)>5);
  bb=bb(cc);

  if 0
    lim600 = 0;
    if lim600
      aa = find(press(bb) >200 & press(bb) <= 600);
    else
      aa = find(press(bb) >0 & press(bb) <= 200);
    end
    myax = 1;
    %miny = 200;
    miny = 0;
    if lim600
      maxy = 600;
    else
      maxy = 200;
    end
  else
    myax=0;
    aa=1:length(bb);
    miny = 0;
    if ii == 1
      if cast == 15
	maxy = 50;
      elseif cast == 17
	maxy = 65;
      elseif cast == 20
	maxy = 65;
      elseif cast == 22
	maxy = 65;
      elseif cast == 24
	maxy = 45;
      end
    elseif ii == 2
      maxy = 40;
      if cast == 20
	maxy = 55;	
      end
       if cast == 22
	maxy = 50;	
      end
       if cast == 24
	maxy = 40;	
      end
   elseif ii == 3
      maxy = 30;
      if cast == 20
	maxy = 45;	
      end
      if cast == 24
	maxy = 40;	
      end
    elseif ii == 4
      maxy = 20;
      if cast == 20
	maxy = 35;	
      elseif cast == 24
	maxy = 30;	
      end
    elseif ii == 5
      if cast == 20
	maxy = 15;
      end
    end
  end
  
  ntics = 10;
  ax1 = axes('Position',[.13 .11 .7750 .775]);
  h1 = line(theta0(bb(aa)),press(bb(aa)));
  set(ax1,'ydir','reverse')
  if myax
    axlim = axis;
    axis([axlim(1:2) miny maxy])
  end
  axis([minth maxth miny maxy])
  grid
  xlabel(['Theta (',setstr(176),' C) (-)'])
  ylabel('Pressure (db)')
  %title(['Cast: ',fn,' - ',dirstr(ii,:),' cast'])
  %keyboard
  hold on
  ax2 = axes('Position',get(ax1,'Position'),'xaxislocation','top', ...
	     'color','none','ydir','reverse');
  h2 = line(s0(bb(aa)),press(bb(aa)),'linestyle','--','parent',ax2);  
  if myax
    axlim = axis;
    axis([axlim(1:2) miny maxy])
  end
  axis([minsal maxsal miny maxy])
  xlabel('Salinity (PSU) (- - -)')
  %title(['Cast: ',fn,' - ',dirstr(ii,:),' cast ',Lat,' ',Lon,' ', ...
	% date,' ',time])
  title(['Cast: ',fn,'  ',Lat,' ',Lon,' ', ...
	 date,' ',time])
  box on
  xlimits = get(ax1,'XLim');
  xinc = (xlimits(2)-xlimits(1))/ntics;
  set(ax1,'XTick',[xlimits(1):xinc:xlimits(2)])
  x2limits = get(ax2,'XLim');
  xinc = (x2limits(2)-x2limits(1))/ntics;
  set(ax2,'XTick',[x2limits(1):xinc:x2limits(2)])

  ax3 = axes('Position',get(ax1,'Position'),'visible','off', ...
	     'color','none','ydir','reverse');
  h3 = line(sigma0(bb(aa))-1000,press(bb(aa)), ...
	    'linestyle','-.','parent',ax3);
  if myax
    axlim = axis;
    axis([axlim(1:2) miny maxy])
  end
  axis([minsig maxsig miny maxy])

  ax3lim = get(ax3,'XLim');
  xinc = (ax3lim(2)-ax3lim(1))/ntics;
  set(ax3,'XTick',[ax3lim(1):xinc:ax3lim(2)])
  axpos = get(ax1,'Position');
  ax4 = axes('position',[axpos(1) .05 axpos(3) .005]);
  set(ax4,'xlim',ax3lim)
  set(ax4,'XTick',[ax3lim(1):xinc:ax3lim(2)],'yticklabel',[])
  txl=xlabel('Sigma (kg/m^3) ( _ . _ . _ )');
  set(txl,'interpreter','none')
  orient tall
  wysiwyg
  %end

  if printit
    print('-dpng',['../process/oc4152',int2str(cast),letter(ii)]);
    % print -Poceanpr2
  end

  if letter == 'd'
    fl_ind = [3404:10110];
  elseif letter == 'c'
    fl_ind = [3009:8944];
  elseif letter == 'b'
    fl_ind = [3562:10715];
  elseif letter == ' '
    fl_ind = [3957:11068];
  end

% do stats of 'flight':

end
