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     ! 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.