%                tide.m

%***************************** frequencies *****************************
fm2=1.40519e-4;
fs2=1.45444e-4;
fn2=1.37880e-4;
fk2=1.45842e-4;

fk1=0.72921e-4;
fo1=0.67598e-4;
fp1=0.72523e-4;
fq1=0.64959e-4;

fmsf=0.053234e-4;
fmm=0.026392e-4;
fssa=0.0038921e-4;

% *************************************************************************

sec=(hour)*3600;

D=day+365*(year-1975)+floor((year-1973)/4);
T=(27392.500528+1.0000000356.*D)/36525;

%***************************** Astronomical angles **************************
h0=279.69668+36000.768930485.*T+3.03e-4.*T.^2;
s0=270.434358+481267.88314137.*T-.001133.*T.^2+1.9e-6.*T.^3;
p0=334.329653+4069.0340329575.*T-.010325.*T.^2-1.2e-5.*T.^3;
h0=rem(h0,360);
s0=rem(s0,360);
p0=rem(p0,360);

%**************************   Greenwich phases of constituents **************
xm2=2*h0-2*s0;
xs2=0;
xn2=2*h0-3*s0+p0;
xk2=2*h0;

xk1=h0+90;
xo1=h0-2*s0-90;
xp1=-h0-90;
xq1=h0-3*s0-90+p0;   %  the +p0 was handwritten in by Spiesberger

xmsf=2*s0;
xmm=s0-p0;
xssa=2*h0;

%******************************************************************************

a=input('type 1 to calculate height ');

if a==1

%************************** Gay Head Data *************************
Am2=53;
As2=6.8;
An2=11.2;
Ak1=7.2;
Ao1=0;
Az0=6.2;
Amm=4.2;
Amsf=3.4;

delm2=33.8;
dels2=72.7;
deln2=331.8;
delk1=115.5;
delo1=196.6;
delmm=135;
delmsf=45;


height=Am2*cos(fm2*sec+(xm2-delm2)*pi/180)+ ..
        As2*cos(fs2*sec-(xs2-dels2)*pi/180)+ ..
        An2*cos(fs2*sec-(xn2-deln2)*pi/180)+ ..
        Ao1*cos(fo1*sec-(xo1-delo1)*pi/180)+ ..
        Ak1*cos(fk1*sec-(xk1-delk1)*pi/180);
height=height+..
        Amm*cos(fmm*sec-(xmm-delmm)*pi/180)+ ..
        Ak1*cos(fmsf*sec-(xmsf-delmsf)*pi/180)+ ..
        Az0;
end
