
% MATLAB routine to process Albatros SCS data after a cruise

% This version modified Jan 1998 assumes .aco files have been copied
% from monk's ~albscs directory under public_html/cruise_name
% or downloaded through the web & renamed to *.new after ridcomma.f
% gets rid of commas.

% Note: old versions assumed  the "*.lab" or ".asp" files are up on the UNIX

% prompt for cruise code ("cc") if it doesn't exist yet
ecc=exist('cc');
if ecc==0;
  cc=input('enter cruise code (for example: al9403): ','s');
end;


% Prompt user with menu operation
phase=menu('What?','load and plot raw data',...
                   'range&delta check data',...
                   'extract and merge data',...
                   'plot merged data time series',...
                   'plot ship track with labels',...
                   'exit');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if phase==1;
%
%  Have the user click on variables desired
%
  figure(1);
  

% Add some buttons to figure one to allow user options
  byat=uicontrol('style','checkbox','string',' air temp','Position',   [ 10 20 150 20]);
  byst=uicontrol('style','checkbox','string',' sea temp','Position',   [160 20 150 20]);
  bys=uicontrol('style','checkbox','string',' salinity','Position',    [300 20 150 20]);
  byf=uicontrol('style','checkbox','string',' fluorescence','Position',[ 10 50 150 20]);
  byw=uicontrol('style','checkbox','string','wind spd&dir','Position',  [160 50 150 20]);
  byspdd=uicontrol('style','checkbox','string','dslshipspd','Position',[ 10  80 150 20]);
  byspde=uicontrol('style','checkbox','string','edoshipspd','Position',[160  80 150 20]);
  byspdg=uicontrol('style','checkbox','string','gpsshipspd','Position',[300  80 150 20]);

%  allow user to choose a navigation method
   gs=uimenu('label','Lat & Lon (default= GPS)');
   uimenu(gs,'label','GPS   ','Callback','p=1;');
   uimenu(gs,'label','Loran ','Callback','p=2;');
   uimenu(gs,'label','None  ','Callback','p=0;');
%  hardcode above choice to GPS (loran not an option)
   p=1;


%
  disp(' ')
  disp('')
  disp(' ***************************************************')
  disp('     Click  the variables you want to work with. ')
  disp(' ')
  disp('     Then click and press return in this window  when done. ')
  disp(' ***************************************************')
  disp(' ')
  disp('')
  pause;

  set(byat,'Visible','off');
  set(byst,'Visible','off');
  set(bys,'Visible','off');
  set(byf,'Visible','off');
  set(byw,'Visible','off');
  set(byspdd,'Visible','off');
  set(byspde,'Visible','off');
  set(byspdg,'Visible','off');
  set(gs,'visible','off');

  vat=get(byat,'Value');
  vst=get(byst,'value');
  vs=get(bys,'Value');
  vf=get(byf,'value');
  vw=get(byw,'Value');
  vspdd=get(byspdd,'value');
  vspde=get(byspde,'Value');
  vspdg=get(byspdg,'Value');


%  Note: the variables names starting with "t" are time (yeardays)
%        and those starting with "y" are the variables themselves.
%        For example, yat, yst, ys, and yf are air temp, sea temp,
%        salinity, and fluorescence, respectively.
%        We are going to load only those asked for.
disp('Loading data');

  if vat==1;

%  load air temperature 
	eval(['load /data5/jmanning/ship/',cc,'/young.new']);
	tat = young(:,2);
	yat = young(:,5);
        plot(tat-1,yat)
        yd2calab(str2num(cc(3:4)))
        title([cc,' raw air temp']);
        ylabel('degC')
        pt=menu('Hardcopy?','yes','no');
        if pt==1;
          eval(['print -dps ',cc,'airt']);
          eval(['!lpr -h ',cc,'airt.ps']);
        end;
  end;

  
  if vst==1;

%   load seatemp
	eval(['load /data5/jmanning/ship/',cc,'/furuno.new']);
	tst = furuno(:,2);
	yst = furuno(:,5);
        plot(tst-1,yst)
        yd2calab(str2num(cc(3:4)))
        title([cc,' raw sea temp']);
        ylabel('degC')
        pt=menu('Hardcopy?','yes','no');
        if pt==1;
          eval(['print -dps ',cc,'sst']);
          eval(['!lpr -h ',cc,'sst.ps']);
        end;
  end;

  if vs==1;

