nutate.c (1754B)
1 #include "astro.h" 2 3 void 4 nutate(void) 5 { 6 7 /* 8 * uses radian, radsec 9 * sets phi, eps, dphi, deps, obliq, gst, tobliq 10 */ 11 12 double msun, mnom, noded, dmoon; 13 14 /* 15 * nutation of the equinoxes is a wobble of the pole 16 * of the earths rotation whose magnitude is about 17 * 9 seconds of arc and whose period is about 18.6 years. 18 * 19 * it depends upon the pull of the sun and moon on the 20 * equatorial bulge of the earth. 21 * 22 * phi and eps are the two angles which specify the 23 * true pole with respect to the mean pole. 24 * 25 * all coeffieients are from Exp. Supp. pp.44-45 26 */ 27 28 mnom = 296.104608 + 13.0649924465*eday + 9.192e-3*capt2 29 + 14.38e-6*capt3; 30 mnom *= radian; 31 msun = 358.475833 + .9856002669*eday - .150e-3*capt2 32 - 3.33e-6*capt3; 33 msun *= radian; 34 noded = 11.250889 + 13.2293504490*eday - 3.211e-3*capt2 35 - 0.33e-6*capt3; 36 noded *= radian; 37 dmoon = 350.737486 + 12.1907491914*eday - 1.436e-3*capt2 38 + 1.89e-6*capt3; 39 dmoon *= radian; 40 node = 259.183275 - .0529539222*eday + 2.078e-3*capt2 41 + 2.22e-6*capt3; 42 node *= radian; 43 44 /* 45 * long period terms 46 */ 47 48 phi = 0.; 49 eps = 0.; 50 dphi = 0.; 51 deps = 0.; 52 53 54 icosadd(nutfp, nutcp); 55 phi = -(17.2327+.01737*capt)*sin(node); 56 phi += sinadd(4, node, noded, dmoon, msun); 57 58 eps = cosadd(4, node, noded, dmoon, msun); 59 60 /* 61 * short period terms 62 */ 63 64 65 dphi = sinadd(4, node, noded, mnom, dmoon); 66 67 deps = cosadd(3, node, noded, mnom); 68 69 phi = (phi+dphi)*radsec; 70 eps = (eps+deps)*radsec; 71 dphi *= radsec; 72 deps *= radsec; 73 74 obliq = 23.452294 - .0130125*capt - 1.64e-6*capt2 75 + 0.503e-6*capt3; 76 obliq *= radian; 77 tobliq = obliq + eps; 78 79 gst = 99.690983 + 360.9856473354*eday + .000387*capt2; 80 gst -= 180.; 81 gst = fmod(gst, 360.); 82 if(gst < 0.) 83 gst += 360.; 84 gst *= radian; 85 gst += phi*cos(obliq); 86 }