Return to nutate.c CVS log

Up to [Research Unix] / researchv10no / cmd / sky

researchv10 Norman

#include "sky.h"

extern struct nuttab
{
	float f[2];
	char c[5];
} nuttab[];

nutate()
{

/*
 *	uses radian, radsec, eday
 *	sets psi, eps, dpsi, deps, obliq, gst, tobliq
 */

	double msun, mnom, noded, dmoon;
	struct nuttab *pp;
	double arg;
	double Eday;
	double capT, capT2, capT3;

/*
 *	nutation of the equinoxes is a irregular wobble of the
 *	pole of the earths rotation whose magnitude is about
 *	9 seconds of arc and whose period is about 18.6 years.
 *
 *	it depends upon the pull of the sun and moon on the
 *	equatorial bulge of the earth.
 *
 *	psi and eps are the two angles which specify the
 *	true pole with respect to the mean pole.
 *
 *	all coeffieients are from Supplement to A.E. 1984
 */

	pp = &nuttab[0];
	Eday = eday - 36525;
	capT = capt-1;
	capT2 = capT*capT;
	capT3 = capT2*capT;

	mnom = 134.96298139 + 13.06499294724*Eday + 8.6972e-3*capT2
		 + 17.8e-6*capT3;
	msun = 357.52772333 + 0.98560028309*Eday - 0.16028e-3*capT2
		 - 3.33e-6*capT3;
	noded = 93.27191028 + 13.22935024060*Eday - 3.6825e-3*capT2
		 + 3.05e-6*capT3;
	dmoon = 297.85036306 + 12.19074911650*Eday - 1.91417e-3*capT2
		 + 5.28e-6*capT3;
	node = 125.04452222 - .05295376484*Eday + 2.07083e-3*capT2
		 + 2.22e-6*capT3;
	mnom *= radian;
	msun *= radian;
	noded *= radian;
	dmoon *= radian;
	node *= radian;

/*
 *	long period terms
 */


	psi = -(17.1996+.01742*capT)*sin(node);
	eps = (9.2025+0.00089*capT)*cos(node);
	for(;;){
		if(pp->f[0]==0.){
			pp++;
			break;
		}
		arg = pp->c[0]*mnom + pp->c[1]*msun + pp->c[2]*noded
			+ pp->c[3]*dmoon + pp->c[4]*node;
		psi += pp->f[0]*sin(arg);
		eps += pp->f[1]*cos(arg);
		pp++;
	}

/*
 *	short period terms
 */

	dpsi = 0.;
	deps = 0.;
	if((flags&APPARENT)==0){
		for(;;){
			if(pp->f[0]==0.){
				pp++;
				break;
			}
			arg = pp->c[0]*mnom + pp->c[1]*msun + pp->c[2]*noded
				+ pp->c[3]*dmoon + pp->c[4]*node;
			dpsi += pp->f[0]*sin(arg);
			deps += pp->f[1]*cos(arg);
			pp++;
		}
	}


	psi = (psi+dpsi)*radsec;
	eps = (eps+deps)*radsec;
	dpsi *= radsec;
	deps *= radsec;

	obliq = 23.43929111 - .013004167*capT - 0.164e-6*capT2
		 + 0.5036e-6*capT3;
	obliq *= radian;
	tobliq = obliq + eps;

	gst = 100.460618375 + 360.98564736629*Eday + 0.3882667e-3*capT2
		 - 0.025e-6*capT3;
	gst -= 180.;
	gst = fmod(gst, 360.);
	if(gst < 0.)
		gst += 360.;
	gst *= radian;
	gst += psi*cos(tobliq);

/*
 * print true obliquity, mean obliquity, mean GST, and true GST.
 */

	if(flags&GOOBIE){
		printf("Mean G. Sid. Time: ");
		prhms(gst - psi*cos(obliq), 4);
		printf("\n");

		printf("App. G. Sid. Time: ");
		prhms(gst,4);
		printf("\n");

		printf("Goobie: %6.3f %6.3f %6.3f %6.3f\n", psi/radsec,
			eps/radsec, dpsi/radsec, deps/radsec);

		printf("Obliquity: ");
		prdms(tobliq,3);
		printf("\n");
	}

}

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.