
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%         PHASE 4  ESTIMATING THE CURRENT & PLOT TRAJECTORY
%                  (with wind driven component optional) 
%
%
if phase==5;
%
% a routine to read tidal current speed & direction 
% for each of the 12 hours after maximum flood from NOAA
% tide tables & rotary currents (Table 5 of the annual book)
%
%  first, load tide data from a particular location
%
  load gb10.dat
%
%  second, load flood tide times from Pollack Rip 
% where three columns refer to day, am flood, and pm flood in local time.
%
  load maxf.dat;
%
%  third, figure the tide spd and dir given the time.
%
% Where we use MATLAB "clock" command
% which passes a six element array: [1992 8 4 10 20 32]
% for [yr, mth, day,hr, min, sec] assumes PC is set to GMT
%
  c=fix(clock);
%
% find the local hour
  lhr= (c(4)*100+c(5)/60)-400;
  if lhr<0;
    lhr=lhr+2400;
    c(3)=c(3)-1;
  end;
%
% pick off max flood times for this day
%
  id=find(maxf(:,1)==c(3));
%
% calculate hours after pollack rip's maximum flood
%
  if lhr<maxf(id,2);
%
%  get the previous day's max flood
%
    mf=maxf(id-1,3);
    lhr=2400+lhr
%
elseif lhr>=maxf(id,2)&lhr<maxf(id,3);
%
% pick the am flood
%
    mf=maxf(id,2);
%
% else pick the pm flood
%
  else mf=maxf(id,3);
  end;
%
% calc difference
%
  lhr=(lhr/100);
  dechr=(lhr-floor(lhr))*100/60.;
  lhr=floor(lhr)+dechr;
%
  mf=(mf/100);
  dechr=(mf-floor(mf))*100/60.;
  mf=floor(mf)+dechr;
  diff=lhr-mf;
%
  if diff>12
    disp('something is wrong with input table data');
  end;
%
% round diff off to nearest integer
%
  trunc=floor(diff);
  remain=diff-trunc;
  if remain>=.5;
    diff=trunc+1;
  else
    diff=trunc;
  end;
%
% Note that "DIFF" is our tidal time in hours!
%
% pick spd and direction for this phase of the tide
%
  id2=find(gb10(:,1)==diff);
%
  disp(' Tide at present:');
  disp(['speed:     ',num2str(gb10(id2,2)),' knots'])
  disp(['direction: ',num2str(gb10(id2,3)),' towards'])
%
  [unow,vnow]=sd2uv(gb10(id2,2),gb10(id2,3));

    esp=exist('strack');
    if esp==0;
      load strack.dat;
      slat=strack(:,2);
      slon=strack(:,3);
    end;


%pasthr=input('How far back in time (hours) do you want?');
pasthr=4;

addwind=menu('Do you want the wind driven component?','yes','no');
  if addwind==1;

% Here is where the optional wind driven component is
% added to the tidal component with simple Ekman theory
% 
%       velocity at surface (cm/s)  = 1.5 * wind (m/s)
%
% Read the wind data which is in the same file as Ship Track
% extracting only the last "pasthr" worth of wind

    now=julian(c(:,2),c(:,3),c(:,1),c(:,4)*100+c(:,5))
    goback=now-pasthr/24;
    idw=find(strack(:,1)>goback);
    if length(idw)>0;
      wind=strack(idw,4:5);
  
% Now calculate ekman velocities
      for i=1:length(wind);
        [ue(i),ve(i)]=ekman(wind(i,1),wind(i,2));
      end;
    end;
  end;


%
% Now plot the arrow
% set the origin of the current vector arrow at the
% most recent ship pos'n unless that is not loaded
% then ask the user for option
%
  if esp==0;
    orig=menu('Where do you want the vector origin?',...
            'most recent ship position',...
            'most recent drifter position',...
            'user specified position');
    if orig==1;
      load strack.dat
      slat=strack(:,2);
      slon=strack(:,3);
      tx=slon(length(slon));
      ty=slat(length(slat));
    end;
%
    if orig==2;
      edp=exist('dtrack');
      if edp==0;
        load dtrack.dat
      end;
      tx=dtrack(length(dtrack),1);
      ty=dtrack(length(dtrack),2);
    end;
%
    if orig==3;
      ty=input('enter latitude in decimal degrees: ');
      tx=input('enter longtitude in decimal degrees: ');
    end;
%
  else;
    tx=slon(length(slon));
    ty=slat(length(slat));
  end;
%
  bound=axis;
  if addwind==0;
    arrow(tx,ty,unow*100+ue(length(ue)),vnow*100+ve(length(ve)),100,-1,1,[bound]);
  else;
    arrow(tx,ty,unow*100,vnow*100,100,-1,1,[bound]);
  end;
  hold on;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  PHASE 5b  - Draw recent current path
%
% send u & v tidal estimates from the last pasthr (=4)
% to the subroutine "uv2path.m"
%
  ik=0;
  for j=diff:-1:diff-pasthr;
    if j<=0;
      j=j+12;
    end;
    ik=ik+1;
    u(ik)=.1 - 0.70 * cos(2*pi*j/12) + 0.11 * sin(2*pi*j/12);
    v(ik)=-0.07 + 0.35 * cos(2*pi*j/12) + 0.81 * sin(2*pi*j/12);
  end;


  if addwind==0

%  add on the wind estimates spline fitted to the same hours
%  as the tide was
    sue=spline(strack(idw,1),ue,[now-pasthr/24:1/24:now])  
    sve=spline(strack(idw,1),ve,[now-pasthr/24:1/24:now])  
    u=u+sue;
    v=v+sve;
  end;

  tinc=-1;
  [ptlat,ptlon]=uv2path(tx,ty,u,v,tinc);
  plot(ptlon,ptlat,'g')
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  PHASE 5c  - Draw future tidal path
%
%%%
  ik=0;
  for j=diff:1:diff+pasthr;
%
     if j>12;
      j=j-12;
    end;
    ik=ik+1;
    uf(ik)=.1 - 0.70 * cos(2*pi*j/12) + 0.11 * sin(2*pi*j/12);
    vf(ik)=-0.07 + 0.35 * cos(2*pi*j/12) + 0.81 * sin(2*pi*j/12);
  end;

  if addwind==0

%  add on the wind estimates spline fitted to the same hours
%  as the tide was
    suef=spline(strack(idw,1),ue,[now:1/24:now+pasthr/24])  
    svef=spline(strack(idw,1),ve,[now:1/24:now+pasthr/24])  
    uf=uf+suef;
    vf=vf+svef;
  end;

%
  tinc=1;
  [ptlat,ptlon]=uv2path(tx,ty,uf,vf,tinc);
  plot(ptlon,ptlat,'c')
end;
