function [V ET c]=POTENTIAL(t,lon,lat)
% USAGE: [V ET c]=POTENTIAL(t,lon,lat)
% Obtain tidal potential based on the Appendixes of Munk and Cartwright (1966)
%
% INPUTS:
% t   [Nt x 1] Matlab datenumbers
% lon [1 x 1]  longitude (degrees) 
% lat [1 x 1]  latitude  (degrees) 
%
% OUTPUTS:
% V   [Nt x 1]  Tidal potential  [m^2/s^2]
% ET  [Nt x 1]  Equilibrium tide [m]
% c   [Nt x n]  Complex time dependence of spherical harmonic c=a+ib [m]
%
% Sam Kelly 21 AUG 2015 (smkelly@d.umn.edu)

% Convert matlab datenum into hours since 1 JAN 1900 00:00:00
t=24*(t-datenum(1900,1,1,0,0,0)); 

if size(t,1)==1; t=t.';end
Nt=length(t);

% Calculate Julian centuries since 31 DEC 1899 12:00:00
T=(t+12)/(24*36525);

% Coefficients
ii=complex(0,1);
deg2rad=pi/180;
M_earth=5.9722*10^24; % Earth's mass
a=6378160; % Earth's mean radius
G=6.67300*10^(-11); % Gravitational constant
g=G*M_earth/a^2; % Local gravity
K_moon=[0 .358378 .005946]; % Constants A2
K_sun=[0 .164577 0.000007]; % Constants A2

% Calculate the ephemeris
h_sun=4.8816280+628.3319509*T+0.0000052*T.^2;
p_sun=4.9082295+0.0300053*T+0.0000079*T.^2;
e_sun=0.01675104-0.00004180*T-0.000000126*T.^2;
omega_sun=23.452*deg2rad;

h_moon=4.7200089+8399.7092745*T+0.0000346*T.^2;
p_moon=5.8351526+71.0180412*T-0.0001801*T.^2;
n_moon=4.5236016-33.7571463*T+0.0000363*T.^3;
e_moon=0.054900;
m_moon=0.074804;
i_moon=5.145*deg2rad;

% Sun's polar angle (Z) and east logintude (L)
l_sun=h_sun+2*e_sun.*sin(h_sun-p_sun)+5/4*e_sun.^2.*sin(2*(h_sun-p_sun));
Z_sun=acos(sin(l_sun).*sin(omega_sun));
si_sun=2*atan(sin(1/2*(pi/2+omega_sun))./sin(1/2*(pi/2-omega_sun)).*tan(1/2*(l_sun+Z_sun-pi/2)));
si_sun(si_sun<0)=si_sun(si_sun<0)+2*pi;
si_G=15*deg2rad*(t-12)+h_sun;
si_G(si_G<0)=si_G(si_G<0)+2*pi;
L_sun=si_sun-si_G;

% Ratio of average to present distance
RR_sun=1+e_sun.*cos(h_sun-p_sun)+e_sun.^2.*cos(2*(h_sun-p_sun));

% Moons's polar angle (Z) and east logintude (L)
omega_moon=acos(cos(omega_sun).*cos(i_moon)-sin(omega_sun).*sin(i_moon).*cos(n_moon));
nu=asin(sin(i_moon).*sin(n_moon)./sin(omega_moon));% omega_sun or omega_moon?
sinAOMEGA=sin(n_moon).*sin(omega_sun)./sin(omega_moon);
cosAOMEGA=cos(n_moon).*cos(nu)+cos(omega_sun).*sin(n_moon).*sin(nu);
AOMEGA=2*atan(sinAOMEGA./(1+cosAOMEGA));
sigma_moon=h_moon-n_moon+AOMEGA;
l_moon=sigma_moon+2*e_moon.*sin(h_moon-p_moon)+5/4*e_moon.^2.*sin(2*(h_moon-p_moon))...
    +m_moon.*e_moon.*(15/4+263/16.*m_moon).*sin(h_moon-2*h_sun+p_moon)+m_moon.^2.*(11/8+59/12.*m_moon+75/16.*e_moon.^2./m_moon)*sin(2*(h_moon-h_sun))...
    +17/8*m_moon.^2.*e_moon.*sin(3*h_moon-2*h_sun-p_moon)+77/16*m_moon.^2*e_sun.*sin(2*h_moon-3*h_sun+p_sun);
