/*  * This file is part of the UCB release of Plan 9. It is subject to the license * terms in the LICENSE file found in the top-level directory of this * distribution and at http://akaros.cs.berkeley.edu/files/Plan9License. No * part of the UCB release of Plan 9, including this file, may be copied, * modified, propagated, or distributed except according to the terms contained * in the LICENSE file. */ #include "astro.h" voidmerc(void){	double pturbl, pturbr;	double lograd;	double dele, enom, vnom, nd, sl;	double q0, v0, t0, j0 , s0;	double lsun, elong, ci, dlong; 	ecc = .20561421 + .00002046*capt - 0.03e-6*capt2;	incl = 7.0028806 + .0018608*capt - 18.3e-6*capt2;	node = 47.145944 + 1.185208*capt + .0001739*capt2;	argp = 75.899697 + 1.555490*capt + .0002947*capt2;	mrad = .3870986;	anom = 102.279381 + 4.0923344364*eday + 6.7e-6*capt2;	motion = 4.0923770233; 	q0 = 102.28  + 4.092334429*eday;	v0 = 212.536 + 1.602126105*eday;	t0 = -1.45  + .985604737*eday;	j0 = 225.36 + .083086735*eday;	s0 = 175.68 + .033455441*eday; 	q0 *= radian;	v0 *= radian;	t0 *= radian;	j0 *= radian;	s0 *= radian; 	incl *= radian;	node *= radian;	argp *= radian;	anom = fmod(anom, 360.)*radian;  	enom = anom + ecc*sin(anom);	do {		dele = (anom - enom + ecc * sin(enom)) /			(1. - ecc*cos(enom));		enom += dele;	} while(fabs(dele) > converge);	vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),			cos(enom/2.));	rad = mrad*(1. - ecc*cos(enom)); 	icosadd(mercfp, merccp);	pturbl =  cosadd(2, q0, -v0);	pturbl += cosadd(2, q0, -t0);	pturbl += cosadd(2, q0, -j0);	pturbl += cosadd(2, q0, -s0); 	pturbr =  cosadd(2, q0, -v0);	pturbr += cosadd(2, q0, -t0);	pturbr += cosadd(2, q0, -j0); /* *	reduce to the ecliptic */ 	lambda = vnom + argp + pturbl*radsec;	nd = lambda - node;	lambda = node + atan2(sin(nd)*cos(incl), cos(nd)); 	sl = sin(incl)*sin(nd);	beta = atan2(sl, pyth(sl)); 	lograd = pturbr*2.30258509;	rad *= 1. + lograd; 	motion *= radian*mrad*mrad/(rad*rad);	semi = 3.34; 	lsun = 99.696678 + 0.9856473354*eday;	lsun *= radian;	elong = lambda - lsun;	ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));	dlong = atan2(pyth(ci), ci)/radian;	mag = -.003 + .01815*dlong + .0001023*dlong*dlong; 	helio();	geo();}