function [utide,vtide] = otps_tidecurrent_function(rlon,rlat,rtime)
% this function is to extract barotropic tidal velocity from OTPS database
% TMD 2.05
% Input are 1D arrays of coordinates (rlon and rlat) and time (rtime; 
% matlab time) with the same length
% Output are 1D arrays of tidal velocity (utide and vtide)
%
% Weifeng Gordon Zhang 
% July 6, 2019 (during Cruise TN168)


otpsdir  = '~/research/matlab/OTPS/TMD2.05/DATA/';
tidecon  = {'m2', 's2',  'n2',  'k2',  'k1',  'o1',  'q1',  'm4'};
tideprd  = [12.42,12.00, 12.66, 11.97, 23.93, 25.82, 26.87, 6.21]; % tidal periods in hours
refdate  = [  2019   01   01    0      0     0];   % GMT
disp(['Tidal Start Time = ' datestr(refdate)])
Ntide    = length(tidecon);

Lp    = length(rtime);
utide = NaN([Lp 1]);
vtide = NaN([Lp 1]);

geamp    = zeros([Lp length(tidecon)]);
gepha    = zeros([Lp length(tidecon)]);
guamp    = zeros([Lp length(tidecon)]);
gupha    = zeros([Lp length(tidecon)]);
gvamp    = zeros([Lp length(tidecon)]);
gvpha    = zeros([Lp length(tidecon)]);

