Return to mars.c CVS log

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

Annotation of researchv10no/cmd/sky/mars.c, revision 1.1

1.1     ! root        1: #include "sky.h"
        !             2: 
        !             3: extern struct marst
        !             4: {
        !             5:        float f[2];
        !             6:        char c[3];
        !             7: } marst[];
        !             8: 
        !             9: mars()
        !            10: {
        !            11:        double pturbl, pturbb, pturbr;
        !            12:        double lograd;
        !            13:        double dele, enom, vnom, nd, sl;
        !            14:        double q0, v0, t0, m0, j0 , s0, u0;
        !            15:        double lsun, elong, ci, dlong;
        !            16:        double planp[8];
        !            17:        struct marst *pp = &marst[0];
        !            18:        double olong;
        !            19:        double temp;
        !            20: 
        !            21: /*
        !            22:  *     The arguments nnd coefficients are taken from
        !            23:  *     Simon Newcomb, Tables of the Heliocentric Motion
        !            24:  *     of Mars
        !            25:  *     A.P.A.E. VI, part 4 (1895).
        !            26:  *
        !            27:  *     Here are the mean orbital elements.
        !            28:  */
        !            29: 
        !            30:        object = "Mars        ";
        !            31:        ecc = .09331290 + .000092064*capt - 0.077e-6*capt2;
        !            32:        incl = 1.850333 - 6.75e-4*capt + 12.61e-6*capt2;
        !            33:        node = 48.786442 + .770992*capt - 1.39e-6*capt2
        !            34:         - 5.33e-6*capt3;
        !            35:        argp = 334.218203 + 1.840758*capt + 1.299e-4*capt2
        !            36:         - 1.19e-6*capt3;
        !            37:        mrad = 1.52368840;
        !            38:        anom = 319.529425 + .5240207666*eday + 1.808e-4*capt2
        !            39:         + 1.19e-6*capt3;
        !            40:        motion = 0.5240711638;
        !            41: 
        !            42:        incl *= radian;
        !            43:        node *= radian;
        !            44:        argp *= radian;
        !            45:        anom = fmod(anom, 360.)*radian;
        !            46:        motion *= radian;
        !            47: 
        !            48: /*
        !            49:  *     Conventional mean anomalies of perturbing planets.
        !            50:  */
        !            51: 
        !            52:        q0 = 102.28  + 4.092334429*eday;
        !            53:        v0 = 212.388 + 1.60211831*eday;
        !            54:        t0 = 358.415  + .98559696*eday;
        !            55:        m0 = 319.530 + .52402078*eday;
        !            56:        j0 = 225.209 + .08307904*eday + 0.332*sin((134.4+38.5*capt)*radian);
        !            57:        s0 = 175.533 + .03344747*eday - 0.808*sin((134.4+38.5*capt)*radian);
        !            58:        u0 = 74.188 + 0.0117193*eday;
        !            59: 
        !            60:        q0 *= radian;
        !            61:        v0 *= radian;
        !            62:        t0 *= radian;
        !            63:        m0 *= radian;
        !            64:        j0 *= radian;
        !            65:        s0 *= radian;
        !            66:        u0 *= radian;
        !            67: 
        !            68:        planp[1] = q0;
        !            69:        planp[2] = v0;
        !            70:        planp[3] = t0;
        !            71:        planp[4] = m0;
        !            72:        planp[5] = j0;
        !            73:        planp[6] = s0;
        !            74:        planp[7] = u0;
        !            75: 
        !            76: /*
        !            77:  *     Computation of long period terms affecting the mean anomaly.
        !            78:  *             4*mars - 7*earth + 3*venus
        !            79:  *             3*jupiter - 8*mars + 4*earth
        !            80:  *             2*jupiter - - 6*mars + 3*earth
        !            81:  *             2*saturn - 2*mars + earth
        !            82:  *             jupiter - 2*mars + earth
        !            83:  *             5*saturn - 2*jupiter
        !            84:  */
        !            85: 
        !            86:        anom = anom
        !            87:         - (37.05 + 13.50*capt)*radsec
        !            88:         + 0.606*radsec*sin((212.87+119.051*capt)*radian)
        !            89:         + 52.490*radsec*sin((47.48+19.771*capt)*radian)
        !            90:         + 0.319*radsec*sin((116.88+773.444*capt)*radian)
        !            91:         + 0.130*radsec*sin((74.00+163.00*capt)*radian)
        !            92:         + 0.009*radsec*sin((325.00+753.67*capt)*radian)
        !            93:         + 0.280*radsec*sin((300.00+40.8*capt)*radian);
        !            94: 
        !            95: /*
        !            96:  *     Computation of elliptic orbit.
        !            97:  */
        !            98: 
        !            99:        enom = anom + ecc*sin(anom);
        !           100:        do {
        !           101:                dele = (anom - enom + ecc * sin(enom)) /
        !           102:                        (1. - ecc*cos(enom));
        !           103:                enom += dele;
        !           104:        } while(fabs(dele) > 1.e-8);
        !           105:        vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
        !           106:                        cos(enom/2.));
        !           107:        rad = mrad*(1. - ecc*cos(enom));
        !           108: 
        !           109: /*
        !           110:  *     Perturbations in longitude.
        !           111:  */
        !           112: 
        !           113:        pturbl = 0.043*sin(2.*anom);
        !           114:        for(;;){
        !           115:                if(pp->f[0]==0.){
        !           116:                        pp++;
        !           117:                        break;
        !           118:                }
        !           119:                pturbl += pp->f[0]*cos(pp->f[1] + pp->c[0]*m0 + pp->c[1]*planp[pp->c[2]]);
        !           120:                pp++;
        !           121:        }
        !           122: 
        !           123: /*
        !           124:  *     Perturbations in latitude.
        !           125:  */
        !           126: 
        !           127:        pturbb = 0.;
        !           128:        for(;;){
        !           129:                if(pp->f[0]==0.){
        !           130:                        pp++;
        !           131:                        break;
        !           132:                }
        !           133:                pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*m0 + pp->c[1]*planp[pp->c[2]]);
        !           134:                pp++;
        !           135:        }
        !           136: 
        !           137: /*
        !           138:  *     Perturbations in log radius vector.
        !           139:  */
        !           140: 
        !           141:        pturbr = 0.;
        !           142:        for(;;){
        !           143:                if(pp->f[0]==0.){
        !           144:                        pp++;
        !           145:                        break;
        !           146:                }
        !           147:                pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*m0 + pp->c[1]*planp[pp->c[2]]);
        !           148:                pp++;
        !           149:        }
        !           150:        pturbr *= 1.e-6;
        !           151: 
        !           152: /*
        !           153:  *     reduce to the ecliptic
        !           154:  */
        !           155: 
        !           156:        olong = vnom + argp + pturbl*radsec;
        !           157:        nd = olong - node;
        !           158:        lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
        !           159: 
        !           160:        sl = sin(incl)*sin(nd);
        !           161:        beta = atan2(sl, sqrt(1.-sl*sl)) + pturbb*radsec;
        !           162: 
        !           163:        lograd = pturbr*2.30258509;
        !           164:        rad *= 1. + lograd;
        !           165: 
        !           166: /*
        !           167:  *     Compute motion for planetary aberration.
        !           168:  */
        !           169: 
        !           170:        temp = motion*mrad*mrad*sqrt(1.-ecc*ecc)/(rad*rad);
        !           171:        ldot = temp*sin(2.*(lambda-node))/sin(2.*(olong-node));
        !           172:        bdot = temp*sin(incl)*cos(lambda-node);
        !           173:        rdot = motion*mrad*ecc*sin(olong-argp)/sqrt(1.-ecc*ecc);
        !           174: 
        !           175: /*
        !           176:  *     Compute magnitude.
        !           177:  */
        !           178: 
        !           179:        lsun = 99.696678 + 0.9856473354*eday;
        !           180:        lsun *= radian;
        !           181:        elong = lambda - lsun;
        !           182:        ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
        !           183:        dlong = atan2(sqrt(1.-ci*ci), ci)/radian;
        !           184:        mag = -1.30 + .01486*dlong;
        !           185: 
        !           186:        semi = 4.68;
        !           187: 
        !           188:        helio();
        !           189:        geo();
        !           190: 
        !           191: }

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