function [usum,vsum,ydm]=moody(comp,const,doi,hmh);
disp('WARNING: 1 Hour Added to output time')

% function [usum,vsum,ydm]=moody([comp],[const],[doi],[hmh]);
%
% .... a  routine to estimate & plot multi-constiuent tide from Moody Constants
%    comp       =1 for e/w
%               =2 for n/s component
%               =3 to just pass back usum & vsum
%    const      =constiuent such that [1 2 3 4 5] refers to M2,N2,S2,K1,O1
%      doi      =start time in the form of [yy mth day hhmm] local
%      hmh      =number of hours 
%    usum & vsum=estimate of tide
%     ydm       = yearday in "yrday0_gmt" format
%
% please note hardcodes
% -JManning 


% default to plot the n/s component
if nargin==0;
  comp=2;
end;

if nargin<2;
  const=[1 2 3 4 5];
end;

if nargin<3;
  doi=input('Enter yr mth day hhmm(local) to start estimate: ');
end;
if nargin<4;
  hmh=input('How many hours?');
end;

lc=length(const);

ws=menu('Which starting date?','6 may 0100 local 1995',...
                               '5 may 0100 local 1998');
jdi=julian(doi(2),doi(3),doi(1),doi(4));

if ws==1;
  disp('Currently hardcoded to start 6 may 0100 local 1995');
  %m2--(Vo+u)
  Vou(1)=18.2+314.22+238.09+28.98;
  %n2--(Vo+u)
  Vou(2)=346.4+186.42+172.77+28.44;
  %s2--(Vo+u)
  Vou(3)=0+0+0+30;
  %k1--(Vo+u)
  Vou(4)=15.6+118.28+4.93+15.04;
  %o1--(Vo+u)
  Vou(5)=359.7+195.94+233.16+13.94;
  % calculate hours past 6 May 1995 @0100 local (or 6 May 1995 0500 GMT)
  syd=125;
  if doi(1)==95;
    hr1=(jdi-(syd+1/24))*24;
  else;
    if doi(1)>=96&doi(2)>=4;leap=1;else;leap=0;end;
    hr1=(365-(syd+1/24))*24+(doi(1)-96)*365*24+leap*24+jdi*24;
  end;
elseif ws==2;
  disp('Currently hardcoded to start 5 may 0100 local 1998');
  %m2--(Vo+u)
  Vou(1)=294.0+314.22+262.47+28.98;
  %n2--(Vo+u)
  Vou(2)=83.0+186.42+210.21+28.44;
  %s2--(Vo+u)
  Vou(3)=0+0+0+30;
  %k1--(Vo+u)
  Vou(4)=6.0+118.28+3.94+15.04;
  %o1--(Vo+u)
  Vou(5)=290.5+195.94+258.53+13.94;
  % calculate hours past 5 May 1998 @0100 local (or 5 May 1998 0500 GMT)
  syd=124;
  hr1=(jdi-(syd+1/24))*24;

end;

  hrphi(1)=28.98;
  hrphi(2)=28.44;
  hrphi(3)=30;
  hrphi(4)=15.04;
  hrphi(5)=13.94;


disp('Currently hardcoded to daylight savings time.');
%wS/15=(360/12.42)*60/15=115.94;
wS15=115.94;

disp('Currently hardcoded to station "K" at 15m');
% calculate phi for top 5 constiuents given the phase
ampu=[42.2 139;8.7 109;5.5 173;5.5 156;1.8 110];
ampv=[49.6 30;11 1;7.2 62;4.8 77;1.7 58];

for i=1:lc;
  phiu(i)=ampu(i,2)-rem(Vou(i),360)-wS15+360;
  phiv(i)=ampv(i,2)-rem(Vou(i),360)-wS15+360;
end;


%Filling in the equation we get:
ik=0;
for i=hr1:hr1+hmh;
  for k=1:lc;
    um(k)=ampu(k,1)*cos(pi*(hrphi(k)*i-phiu(k))/180)/100;
    vm(k)=ampv(k,1)*cos(pi*(hrphi(k)*i-phiv(k))/180)/100;
  end;
  ik=ik+1;
  ydm(ik)=jdi+ik/24-1/24+4/24+1/24;
  usum(ik)=sum(um);
  vsum(ik)=sum(vm);
end;

if comp==1;plot(ydm,usum,'linewidth',2);
elseif comp==2;plot(ydm,vsum);end;
%yd2calab(doi(1))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  EXTRA OVERLAYS %%%%%%%%%%%%%%%
% plot Candela's estimate
%hold on;
%candelah=hr1+(365*24)*4+24+(syd+1/24)*24+4:hr1+(365*24)*4+24+(syd+1/24)*24+hmh+4;
%[uc,vc]=detide(-67.567*ones(length(ydm),1),40.20667*ones(length(ydm),1),candelah);
%[uc,vc]=detide(-67.566*ones(length(ydm),1),41.066*ones(length(ydm),1),candelah);
%plot(ydm,vc,'r','linewidth',2);

% plot QUODDY's 
%load /data5/jmanning/model/bank150/may98/r71/ts.dat;
%plot(ts(:,6),ts(:,10),'g');
