Return to sat.c CVS log

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

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

1.1       root        1: #include "sky.h"
                      2: 
                      3: sat()
                      4: {
                      5:        double pturbl, pturbb, pturbr;
                      6:        double lograd;
                      7:        double dele, enom, vnom, nd, sl;
                      8:        double q0, v0, t0, j0, s0;
                      9: 
                     10:        double capj, capn, eye, comg, omg;
                     11:        double sb, su, cu, u, b, up;
                     12:        double sd, ca, sa;
                     13: 
                     14:        object = "Saturn      ";
                     15: 
                     16:        ecc = .0558900 - .000347*capt;
                     17:        incl = 2.49256 - .0044*capt;
                     18:        node = 112.78364 + .87306*capt;
                     19:        argp = 91.08897 + 1.95917*capt;
                     20:        mrad = 9.538843;
                     21:        anom = 175.47630 + .03345972*eday - .56527*capt;
                     22:        motion = 120.4550/3600.;
                     23: 
                     24:        q0 = 102.28  + 4.092334429*eday;
                     25:        q0 *= radian;
                     26:        v0 = 212.536 + 1.602126105*eday;
                     27:        v0 *= radian;
                     28:        t0 = -1.45  + .985604737*eday;
                     29:        t0 *= radian;
                     30:        j0 = 225.36 + .083086735*eday;
                     31:        j0 *= radian;
                     32:        s0 = 175.68 + .033455441*eday;
                     33:        s0 *= radian;
                     34: 
                     35:        anom = anom;
                     36:        incl *= radian;
                     37:        node *= radian;
                     38:        argp *= radian;
                     39:        anom = fmod(anom,360.)*radian;
                     40: 
                     41:        enom = anom + ecc*sin(anom);
                     42:        do {
                     43:                dele = (anom - enom + ecc * sin(enom)) /
                     44:                        (1. - ecc*cos(enom));
                     45:                enom += dele;
                     46:        } while(fabs(dele) > 1.e-8);
                     47:        vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
                     48:                cos(enom/2.));
                     49:        rad = mrad*(1. - ecc*cos(enom));
                     50: 
                     51:        lambda = vnom + argp;
                     52: 
                     53:        pturbl = 0.;
                     54: 
                     55:        lambda += pturbl*radsec;
                     56: 
                     57:        pturbb = 0.;
                     58: 
                     59:        pturbr = 0.;
                     60: 
                     61: /*
                     62:  *     reduce to the ecliptic
                     63:  */
                     64: 
                     65:        nd = lambda - node;
                     66:        lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
                     67: 
                     68:        sl = sin(incl)*sin(nd) + pturbb*radsec;
                     69:        beta = atan2(sl, sqrt(1.-sl*sl));
                     70: 
                     71:        lograd = pturbr*2.30258509;
                     72:        rad *= 1. + lograd;
                     73: 
                     74: 
                     75:        lambda -= 1185.*radsec;
                     76:        beta -= 51.*radsec;
                     77: 
                     78:        motion *= radian*mrad*mrad/(rad*rad);
                     79:        semi = 83.33;
                     80: 
                     81: /*
                     82:  *     here begins the computation of magnitude
                     83:  *     first find the geocentric equatorial coordinates of Saturn
                     84:  */
                     85: 
                     86:        sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
                     87:                sin(beta)*cos(obliq));
                     88:        sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
                     89:                sin(beta)*sin(obliq));
                     90:        ca = rad*cos(beta)*cos(lambda);
                     91:        sd += zms;
                     92:        sa += yms;
                     93:        ca += xms;
                     94:        alpha = atan2(sa,ca);
                     95:        delta = atan2(sd,sqrt(sa*sa+ca*ca));
                     96: 
                     97: /*
                     98:  *     here are the necessary elements of Saturn's rings
                     99:  *     cf. Exp. Supp. p. 363ff.
                    100:  */
                    101: 
                    102:        capj = 6.9056 - 0.4322*capt;
                    103:        capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
                    104:        eye = 28.0743 - 0.0128*capt;
                    105:        comg = 168.1179 + 1.3936*capt;
                    106:        omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
                    107: 
                    108:        capj *= radian;
                    109:        capn *= radian;
                    110:        eye *= radian;
                    111:        comg *= radian;
                    112:        omg *= radian;
                    113: 
                    114: /*
                    115:  *     now find saturnicentric ring-plane coords of the earth
                    116:  */
                    117: 
                    118:        sb = sin(capj)*cos(delta)*sin(alpha-capn) -
                    119:                cos(capj)*sin(delta);
                    120:        su = cos(capj)*cos(delta)*sin(alpha-capn) +
                    121:                sin(capj)*sin(delta);
                    122:        cu = cos(delta)*cos(alpha-capn);
                    123:        u = atan2(su,cu);
                    124:        b = atan2(sb,sqrt(su*su+cu*cu));
                    125: 
                    126: /*
                    127:  *     and then the saturnicentric ring-plane coords of the sun
                    128:  */
                    129: 
                    130:        su = sin(eye)*sin(beta) +
                    131:                cos(eye)*cos(beta)*sin(lambda-comg);
                    132:        cu = cos(beta)*cos(lambda-comg);
                    133:        up = atan2(su,cu);
                    134: 
                    135: /*
                    136:  *     at last, the magnitude
                    137:  */
                    138: 
                    139: 
                    140:        sb = sin(b);
                    141:        mag = -8.68 +2.52*fabs(up+omg-u)-
                    142:                2.60*fabs(sb) + 1.25*(sb*sb);
                    143: 
                    144:        helio();
                    145:        geo();
                    146: }

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