%     load thermosaliniograph salinity data
	eval(['load /data5/jmanning/ship/',cc,'/tsg.new']);
	ts = tsg(:,2);
	ys = tsg(:,7);
        plot(ts-1,ys)
        yd2calab(str2num(cc(3:4)))
        ylabel('Salinity (PSU)')
        title([cc,' raw salt']);
        pt=menu('Hardcopy?','yes','no');
        if pt==1;
          eval(['print -dps ',cc,'sal']);
          eval(['!lpr -h ',cc,'sal.ps']);
        end;
  end;


  if vf==1;

%  load fluorescence
	eval(['load /data5/jmanning/ship/',cc,'/fluoro.new']);
	tf = fluoro(:,2);
	yf = fluoro(:,5);
        plot(tf-1,yf)
        ylabel('volts')
        yd2calab(str2num(cc(3:4)))
        title([cc,' raw fluorescence']);
        pt=menu('Hardcopy?','yes','no');
        if pt==1;
          eval(['print -dps ',cc,'flur']);
          eval(['!lpr -h ',cc,'flur.ps']);
        end;
  end;


  if vw==1;

%  load windspeed and direction  and convert to U&V
	eval(['load /data5/jmanning/ship/',cc,'/youngw.new']);
	tw = youngw(:,2);
        disp('calculating u & v wind')
	[ywu,ywv] = sd2uv(youngw(:,7),youngw(:,8));
        subplot(2,1,1)
          plot(tw-1,ywu)
          ylabel('Westward Wind (Knots)')
          set(gca,'xticklabels','')
        subplot(2,1,2)
          plot(tw-1,ywv)
          ylabel('Southward Wind (Knots)')
          yd2calab(str2num(cc(3:4)))

        pt=menu('Hardcopy?','yes','no');
        if pt==1;
          eval(['print -dps ',cc,'wind']);
          eval(['!lpr -h ',cc,'wind.ps']);
        end;
        subplot(1,1,1)
  end;


  if vspdd==1;

%  load ship speed (dsl)
	eval(['load /data5/jmanning/ship/',cc,'/dslsort.all']);
	tspdd=dslsort(:,2);
	yspdd = dslsort(:,5);
        plot(tspdd-1,yspdd)
  end;


  if vspdd==1;

%  load ship speed (edo)
	eval(['load /data5/jmanning/ship/',cc,'/edosort.all']);
	tspde=edosort(:,2);
	yspde=edosort(:,5);
        plot(tspde-1,yspde)
  end;

  if vspdg==1;

%  load ship speed (gps)
	eval(['load /data5/jmanning/ship/',cc,'/gpsort.all']);
	tspdg=gpsort(:,2);
	yspdg=gpsort(:,5);
        plot(tspdg-1,yspdg)
  end;


  if p==1;

% Load GPS position
     wc=2;
     if wc==1;
        fil=['/data5/jmanning/ship/',cc,'/gps1.new'];
        fid=fopen(fil);
%        eval(['fid=fopen(''/data5/jmanning/ship/',cc,'/this.dat''']);
%        fid=fopen('/data5/jmanning/ship/al9707/gps1.new');
        tgps=[0];ygps=[0 0];
        
        lin=0;
        while lin~=-1
%         just get once per minute
          for k=1:23;fgetl(fid);end
          lin=fgetl(fid);
          id=find(lin==',');
          tgps=[tgps str2num(lin(id(1)+1:id(2)-1))];
          ygps=[ygps; [str2num(lin(id(5)+1:id(6)-1)) str2num(lin(id(6)+1:id(7)-1))]];
        end;
        tgps=tgps(2:length(tgps));
        ygps=ygps(2:length(ygps),:);

     else;          
        eval(['load /data5/jmanning/ship/',cc,'/gps1.new']);
	tgps = gps1(:,3)+gps1(:,4);
	ygps = gps1(:,6:7);
     end;
        plot(ygps(:,2),ygps(:,1),'r','linewidth',3)
        lllabel
        title(cc)
        pt=menu('Hardcopy?','yes','no');
        if pt==1;
          eval(['print -dps ',cc,'shptrk']);
          eval(['!lpr -h ',cc,'shptrk.ps']);
        end;
  elseif p==2;
 
% Load Loran
        disp('Not an option yet');
  end;

end;



%%%%%%%%%%%%%%%%%%%%%%%%%  PHASE 2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
if phase==2;
%
%                        CLEAN THE DATA
%
% first, set the range and deltas for each variable
% where the first two elements are the range and the third is the delta.
% Note that we check the "time" of each variable as well as the
% the variable itself. Note all times have the same range and delta ("trd").
% We use a function called "rngdelchk.m" to do most of the work.

