%                  tidenew.m
%             need day, year

thetah=280.190;
thetap1=281.221;
thetas=277.026;
thetap=334.384;
thetan=259.156;
thetat=0;

fh=.9856473;
fp1=.0000471;
fs=13.1763968;
fp=0.1114040;
fn=-.0529539;
ft=360;

timezone=5;

name=[' 1  2   3  4  5  6  7  8   9 10 11 12 13  14'
      'z0 mm msf q1 o1 p1 j1 k1 mu2 n2 m2 l2 s2  k2'];
%    s  h  p  p1 90  eps nu  nu'  ind  type  node corr
 nn=[0  0  0  0   0  0     0   0   1    0     1 % z0
     1  0 -1  0   0  0     0   0   2    0     1 % mm
     2 -2  0  0   0  0     0   0   3    0     1 % msf
     -3 1  1  0  -1  0    -2  -1   6    1     3 % q1
     -2 1  0  0  -1  0    -2  -1   7    1     3 % o1
     0 -1  0  0  -1  0     0   0   9    1    10 % p1
     0  1  0  0   1  0    0  -1    10   1     9 % k1
     1  1  -1 0  +1  0    -1   0   11   1     4 %j1
    -4 +4 +0  +0  +0 +2   -2  +0   15   2     6 % mu2
    -3 +2 +1  +0  0  +2   -2  +0   16   2     6 % n2
    -2 +2 +0  +0  +0 +2   -2  +0   17   2     6 % m2
    -1 +2 -1  +0  2  +2   -2  +0   18   2     6 % l2 *not complete
     0  0  0  0   0  0     0   0   19   2    10 % s2
    -0 +2 +0  +0  +0 +0   -0  +2   20   2     6 % k2 *not complete
];
nconst=size(nn);


D=day-1 +365*year+floor((year-1)/4);

phaseh=rem(thetah+D*fh,360);
phasep1=rem(thetap1+D*fp1,360);
phases=rem(thetas+D*fs,360);
phasep=rem(thetap+D*fp,360);
phasen=rem(thetan+D*fn,360);


%...............................nodal phase corrections..............
eps=11.76*sin(phasen*pi/180);
nu=12.87*sin(phasen*pi/180);
nuprime=8.85*sin(phasen*pi/180);

% nodal amplitude corrections

I=(23.155+4.845*cos(phasen*pi/180))*pi/180;
fa(1)=(2/3-(sin(I))^2)/0.5021;
fa(2)=(sin(I))^2/0.1578;
fa(3)=sin(I)*(cos(0.5*I))^2/0.38;
fa(4)=sin(2*I)/0.7214;
fa(5)=sin(I)*(sin(0.5*I))^2/0.0164;
fa(6)=(cos(0.5*I))^4/0.9154;
fa(7)=(sin(I))^2/0.1565;
fa(8)=(sin(0.5*I))^4/0.0017;
fa(9)=(0.2523*(sin(2*I))^2+0.1689*sin(2*I)*cos(nu)+0.0283)^.5/0.5305;
fa(10)=1;
for i=1:nconst;
phase(i)=rem(..
        nn(i,1)*phases+nn(i,2)*phaseh+nn(i,3)*phasep+nn(i,4)*phasep1+..
        nn(i,5)*90 +  nn(i,6)*eps+nn(i,7)*nu+nn(i,8)*nuprime ,360);
if phase(i)<0 
       phase(i)=phase(i)+360;
end
freq(i)=(nn(i,1)*fs+nn(i,2)*fh+nn(i,3)*fp+nn(i,4)*fp1+ft*nn(i,10))/24;
end
name
year
phase


jd=day:1/24:day+nday;   
t=(0:length(jd)-1)+timezone;     %  hours in greenwich time rel. to start
ind=nn(:,9);
u=zeros(t);
for i=1:nconst;
%for i=[1 10 11 13];    %   just for m2
%     I don't know if I have gotten the major and minor axes squared away

%    major axis
u=u+fa(nn(i,11))* ghdeep(ind(i),3)*cos(ghdeep(ind(i),7)*pi/180)* ..
                  cos((freq(i)*t+phase(i)-ghdeep(ind(i),9))*pi/180);
%   minor axis
u=u-fa(nn(i,11))*  ghdeep(ind(i),5)*sin(ghdeep(ind(i),7)*pi/180)* ..
                   sin((freq(i)*t+phase(i)-ghdeep(ind(i),9))*pi/180) ;

end

