function [z,u,v,depthi,loni,lati,jdm]=tide_track(model_name,tides_to_use,track)
% TIDE_TRACK  interpolates from ADCIRC model onto shiptrack
% to find the tidal constituents along the shiptrack, and then
% calculate the tidal elevations along the shiptrack.
% Usage: [z,u,v,depthi,loni,lati,jdm]=tide_track(model_name,tides_to_use,track)
% 
%  Inputs: model_name = char string, one of: ['ec_1995','ec_2001','enpac_2003']
%          tides_to_use = char string array, name constituents, e.g. ['N2';'K1';'N2';'M2'];
%          track  = 8 column array containing gregorian time, lon and lat
%                  [yyyy(:) mo(:) da(:) hr(:)  mi(:)  sc(:)  lon(:) lat(:)] 
%Example:  
%   track=[ 2000 11 1 0 0 0  -73.8333  40.4667    % ambrose
%           2000 11 2 0 0 0  -74.0167  40.4667];  % sandyhook
%   tides_to_use = ['Z0';'O1';'K1';'N2';'M2';'S2'];
%   model_name='ec_2001';
%   [z,u,v,depthi,loni,lati,jdm]=tide_track(model_name,tides_to_use,track)

% Rich Signell  
% v1.0  Dec  6, 2000
% v1.1  May 13, 2005  modified to read ec2001 data
% v1.2  Nov 11, 2005  modified to read enpac2003 data

% ENPAC example
% track=[ 2002 11 1 0 0 0 -118.371 33.7165     
%         2002 11 1 0 0 0 -118.367 33.7233 ];
% tides_to_use = ['K1';'O1';'P1';'Q1';'M2';'S2';'N2';'K2'];
% model_name='enpac_2003'
% [z,u,v,depthi,loni,lati,jdm]=tide_track(model_name,tides_to_use,track)

[npoints,ncols]=size(track);

g=track(:,[1:6]);  % UTC time in columns 1:6
jdm=datenum(g);    
loni=track(:,7);   % lon in 7th column
lati=track(:,8);   % lat in 8th column

%Load triangle-based tidal constituent data 

%model_name='ec_1995';
%model_name='ec_2001';
%model_name='enpac_2003';

%tides_to_use=['Z0';'O1';'K1';'N2';'M2';'S2'];

switch model_name
case 'ec_1995'
  load adcirc_ec1995_fix.mat tri lon lat u v elev depth periods freq names
%Specify which tidal constituents to use.  Look at "names" in the .mat file to
%see what is available.  Program will halt if these are not found.
% ec95d tides contain ['Z0';'O1';'K1';'N2';'M2';'S2';'M4';'M6'];

case 'ec_2001'
  load adcirc_ec2001v2e_fix.mat tri lon lat u v elev depth periods freq names
%Specify which tidal constituents to use.  Look at "names" in the .mat file to
%see what is available.  Program will halt if these are not found.
% ec2001v2e tides contain ['Z0';'K1';'O1';'Q1';'M2';'S2';'N2';'K2';'M4';'M6'];

case 'enpac_2003'
  load adcirc_enpac2003_tide.mat tri lon lat u v elev depth periods freq names
%Specify which tidal constituents to use.  Look at "names" in the .mat file to
%see what is available.  Program will halt if these are not found.
% enpac2003 tides contain ['Z0';'K1';'O1';'P1';'Q1';'M2';'S2';'N2';'K2';'M4';'M6'];
otherwise
      disp('Error:invalid syntax');  help tide_track; return
end


% make sure that tidal data is double 
tri=double(tri);
lon=double(lon);
lat=double(lat);
u=double(u);
v=double(v);
elev=double(elev);
depth=double(depth);

[Ntide,ncol]=size(tides_to_use);
for k=1:Ntide
  imatch=strmatch(tides_to_use(k,:),names);
  if isempty(imatch),
    disp(['Error: Did not find constituent ' tides_to_use(k,:) ' in .mat file']);
    return
  end
  indx(k)=imatch;
end
periodsi=periods(indx);
%---------------------------------------------------------------------
%  Read in application grid data.
%---------------------------------------------------------------------


% Calculate V,U,F arguments using Rich P's T_Tide.  
% These are the astronomical offsets
% that allow you to go from  amplitude and Greenwich phase to 
% real tidal predictions for a certain time period.  See the
% T_TIDE toolkit and the T_VUF.M program for more info.

a=t_getconsts;

depthi=triterp(tri,lon,lat,depth,loni,lati);

%  Interpolate elevation and velocity from ADCIRC model data using
%  "triterp", linear-interpolation from triangulated grid.
for kk=1:Ntide,
  k=indx(kk);
  disp(['Interpolating ' names(k,:)]);
  format loose

  ei(:,kk)=triterp(tri,lon,lat,elev(:,k),loni,lati);
  ui(:,kk)=triterp(tri,lon,lat,u(:,k),loni,lati);
  vi(:,kk)=triterp(tri,lon,lat,v(:,k),loni,lati);

% find T_TIDE constituent names (a.name) that match 
% selected ADCIRC constituent names (names(tides_to_use))
  iconst(kk)=strmatch(names(k,:), a.name);
end
%***********************************************************************
% This is the call to t_vuf that
% will correct the phase to be at the user specified time.  Also, the amplitude
% is corrected for nodal adjustment.

omega=2*pi*a.freq(iconst);   %tidal frequencies in radians/hour

% Reference latitude for 3rd order satellites (degrees) is 
% set to 55.  You don't need to adjust this to your local latitude
% It could also be set to NaN as in Xtide, with very little effect.  
% See T_VUF for more info.

reflat=55;  
%[vv,uu,ff]=t_vuf(jdm(1),iconst,reflat); %T_TIDE V1.1 syntax
[vv,uu,ff]=t_vuf('nodal',jdm(1),iconst,reflat); %T_TIDE V1.2 syntax

vv=vv*2*pi;  % convert v to phase in radians
uu=uu*2*pi;  % convert u to phase in radians

thours=(jdm-jdm(1))*24;
nt=length(thours);
z=zeros(nt,1);
u=zeros(nt,1);
v=zeros(nt,1);
for i=1:Ntide
    z=z+ff(i).*abs(ei(:,i)).*cos(vv(i)+thours(:)*omega(i)+uu(i)-angle(ei(:,i)));
    u=u+ff(i).*abs(ui(:,i)).*cos(vv(i)+thours(:)*omega(i)+uu(i)-angle(ui(:,i)));
    v=v+ff(i).*abs(vi(:,i)).*cos(vv(i)+thours(:)*omega(i)+uu(i)-angle(vi(:,i)));
end