function [astro,ader] = astron(jd)
% [astro,ader] = astron(jd)
% computes the astronomical variables a=[s,h,p,np,pp] (deg) and 
% their time derivatives ader=[ds,dh,dp,dnp,dpp] (deg/day) at the 
% time jd (given in julian days) relative to the reference time 
% origin 12:00 UT Dec 31, 1899.
%
%  INPUT:   jd=julian day
%
%  OUTPUT:  astro=[tau s h p np pp]
%           ader=[dtau ds dh dp dnp dpp]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 12/29/98; 1/11/98
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute number of days from 12:00 UT Dec 31, 1899
d1=jd-julian(1899,12,31,12);
d2=0.0001.*d1;

% compute astronomical constants at time d1
coef=[1.0,d1,d2.*d2,d2.^3];

sc=[  270.434164,13.1763965268,-0.000085,  0.000000039];
hc=[  279.696678, 0.9856473354, 0.00002267,0.000000000];
pc=[  334.329556, 0.1114040803,-0.0007739,-0.00000026];
npc=[-259.183275, 0.0529539222,-0.0001557,-0.00000005];
ppc=[ 281.220833, 0.0000470684, 0.0000339, 0.00000007];

astro=[sc;hc;pc;npc;ppc]*coef'./360.0;
astro=astro-floor(astro);% units cycles so just keep fraction of cycle

% compute lunar time tau, based on fractional part of solar day
hh=jd-floor(jd);
tau=hh+astro(2)-astro(1);
astro=[tau;astro];

% compute rates of change
dcoef=[1.0,2.0e-8.*d1,3.0e-12.*d1.*d1];

dsc=[13.1763965268,-0.000085  , 0.000000039];
dhc=[ 0.9856473354, 0.00002267, 0.00       ];
dpc=[ 0.1114040803,-0.0007739 ,-0.00000026 ];
dnpc=[0.0529539222,-0.0001557 ,-0.00000005 ];
dppc=[0.0000470684, 0.0000339 , 0.00000007 ];

ader=[dsc;dhc;dpc;dnpc;dppc]*dcoef'./(360.0./365.0);
dtau=365.0+ader(2)-ader(1);
ader=[dtau;ader];

