nutate.m

球面天文学（岁差、章动）和天体力学（行星星历表）的原始计算公式、算法和程序。matlab编译通过。主要包括：公历/儒略历转换为儒略日

function an=nutate(J)
%Nutation in longitude and obliquity
% input: julian date
% return: Nutation in longitude and obliquity (in radians and ")

header;

% The answers are posted here by nutlo():
  
% Each term in the expansion has a trigonometric
% argument given by
%   W = i*MM + j*MS + k*FF + l*DD + m*OM
% where the variables are defined below.
% The nutation in longitude is a sum of terms of the
% form (a + bT) * sin(W). The terms for nutation in obliquity
% are of the form (c + dT) * cos(W).  The coefficients
% are arranged in the tabulation as follows:
% 
% Coefficient:
% i  j  k  l  m      a      b      c     d
% 0, 0, 0, 0, 1, -171996, -1742, 92025, 89,
% The first line of the table, above, is done separately
% since two of the values do not fit into 16 bit integers.
% The values a and c are arc seconds times 10000.  b and d
% are arc seconds per Julian century times 100000.  i through m
% are integers.  See the program for interpretation of MM, MS,
% etc., which are mean orbital elements of the Sun and Moon.
%
% If terms with coefficient less than X are omitted, the peak
% errors will be:
%
%   omit	error,		  omit	error,
%   a <	longitude	  c <	obliquity
% .0005"	.0100"		.0008"	.0094"
% .0046	.0492		.0095	.0481
% .0123	.0880		.0224	.0905
% .0386	.1808		.0895	.1129

nt=[...   %105*9
 0, 0, 0, 0, 2, 2062, 2,-895, 5,
-2, 0, 2, 0, 1, 46, 0,-24, 0,
 2, 0,-2, 0, 0, 11, 0, 0, 0,
-2, 0, 2, 0, 2,-3, 0, 1, 0,
 1,-1, 0,-1, 0,-3, 0, 0, 0,
 0,-2, 2,-2, 1,-2, 0, 1, 0,
 2, 0,-2, 0, 1, 1, 0, 0, 0,
 0, 0, 2,-2, 2,-13187,-16, 5736,-31,
 0, 1, 0, 0, 0, 1426,-34, 54,-1,
 0, 1, 2,-2, 2,-517, 12, 224,-6,
 0,-1, 2,-2, 2, 217,-5,-95, 3,
 0, 0, 2,-2, 1, 129, 1,-70, 0,
 2, 0, 0,-2, 0, 48, 0, 1, 0,
 0, 0, 2,-2, 0,-22, 0, 0, 0,
 0, 2, 0, 0, 0, 17,-1, 0, 0,
 0, 1, 0, 0, 1,-15, 0, 9, 0,
 0, 2, 2,-2, 2,-16, 1, 7, 0,
 0,-1, 0, 0, 1,-12, 0, 6, 0,
-2, 0, 0, 2, 1,-6, 0, 3, 0,
 0,-1, 2,-2, 1,-5, 0, 3, 0,
 2, 0, 0,-2, 1, 4, 0,-2, 0,
 0, 1, 2,-2, 1, 4, 0,-2, 0,
 1, 0, 0,-1, 0,-4, 0, 0, 0,
 2, 1, 0,-2, 0, 1, 0, 0, 0,
 0, 0,-2, 2, 1, 1, 0, 0, 0,
 0, 1,-2, 2, 0,-1, 0, 0, 0,
 0, 1, 0, 0, 2, 1, 0, 0, 0,
-1, 0, 0, 1, 1, 1, 0, 0, 0,
 0, 1, 2,-2, 0,-1, 0, 0, 0,
 0, 0, 2, 0, 2,-2274,-2, 977,-5,
 1, 0, 0, 0, 0, 712, 1,-7, 0,
 0, 0, 2, 0, 1,-386,-4, 200, 0,
 1, 0, 2, 0, 2,-301, 0, 129,-1,
 1, 0, 0,-2, 0,-158, 0,-1, 0,
-1, 0, 2, 0, 2, 123, 0,-53, 0,
 0, 0, 0, 2, 0, 63, 0,-2, 0,
 1, 0, 0, 0, 1, 63, 1,-33, 0,
-1, 0, 0, 0, 1,-58,-1, 32, 0,
-1, 0, 2, 2, 2,-59, 0, 26, 0,
 1, 0, 2, 0, 1,-51, 0, 27, 0,
 0, 0, 2, 2, 2,-38, 0, 16, 0,
 2, 0, 0, 0, 0, 29, 0,-1, 0,
 1, 0, 2,-2, 2, 29, 0,-12, 0,
 2, 0, 2, 0, 2,-31, 0, 13, 0,
 0, 0, 2, 0, 0, 26, 0,-1, 0,
-1, 0, 2, 0, 1, 21, 0,-10, 0,
-1, 0, 0, 2, 1, 16, 0,-8, 0,
 1, 0, 0,-2, 1,-13, 0, 7, 0,
-1, 0, 2, 2, 1,-10, 0, 5, 0,
 1, 1, 0,-2, 0,-7, 0, 0, 0,
 0, 1, 2, 0, 2, 7, 0,-3, 0,
 0,-1, 2, 0, 2,-7, 0, 3, 0,
 1, 0, 2, 2, 2,-8, 0, 3, 0,
 1, 0, 0, 2, 0, 6, 0, 0, 0,
 2, 0, 2,-2, 2, 6, 0,-3, 0,
 0, 0, 0, 2, 1,-6, 0, 3, 0,
 0, 0, 2, 2, 1,-7, 0, 3, 0,
 1, 0, 2,-2, 1, 6, 0,-3, 0,
 0, 0, 0,-2, 1,-5, 0, 3, 0,
 1,-1, 0, 0, 0, 5, 0, 0, 0,
 2, 0, 2, 0, 1,-5, 0, 3, 0, 
 0, 1, 0,-2, 0,-4, 0, 0, 0,
 1, 0,-2, 0, 0, 4, 0, 0, 0,
 0, 0, 0, 1, 0,-4, 0, 0, 0,
 1, 1, 0, 0, 0,-3, 0, 0, 0,
 1, 0, 2, 0, 0, 3, 0, 0, 0,
 1,-1, 2, 0, 2,-3, 0, 1, 0,
-1,-1, 2, 2, 2,-3, 0, 1, 0,
-2, 0, 0, 0, 1,-2, 0, 1, 0,
 3, 0, 2, 0, 2,-3, 0, 1, 0,
 0,-1, 2, 2, 2,-3, 0, 1, 0,
 1, 1, 2, 0, 2, 2, 0,-1, 0,
-1, 0, 2,-2, 1,-2, 0, 1, 0,
 2, 0, 0, 0, 1, 2, 0,-1, 0,
 1, 0, 0, 0, 2,-2, 0, 1, 0,
 3, 0, 0, 0, 0, 2, 0, 0, 0,
 0, 0, 2, 1, 2, 2, 0,-1, 0,
-1, 0, 0, 0, 2, 1, 0,-1, 0,
 1, 0, 0,-4, 0,-1, 0, 0, 0,
-2, 0, 2, 2, 2, 1, 0,-1, 0,
-1, 0, 2, 4, 2,-2, 0, 1, 0,
 2, 0, 0,-4, 0,-1, 0, 0, 0,
 1, 1, 2,-2, 2, 1, 0,-1, 0,
 1, 0, 2, 2, 1,-1, 0, 1, 0,
-2, 0, 2, 4, 2,-1, 0, 1, 0,
-1, 0, 4, 0, 2, 1, 0, 0, 0,
 1,-1, 0,-2, 0, 1, 0, 0, 0,
 2, 0, 2,-2, 1, 1, 0,-1, 0,
 2, 0, 2, 2, 2,-1, 0, 0, 0,
 1, 0, 0, 2, 1,-1, 0, 0, 0,
 0, 0, 4,-2, 2, 1, 0, 0, 0,
 3, 0, 2,-2, 2, 1, 0, 0, 0,
 1, 0, 2,-2, 0,-1, 0, 0, 0,
 0, 1, 2, 0, 1, 1, 0, 0, 0,
-1,-1, 0, 2, 1, 1, 0, 0, 0,
 0, 0,-2, 0, 1,-1, 0, 0, 0,
 0, 0, 2,-1, 2,-1, 0, 0, 0,
 0, 1, 0, 2, 0,-1, 0, 0, 0,
 1, 0,-2,-2, 0,-1, 0, 0, 0,
 0,-1, 2, 0, 1,-1, 0, 0, 0,
 1, 1, 0,-2, 1,-1, 0, 0, 0,
 1, 0,-2, 2, 0,-1, 0, 0, 0,
 2, 0, 0, 2, 0, 1, 0, 0, 0,
 0, 0, 2, 4, 2,-1, 0, 0, 0,
 0, 1, 0, 1, 0, 1, 0, 0, 0];

base=1296000;

jdnut = J;

% Julian centuries from 2000 January 1.5, barycentric dynamical time

T = (J-J2000)/36525.0;
T2 = T * T;
T10 = T / 10.0;

% Fundamental arguments in the FK5 reference system.   

% longitude of the mean ascending node of the lunar orbit on the ecliptic, measured from the mean equinox of date
  
OM = (mod(-6962890.539 * T + 450160.280,base) + (0.008 * T + 7.455) * T2) * STR;

% mean longitude of the Sun minus the mean longitude of the Sun's perigee
  
MS = (mod(129596581.224 * T + 1287099.804,base) - (0.012 * T + 0.577) * T2)  * STR;

% mean longitude of the Moon minus the mean longitude of the Moon's perigee
  
MM = (mod(1717915922.633 * T + 485866.733,base) + (0.064 * T + 31.310) * T2)  * STR;

% mean longitude of the Moon minus the mean longitude of the Moon's node
  
FF = (mod(1739527263.137 * T + 335778.877,base) + (0.011 * T - 13.257) * T2)   * STR;

% mean elongation of the Moon from the Sun.
  
DD = (mod(1602961601.328 * T + 1072261.307,base) + (0.019 * T - 6.891) * T2)   * STR;

% Calculate sin( i*MM ), etc. for needed multiple angles
  
for i=1:3
   ss(1,i)=sin(i*MM);
   cc(1,i)=cos(i*MM);
end
 
for i=1:2
   ss(2,i)=sin(i*MS);
   cc(2,i)=cos(i*MS);
end

for i=1:4
   ss(3,i)=sin(i*FF);
   cc(3,i)=cos(i*FF);
end
  
for i=1:4
   ss(4,i)=sin(i*DD);
   cc(4,i)=cos(i*DD);
end

for i=1:2
   ss(5,i)=sin(i*OM);
   cc(5,i)=cos(i*OM);
end

C = 0.0;
D = 0.0;
p = nt; % point to start of table  
ip=1;
for i=1:105
    ip=1;
% argument of sine and cosine  
	k1 = 0;
	cv = 0.0;
	sv = 0.0;
	for m=1:5
		   j = p(i,ip);
         ip=ip+1;
		if( j~=0 )
			k =abs( j);
         su = ss(m,k); % sin(k*angle)  
			if( j < 0 )
            su = -su;
         end
         cu = cc(m,k);
			if( k1 == 0 )
				 % set first angle  
				sv = su;
				cv = cu;
				k1 = 1;
			else
				 % combine angles  
				sw = su*cv + cu*sv;
				cv = cu*cv - su*sv;
				sv = sw;
			end
    end
 end % end for
 
 % longitude coefficient  
f=p(i,6);
k=p(i,7);
if k~=0
   f=f+T10*k;
end

   %f  = *p++;
%	if( (k = *p++) != 0 )
%		f += T10 * k;

% obliquity coefficient  
g=p(i,8);
k=p(i,9);
if k~=0
   f=f+T10*k;
end

%	g = p++;
%	if( (k = *p++) != 0 )
%		g += T10 * k;

% accumulate the terms  
	C =C+ f * sv;
	D =D+ g * cv;
end % end for i

% first terms, not in table:  
C =C+ (-1742.*T10 - 171996.)*ss(5,1);	% sin(OM)  
D =D+ (   89.*T10 +  92025.)*cc(5,1);	% cos(OM)  

fprintf( 1,'\n\nnutation: in longitude %.3f\", in obliquity %.3f\"\n', .0001*C, D/10000 );
 
% Save answers, expressed in radians  
nutl = 0.0001 * STR * C;
nuto = 0.0001 * D; % in "



% Nutation using nutation in longitude and obliquity from nutlo()and obliquity of the ecliptic from epsiln()
% both calculated for Julian date J.

% * p[] = equatorial rectangular position vector of object for mean ecliptic and equinox of date.
  
eps=epsiln(J); % obliquity of date  

f = eps + nuto;

an=[nutl,f];