Return to merc.c CVS log

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

Annotation of researchv10dc/cmd/sky/merc.c, revision 1.1

1.1     ! root        1: #include "sky.h"
        !             2: 
        !             3: extern struct merctab
        !             4: {
        !             5:        float f[2];
        !             6:        char c[3];
        !             7: } merctab[];
        !             8: 
        !             9: merc()
        !            10: {
        !            11:        double pturbl, pturbb, pturbr;
        !            12:        double lograd;
        !            13:        double dele, enom, vnom, nd, sl;
        !            14:        double q0, v0, t0, m0, j0 , s0;
        !            15:        double lsun, elong, ci, dlong;
        !            16:        double planp[7];
        !            17:        struct merctab *pp = &merctab[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 Mercury
        !            25:  *     A.P.A.E. VI, part 2 (1895).
        !            26:  *
        !            27:  *     Here are the mean orbital elements.
        !            28:  */
        !            29: 
        !            30:        object = "Mercury     ";
        !            31:        ecc = .20561421 + .00002046*capt - 0.03e-6*capt2;
        !            32:        ecc -= .00000250; /*empirical*/
        !            33:        ecc -= .00000003; /*empirical*/
        !            34:        incl = 7.0028806 + .0018608*capt - 18.3e-6*capt2;
        !            35:        incl += .00025;/* empirical */
        !            36:        node = 47.145944 + 1.185208*capt + .0001739*capt2;
        !            37:        node -= 2.2/3600.; /* empirical */
        !            38:        argp = 75.899697 + 1.555490*capt + .0002947*capt2;
        !            39:        argp += 0.06/3600. + 1.89/3600.; /* empirical */
        !            40:        mrad = .3870986;
        !            41:        anom = 102.279381 + 4.0923344364*eday + 6.7e-6*capt2;
        !            42:        anom += 2.07/3600. - 1.89/3600.; /* empirical */
        !            43:        motion = 4.0923770233;
        !            44: 
        !            45:        incl *= radian;
        !            46:        node *= radian;
        !            47:        argp *= radian;
        !            48:        anom = fmod(anom, 360.)*radian;
        !            49:        motion *= radian;
        !            50: 
        !            51: /*
        !            52:  *     Conventional mean anomalies of perturbing planets.
        !            53:  */
        !            54: 
        !            55:        q0 = 102.28  + 4.092334429*eday;
        !            56:        v0 = 212.536 + 1.602126105*eday;
        !            57:        t0 = -1.45  + .985604737*eday;
        !            58:        m0 = 319.66 + .524028480*eday;
        !            59:        j0 = 225.36 + .083086735*eday;
        !            60:        s0 = 175.68 + .033455441*eday;
        !            61: 
        !            62:        q0 *= radian;
        !            63:        v0 *= radian;
        !            64:        t0 *= radian;
        !            65:        m0 *= radian;
        !            66:        j0 *= radian;
        !            67:        s0 *= radian;
        !            68: 
        !            69:        planp[1] = q0;
        !            70:        planp[2] = v0;
        !            71:        planp[3] = t0;
        !            72:        planp[4] = m0;
        !            73:        planp[5] = j0;
        !            74:        planp[6] = s0;
        !            75: 
        !            76: /*
        !            77:  *     Computation of long period terms affecting the mean anomaly.
        !            78:  */
        !            79: 
        !            80:        anom = anom;
        !            81: 
        !            82: /*
        !            83:  *     Computation of elliptic orbit.
        !            84:  */
        !            85: 
        !            86:        enom = anom + ecc*sin(anom);
        !            87:        do {
        !            88:                dele = (anom - enom + ecc * sin(enom)) /
        !            89:                        (1. - ecc*cos(enom));
        !            90:                enom += dele;
        !            91:        } while(fabs(dele) > 1.e-8);
        !            92:        vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
        !            93:                        cos(enom/2.));
        !            94:        rad = mrad*(1. - ecc*cos(enom));
        !            95: 
        !            96: /*
        !            97:  *     Perturbations in longitude.
        !            98:  */
        !            99: 
        !           100:        pturbl = 0.;
        !           101:        for(;;){
        !           102:                if(pp->f[0]==0.){
        !           103:                        pp++;
        !           104:                        break;
        !           105:                }
        !           106:                pturbl += pp->f[0]*cos(pp->f[1] + pp->c[0]*q0 + pp->c[1]*planp[pp->c[2]]);
        !           107:                pp++;
        !           108:        }
        !           109: 
        !           110: /*
        !           111:  *     Perturbations in latitude.
        !           112:  */
        !           113: 
        !           114:        pturbb = 0.;
        !           115:        for(;;){
        !           116:                if(pp->f[0]==0.){
        !           117:                        pp++;
        !           118:                        break;
        !           119:                }
        !           120:                pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*q0 + pp->c[1]*planp[pp->c[2]]);
        !           121:                pp++;
        !           122:        }
        !           123: 
        !           124: /*
        !           125:  *     Perturbations in log radius vector.
        !           126:  */
        !           127: 
        !           128:        pturbr = 0.;
        !           129:        for(;;){
        !           130:                if(pp->f[0]==0.){
        !           131:                        pp++;
        !           132:                        break;
        !           133:                }
        !           134:                pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*q0 + pp->c[1]*planp[pp->c[2]]);
        !           135:                pp++;
        !           136:        }
        !           137:        pturbr *= 1.e-6;
        !           138: 
        !           139: /*
        !           140:  *     reduce to the ecliptic
        !           141:  */
        !           142: 
        !           143:        olong = vnom + argp + pturbl*radsec;
        !           144:        nd = olong - node;
        !           145:        lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
        !           146: 
        !           147:        sl = sin(incl)*sin(nd);
        !           148:        beta = atan2(sl, sqrt(1.-sl*sl)) + pturbb*radsec;
        !           149: 
        !           150:        lograd = pturbr*2.30258509;
        !           151:        rad *= 1. + lograd;
        !           152: 
        !           153: /*
        !           154:  *     Compute motion for planetary aberration.
        !           155:  */
        !           156: 
        !           157:        temp = motion*mrad*mrad*sqrt(1.-ecc*ecc)/(rad*rad);
        !           158:        ldot = temp*sin(2.*(lambda-node))/sin(2.*(olong-node));
        !           159:        bdot = temp*sin(incl)*cos(lambda-node);
        !           160:        rdot = motion*mrad*ecc*sin(olong-argp)/sqrt(1.-ecc*ecc);
        !           161: 
        !           162: /*
        !           163:  *     Compute magnitude.
        !           164:  */
        !           165: 
        !           166:        lsun = 99.696678 + 0.9856473354*eday;
        !           167:        lsun *= radian;
        !           168:        elong = lambda - lsun;
        !           169:        ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
        !           170:        dlong = atan2(sqrt(1.-ci*ci), ci)/radian;
        !           171:        mag = -.003 + .01815*dlong + .0001023*dlong*dlong;
        !           172: 
        !           173:        semi = 3.34;
        !           174: 
        !           175:        helio();
        !           176:        geo();
        !           177: 
        !           178: }

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