1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495 |
- /*
- * 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"
- void
- nutate(void)
- {
- /*
- * uses radian, radsec
- * sets phi, eps, dphi, deps, obliq, gst, tobliq
- */
- double msun, mnom, noded, dmoon;
- /*
- * nutation of the equinoxes is a 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.
- *
- * phi and eps are the two angles which specify the
- * true pole with respect to the mean pole.
- *
- * all coeffieients are from Exp. Supp. pp.44-45
- */
- mnom = 296.104608 + 13.0649924465*eday + 9.192e-3*capt2
- + 14.38e-6*capt3;
- mnom *= radian;
- msun = 358.475833 + .9856002669*eday - .150e-3*capt2
- - 3.33e-6*capt3;
- msun *= radian;
- noded = 11.250889 + 13.2293504490*eday - 3.211e-3*capt2
- - 0.33e-6*capt3;
- noded *= radian;
- dmoon = 350.737486 + 12.1907491914*eday - 1.436e-3*capt2
- + 1.89e-6*capt3;
- dmoon *= radian;
- node = 259.183275 - .0529539222*eday + 2.078e-3*capt2
- + 2.22e-6*capt3;
- node *= radian;
- /*
- * long period terms
- */
- phi = 0.;
- eps = 0.;
- dphi = 0.;
- deps = 0.;
- icosadd(nutfp, nutcp);
- phi = -(17.2327+.01737*capt)*sin(node);
- phi += sinadd(4, node, noded, dmoon, msun);
- eps = cosadd(4, node, noded, dmoon, msun);
- /*
- * short period terms
- */
- dphi = sinadd(4, node, noded, mnom, dmoon);
- deps = cosadd(3, node, noded, mnom);
- phi = (phi+dphi)*radsec;
- eps = (eps+deps)*radsec;
- dphi *= radsec;
- deps *= radsec;
- obliq = 23.452294 - .0130125*capt - 1.64e-6*capt2
- + 0.503e-6*capt3;
- obliq *= radian;
- tobliq = obliq + eps;
- gst = 99.690983 + 360.9856473354*eday + .000387*capt2;
- gst -= 180.;
- gst = fmod(gst, 360.);
- if(gst < 0.)
- gst += 360.;
- gst *= radian;
- gst += phi*cos(obliq);
- }
|