text: /* nutation (in IAU (1980) expression) and abberation; stern
   * on an HP PA processor, this reproduces the Almanac nutation values
   * (given to 0.001") EXACTLY over 750 days (1995 and 1996)
   */
  #include <stdio.h>
  #include <math.h>
  
  #include "astro.h"
  
  #define NUT_SCALE	1e4
  #define NUT_SERIES	106
  #define NUT_MAXMUL	4
  #define SECPERCIRC	(3600.*360.)
  
  /* Delaunay arguments, in arc seconds; they differ slightly from ELP82B */
  static double delaunay[5][4] = {
      {485866.733,  1717915922.633, 31.310,  0.064}, /* M', moon mean anom */
      {1287099.804, 129596581.224,  -0.577, -0.012}, /* M, sun mean anom */
      {335778.877,  1739527263.137, -13.257, 0.011}, /* F, moon arg lat */
      {1072261.307, 1602961601.328, -6.891,  0.019}, /* D, elong moon sun */
      {450160.280,  -6962890.539,   7.455,   0.008}, /* Om, moon l asc node */
  };
  
  /* multipliers for Delaunay arguments */
  static short multarg[NUT_SERIES][5] = {
  	/* bounds:  -2..3, -2..2, -2/0/2/4, -4..4, 0..2 */
      {0, 0, 0, 0, 1},
      {0, 0, 0, 0, 2},
      {-2, 0, 2, 0, 1},
      {2, 0, -2, 0, 0},
      {-2, 0, 2, 0, 2},
      {1, -1, 0, -1, 0},
      {0, -2, 2, -2, 1},
      {2, 0, -2, 0, 1},
      {0, 0, 2, -2, 2},
      {0, 1, 0, 0, 0},
      {0, 1, 2, -2, 2},
      {0, -1, 2, -2, 2},
      {0, 0, 2, -2, 1},
      {2, 0, 0, -2, 0},
      {0, 0, 2, -2, 0},
      {0, 2, 0, 0, 0},
      {0, 1, 0, 0, 1},
      {0, 2, 2, -2, 2},
      {0, -1, 0, 0, 1},
      {-2, 0, 0, 2, 1},
      {0, -1, 2, -2, 1},
      {2, 0, 0, -2, 1},
      {0, 1, 2, -2, 1},
      {1, 0, 0, -1, 0},
      {2, 1, 0, -2, 0},
      {0, 0, -2, 2, 1},
      {0, 1, -2, 2, 0},
      {0, 1, 0, 0, 2},
      {-1, 0, 0, 1, 1},
      {0, 1, 2, -2, 0},
      {0, 0, 2, 0, 2},
      {1, 0, 0, 0, 0},
      {0, 0, 2, 0, 1},
      {1, 0, 2, 0, 2},
      {1, 0, 0, -2, 0},
      {-1, 0, 2, 0, 2},
      {0, 0, 0, 2, 0},
      {1, 0, 0, 0, 1},
      {-1, 0, 0, 0, 1},
      {-1, 0, 2, 2, 2},
      {1, 0, 2, 0, 1},
      {0, 0, 2, 2, 2},
      {2, 0, 0, 0, 0},
      {1, 0, 2, -2, 2},
      {2, 0, 2, 0, 2},
      {0, 0, 2, 0, 0},
      {-1, 0, 2, 0, 1},
      {-1, 0, 0, 2, 1},
      {1, 0, 0, -2, 1},
      {-1, 0, 2, 2, 1},
      {1, 1, 0, -2, 0},
      {0, 1, 2, 0, 2},
      {0, -1, 2, 0, 2},
      {1, 0, 2, 2, 2},
      {1, 0, 0, 2, 0},
      {2, 0, 2, -2, 2},
      {0, 0, 0, 2, 1},
      {0, 0, 2, 2, 1},
      {1, 0, 2, -2, 1},
      {0, 0, 0, -2, 1},
      {1, -1, 0, 0, 0},
      {2, 0, 2, 0, 1},
      {0, 1, 0, -2, 0},
      {1, 0, -2, 0, 0},
      {0, 0, 0, 1, 0},
      {1, 1, 0, 0, 0},
      {1, 0, 2, 0, 0},
      {1, -1, 2, 0, 2},
      {-1, -1, 2, 2, 2},
      {-2, 0, 0, 0, 1},
      {3, 0, 2, 0, 2},
      {0, -1, 2, 2, 2},
      {1, 1, 2, 0, 2},
      {-1, 0, 2, -2, 1},
      {2, 0, 0, 0, 1},
      {1, 0, 0, 0, 2},
      {3, 0, 0, 0, 0},
      {0, 0, 2, 1, 2},
      {-1, 0, 0, 0, 2},
      {1, 0, 0, -4, 0},
      {-2, 0, 2, 2, 2},
      {-1, 0, 2, 4, 2},
      {2, 0, 0, -4, 0},
      {1, 1, 2, -2, 2},
      {1, 0, 2, 2, 1},
      {-2, 0, 2, 4, 2},
      {-1, 0, 4, 0, 2},
      {1, -1, 0, -2, 0},
      {2, 0, 2, -2, 1},
      {2, 0, 2, 2, 2},
      {1, 0, 0, 2, 1},
      {0, 0, 4, -2, 2},
      {3, 0, 2, -2, 2},
      {1, 0, 2, -2, 0},
      {0, 1, 2, 0, 1},
      {-1, -1, 0, 2, 1},
      {0, 0, -2, 0, 1},
      {0, 0, 2, -1, 2},
      {0, 1, 0, 2, 0},
      {1, 0, -2, -2, 0},
      {0, -1, 2, 0, 1},
      {1, 1, 0, -2, 1},
      {1, 0, -2, 2, 0},
      {2, 0, 0, 2, 0},
      {0, 0, 2, 4, 2},
      {0, 1, 0, 1, 0}
  };
  
  /* amplitudes which  have secular terms; in 1/NUT_SCALE arc seconds
   * {index, constant dPSI, T/10 in dPSI, constant in dEPS, T/10 in dEPS}
   */
  static long ampsecul[][5] = {
      {0  ,-171996 ,-1742 ,92025 ,89},
      {1  ,2062    ,2     ,-895  ,5},
      {8  ,-13187  ,-16   ,5736  ,-31},
      {9  ,1426    ,-34   ,54    ,-1},
      {10 ,-517    ,12    ,224   ,-6},
      {11 ,217     ,-5    ,-95   ,3},
      {12 ,129     ,1     ,-70   ,0},
      {15 ,17      ,-1    ,0     ,0},
      {17 ,-16     ,1     ,7     ,0},
      {30 ,-2274   ,-2    ,977   ,-5},
      {31 ,712     ,1     ,-7    ,0},
      {32 ,-386    ,-4    ,200   ,0},
      {33 ,-301    ,0     ,129   ,-1},
      {37 ,63      ,1     ,-33   ,0},
      {38 ,-58     ,-1    ,32    ,0},
      /* termination */  { -1, }
  };
  
  /* amplitudes which only have constant terms; same unit as above
   * {dPSI, dEPS}
   * indexes which are already in ampsecul[][] are zeroed
   */
  static short ampconst[NUT_SERIES][2] = {
      {0,0},
      {0,0},
      {46,-24},
      {11,0},
      {-3,1},
      {-3,0},
      {-2,1},
      {1,0},
      {0,0},
      {0,0},
      {0,0},
      {0,0},
      {0,0},
      {48,1},
      {-22,0},
      {0,0},
      {-15,9},
      {0,0},
      {-12,6},
      {-6,3},
      {-5,3},
      {4,-2},
      {4,-2},
      {-4,0},
      {1,0},
      {1,0},
      {-1,0},
      {1,0},
      {1,0},
      {-1,0},
      {0,0},
      {0,0},
      {0,0},
      {0,0},
      {-158,-1},
      {123,-53},
      {63,-2},
      {0,0},
      {0,0},
      {-59,26},
      {-51,27},
      {-38,16},
      {29,-1},
      {29,-12},
      {-31,13},
      {26,-1},
      {21,-10},
      {16,-8},
      {-13,7},
      {-10,5},
      {-7,0},
      {7,-3},
      {-7,3},
      {-8,3},
      {6,0},
      {6,-3},
      {-6,3},
      {-7,3},
      {6,-3},
      {-5,3},
      {5,0},
      {-5,3},
      {-4,0},
      {4,0},
      {-4,0},
      {-3,0},
      {3,0},
      {-3,1},
      {-3,1},
      {-2,1},
      {-3,1},
      {-3,1},
      {2,-1},
      {-2,1},
      {2,-1},
      {-2,1},
      {2,0},
      {2,-1},
      {1,-1},
      {-1,0},
      {1,-1},
      {-2,1},
      {-1,0},
      {1,-1},
      {-1,1},
      {-1,1},
      {1,0},
      {1,0},
      {1,-1},
      {-1,0},
      {-1,0},
      {1,0},
      {1,0},
      {-1,0},
      {1,0},
      {1,0},
      {-1,0},
      {-1,0},
      {-1,0},
      {-1,0},
      {-1,0},
      {-1,0},
      {-1,0},
      {1,0},
      {-1,0},
      {1,0}
  };
  
  /* given the modified JD, mj, find the nutation in obliquity, *deps, and
   * the nutation in longitude, *dpsi, each in radians.
   */
  void
  nutation (
  double mj,
  double *deps,	/* on input:  precision parameter in arc seconds */
  double *dpsi)
  {
  	static double lastmj = -10000, lastdeps, lastdpsi;
  	double T, T2, T3, T10;			/* jul cent since J2000 */
  	double prec;				/* series precis in arc sec */
  	int i, isecul;				/* index in term table */
  	static double delcache[5][2*NUT_MAXMUL+1];
  			/* cache for multiples of delaunay args
  			 * [M',M,F,D,Om][-min*x, .. , 0, .., max*x]
  			 * make static to have unfilled fields cleared on init
  			 */
  
  	if (mj == lastmj) {
  	    *deps = lastdeps;
  	    *dpsi = lastdpsi;
  	    return;
  	}
  
  	prec = 0.0;
  
  #if 0	/* this is if deps should contain a precision value */
  	prec =* deps;
  	if (prec < 0.0 || prec > 1.0)	/* accept only sane value */
  		prec = 1.0;
  #endif
  
  	/* augment for abundance of small terms */
  	prec *= NUT_SCALE/10;
  
  	T = (mj - J2000)/36525.;
  	T2 = T * T;
  	T3 = T2 * T;
  	T10 = T/10.;
  
  	/* calculate delaunay args and place in cache */
  	for (i = 0; i < 5; ++i) {
  	    double x;
  	    short j;
  
  	    x = delaunay[i][0] +
  		delaunay[i][1] * T +
  		delaunay[i][2] * T2 +
  		delaunay[i][3] * T3;
  
  	    /* convert to radians */
  	    x /= SECPERCIRC;
  	    x -= floor(x);
  	    x *= 2.*PI;
  
  	    /* fill cache table */
  	    for (j = 0; j <= 2*NUT_MAXMUL; ++j)
  		delcache[i][j] = (j - NUT_MAXMUL) * x;
  	}
  
  	/* find dpsi and deps */
  	lastdpsi = lastdeps = 0.;
  	for (i = isecul = 0; i < NUT_SERIES ; ++i) {
  	    double arg = 0., ampsin, ampcos;
  	    short j;
  
  	    if (ampconst[i][0] || ampconst[i][1]) {
  		/* take non-secular terms from simple array */
  		ampsin = ampconst[i][0];
  		ampcos = ampconst[i][1];
  	    } else {
  		/* secular terms from different array */
  		ampsin = ampsecul[isecul][1] + ampsecul[isecul][2] * T10;
  		ampcos = ampsecul[isecul][3] + ampsecul[isecul][4] * T10;
  		++isecul;
  	    }
  
  	    for (j = 0; j < 5; ++j)
  		arg += delcache[j][NUT_MAXMUL + multarg[i][j]];
  
  	    if (fabs(ampsin) >= prec)
  		lastdpsi += ampsin * sin(arg);
  
  	    if (fabs(ampcos) >= prec)
  		lastdeps += ampcos * cos(arg);
  
  	}
  
  	/* convert to radians.
  	 */
  	lastdpsi = degrad(lastdpsi/3600./NUT_SCALE);
  	lastdeps = degrad(lastdeps/3600./NUT_SCALE);
  
  	lastmj = mj;
  	*deps = lastdeps;
  	*dpsi = lastdpsi;
  }
  
  /* given the modified JD, mj, correct, IN PLACE, the right ascension *ra
   * and declination *dec (both in radians) for nutation.
   */
  void
  nut_eq (double mj, double *ra, double *dec)
  {
  	static double lastmj = -10000;
  	static double a[3][3];		/* rotation matrix */
  	double xold, yold, zold, x, y, z;
  
  	if (mj != lastmj) {
  	    double epsilon, dpsi, deps;
  	    double se, ce, sp, cp, sede, cede;
  
  	    obliquity(mj, &epsilon);
  	    nutation(mj, &deps, &dpsi);
  
  	    /* the rotation matrix a applies the nutation correction to
  	     * a vector of equatoreal coordinates Xeq to Xeq' by 3 subsequent
  	     * rotations:  R1 - from equatoreal to ecliptic system by
  	     * rotation of angle epsilon about x, R2 - rotate ecliptic
  	     * system by -dpsi about its z, R3 - from ecliptic to equatoreal
  	     * by rotation of angle -(epsilon + deps)
  	     *
  	     *	Xeq' = A * Xeq = R3 * R2 * R1 * Xeq
  	     * 
  	     *		[ 1       0          0    ]
  	     * R1 =	[ 0   cos(eps)   sin(eps) ]
  	     *		[ 0  - sin(eps)  cos(eps) ]
  	     * 
  	     *		[ cos(dpsi)  - sin(dpsi)  0 ]
  	     * R2 =	[ sin(dpsi)   cos(dpsi)   0 ]
  	     *		[      0           0      1 ]
  	     * 
  	     *		[ 1         0                 0         ]
  	     * R3 =	[ 0  cos(eps + deps)  - sin(eps + deps) ]
  	     *		[ 0  sin(eps + deps)   cos(eps + deps)  ]
  	     * 
  	     * for efficiency, here is a explicitely:
  	     */
  	    
  	    se = sin(epsilon);
  	    ce = cos(epsilon);
  	    sp = sin(dpsi);
  	    cp = cos(dpsi);
  	    sede = sin(epsilon + deps);
  	    cede = cos(epsilon + deps);
  
  	    a[0][0] = cp;
  	    a[0][1] = -sp*ce;
  	    a[0][2] = -sp*se;
  
  	    a[1][0] = cede*sp;
  	    a[1][1] = cede*cp*ce+sede*se;
  	    a[1][2] = cede*cp*se-sede*ce;
  
  	    a[2][0] = sede*sp;
  	    a[2][1] = sede*cp*ce-cede*se;
  	    a[2][2] = sede*cp*se+cede*ce;
  
  	    lastmj = mj;
  	}
  
  	sphcart(*ra, *dec, 1.0, &xold, &yold, &zold);
  	x = a[0][0] * xold + a[0][1] * yold + a[0][2] * zold;
  	y = a[1][0] * xold + a[1][1] * yold + a[1][2] * zold;
  	z = a[2][0] * xold + a[2][1] * yold + a[2][2] * zold;
  	cartsph(x, y, z, ra, dec, &zold);	/* radius should be 1.0 */
  	if (*ra < 0.) *ra += 2.*PI;		/* make positive for display */
  }
  
  