Return to sun.c CVS log

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

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

1.1     ! root        1: #include "sky.h"
        !             2: 
        !             3: extern struct suntab
        !             4: {
        !             5:        float f[2];
        !             6:        char c[3];
        !             7: } suntab[];
        !             8: /*
        !             9:  *     WARNING: Requires sign extension of characters.
        !            10:  */
        !            11: 
        !            12: sun()
        !            13: {
        !            14:        double mmerc, mven, merth, mmars, mjup, msat;
        !            15:        double dmoon, mmoon, gmoon;
        !            16:        double pturbb, pturbl, pturbr, lograd;
        !            17:        double planp[7];
        !            18:        struct suntab *pp;
        !            19:        double Eday, capT, capT2, capT3;
        !            20:        pp = &suntab[0];
        !            21:        Eday = eday - 36525.;
        !            22:        capT = Eday/36525.;
        !            23:        capT2 = capT*capT;
        !            24:        capT3 = capT2*capT;
        !            25: 
        !            26: /*
        !            27:  *     The arugments and coefficients are taken from
        !            28:  *     Simon Newcomb, "Tables of the Motion of the Earth upon
        !            29:  *     its Axis and About the Sun", A.P.A.E. VI (1895)
        !            30:  */
        !            31: 
        !            32:        object = "sun         ";
        !            33:        incl = 0.;
        !            34:        ecc = .01670911 - 4.205e-5 * capT - 1.26e-7*capT2;
        !            35:        node = 0.;
        !            36:        argp = 282.938352 + .0000470774*Eday + .000462*capT2
        !            37:                 + .000003*capT3;
        !            38:        mrad = 1.0;
        !            39:        anom = 357.527723 + .985600283*Eday - .000159*capT2
        !            40:                 - .000003*capT3;
        !            41:        motion = .9856473354;
        !            42: 
        !            43:        argp -= 0.00085; /*empirical*/
        !            44:        anom += 0.00085; /*empirical*/
        !            45: 
        !            46:        dmoon = 350.737681+12.1907491914*eday-.001436*capt2;
        !            47:        gmoon = 11.250889 + 13.2293504490*eday - .003212*capt2;
        !            48:        mmoon = 296.104608 + 13.0649924465*eday + 9.192e-3*capt2;
        !            49:        mmerc = 102.19 + 4.092329959*eday + .19;
        !            50:        mven  = 212.448 + 1.602121635*eday + .17;
        !            51:        merth = 358.476 + 0.985600267*eday + .19;
        !            52:        mmars = 319.590 + .524024095*eday + .11;
        !            53:        mjup  = 225.269 + .083082362*eday + .34;
        !            54:        msat  = 175.593 + .033450794*eday + .24;
        !            55: /*     muran = 72.248 + .011722663*eday + .19;*/
        !            56: 
        !            57:        dmoon = fmod(dmoon, 360.)*radian;
        !            58:        gmoon = fmod(gmoon, 360.)*radian;
        !            59:        mmoon = fmod(mmoon, 360.)*radian;
        !            60:        mmerc *= radian;
        !            61:        mven  *= radian;
        !            62:        merth *= radian;
        !            63:        mmars *= radian;
        !            64:        mjup *= radian;
        !            65:        msat  *= radian;
        !            66: 
        !            67: 
        !            68: /*
        !            69:        long period terms in the mean anomaly - the arguments
        !            70:        of the terms are:
        !            71:                4*mars - 7*earth + 3*venus
        !            72:                3*jup - 8*mars + 4*earth
        !            73:                13*earth - 8*venus
        !            74:                8*earth - 15*mars
        !            75: */
        !            76: 
        !            77:        anom = anom
        !            78:          + 0.266/3600.*sin((31.8 + 119.0*capt)*radian)
        !            79:          + 6.40 /3600.*sin((231.19 + 20.20*capt)*radian)
        !            80:          + 1.870/3600.*sin((57.24+150.27*capt)*radian);
        !            81: 
        !            82:        incl *= radian;
        !            83:        node *= radian;
        !            84:        argp *= radian;
        !            85:        anom = fmod(anom, 360.)*radian;
        !            86: 
        !            87: /*
        !            88:  *     computation of elliptic orbit
        !            89:  */
        !            90: 
        !            91:        lambda = anom + argp;
        !            92: 
        !            93:        pturbl = (2.*ecc - ecc*ecc*ecc/4.)*sin(anom)
        !            94:                + (ecc*ecc*5./4. - ecc*ecc*ecc*ecc*11./24.)*sin(2.*anom)
        !            95:                + (ecc*ecc*ecc*13./12.)*sin(3.*anom)
        !            96:                + (ecc*ecc*ecc*ecc*103./96.)*sin(4.*anom);
        !            97: 
        !            98:        pturbl -= 0.27*radsec*sin(anom);
        !            99: 
        !           100:        lambda += pturbl;
        !           101: 
        !           102:        beta = 0.;
        !           103: 
        !           104:        lograd = (ecc*ecc/4. + ecc*ecc*ecc*ecc/32.)
        !           105:                 + (-ecc + ecc*ecc*ecc*3./8.)*cos(anom)
        !           106:                 + (-ecc*ecc*3./4. + ecc*ecc*ecc*ecc*11./24.)*cos(2.*anom)
        !           107:                 + (-ecc*ecc*ecc*17./24.)*cos(3.*anom)
        !           108:                 + (-ecc*ecc*ecc*ecc*71./96.)*cos(4.*anom);
        !           109: 
        !           110:        lograd += 0.00000035; /*empirical*/
        !           111:        lograd += 0.00000070*cos(anom); /*empirical*/
        !           112: 
        !           113:        lograd = lograd/2.30258509;
        !           114: 
        !           115: /*
        !           116:  *     Computation of perturbations to elliptic orbit.
        !           117:  */
        !           118: 
        !           119:        planp[1] = mmerc;
        !           120:        planp[2] = mven;
        !           121:        planp[3] = merth;
        !           122:        planp[4] = mmars;
        !           123:        planp[5] = mjup;
        !           124:        planp[6] = msat;
        !           125: 
        !           126:        pturbl = 0.;
        !           127:        for(;;){
        !           128:                if(pp->f[0]==0.){
        !           129:                        pp++;
        !           130:                        break;
        !           131:                }
        !           132:                pturbl += pp->f[0]*cos(pp->f[1] + pp->c[0]*merth + pp->c[1]*planp[pp->c[2]]);
        !           133:                pp++;
        !           134:        }
        !           135: 
        !           136:        pturbl = pturbl
        !           137:          + 6.454*sin(dmoon)
        !           138:          + 0.015*sin(dmoon) /* empirical */
        !           139:          + 0.013*sin(3.*dmoon)
        !           140:          + 0.177*sin(dmoon + mmoon)
        !           141:          - 0.424*sin(dmoon - mmoon)
        !           142:          + 0.039*sin(3.*dmoon - mmoon)
        !           143:          - 0.064*sin(dmoon + merth)
        !           144:          + 0.172*sin(dmoon - merth)
        !           145: /*
        !           146:          - 0.013*sin(dmoon-mmoon-merth)
        !           147:          - 0.013*sin(2.*gmoon-2.*dmoon)
        !           148: */
        !           149:        ;
        !           150: 
        !           151:        pturbl *= radsec;
        !           152: 
        !           153: 
        !           154:        pturbb = 0.;
        !           155:        for(;;){
        !           156:                if(pp->f[0]==0.){
        !           157:                        pp++;
        !           158:                        break;
        !           159:                }
        !           160:                pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*merth + pp->c[1]*planp[pp->c[2]]);
        !           161:                pp++;
        !           162:        }
        !           163: 
        !           164: 
        !           165:        pturbb = pturbb
        !           166:          + 0.576*sin(gmoon)
        !           167:          + 0.016*sin(gmoon + mmoon)
        !           168:          - 0.047*sin(gmoon - mmoon)
        !           169:          + 0.021*sin(2.*dmoon - gmoon)
        !           170: /*
        !           171:          + 0.005*sin(2.*dmoon-gmoon-mmoon)
        !           172:          + 0.005*sin(gmoon+merth)
        !           173:          + 0.005*sin(gmoon-merth)
        !           174: */
        !           175:        ;
        !           176: 
        !           177:        pturbb *= radsec;
        !           178: 
        !           179: 
        !           180:        pturbr = 0.;
        !           181:        for(;;){
        !           182:                if(pp->f[0]==0.){
        !           183:                        pp++;
        !           184:                        break;
        !           185:                }
        !           186:                pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*merth + pp->c[1]*planp[pp->c[2]]);
        !           187:                pp++;
        !           188:        }
        !           189: 
        !           190: 
        !           191:        pturbr = pturbr
        !           192:          +13.360e-6*cos(dmoon)
        !           193:          + 0.030e-6*cos(3.*dmoon)
        !           194:          + 0.370e-6*cos(dmoon + mmoon)
        !           195:          - 1.330e-6*cos(dmoon - mmoon)
        !           196:          + 0.080e-6*cos(3.*dmoon - mmoon)
        !           197:          - 0.140e-6*cos(dmoon + merth)
        !           198:          + 0.360e-6*cos(dmoon - merth)
        !           199: /*
        !           200:          - 0.030e-6*cos(dmoon-mmoon-merth)
        !           201:          + 0.030e-6*cos(2.*gmoon-2.*dmoon)
        !           202: */
        !           203:        ;
        !           204:        lambda += pturbl;
        !           205:        if(lambda > 2.*pi)
        !           206:                lambda -= 2.*pi;
        !           207: 
        !           208:        beta += pturbb;
        !           209: 
        !           210:        lograd = (lograd+pturbr) * 2.30258509;
        !           211:        rad = 1. + lograd * (1. + lograd * (.5 + lograd/6.));
        !           212: 
        !           213:        motion *= radian*mrad*mrad/(rad*rad);
        !           214: 
        !           215:        semi = 961.182;
        !           216:        if(flags & OCCULT)
        !           217:                semi = 959.63;
        !           218:        mag = -26.5;
        !           219: }

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