%
  yatrd=[-50 50 1];
  ystrd=[-20 40 2];
  ysrd=[20 38 1];
  yfrd=[0 200 1];
  ywrd=[-100 100 1];
  yspdd=[0 20 1];
  yspde=[0 20 1];
  yspdg=[0 20 1];
  ygpsrd1=[38 44 .1];
  ygpsrd2=[-71 -65 .1];

%
% keep track of the minimum # of pts
%
  minpts=1000000;
%
  if vat==1;
    disp('cleaning air temp')
    ctat=rngdelchk(tat,[0 366 mean(abs(diff(tat)))+2*std(diff(tat))],0);
    cyat=rngdelchk(yat,[yatrd(1:2) mean(abs(diff(yat)))+2*std(diff(yat))],0);
    if length(ctat)<minpts; minpts=length(ctat); end;
  end;
%
  if vst==1;
    disp('cleaning sea temp')
    ctst=rngdelchk(tst,[0 366 mean(abs(diff(tst)))+2*std(diff(tst))],0);
    cyst=rngdelchk(yst,[ystrd(1:2) mean(abs(diff(yst)))+2*std(diff(yst))],0);
    if length(ctst)<minpts; minpts=length(ctst); end;
  end;
%
  if vs==1;
    disp('cleaning salinity')
    cts=rngdelchk(ts,[0 366 mean(abs(diff(ts)))+2*std(diff(ts))],0);
    cys=rngdelchk(ys,[ysrd(1:2) mean(abs(diff(ys)))+2*std(diff(ys))],0);
    if length(cts)<minpts; minpts=length(cts); end;
  end;
%
  if vf==1;
    disp('cleaning fluorescene')
    ctf=rngdelchk(tf,[0 366 mean(abs(diff(tf)))+2*std(diff(tf))],0);
    cyf=rngdelchk(yf,[yfrd(1:2) mean(abs(diff(yf)))+2*std(diff(yf))],0);
    if length(ctf)<minpts; minpts=length(ctf); end;
  end;
%
  if vw==1;
    disp('cleaning wind')
    ctw=rngdelchk(tw,[0 366 mean(abs(diff(tw)))+2*std(diff(tw))],0);
    cywu=rngdelchk(ywu,[ywrd(1:2) mean(abs(diff(ywu)))+2*std(diff(ywu))],0);
    cywv=rngdelchk(ywv,[ywrd(1:2) mean(abs(diff(ywv)))+2*std(diff(ywv))],0);
    if length(ctw)<minpts; minpts=length(ctw); end;
  end;
%
  if p==1;
    disp('cleaning GPS')
    ctgps=rngdelchk(tgps,[0 366 mean(abs(diff(tgps)))+2*std(diff(tgps))],0);
    cygps1=rngdelchk(ygps(:,1),[ygpsrd1(1:2) mean(abs(diff(ygps(:,1))))+2*std(diff(ygps(:,1)))],0);
    cygps2=rngdelchk(ygps(:,2),[ygpsrd2(1:2) mean(abs(diff(ygps(:,2))))+2*std(diff(ygps(:,2)))],0);
    cygps=[cygps1 cygps2];
    clear cygps1 cygps2;
    if length(ctgps)<minpts; minpts=length(ctgps); end;
  end;
%
end; 
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% PHASE 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if phase==3;
%
%                    EXTRACT and  MERGE DATA
%
%
% Now we begin the actual merge by:
%
%     Fit all data to certain time interval
%
  tint=input(' Enter desired rate of subsample in minutes:');
%
% convert this to fraction of a yearday
%
  fyd=(tint/60)/24;
%
%

%  mth(1)=menu('Start Month?','1','2','3','4','5','6','7','8','9','10','11','12');
  mth(1)=input('start month?');
  day(1)=input('start day?');
  gmt(1)=input('start hour and minute in GMT time (example 1345)');
  mth(2)=input('End  month?');
  day(2)=input('End day?');
  gmt(2)=input('End hour and minute (example 1345)');
  jd(1)=julian(mth(1),day(1),98,gmt(1));
  jd(2)=julian(mth(2),day(2),98,gmt(2));

%
% Now, generate a time vector desired 
% Note: We add one to be consistent with AlbIv convention
%

  timevec=[jd(1):fyd:jd(2)]+1;
%
  y=timevec';
%
% user chooses method of merging
  meth=menu('What method to fit data to regular interval?',...
            'spline',...
            'linear interp');

  disp(' ');
  disp('Hold on ...  fit in progress.')

  hold off;
