function [v,u,f]=vuf(stime,nobsu,konx,xlat);
% [v,u,f]=vuf(stime,nobsu,konx,xlat)
% calculates the astronomical phase argument, nodal modulation 
% amplitude and phase corrections.
%
%  INPUTS:   stime=start time [year,month,day,hour]
%            nobsu=number of points analyzed (odd)
%            konx=number of constituents
%            xlat=latitude of observation (decimal degree)
%
%  OUTPUTS:  v=astronomical phase argument
%            u=nodal modulation phase
%            f=amplitude corrections
%            vux=v+u (total phase correction)

% CALLS: opnvuf.m - loads Doodson numbers, satellite and 
%                   shallow water constituent information
%        astron.m - calculates astronomical arguments
%        julian.m - julian day from gregorian date
%        find_const.m - finds match of constituent name to
%                       list of constituents
%
% DATA:  all called by opnvuf.m
%        doodson.dat - Doodson numbers (6) for 45 major constituents
%        satellite.dat - satellite constituent information 
%        shallow.dat - shallow water constituent information

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1/11/99
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

opnvuf;

slat=sin(pi.*xlat./180);

% determine julian day of mid-point of time series
syr=stime(1);smon=stime(2);sday=stime(3);shr=stime(4);
jdmid=julian(syr,smon,sday,shr)+floor(nobsu./2)./24.0;

% calculate astronomical arguments at mid-point of data time series
[astro,ader]=astron(jdmid);

% calculate correction to get phase relative to Greenwich
vdbl=doodson*astro+semi;
v=rem(vdbl,2);

% satellite amplitude ratio adjustment for latitude made
rr=ee; % no amplitude correction
j=find(ir==1); % latitude correction for diurnal constituents
rr(j)=ee(j).*0.36309.*(1.0-5.0.*slat.*slat)./slat;
j=find(ir==2); % latitude correction for semi-diurnal constituents
rr(j)=ee(j).*2.59808.*slat;

% calculate nodal amplitude and phase corrections
ast3=astro(4:6);
uudbl=del*ast3+ph;
uu=uudbl-round(uudbl);
crr=(rr.*cos(2.*pi.*uu));
srr=(rr.*sin(2.*pi.*uu));
jbase=0;
for k=1:ntidal
 j1=jbase+1;
 jl=jbase+nj(k);
 if (jl>0)
   j=[j1:jl];
   sumc(k)=1.0+sum(crr(j));
   sums(k)=sum(srr(j));
  else
   sumc(k)=1;
   sums(k)=0;
  end
 jbase=jl;
end
sumc=sumc.';sums=sums.';
f=sqrt(sumc.*sumc+sums.*sums);
u=atan2(sums,sumc)./(2.*pi);
vu=u+v;

% compute amplitude and phase corrections for shallow water constituents 
kon=[kon(:,1:4);konsh(:,1:4)];
nsh=length(konsh);
fsh=ones(nsh,1);
ush=zeros(nsh,1);
jr=1:4;
for k=1:nsh;
 for j=1:4;
 jg=jr+(j-1).*5;
 jsh=find_const(konco(k,jg),kon(:,1:4));
   if (length(jsh)>0);
   fsh(k)=fsh(k).*f(jsh).^abs(coef(k,j));
   ush(k)=ush(k)+vu(jsh).*coef(k,j);
   end
 end
end

f=[f;fsh];
u=[u;ush];
v=[v;zeros(nsh,1)];

% get constituents to be used
% CONSTITUENTS ARE NOT IN SAME ORDER SO NEED TO STEP THROUGH TO FIND
% APPROPRIATE ONES
nc=length(konx);
for k=1:nc
 j=find_const(konx(k,:),kon);
 fx(k)=f(j);
 ux(k)=u(j);
 vx(k)=v(j);
 end
f=fx.';v=vx.';u=ux.';
