Return to nutate.c CVS log

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

Annotation of researchv10no/cmd/sky/nutate.c, revision 1.1.1.1

1.1       root        1: #include "sky.h"
                      2: 
                      3: extern struct nuttab
                      4: {
                      5:        float f[2];
                      6:        char c[5];
                      7: } nuttab[];
                      8: 
                      9: nutate()
                     10: {
                     11: 
                     12: /*
                     13:  *     uses radian, radsec, eday
                     14:  *     sets psi, eps, dpsi, deps, obliq, gst, tobliq
                     15:  */
                     16: 
                     17:        double msun, mnom, noded, dmoon;
                     18:        struct nuttab *pp;
                     19:        double arg;
                     20:        double Eday;
                     21:        double capT, capT2, capT3;
                     22: 
                     23: /*
                     24:  *     nutation of the equinoxes is a irregular wobble of the
                     25:  *     pole of the earths rotation whose magnitude is about
                     26:  *     9 seconds of arc and whose period is about 18.6 years.
                     27:  *
                     28:  *     it depends upon the pull of the sun and moon on the
                     29:  *     equatorial bulge of the earth.
                     30:  *
                     31:  *     psi and eps are the two angles which specify the
                     32:  *     true pole with respect to the mean pole.
                     33:  *
                     34:  *     all coeffieients are from Supplement to A.E. 1984
                     35:  */
                     36: 
                     37:        pp = &nuttab[0];
                     38:        Eday = eday - 36525;
                     39:        capT = capt-1;
                     40:        capT2 = capT*capT;
                     41:        capT3 = capT2*capT;
                     42: 
                     43:        mnom = 134.96298139 + 13.06499294724*Eday + 8.6972e-3*capT2
                     44:                 + 17.8e-6*capT3;
                     45:        msun = 357.52772333 + 0.98560028309*Eday - 0.16028e-3*capT2
                     46:                 - 3.33e-6*capT3;
                     47:        noded = 93.27191028 + 13.22935024060*Eday - 3.6825e-3*capT2
                     48:                 + 3.05e-6*capT3;
                     49:        dmoon = 297.85036306 + 12.19074911650*Eday - 1.91417e-3*capT2
                     50:                 + 5.28e-6*capT3;
                     51:        node = 125.04452222 - .05295376484*Eday + 2.07083e-3*capT2
                     52:                 + 2.22e-6*capT3;
                     53:        mnom *= radian;
                     54:        msun *= radian;
                     55:        noded *= radian;
                     56:        dmoon *= radian;
                     57:        node *= radian;
                     58: 
                     59: /*
                     60:  *     long period terms
                     61:  */
                     62: 
                     63: 
                     64:        psi = -(17.1996+.01742*capT)*sin(node);
                     65:        eps = (9.2025+0.00089*capT)*cos(node);
                     66:        for(;;){
                     67:                if(pp->f[0]==0.){
                     68:                        pp++;
                     69:                        break;
                     70:                }
                     71:                arg = pp->c[0]*mnom + pp->c[1]*msun + pp->c[2]*noded
                     72:                        + pp->c[3]*dmoon + pp->c[4]*node;
                     73:                psi += pp->f[0]*sin(arg);
                     74:                eps += pp->f[1]*cos(arg);
                     75:                pp++;
                     76:        }
                     77: 
                     78: /*
                     79:  *     short period terms
                     80:  */
                     81: 
                     82:        dpsi = 0.;
                     83:        deps = 0.;
                     84:        if((flags&APPARENT)==0){
                     85:                for(;;){
                     86:                        if(pp->f[0]==0.){
                     87:                                pp++;
                     88:                                break;
                     89:                        }
                     90:                        arg = pp->c[0]*mnom + pp->c[1]*msun + pp->c[2]*noded
                     91:                                + pp->c[3]*dmoon + pp->c[4]*node;
                     92:                        dpsi += pp->f[0]*sin(arg);
                     93:                        deps += pp->f[1]*cos(arg);
                     94:                        pp++;
                     95:                }
                     96:        }
                     97: 
                     98: 
                     99:        psi = (psi+dpsi)*radsec;
                    100:        eps = (eps+deps)*radsec;
                    101:        dpsi *= radsec;
                    102:        deps *= radsec;
                    103: 
                    104:        obliq = 23.43929111 - .013004167*capT - 0.164e-6*capT2
                    105:                 + 0.5036e-6*capT3;
                    106:        obliq *= radian;
                    107:        tobliq = obliq + eps;
                    108: 
                    109:        gst = 100.460618375 + 360.98564736629*Eday + 0.3882667e-3*capT2
                    110:                 - 0.025e-6*capT3;
                    111:        gst -= 180.;
                    112:        gst = fmod(gst, 360.);
                    113:        if(gst < 0.)
                    114:                gst += 360.;
                    115:        gst *= radian;
                    116:        gst += psi*cos(tobliq);
                    117: 
                    118: /*
                    119:  * print true obliquity, mean obliquity, mean GST, and true GST.
                    120:  */
                    121: 
                    122:        if(flags&GOOBIE){
                    123:                printf("Mean G. Sid. Time: ");
                    124:                prhms(gst - psi*cos(obliq), 4);
                    125:                printf("\n");
                    126: 
                    127:                printf("App. G. Sid. Time: ");
                    128:                prhms(gst,4);
                    129:                printf("\n");
                    130: 
                    131:                printf("Goobie: %6.3f %6.3f %6.3f %6.3f\n", psi/radsec,
                    132:                        eps/radsec, dpsi/radsec, deps/radsec);
                    133: 
                    134:                printf("Obliquity: ");
                    135:                prdms(tobliq,3);
                    136:                printf("\n");
                    137:        }
                    138: 
                    139: }

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