%
%
% Begin by fitting the lat & lon
%
  if p==1;

%   get rid of non-distinct time records 
    id=find(diff(ctgps)>0);

    if meth==1;

      id2=find(isnan(cygps(id,1))~=1);
      id3=find(isnan(cygps(id,2))~=1);

% calculate a) subsample rate ("ssr") depending on the size of the file and
%           b) 10pt time window ("win") in terms of hours
      ssr=ceil(length(id2)/minpts);
      win=mean(diff(ctgps(id(id2))))*240;

      [stgps,sygps1,npts]=clean(ctgps(id(id2(1:ssr:length(id2)))),cygps(id(id2(1:ssr:length(id2))),1),win*ssr,timevec);
      [stgps,sygps2,npts]=clean(ctgps(id(id3(1:ssr:length(id3)))),cygps(id(id3(1:ssr:length(id3))),2),win*ssr,timevec);

    elseif meth==2;
      idbad=find(diff(ctgps(id))<0);
      if length(idbad)>0;
        disp(['Note: ctgps time vector still has ',numstr(length(idbad)),' time reversals']);
        keyboard;
      end;
      ssr=ceil(length(id)/minpts);
      sygps1=interp1(ctgps(id(1:ssr:length(id))),cygps(id(1:ssr:length(id)),1),timevec);
      sygps2=interp1(ctgps(id(1:ssr:length(id))),cygps(id(1:ssr:length(id)),2),timevec);
    end;

    subplot(1,1,1);
    y=[y sygps1 sygps2];
    disp('Lat/Lon done...')
  end;
%
% Now perform the same operations on all the variables
%
  if vat==1;

%   get rid of non-distinct time records 
    id=find(diff(ctat)>0);

    if meth==1;
      id2=find(diff(ctat(id))~=0);
      ssr=ceil(length(id2)/minpts);
      win=meannnan(diff(ctat(id(id2))))*240;
      [stat,syat,npts]=clean(ctat(id(id2(1:ssr:length(id2)))),cyat(id(id2(1:ssr:length(id2)))),win*ssr,timevec);

    elseif meth==2;
      idbad=find(diff(ctat(id))<0);
      if length(idbad)>0;
        disp('ctat time vector has problems that need correction');
        keyboard;
      end;
      syat=interp1(ctat(id(1:ssr:length(id))),cyat(id(1:ssr:length(id))),timevec);
    end;

    subplot(1,1,1);
    y=[y syat];
    disp('airt done...')
  end;
%
%
  if vst==1;
    id=find(diff(ctst)>0);
    if meth==1;
      id2=find(isnan(cyst(id))~=1);
      win=meannnan(diff(ctst(id(id2))))*240;
      ssr=ceil(length(id2)/minpts);
      [stst,syst,npts]=clean(ctst(id(id2(1:ssr:length(id2)))),cyst(id(id2(1:ssr:length(id2)))),win*ssr,timevec);

    elseif meth==2;
      idbad=find(diff(ctst(id))<0);
      if length(idbad)>0;
        disp('ctst time vector has problems that need correction');
        keyboard;
      end;
      syst=interp1(ctst(id(1:ssr:length(id))),cyst(id(1:ssr:length(id))),timevec);
    end;

    y=[y syst];
    disp('sst done...')
  end;
%
%
  if vs==1;
    id=find(diff(cts)>0);
    if meth==1;

      id2=find(isnan(cts(id))~=1);
      win=meannnan(diff(cts(id(id2))))*240;
      ssr=ceil(length(id2)/minpts);
      [sts,sys]=clean(cts(id(id2(1:ssr:length(id2)))),cys(id(id2(1:ssr:length(id2)))),win*ssr,timevec);

    elseif meth==2;
      idbad=find(diff(cts(id))<0);
      if length(idbad)>0;
        disp('cts time vector has problems that need correction');
        keyboard;
      end;
      sys=interp1(cts(id(1:ssr:length(id))),cys(id(1:ssr:length(id))),timevec);
    end;

    y=[y sys];
    disp('salt done...')
  end;
