function [GT,DATE]=tides(COLAT,LONG,YEAR,MONTH,DAY,HOUR,MINUTE,DUT)
%
%[GT,DATE]=tides(COLAT,LONG,YEAR,MONTH,DAY,HOUR,MINUTE,DUT)
%     or
%[GT,DATE]=tides(COLAT,LONG,YEAR*ones(N,1),ones(N,1)*MONTH,ones(N,1)*DAY,[0:N-1]',zeros(N,1),DUT)
%  for N consecutive hourly samples 
%
%  {COLAT} colatitude in degrees
%  {LONG}  longitude in degrees
%  {DAY} {MONTH} {YEAR} {HOUR} {MINUTE}  vectors of date and local time
%  {DUT} hours to be added to local time to get UT
%  {GT}  vector of tidal predictions in mgal 
%  {DATE} vector of time in hours from 0 hrs UT on the first day
%
% 
%  ************************************************************************
%  TIDES IS A GRAVITY TIDE PREDICTION PROGRAM BASED ON THE XI QINWEN SERIES
%  OF 3070 TIDAL CONSTITUENTS. IT IS DERIVED FROM THE FORTRAN PROGRAM GWAVE
%  TIDES USES ONLY 130 TERMS OF THE 3070 IN THE FULL XI CATALOGUE, SO THE 
%  ACCELERATION OF THE TIDE RAISING POTENTIAL IS ONLY GOOD TO ABOUT a 
%  MICROGAL. NO ALLOWANCE IS MADE FOR THE ELLIPSOIDAL EARTH, OR THE FREE
%  CORE NUTATION, WHICH INCURS FURTHER ERRORS OF ABOUT A MICROGAL. 
%  LASTLY, THERE IS NO OCEAN TIDE LOADING, WHICH AT COASTAL STATIONS MIGHT
%  MEAN A FURTHER ERROR OF TEN MICROGAL OR MORE.
%  THUS, THE PREDICTIONS OF THIS PROGRAM ARE PROBABLY ACCURATE TO ABOUT 
%  TEN MICROGAL OR BETTER AT MOST STATIONS. 

%***********************************************************************************
%
%                              Jim Merriam
%                       University of Saskatchewan
%                          jim.merriam@usask.ca
%                              (Feb, 1997)  
%
%***********************************************************************************
%
DEGRAD=pi/180.0;
disp('                     - TIDES - ')
disp('    COMPUTES TIDAL GRAVITY FROM THE XI QINWEN SERIES ')
disp('    USING 130 TERMS (IN xitides.mat)). NO ALLOWANCE ')
disp('    IS MADE FOR ELLIPSOIDAL FIGURE OF EARTH, THE FREE')
disp('    CORE NUTATION, OR OCEAN TIDE LOADING. PREDICTIONS')
disp('    ARE ACCURATE TO ABOUT TEN MICROGAL AT MOST STATIONS')
disp(' ')

DAY1=DAY(1); % day month and year of first prediction
MTH1=MONTH(1);
YR1=YEAR(1);
TIME=HOUR+MINUTE/60.0 + DUT;  % UT time of observations
disp(' ')
LONG=LONG*DEGRAD;  % CONVERT LAT AND LONG OF BASE TO RADIANS
COLAT=COLAT*DEGRAD;
% CONVERT COLAT TO GEOCENTRIC
COLAT=COLAT+0.003352811*sin(2.0*COLAT);
%  ***********************************************************************
%  CALCULATE TD1, THE TIME IN DAYS FROM 12 HRS JAN 1 2000 TO O HRS ON THE
%  THE FIRST DAY FOR WHICH A PREDICTION IS DESIRED, ALLOWING FOR LEAP YEARS,
%  AND FOUR CENTURY EXCEPTIONS.
%  ***********************************************************************
TD1=round(365.25*YR1)+round(30.6001*(MTH1+1))+round(YR1/400)-round(YR1/100);
  if(MTH1<=2)
  TD1=round(365.25*(YR1-1))+round(30.6001*(MTH1+13))+round(YR1/400)-round(YR1/100);
  end 
TD1=TD1+1720996.50-2451545.00+DAY1;
DOOD=2.6335838E+04;  % DOODSON'S CONSTANT
TJ=TD1/36525.0;  %  CALCULATE TJ, THE TIME IN JULIAN CENTURIES FROM J2000.0

%  CALCULATE ETUT, THE DIFFERENCE BETWEEN ET. AND UT, WHICH MUST BE
%  ADDED TO TJ TO PRODUCE A PREDICTION FOR UT. ETUT IS CALCULATED AS
%  A LINEAR RATE OF 1 SEC/YR FROM JAN 1 1991.
ETUT=(TJ*100.0+9.0)+32.1840+26.0;
TJ=TJ+ETUT/(3600.0*24.0*36525.0);
%  CALCULATE THE RIGHT ASCENSION OF THE FICTITIOUS MEAN SUN IN UT
R=18.69737460+2400.05133690*TJ+0.00002586220*TJ^2 -1.7222*10^(-9)*TJ^3;
R=R*15.0;
%
%  CALCULATE THE FUNDAMENTAL ARGUMENTS OF THE TIDAL CONSTITUENTS, 
%  INCLUDING ADDITIVE TERMS DUE TO PLANETARY INTERACTIONS.
% 
%  XN - THE NEGATIVE OF THE LONGITUDE OF THE MEAN ASCENDING LUNAR NODE.
%  S  - THE GEOMETRIC MEAN LONGITUDE OF THE MOON.
%  P1 - THE MEAN LONGITUDE OF SOLAR PERIGEE.
%  P  - THE MEAN LONGITUDE OF LUNAR PERIGEE.
%  H  - THE GEOMETRIC MEAN LONGITUDE OF THE SUN
%
XN=-125.044520+1934.136260*TJ-0.002070*TJ^2-0.0000020*TJ^3;
XN=XN-0.026660*sin(-XN*DEGRAD) ...
 -0.004330*sin((-XN+272.75-2.30*TJ)*DEGRAD) ...
 -0.000420*sin((188.83-1934.140*(TJ+1))*DEGRAD) ...
 -0.000170*sin((51.19+20.20*(TJ+1))*DEGRAD) ...
 -0.000050*sin((193.44+20.20*(TJ+1))*DEGRAD);
 
S=218.316430+481267.881280*TJ-0.001610*TJ^2+0.0000050*TJ^3 ...
 +0.003960*sin((60.57-132.87*TJ)*DEGRAD)+0.002020*sin(-XN*DEGRAD) ...
 +0.000233*sin((51.19+20.20*(TJ+1))*DEGRAD)+0.000086*sin((84.11+16.2*(TJ+1))*DEGRAD) ...
 +0.000078*sin((174.23-1936.44*(TJ+1))*DEGRAD)+0.000065*sin((304.33-150.6790*(TJ+1))*DEGRAD) ...
 +0.000035*sin((235.96-1034.100*(TJ+1))*DEGRAD)+0.000030*sin((145.27-287.550*(TJ+1))*DEGRAD) ...
 +0.000011*sin((230.91-6.08*(TJ+1))*DEGRAD)+0.000011*sin((211.50+119.0*(TJ+1))*DEGRAD) ...
 +0.000021*sin((76.93+916.720*(TJ+1))*DEGRAD)+0.000008*sin((231.63+4311.540*(TJ+1))*DEGRAD) ...
 +0.000009*sin((275.53-88.34*(TJ+1))*DEGRAD)+0.000015*sin((276.67-723.020*(TJ+1))*DEGRAD) ...
 +0.000003*sin((113.15+619.690*(TJ+1))*DEGRAD)+0.000004*sin((269.69-507.650*(TJ+1))*DEGRAD) ...
 +0.000004*sin((188.12-658.330*(TJ+1))*DEGRAD)+0.000007*sin((46.75-890.440*(TJ+1))*DEGRAD) ...
 +0.000005*sin((10.59+345.700*(TJ+1))*DEGRAD);

P1=282.938350+1.719460*TJ+0.000460*TJ^2+0.0000030*TJ^2;
P=83.353450+4069.013880*TJ-0.010310*TJ^2-0.000010*TJ^3 ...
 -0.000580*sin((71.40+20.20*TJ)*DEGRAD)...
 -0.000580*sin(-XN*DEGRAD)...
 -0.000230*sin((174.23-1936.44*(TJ+1))*DEGRAD)...
 -0.000160*sin((304.33-150.68*(TJ+1))*DEGRAD)...
 -0.000030*sin((193.44-132.87*(TJ+1))*DEGRAD)...
 -0.000030*sin((211.50+119.00*(TJ+1))*DEGRAD);

H=280.466070+36000.769800*TJ+0.000300*TJ^2;
%  ADDITIVE TERMS IN H
H=H+0.001780*sin((251.39+20.20*TJ)*DEGRAD)...
 +(0.00051830-0.00000440*TJ)*sin((207.51+150.27*TJ)*DEGRAD)...
 +0.0000738*sin((31.8+119.0*(TJ+1))*DEGRAD)...
 +0.0000561202*sin((315.6+893.3*(TJ+1))*DEGRAD)...
 -0.002000*sin((67.20+32964.47*TJ)*DEGRAD)...
 -0.001540*sin((16.85-45036.89*TJ)*DEGRAD)...
 +0.001340*sin((81.51+22518.44*TJ)*DEGRAD)...
 -0.004790*(TJ+1)*(1.0+0.003*(TJ+1))*sin((H-P1)*DEGRAD) ...  % SECULAR TERM IN ECCENTRICITY
 +0.001790*sin((S-H)*DEGRAD);  % LUNAR INEQUALITY
%  CALCULATE XNUT, THE NUTATION, AND ADD TO ALL THE FUNDAMENTAL ARGUMENTS.

XNUT=-0.00478*sin(-XN*DEGRAD)-0.00037*sin(2.0*H*DEGRAD);
S=S+XNUT;
H=H+XNUT;
P=P+XNUT;
P1=P1+XNUT;
XN=XN-XNUT;
EQEQ=XNUT*cos(23.43939*DEGRAD);
R=R+EQEQ;
%  **********************************************************************
%  READ IN THE DATA IN XITIDES.MAT THESE ARE;
%
%  ITAU - AN INTEGER MULTIPLE OF MEAN LUNAR TIME.
% NS,NH,NP,NN,NP1 - INTEGER MULTIPLES OF THE ASSOCIATED FUNDAMENTAL
%                   ARGUMENTS.
% N,M - DEGREE AND ORDER OF THE CONSTITUENT.
% OMEGA - THE FREQUENCY OF THE CONSTITUEENT IN DEGREES PER SOLAR HOUR.
% IAMP - AN INTEGER BETWEEN +/-1 AND 908175 SPECIFYING THE RELATIVE
%        AMPLITUDE OF THE CONSTITUENTS. FOR EXAMPLE, THE ENTRY FOR THE
%  M2 WAVE (THE 465 TH ENTRY IN THE FILE) IS
%
%  2  0  0  0  0  0  2 2 28.98410424 908175
%
%  xtides is a reduced xi data set. There are only 130 terms.
%  *************************************************************************
%
load xitides
ITAU=xitides(:,1);
NS=xitides(:,2);
NH=xitides(:,3);
NP=xitides(:,4);
NN=xitides(:,5);
NP1=xitides(:,6);
N=xitides(:,7);
M=xitides(:,8);
OMEGA=xitides(:,9);
IAMP=xitides(:,10);
FREQ=OMEGA*DEGRAD; % frequency in rad/hr

%  CALCULATE THE COEFFICIENT OF THE GEOGRAPHIC DEPENDANCE OF TIDAL GRAVITY
%  FOR EACH DEGREE (N) AND ORDER (M). (SEE NOTE IN SUBROUTINE GFACT)
XL=ones(size(N));

I=find(N==2&M==0);,XL(I)=-1.0*ones(size(I));
I=find(N==2&M==1);,XL(I)=2.0*ones(size(I))/3.0;
I=find(N==2&M==2);,XL(I)=1.0*ones(size(I))/3.0;
I=find(N==3&M==0);,XL(I)=-1.11803*2.0*ones(size(I));
I=find(N==3&M==1);,XL(I)=-.726180*2.0*ones(size(I))/3.0;
I=find(N==3&M==2);,XL(I)=2.598080*ones(size(I))/15.0;
I=find(N==3&M==3);,XL(I)=ones(size(I))/15.0;
%I=find(N==4&M==0);,XL(I)=0.1250000*8.0*ones(size(I));
%I=find(N==4&M==1);,XL(I)=-0.473460*4.0*ones(size(I))/5.0;
%I=find(N==4&M==2);,XL(I)=-0.777780*2.0*ones(size(I))/15.0;
%I=find(N==4&M==3);,XL(I)=3.079200*ones(size(I))/105.0;
%I=find(N==4&M==4);,XL(I)=1.0*ones(size(I))/105.0;


%  **********************  START OF GFACT  ******************************
% ADAPTED FROM THE FORTRAN PROGRAM GFACT
% GFACT RETURNS  DELTAE, AND PNM. PNM IS THE DEGREE N ORDER M
% ASSOCIATED LEGENDRE FUNCTION, AND DELTAE IS THE ELASTIC GRAVIMETRIC 
% FACTOR. NO ADJUSTMENT IS DONE FOR THE FCN, OR ELLIPSOIDAL SHAPE.
% ERRORS FROM THESE OMISSIONS AMOUNT TO ABOUT TWO MICROGAL.
%
% TIDAL GRAVITY IS PNM*DELTAE*[]
% WHERE [] IS THE QUANTITY CALCULATED BY GTIDE AND GWAVE.
%
% ***********************************************************************
%
%                              FEB 1997
%
% ***********************************************************************
%
% THIS PROGRAM WAS WRITTEN BY J. MERRIAM AT THE UNIVERSITY OF SASKATCHEWAN
%
% PROBLEMS WITH THIS PROGRAM SHOULD BE ADDRESSED TO
%
%                       MERRIAM@GEOID.USASK.CA
%
% *************************************************************************
%
  DELTAO=[1.1563 1.1533 1.1589 1.0719 1.035 1.022];
%  DELTAO(1-6) ARE DEGREE 2,0 2,1 2,2 3 4 5 GRAVIMETRIC FACTORS
%  THETA - COLATITUDE 
%
  PNM=ones(size(N));  % create vector for legendre functions
  DELTAE=ones(size(N)); % create vector for gravimetric factors
  I=find(N==2&M==0);  % degree 2 zonal 
  DELTAE(I)=DELTAO(1)*ones(size(I)); 
  PNM(I)=ones(size(I))*(3.0*cos(COLAT)^2-1)/2.0;
  I=find(N==2&M==1);  %  degree 2 tesserals (diurnals)
  DELTAE(I)=ones(size(I))*DELTAO(2);
  PNM(I)=ones(size(I))*3.0*cos(COLAT)*sin(COLAT);
  I=find(N==2&M==2);  %  degree 2 sectorals
  DELTAE(I)=ones(size(I))*(DELTAO(3)-0.00100*sqrt(3.0/4.0)*(7.0*cos(COLAT)^2-1.0));
  PNM(I)=ones(size(I))*3.0*sin(COLAT)^2;
  I=find(N==3); % degree 3 
  DELTAE(I)=ones(size(I))*DELTAO(4);
% degree 4 and 5 are omitted
  I=find(N==3&M==0);,PNM(I)=ones(size(I))*(5.0*cos(COLAT)^3-3*cos(COLAT))/2.0;
  I=find(N==3&M==1);,PNM(I)=ones(size(I))*3.0*(5.0*cos(COLAT)^2-1)*sin(COLAT)/2.0;
  I=find(N==3&M==2);,PNM(I)=ones(size(I))*15.0*cos(COLAT)*sin(COLAT)^2;
  I=find(N==3&M==3);,PNM(I)=ones(size(I))*15.0*sin(COLAT)^3;
  DELTA=PNM.*DELTAE;
%
% ******************  END OF GFACT  ****************************************************************
%
% CALCULATE THE AMPLITUDE OF EACH SUMMED CONSTITUENT (IN GALS), AND
% ITS PHASE AT 0 HRS UT ON THE FIRST DAY FOR WHICH A PREDICTION IS DESIRED.

AMP=DELTA.*IAMP.*XL*DOOD.*N/(6.378137*10^(8));
PHA=NH*H+R*ITAU+(NS-ITAU)*S+NP*P+NN*XN+NP1*P1+M*LONG*180.0/pi;


I=find(N==2&M==1);
PHA(I)=PHA(I)-90.0;    %IF N+M is odd subtract 90 from phase
I=find(N==3&M==0);
PHA(I)=PHA(I)-90.0;
I=find(N==3&M==2);
PHA(I)=PHA(I)-90.0;
PHA=PHA*DEGRAD;
%
% IF THE AMPLITUDE OF THE WAVE IS NEGATIVE, CHANGE ITS
% SIGN AND SUBTRACT PI FROM THE PHASE.
%
I=find(AMP<0);  
AMP(I)=-AMP(I);
PHA(I)=PHA(I)-pi;
PHASE=atan2(sin(PHA),cos(PHA));

DATE=ones(size(HOUR));

  for I=1:length(HOUR)  % cycle through observations
      % calculate TDX the time in days from 12 hrs Jan 1 2000
      % to 0 hrs on the day of the observation allowing for leap
      % years and four century exceptions.
  TDX=round(365.25*YEAR(I))+round(30.6001*(MONTH(I)+1))+round(YEAR(I)/400)-round(YEAR(I)/100);
    if(MONTH(I)<=2)TDX= ...
    round(365.25*(YEAR(I)-1))+round(30.6001*(MONTH(I)+13))+round(YEAR(I)/400)-round(YEAR(I)/100);
    end
  TDX=TDX+1720996.50-2451545.00+DAY(I);
  DATE(I)=(TDX-TD1)*24.0+TIME(I); % time in hours from 0 hrs on the first day
  ARG=(FREQ*DATE(I)+PHASE);    % phase arg of tide
%  SUM OVER ALL CONSTITUENTS 
  GT(I)=sum(-AMP.*cos(ARG)/1000);
  end