% to get the tidal amplitude
for it = 1:length(tidecon)
    % tidal elevation
    tfname = [otpsdir 'h_' tidecon{it} '_tpxo9_atlas_30.nc'];
    if it == 1
        lonz   = ncread(tfname,'lon_z');
        latz   = ncread(tfname,'lat_z');
        tiind  = find(lonz<(min(rlon)+360));
        siind  = tiind(end);
        tiind  = find(lonz>(max(rlon)+360));
        eiind  = tiind(1);       
        tjind  = find(latz < min(rlat));
        sjind  = tjind(end);
        tjind  = find(latz > max(rlat));
        ejind  = tjind(1);
        lonz1  = lonz(siind:eiind);
        latz1  = latz(sjind:ejind);
        lonz2  = repmat(lonz1(:)',[length(latz1) 1]);
        latz2  = repmat(latz1(:),[1 length(lonz1)]);
        ilen   = eiind-siind+1;
        jlen   = ejind-sjind+1;
    end
    hRe  = double(ncread(tfname,'hRe',[sjind siind],[jlen ilen]))./1000;
    hIm  = double(ncread(tfname,'hIm',[sjind siind],[jlen ilen]))./1000;
    eamp = abs(hRe+sqrt(-1).*hIm);
    epha = atan2(-hIm,hRe)./pi*180;
    Fa   = scatteredInterpolant(lonz2(:),latz2(:),eamp(:));
    geamp(:,it) = Fa(rlon(:)+360,rlat(:));
    Fp   = scatteredInterpolant(lonz2(:),latz2(:),epha(:));
    gepha(:,it) = Fp(rlon(:)+360,rlat(:));
  
    % tidal current
    ufname = [otpsdir 'u_' tidecon{it} '_tpxo9_atlas_30.nc'];
    if it == 1 % coordinate
        lonu   = ncread(ufname,'lon_u');
        latu   = ncread(ufname,'lat_u');
        tiind  = find(lonu<(min(rlon)+360));
        usiind = tiind(end);
        tiind  = find(lonu>(max(rlon)+360));
        ueiind = tiind(1);       
        tjind  = find(latu < min(rlat));
        usjind = tjind(end);
        tjind  = find(latu > max(rlat));
        uejind = tjind(1);
        lonu1  = lonu(usiind:ueiind);
        latu1  = latu(usjind:uejind);
        lonu2  = repmat(lonu1(:)',[length(latu1) 1]);
        latu2  = repmat(latu1(:),[1 length(lonu1)]);
        uilen  = ueiind-usiind+1;
        ujlen  = uejind-usjind+1;
        
        % the y-distance for computing velocity from transport
        for ii = 1:length(lonu1)
            for jj = 2:length(latu1)
                tudist(jj,ii) = sw_dist([latu1(jj-1) latu1(jj)],[lonu1(ii) lonu1(ii)],'km').*1000;
            end
            tudist(1,ii)      = sw_dist([latu(usjind-1) latu1(1)],[lonu1(ii) lonu1(ii)],'km').*1000;
            tudist(jj+1,ii)   = sw_dist([latu1(jj) latu(uejind+1)],[lonu1(ii) lonu1(ii)],'km').*1000;                   
        end
        udist = (tudist(1:end-1,:)+tudist(2:end,:))./2;
        
        lonv   = ncread(ufname,'lon_v');
        latv   = ncread(ufname,'lat_v');
        tiind  = find(lonv<(min(rlon)+360));
        vsiind = tiind(end);
        tiind  = find(lonv>(max(rlon)+360));
        veiind = tiind(1);       
        tjind  = find(latv < min(rlat));
        vsjind = tjind(end);
        tjind  = find(latv > max(rlat));
        vejind = tjind(1);
        lonv1  = lonv(vsiind:veiind);
        latv1  = latv(vsjind:vejind);
        lonv2  = repmat(lonv1(:)',[length(latv1) 1]);
        latv2  = repmat(latv1(:),[1 length(lonv1)]);
        vilen  = veiind-vsiind+1;
        vjlen  = vejind-vsjind+1;
        % the x-distance for computing velocity from transport
        for jj = 1:length(latv1)
            for ii = 2:length(lonv1)
                tvdist(jj,ii) = sw_dist([latv1(jj) latv1(jj)],[lonv1(ii-1) lonv1(ii)],'km').*1000;
            end
            tvdist(jj,1)      = sw_dist([latv1(jj) latv1(jj)],[lonv(vsiind-1) lonv1(1)],'km').*1000;
            tvdist(jj,ii+1)   = sw_dist([latv1(jj) latv1(jj)],[lonv1(ii) lonv(veiind+1)],'km').*1000;                   
        end
        vdist = (tvdist(:,1:end-1)+tvdist(:,2:end))./2;
    end
    % u component
    uRe  = double(ncread(ufname,'uRe',[usjind usiind],[ujlen uilen]))./1000;
    uIm  = double(ncread(ufname,'uIm',[usjind usiind],[ujlen uilen]))./1000;
    uamp = abs(uRe+sqrt(-1).*uIm)./udist; % divided by distance to get velocity
    upha = atan2(-uIm,uRe)./pi*180;
    Fa   = scatteredInterpolant(lonu2(:),latu2(:),uamp(:));
    guamp(:,it) = Fa(rlon(:)+360,rlat(:));
    Fp   = scatteredInterpolant(lonu2(:),latu2(:),upha(:));
    gupha(:,it) = Fp(rlon(:)+360,rlat(:));

    % v component
    vRe  = double(ncread(ufname,'vRe',[vsjind vsiind],[vjlen vilen]))./1000;
    vIm  = double(ncread(ufname,'vIm',[vsjind vsiind],[vjlen vilen]))./1000;
    vamp = abs(vRe+sqrt(-1).*vIm)./vdist;  % divided by distance to get velocity
    vpha = atan2(-vIm,vRe)./pi*180;
    Fa   = scatteredInterpolant(lonv2(:),latv2(:),vamp(:));
    gvamp(:,it) = Fa(rlon(:)+360,rlat(:));
    Fp   = scatteredInterpolant(lonv2(:),latv2(:),vpha(:));
    gvpha(:,it) = Fp(rlon(:)+360,rlat(:));
end

% 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.
a = t_getconsts;
for it = 1:length(tideprd)
    iconst(it) = strmatch(upper(tidecon{it}), a.name);
end
reflat     = 55;
[vv,uu,ff] = t_vuf('nodal',datenum(refdate),iconst,reflat);
%vv and uu are returned in cycles, so * by 360 to get degrees or * by 2 pi to get radians
vv  = vv*360;  % convert vv to phase in degrees
uu  = uu*360;  % convert uu to phase in degrees

%adjust phase and ampl of tide to account for user specified time
ei_phase     = ones(Lp,Ntide).*NaN;                       % elevation
ei_amplitude = ones(Lp,Ntide).*NaN;                       % elevation
ui_phase     = ones(Lp,Ntide).*NaN;                       % u vel
ui_amplitude = ones(Lp,Ntide).*NaN;                       %   u vel
vi_phase     = ones(Lp,Ntide).*NaN;                       % v vel
vi_amplitude = ones(Lp,Ntide).*NaN;                       %   v vel

for k=1:length(tideprd)
    %       ei_phase(:,:,k) = angle(ei(:,:,k)).*180./pi - vv(k) - uu(k);   % degrees
    ei_phase(:,k) = gepha(:,k) - vv(k) - uu(k);   % degrees
    ei_amplitude(:,k) = abs(geamp(:,k)) .* ff(k);
    
    %       ui_phase(:,:,k) = angle(ui(:,:,k)).*180./pi - vv(k) - uu(k);   % degrees
    ui_phase(:,k) = gupha(:,k) - vv(k) - uu(k);   % degrees
    ui_amplitude(:,k) = abs(guamp(:,k)) .* ff(k);
    
    %       vi_phase(:,:,k) = angle(vi(:,:,k)).*180./pi - vv(k) - uu(k);   % degrees
    vi_phase(:,k) = gvpha(:,k) - vv(k) - uu(k);   % degrees
    vi_amplitude(:,k) = abs(gvamp(:,k)) .* ff(k);
end
%Tide.period    = periods(indx);
Tide.period    = tideprd;
Tide.component = a.name(iconst,:);
Tide.Ephase    = ei_phase(:,:,:);
Tide.Eamp      = ei_amplitude(:,:,:);

% convert velocity amplitude/phase info to ellipse info
[major,eccentricity,inclination,phase]=ap2ep(ui_amplitude,ui_phase,vi_amplitude,vi_phase);
Tide.Cmax   = major;
Tide.Cmin   = major.*eccentricity;
Tide.Cangle = inclination;
Tide.Cphase = phase;

% compute the tidal velocity following set_tides.F in ROMS
Tide.Cangle = Tide.Cangle.*pi./180;
Tide.Cphase = Tide.Cphase.*pi./180;

cff = 2.0.*pi.*(rtime - datenum(refdate)).*86400;
for itide  = 1:length(tidecon)
    omega  = cff ./ (Tide.period(itide).*3600);
    angle  = Tide.Cangle(:,itide);
    Cangle = cos(angle);
    Sangle = sin(angle);
    phase  = omega - Tide.Cphase(:,itide);
    Cphase = cos(phase);
    Sphase = sin(phase);
    utide(:,itide) =   Tide.Cmax(:,itide) .* Cangle .* Cphase ...
        - Tide.Cmin(:,itide) .* Sangle .* Sphase;
    vtide(:,itide) =   Tide.Cmax(:,itide) .* Sangle .* Cphase ...
        + Tide.Cmin(:,itide) .* Cangle .* Sphase;
end

utide = sum(utide,2)';
vtide = sum(vtide,2)';