%
%
  if vf==1;
    id=find(diff(ctf)>0);
    if meth==1;
      id2=find(isnan(cyf(id))~=1);
      win=meannnan(diff(ctf(id(id2))))*240;
      ssr=ceil(length(id2)/minpts);
      [stf,syf]=clean(ctf(id(id2(1:ssr:length(id2)))),cyf(id(id2(1:ssr:length(id2)))),win*ssr,timevec);
      disp('please wait...merging fluorescence complete');

    elseif meth==2;
      idbad=find(diff(ctf(id))<0);
      if length(idbad)>0;
        disp('ctf time vector has problems that need correction');
        keyboard;
      end;
      syf=interp1(ctf(id),cyf(id),timevec);
    end;


    y=[y syf];
    disp('flur done...')

  end;
%
%
  if vw==1;
    id=find(diff(ctw)>0);
    if meth==1;
      id1=find(isnan(ctw(id))~=1);
      id2=find(isnan(cywu(id(id1)))~=1);
      id3=find(isnan(cywv(id(id1)))~=1);
      win=meannnan(diff(ctw(id(id1))))*240;
      ssr=ceil(length(id2)/minpts);

      [stw,sywu]=clean(ctw(id(id1(id2(1:ssr:length(id2))))),cywu(id(id1(id2(1:ssr:length(id2))))),win*ssr,timevec);

      [stw,sywv]=clean(ctw(id(id1(id3(1:ssr:length(id3))))),cywv(id(id1(id3(1:ssr:length(id3))))),win*ssr,timevec);
      disp('please wait...merging wind complete');

    elseif meth==2;
      idbad=find(diff(ctw(id))<0);
      if length(idbad)>0;
        disp('ctw time vector has problems that need correction');
        keyboard;
      end;
      sywu=interp1(ctw(id(1:ssr:length(id))),cywu(id(1:ssr:length(id))),timevec);
      sywv=interp1(ctw(id(1:ssr:length(id))),cywv(id(1:ssr:length(id))),timevec);
    end;


     y=[y sywu sywv];
     disp('wind done...')

  end;
%
%
  if vspdd==1;
    id=find(diff(tspdd)>0);
    if meth==1;
      id2=find(isnan(cyspdd(id))~=1);
      win=meannnan(diff(ctspdd))*240;
      ssr=ceil(length(id2)/minpts);
      [stspdd,syspdd]=clean(ctspdd(id(id2(1:ssr:length(id2)))),cyspdd(id(id2(1:ssr:length(id2)))),win*ssr,timevec);
      disp('please wait...');
    end;

    y=[y syspdd];
  end;
%
  if vspdg==1;
    id=find(diff(tspdg)~=0);
    id2=find(isnan(cyspdg(id))~=1);
    win=meannnan(diff(ctspdg))*240;
    ssr=ceil(length(id2)/minpts);
    [stspdg,syspdg]=clean(ctspdg(id(id2(1:ssr:length(id2)))),cyspdg(id(id2(1:ssr:length(id2)))),win*ssr,timevec);
    disp('please wait...');
    y=[y syspdg];
  end;
%
  if vspde==1;
    id=find(diff(tspde)~=0);
    id2=find(isnan(cyspde(id))~=1);
    ssr=ceil(length(id2)/minpts);
    win=meannnan(diff(ctspde))*240;
    [stspde,syspde]=clean(ctspde(id(id2(1:ssr:length(id2)))),cyspde(id(id2(1:ssr:length(id2)))),win*ssr,timevec);
    disp('please wait...');
    y=[y syspde];
  end;
%
%
% Now, we write out our master file in ascii using a m-file
% that does such a thing for any matrix....
%
%  fploopy(y,'/data5/jmanning/ship/al9801/ship.dat');
   eval(['save /data5/jmanning/ship/',cc,'/ship.dat y -ascii']);
%
% plot
%
   [nr,nc]=size(y);
    plot(y(:,1),y(:,2:nc))
%
end;

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%  PHASE 4 %%%%%%%%%%%%%%%%%%%%%%%%%
%
if phase==4;
%
%                 Plot Cleaned DATA 
%
% find the min & max times
%
  startt=floor(min(y(:,1)));

%HARDCODE!!!!!!!!!!!
%  startt=8;


  stopt=ceil(max(y(:,1)));
%  stopt=19;
  plus1=startt+1;
%
% Generate a stack of plots
%
% First, the temperatures
%
  subplot(1,1,1)
  subplot(5,1,1)