Z_moon=acos(sin(l_moon).*sin(omega_moon));
si_moon=2*atan(sin(1/2*(pi/2+omega_moon))./sin(1/2*(pi/2-omega_moon)).*tan(1/2*(l_moon+Z_moon-pi/2)));
si_moon(si_moon<0)=si_moon(si_moon<0)+2*pi;
L_moon=nu+si_moon-si_G;

% Ratio of average to present distance
RR_moon=1+(1+1/6*m_moon.^2).^(-1).*(e_moon.*cos(h_moon-p_moon)+e_moon.^2.*cos(2*(h_moon-p_moon))...
        +m_moon.*e_moon.*(15/8+329/64.*m_moon).*cos(h_moon-2*h_sun+p_moon)+m_moon.^2.*(1+19/6.*m_moon+15/4.*e_moon.^2./m_moon)*cos(2*(h_moon-h_sun))...
        +33/16*m_moon.^2.*e_moon.*cos(3*h_moon-2*h_sun-p_moon)+7/2*m_moon.^2*e_sun.*cos(2*h_moon-3*h_sun+p_sun));


% Calculate the potential
% Initiate output array V
V=zeros([Nt 1]);

% Second spherical harmonic n=3 is much smaller 
n=2;

% Initiate timeseries array c=a+bi
c=zeros([Nt n]);

% Loop through tidal species
for m=1:n % change to 0:n to include permanent potential
    
    % Sun's sperhical harmonic
    theta=Z_sun;
    lambda=L_sun;
    P=legendre(n,cos(theta));
    P=P(m+1,:).';
    if mod(m,2)==1; P=-P; end
    Y_sun=(-1)^m*((2*n+1)/(4*pi))^(1/2)*(factorial(n-m)/factorial(n+m))^(1/2)*P.*exp(ii*m*lambda);
    
    % Moon's sperhical harmonic    
    theta=Z_moon;
    lambda=L_moon;
    P=legendre(n,cos(theta));
    P=P(m+1,:).';
    if mod(m,2)==1; P=-P; end
    Y_moon=(-1)^m*((2*n+1)/(4*pi))^(1/2)*(factorial(n-m)/factorial(n+m))^(1/2)*P.*exp(ii*m*lambda);
        
    % Station spherical harmonic
    theta=(90-lat)*deg2rad; % colatitude
    lambda=lon*deg2rad; % longitude. Corrected for consistency with Cartwright 1973 code

    P=legendre(n,cos(theta));
    P=P(m+1);
    if mod(m,2)==1; P=-P; end
    Y_station=(-1)^m*((2*n+1)/(4*pi))^(1/2)*(factorial(n-m)/factorial(n+m))^(1/2)*P*exp(ii*m*lambda);
    
    % Legendre polynomial
    if m==0
        Pn_sun=4*pi/(2*n+1)*real(conj(Y_sun).*Y_station);
        Pn_moon=4*pi/(2*n+1)*real(conj(Y_moon).*Y_station);
    else
        Pn_sun=4*pi/(2*n+1)*2*real(conj(Y_sun).*Y_station);
        Pn_moon=4*pi/(2*n+1)*2*real(conj(Y_moon).*Y_station);
    end
        
    % Calculate the potential
	V=V+g*RR_sun.^(n+1)*K_sun(n).*Pn_sun+g*RR_moon.^(n+1)*K_moon(n).*Pn_moon;   

	% Calculate time coefficients
	c(:,m)=g.*RR_sun.^(n+1).*K_sun(n).*4.*pi./(2*n+1).*2.*Y_sun+g.*RR_moon.^(n+1).*K_moon(n).*4.*pi./(2*n+1).*2.*Y_moon;	
end

% Calculate the equilibrium tide
ET=V/g;

return