%  [ax,pl,ylab]=plotyy(y(:,1),y(:,4),y(:,1),y(:,5),'SST (degC)',...
%  [0 0 0 min(y(:,1)) max(y(:,1)) min(y(:,4)) max(y(:,4)) min(y(:,5)) max(y(:,5))]);
  [ax,pl,ylab]=plotyy(y(:,1),y(:,4),y(:,1),y(:,5),'',...
  [0 0 0 startt stopt 10 25 8 26]);;
  set(pl(1),'linewidth',1,'color','g')
  set(pl(2),'linewidth',1,'color','r')
  set(ax(1),'ycolor','g')
  set(ax(2),'ycolor','r')
  text(plus1-.5,10,'SST','color','r')
  set(ylab,'color','g');
  ylabel('AirT (degC)');
  title(['GLOBEC cruise ',cc,' June 1998'])
  set(gca,'xticklabels',[' ']);
  set(gca,'position',[.13 .8 .775 .12]);
%
% Next, the salinity
%
  subplot(5,1,2)
  legend off;
  plot(y(:,1),y(:,6),'g');
  ylabel('PSU');
  text(plus1-.8,32.,' Salinity')
  axis([startt stopt 31 32.5]);
  set(gca,'xticklabels',[' ']);
  set(gca,'position',[.13 .66 .775 .12]);
%
% Next, the fluorescsence
%
  subplot(5,1,3)
  plot(y(:,1),y(:,7),'g');
  ylabel('Volts');
  text(plus1,5,' Fluorescence')
  axis([startt stopt 0 60]);
  set(gca,'xticklabels',[' ']);
  set(gca,'position',[.13 .52 .775 .12]);
%
%
% Next, the wind speed
%
  for ii=1:length(y);
    ws(ii)=sqrt(y(ii,8)^2+y(ii,9)^2);
  end;
  subplot(5,1,4)
  plot(y(:,1),ws,'g')
  ylabel('Knots')
  text(plus1,3,' Wind Speed')
  axis([startt stopt 0 30]);
%  set(gca,'xticklabels',[' ']);
  xlabel('Yearday 1998')
  set(gca,'position',[.13 .3 .775 .20]);

%  stkbegin;
%  stkdaxis(startt,stopt,1,94,1);
  set(gcf,'paperposition',[.25 1 8 9.5]);
%
%
%   disp(' You are now in keyboard mode. Type return to continue.')
%
% Next, the ship speeds
%
%  subplot(1,1,1)
%
%  clf;
%  hold
%  plot(stspdd-1,syspdd)
%  ylabel('Knots')
%
%  text(plus1,30,' ship Speed')
%  axis([startt stopt 0 40]);
%  set(gca,'xticklabels',[' ']);
%  set(gca,'position',[.13 .3 .775 .20]);
%
% Finally, the wind wick vector (needs work)
%
% &&&&&&&&&&&&&&&&&&&&&&&
%   see work done for 93 data in gbs93/noaabuoy/scsp.m
%%%%%%%%%%%%%%%%%%%%%%%
%
%  subplot(5,1,5)
%  [u,v]=sd2uv(yws(1:100:length(yws)),sywd(1:100:length(ywd)));
%  y=[[0;u';0] [0;v';0]];
%  stkcombo(94,[startt;tw(1:100:length(stw));stopt],sy,1,1,3,' ')
%  stkdaxis(startt,stopt,1,94,1)
%  text(plus1-1,10,' Wind Vector')
%
% special placement of last plot depends on the stick size at ends
% 
%  set(gca,'position',[-.18 .06 1.09 .2]);
  set(gcf,'paperposition',[.25 1 8 9.5]);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%  PHASE 5 Ship Track w/ Labels %%%%%%%
if phase==5;
 
  load /data5/jmanning/ship/al9808/ship.dat
  y=ship;
%  y=y(1:60:length(y),:);

  % divide the time in quarters
  inc=(max(y(:,1))-min(y(:,1)))/4;
  for i=1:4;
    subplot(2,2,i);
    georbat;
    id=find(y(:,1)>=min(y(:,1))+i*inc-inc&y(:,1)<=min(y(:,1))+i*inc);
    pltdrft(y(id,2),y(id,3),i,28800,98,y(id,1));
    axis([minnnan(y(id,3)) maxnnan(y(id,3)) minnnan(y(id,2)) maxnnan(y(id,2))]);
%    if i<3;hgt=.55;else;hgt=.1;end;
%    if i==1|i==3;lft=.1;else;lft=.55;end;
%    set(gca,'position',[lft hgt .4 .4]);
    lllabel;
    subpos(2,2,i);
    if i==1; title('Cruise Track AL9805 in Four Segments','fontsize',14);end;
    if i==1|i==4;axis([-68.5 -67.25 40.25 41.5]);end;
    if i==2|i==3;lllabel(5,5,40.5,41.25,-68,-67);end;
 
  end;
end;
