Return to moon.c CVS log

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

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

1.1     ! root        1: #include "sky.h"
        !             2: 
        !             3: /*
        !             4:  *     References:
        !             5:  *     Brown,
        !             6:  *     Improved Lunar Ephemeris
        !             7:  *     Supplement to A.E. 1968
        !             8:  *     Transformation corrections.
        !             9:  */
        !            10: 
        !            11: double k1, k2, k3, k4;
        !            12: double mnom, msun, noded, dmoon;
        !            13: double sinx();
        !            14: double cosx();
        !            15: 
        !            16: extern struct  moontab
        !            17: {
        !            18:        float   f;
        !            19:        char    c[4];
        !            20: } moontab[];
        !            21: 
        !            22: extern struct moont1
        !            23: {
        !            24:        float f1[2];
        !            25:        char c1[7];
        !            26: } moont1[];
        !            27: 
        !            28: moon()
        !            29: {
        !            30:        register struct moontab *mp;
        !            31:        register struct moont1 *mp1;
        !            32:        double dlong, lsun, psun;
        !            33:        double eccm, eccs, chp, cpe;
        !            34:        double q0, v0, t0, m0, j0, sn0, l0;
        !            35:        double arg1, arg2, arg3, arg4, arg5, arg6, arg7;
        !            36:        double arg8, arg9, arg10, arg11, arg12, arg13;
        !            37:        double arg14, arg15, arg16, arg17;
        !            38:        double arg18, arg19, arg20, arg21, arg22;
        !            39:        double dgamma, k5, k6;
        !            40:        double lterms, sterms, cterms, nterms, pterms, spterms;
        !            41:        double gamma1, gamma2, gamma3, arglat;
        !            42:        double xmp, ymp, zmp;
        !            43:        double temp1, temp2;
        !            44:        double shsd;
        !            45:        double obl2;
        !            46:        double planp[7];
        !            47: 
        !            48:        object = "moon        ";
        !            49: 
        !            50: /*
        !            51:  *     the fundamental elements - all referred to the epoch of
        !            52:  *     Jan 0.5, 1900 and to the mean equinox of date.
        !            53:  */
        !            54: 
        !            55:        dlong = 270.434164 + 13.1763965268*eday - .001133*capt2
        !            56:                 + 1.9e-6*capt3;
        !            57:        dlong -= .000086; /* empirical*/
        !            58:        argp = 334.329556 + .1114040803*eday - .010325*capt2
        !            59:                 - 12.5e-6*capt3;
        !            60:        node = 259.183275 - .0529539222*eday + .002078*capt2
        !            61:                 + 2.2e-6*capt3;
        !            62:        node += .000100; /* empirical */
        !            63:        lsun = 279.696678 + .9856473354*eday + .000303*capt2;
        !            64:        psun = 281.220844 + .0000470684*eday + .000453*capt2
        !            65:                 + 3.3e-6*capt3;
        !            66: 
        !            67:        dlong = fmod(dlong, 360.);
        !            68:        argp = fmod(argp, 360.);
        !            69:        node = fmod(node, 360.);
        !            70:        lsun = fmod(lsun, 360.);
        !            71:        psun = fmod(psun, 360.);
        !            72: 
        !            73:        eccm = 22639.550;
        !            74:        eccs = .01675104 - .00004180*capt;
        !            75:        incl = 18461.400;
        !            76:        cpe = 124.986;
        !            77:        chp = 3422.451;
        !            78: 
        !            79: /*
        !            80:  *     some subsidiary elements - they are all longitudes
        !            81:  *     and they are referred to the epoch 1/0.5 1900 and
        !            82:  *     to the fixed mean equinox of 1850.0.
        !            83:  */
        !            84: 
        !            85:        q0 = 177.481153 + 4.0923388020*eday;
        !            86:        v0 = 342.069128 + 1.6021304820*eday;
        !            87:        t0 =  98.998753 + 0.9856091138*eday;
        !            88:        m0 = 293.049675 + 0.5240329445*eday;
        !            89:        j0 = 237.352319 + 0.0830912295*eday;
        !            90:        sn0 = 265.869508 + 0.03345974279*eday;
        !            91:        l0 = 269.736239 +13.1763583100*eday;
        !            92: 
        !            93: /*
        !            94:  *     the following are periodic corrections to the
        !            95:  *     fundamental elements and constants.
        !            96:  *     arg4 is the "Great Venus Inequality".
        !            97:  */
        !            98: 
        !            99:        arg1 = 41.1 + 20.2*(capt+.5);
        !           100:        arg2 = dlong - argp + 33. + 3.*t0 - 10.*v0 - 2.6*(capt+.5);
        !           101:        arg3 = dlong + 3.*argp - 4.*lsun + 67. - 23.*t0 + 25.*m0;
        !           102:        arg4 = dlong - argp + 151.1 + 16.*t0 - 18.*v0 - (capt+.5);
        !           103:        arg5 = node;
        !           104:        arg6 = node + 276.2 - 2.3*(capt+.5);
        !           105:        arg7 = 152. + 119.*(capt+0.5);
        !           106:        arg8 = dlong + argp - 2.*lsun + 105. + t0 - 3.*q0;
        !           107:        arg9 = 313.9 + 13.*t0 - 8.*v0;
        !           108:        arg10 = dlong - argp + 112.0 + 29.*t0 - 26.*v0;
        !           109:        arg11 = dlong - argp + 21.*t0 - 21.*v0;
        !           110:        arg12 = dlong + argp - 2.*lsun + 273. + 21.*t0 - 20.*v0;
        !           111:        arg13 = dlong + argp - 2.*lsun + 303. + 8.*t0 - 12.*v0;
        !           112:        arg14 = 2.*lsun - 2.*node + 270. + 6.*t0 - 5.*v0;
        !           113:        arg15 = dlong + 2.*lsun - 3.*argp + 24.*t0 - 24.*v0;
        !           114:        arg16 = dlong - lsun + 262. + 12.*t0 - 15.*v0;
        !           115:        arg17 = dlong - lsun + 190. + 25.*t0 - 23.*v0;
        !           116:        arg18 = 43. - 8.*t0 + 15.*m0;
        !           117:        arg19 = node - lsun + 165. + 2.*m0;
        !           118:        arg20 = node + 290.1 - 0.9*(capt+.5);
        !           119:        arg21 = 115. + 38.5*(capt+.5);
        !           120:        arg22 = 216.0 - 8.*t0 + 15.*m0;
        !           121:        arg1 *= radian;
        !           122:        arg2 *= radian;
        !           123:        arg3 *= radian;
        !           124:        arg4 *= radian;
        !           125:        arg5 *= radian;
        !           126:        arg6 *= radian;
        !           127:        arg7 *= radian;
        !           128:        arg8 *= radian;
        !           129:        arg9 *= radian;
        !           130:        arg10 *= radian;
        !           131:        arg11 *= radian;
        !           132:        arg12 *= radian;
        !           133:        arg13 *= radian;
        !           134:        arg14 *= radian;
        !           135:        arg15 *= radian;
        !           136:        arg16 *= radian;
        !           137:        arg17 *= radian;
        !           138:        arg18 *= radian;
        !           139:        arg19 *= radian;
        !           140:        arg20 *= radian;
        !           141:        arg21 *= radian;
        !           142:        arg22 *= radian;
        !           143: 
        !           144:        dlong += (
        !           145:                   0.84 *sin(arg1)
        !           146:                 +  0.31 *sin(arg2)
        !           147:                 +  0.04 *sin(arg3)
        !           148:                 + 14.27 *sin(arg4)
        !           149:                 +  7.261*sin(arg5)
        !           150:                 +  0.282*sin(arg6)
        !           151:                 +  0.04 *sin(arg7)
        !           152:                 +  0.075*sin(arg8)
        !           153:                 +  0.237*sin(arg9)
        !           154:                 +  0.108*sin(arg10)
        !           155:                 +  0.030*sin(arg11)
        !           156:                 +  0.126*sin(arg12)
        !           157:                 +  0.033*sin(arg13)
        !           158:                 +  0.054*sin(arg14)
        !           159:                 +  0.010*sin(arg15)
        !           160:                 +  0.013*sin(arg16)
        !           161:                 +  0.013*sin(arg17)
        !           162:                 +  0.026*sin(arg18)
        !           163:                 + 0.017*sin(arg19)
        !           164:                 )/3600.;
        !           165: 
        !           166:        argp += (
        !           167:                 - 2.10 *sin(arg1)
        !           168:                 -  0.118*sin(arg4)
        !           169:                 -  2.076*sin(arg5)
        !           170:                 -  0.840*sin(arg6)
        !           171:                  - 0.100*sin(arg7)
        !           172:                 -  0.593*sin(arg9)
        !           173:                 -  0.065*sin(arg18)
        !           174:                )/3600.;
        !           175: 
        !           176:        node += (
        !           177:                   0.63*sin(arg1)
        !           178:                 +  0.17*sin(arg4)
        !           179:                 + 95.96*sin(arg5)
        !           180:                 + 15.58*sin(arg6)
        !           181:                 +  1.86*sin(arg20)
        !           182:                 )/3600.;
        !           183: 
        !           184:        t0 += (
        !           185:                -6.40*sin(arg1)
        !           186:                -0.27*sin(arg7)
        !           187:                -1.89*sin(arg9)
        !           188:                +0.20*sin(arg22)
        !           189:                 )/3600.;
        !           190: 
        !           191:        lsun += (
        !           192:                -6.40*sin(arg1)
        !           193:                -0.27*sin(arg7)
        !           194:                -1.89*sin(arg9)
        !           195:                +0.20*sin(arg22)
        !           196:                )/3600.;
        !           197: 
        !           198:        dgamma = -  4.318*cos(arg5)
        !           199:                 -  0.698*cos(arg6)
        !           200:                 -  0.083*cos(arg20);
        !           201: 
        !           202:        j0 +=
        !           203:                   0.33*sin(arg21);
        !           204: 
        !           205:        sn0 +=
        !           206:                 -  0.83*sin(arg21);
        !           207: 
        !           208: /*
        !           209:  *     the following factors account for the fact that the
        !           210:  *     eccentricity, solar eccentricity, inclination and
        !           211:  *     parallax used by Brown to make up his coefficients
        !           212:  *     are both wrong and out of date.  Brown did the same
        !           213:  *     thing in a different way.
        !           214:  */
        !           215: 
        !           216:        k1 = eccm/22639.500;
        !           217:        k2 = eccs/.01675104;
        !           218:        k3 = 1. + 2.708e-6 + .000108008*dgamma;
        !           219:        k4 = cpe/125.154;
        !           220:        k5 = chp/3422.700;
        !           221: 
        !           222: /*
        !           223:  *     the principal arguments that are used to compute
        !           224:  *     perturbations are the following differences of the
        !           225:  *     fundamental elements.
        !           226:  */
        !           227: 
        !           228:        mnom = dlong - argp;
        !           229:        msun = lsun - psun;
        !           230:        noded = dlong - node;
        !           231:        dmoon = dlong - lsun;
        !           232: 
        !           233: /*
        !           234:  *     solar terms in longitude
        !           235:  */
        !           236: 
        !           237:        lterms = 0.0;
        !           238:        mp = moontab;
        !           239:        for(;;) {
        !           240:                if(mp->f == 0.0){
        !           241:                        mp++;
        !           242:                        break;
        !           243:                }
        !           244:                lterms += sinx(mp->f,
        !           245:                        mp->c[0], mp->c[1],
        !           246:                        mp->c[2], mp->c[3], 0.0);
        !           247:                mp++;
        !           248:        }
        !           249: 
        !           250:        lterms +=
        !           251:                (294.e-9*eday - 9193./1296000.*dgamma + .0013)*sinx(1.0,0,0,0,2,0.)
        !           252:                +(-220.e-9*eday + 5282./1296000.*dgamma)*sinx(1.0,1,0,0,-2,0.);
        !           253: 
        !           254:        planp[1] = q0;
        !           255:        planp[2] = v0;
        !           256:        planp[3] = t0;
        !           257:        planp[4] = m0;
        !           258:        planp[5] = j0;
        !           259:        planp[6] = sn0;
        !           260: 
        !           261: /*
        !           262:  *     planetary terms in longitude
        !           263:  */
        !           264: 
        !           265:        mp1 = moont1;
        !           266:        for(;;){
        !           267:                if(mp1->f1[0] == 0.){
        !           268:                        mp1++;
        !           269:                        break;
        !           270:                }
        !           271:                lterms += sinx(mp1->f1[0], mp1->c1[0], mp1->c1[1],
        !           272:                        mp1->c1[2], mp1->c1[3],
        !           273:                        mp1->c1[4]*t0+mp1->c1[5]*planp[mp1->c1[6]]+mp1->f1[1]);
        !           274:                mp1++;
        !           275:        }
        !           276: 
        !           277:        lterms += sinx(-.189, 0,0,0,0, node) /*IAU*/
        !           278:                + sinx(-.013, -1,0,0,0, node) /*IAU*/
        !           279:                + sinx(-.013, 1,0,0,0, node); /*IAU*/
        !           280:        lterms += sinx(0.019, 0,0,0,0, 5.*t0-3.*v0+node+216.0);
        !           281:        lterms += sinx(0.016, 0,0,0,0, -5.*t0+3.*v0+node+075.0);
        !           282:        lterms += sinx(-.038, 0,0,0,0, 2.*node);
        !           283: 
        !           284: /*
        !           285:  *     solar terms in latitude
        !           286:  */
        !           287: 
        !           288:        sterms = 0.0;
        !           289:        for(;;) {
        !           290:                if(mp->f == 0.0){
        !           291:                        mp++;
        !           292:                        break;
        !           293:                }
        !           294:                sterms += sinx(mp->f,
        !           295:                        mp->c[0], mp->c[1],
        !           296:                        mp->c[2], mp->c[3], 0.0);
        !           297:                mp++;
        !           298:        }
        !           299: 
        !           300:        cterms = 0.0;
        !           301:        for(;;) {
        !           302:                if(mp->f == 0.0){
        !           303:                        mp++;
        !           304:                        break;
        !           305:                }
        !           306:                cterms += cosx(mp->f,
        !           307:                        mp->c[0], mp->c[1],
        !           308:                        mp->c[2], mp->c[3], 0.0);
        !           309:                mp++;
        !           310:        }
        !           311: 
        !           312:        nterms = 0.0;
        !           313:        for(;;) {
        !           314:                if(mp->f == 0.0){
        !           315:                        mp++;
        !           316:                        break;
        !           317:                }
        !           318:                nterms += sinx(mp->f,
        !           319:                        mp->c[0], mp->c[1],
        !           320:                        mp->c[2], mp->c[3], 0.0);
        !           321:                mp++;
        !           322:        }
        !           323: 
        !           324: /*
        !           325:  *     planetary terms in latitude
        !           326:  */
        !           327: 
        !           328:        pterms = 0.;
        !           329:        for(;;){
        !           330:                if(mp1->f1[0] == 0.){
        !           331:                        mp1++;
        !           332:                        break;
        !           333:                }
        !           334:                pterms += sinx(mp1->f1[0], mp1->c1[0], mp1->c1[1],
        !           335:                        mp1->c1[2], mp1->c1[3],
        !           336:                        mp1->c1[4]*t0+mp1->c1[5]*planp[mp1->c1[6]]+mp1->f1[1]);
        !           337:                mp1++;
        !           338:        }
        !           339: 
        !           340:        pterms +=
        !           341:                  sinx(0.014, 0,0,0,0, -2.*t0+2.*v0+l0+285.0)
        !           342:                 + sinx(0.027, 0,0,0,0, -1.*t0+1.*v0+l0+285.0)
        !           343:                 + sinx(0.015, 0,0,0,0, t0-v0+l0+105.0)
        !           344:                 + sinx(0.077, 0,0,0,0, 5.*t0-3.*v0+l0+215.6)
        !           345:                 + sinx(0.025, 0,0,0,0, -6.*t0+4.*v0+l0+255.0)
        !           346:                 + sinx(0.074, 0,0,0,0, -5.*t0+3.*v0+l0+051.6)
        !           347:                 + sinx(0.018, 0,0,0,0, -4.*t0+2.*v0+l0+075.0)
        !           348:                 + sinx(0.010, 0,0,0,0, -3.*t0+v0+l0+075.0)
        !           349:                 + sinx(0.030, 0,0,0,0, 8.*t0-5.*v0+l0+125.0);
        !           350:        pterms +=
        !           351:                 sinx(0.010, 0,0,0,0, -t0+2.*m0+l0+345.0)
        !           352:                 + sinx(0.035, 0,0,0,0, 2.*j0+l0+168.0)
        !           353:                 + sinx(0.018, 0,0,0,0, -2.*j0+l0+024.0)
        !           354:                 + sinx(-.017, 0,0,0,0, l0)
        !           355:                 + sinx(0.083, 0,0,1,0, 2.*node)
        !           356:                 + sinx(0.215, 0,0,0,0, dlong) /*IAU*/
        !           357:                 + sinx(-.012, -1,0,0,0, dlong) /*IAU*/
        !           358:                 + sinx(0.011, 1,0,0,0, dlong); /*IAU*/
        !           359: 
        !           360: /*
        !           361:  *     solar terms in parallax
        !           362:  */
        !           363: 
        !           364:        spterms = 3422.700;
        !           365:        for(;;) {
        !           366:                if(mp->f == 0.0){
        !           367:                        mp++;
        !           368:                        break;
        !           369:                }
        !           370:                spterms += cosx(mp->f,
        !           371:                        mp->c[0], mp->c[1],
        !           372:                        mp->c[2], mp->c[3], 0.0);
        !           373:                mp++;
        !           374:        }
        !           375: 
        !           376: /*
        !           377:  *     planetary terms in parallax
        !           378:  */
        !           379: 
        !           380:        for(;;){
        !           381:                if(mp1->f1[0] == 0.){
        !           382:                        mp1++;
        !           383:                        break;
        !           384:                }
        !           385:                spterms += cosx(mp1->f1[0], mp1->c1[0], mp1->c1[1],
        !           386:                        mp1->c1[2], mp1->c1[3],
        !           387:                        mp1->c1[4]*t0+mp1->c1[5]*planp[mp1->c1[6]]+mp1->f1[1]);
        !           388:                mp1++;
        !           389:        }
        !           390: 
        !           391: /*
        !           392:  *     computation of longitude
        !           393:  */
        !           394: 
        !           395:        lambda = (dlong + lterms/3600.)*radian;
        !           396: 
        !           397: /*
        !           398:  *     computation of latitude
        !           399:  */
        !           400: 
        !           401:        arglat = (noded + sterms/3600.)*radian;
        !           402:        gamma1 = 18519.700 * k3;
        !           403:        gamma2 = -6.241 * k3*k3*k3;
        !           404:        gamma3 = 0.004 * k3*k3*k3*k3*k3;
        !           405: 
        !           406:        k6 = (gamma1 + cterms) / gamma1;
        !           407: 
        !           408:        beta = k6 * (gamma1*sin(arglat) + gamma2*sin(3.*arglat)
        !           409:                 + gamma3*sin(5.*arglat) + nterms)
        !           410:                 + pterms;
        !           411:        if(flags & OCCULT)
        !           412:                beta -= 0.6;
        !           413:        beta *= radsec;
        !           414: 
        !           415: /*
        !           416:  *     computation of parallax
        !           417:  */
        !           418: 
        !           419:        spterms = k5 * spterms *radsec;
        !           420:        hp = spterms + (spterms*spterms*spterms)/6.;
        !           421: 
        !           422:        rad = hp/radsec;
        !           423:        georad = 1.;
        !           424:        semi = .0799 + .272453*(hp/radsec);
        !           425:        if(dmoon < 0.)
        !           426:                dmoon += 360.;
        !           427:        mag = dmoon/360.;
        !           428: 
        !           429: /*
        !           430:  *     change to equatorial coordinates
        !           431:  */
        !           432: 
        !           433:        lambda += psi;
        !           434:        obl2 = obliq + eps;
        !           435:        xmp = georad*cos(lambda)*cos(beta);
        !           436:        ymp = georad*(sin(lambda)*cos(beta)*cos(obl2) - sin(obl2)*sin(beta));
        !           437:        zmp = georad*(sin(lambda)*cos(beta)*sin(obl2) + cos(obl2)*sin(beta));
        !           438: 
        !           439:        alpha = atan2(ymp, xmp);
        !           440:        delta = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
        !           441: 
        !           442: /*
        !           443:  *c    lunar eclipse computation
        !           444:  */
        !           445: 
        !           446:        shsd = 1.0183*hp/radsec - 969.85/rps;
        !           447:        temp1 = sin(shdecl)*sin(delta) + cos(shdecl)*cos(delta)
        !           448:         *cos(shra - alpha);
        !           449:        temp2 = atan2(sqrt(1.-temp1*temp1),temp1)/radsec;
        !           450:        temp2 = semi + shsd - temp2;
        !           451:        temp2 = temp2/(2.*semi);
        !           452:        if(temp2 >= 0.){
        !           453:                if(temp2 < 1.)
        !           454:                        printf("Partial Lunar Eclipse: Mag. = %.4f\n", temp2);
        !           455:                else
        !           456:                        printf("Total Lunar Eclipse: Mag. = %.4f\n", temp2);
        !           457:        }
        !           458: 
        !           459: 
        !           460:        geolam = lambda;
        !           461:        geobet = beta;
        !           462: 
        !           463:        geo();
        !           464: 
        !           465: /*
        !           466:  *     solar eclipse computation
        !           467: */
        !           468: 
        !           469:        if(!((flags&GEO)||(flags&HELIO))){
        !           470:                temp1 = sin(sundec)*sin(decl2) + cos(sundec)*cos(decl2)
        !           471:                 *cos(sunra-ra);
        !           472:                temp1 = atan2(sqrt(1.-temp1*temp1),temp1)/radsec;
        !           473:                if(temp1 <= (semi2+sunsd)){
        !           474:                        temp2 = (semi2+sunsd-temp1)/(2.*sunsd);
        !           475:                        if(temp1 >= fabs(sunsd-semi2))
        !           476:                                printf("Partial Solar Eclipse: Mag. = %.4f\n", temp2);
        !           477:                        else if(temp2 > 1.)
        !           478:                                printf("Total Solar Eclipse: Mag. = %.4f\n", temp2);
        !           479:                        else
        !           480:                                printf("Annular Solar Eclipse: Mag. = %.4f\n", temp2);
        !           481:                }
        !           482:        }
        !           483: 
        !           484: /*
        !           485:  *     constants for occultation computations
        !           486:  */
        !           487: 
        !           488:        moonra = ra;
        !           489:        moonde = decl2;
        !           490:        moonsd = semi2;
        !           491: 
        !           492: }
        !           493: 
        !           494: double
        !           495: sinx(coef,i,j,k,m,angle)
        !           496: double coef, angle;
        !           497: {
        !           498:        double x;
        !           499: 
        !           500:        x = i*mnom + j*msun + k*noded + m*dmoon + angle;
        !           501:        x = coef*sin(x*radian);
        !           502:        if(i < 0)
        !           503:                i = -i;
        !           504:        for(; i>0; i--)
        !           505:                x *= k1;
        !           506:        if(j < 0)
        !           507:                j = -j;
        !           508:        for(; j>0; j--)
        !           509:                x *= k2;
        !           510:        if(k < 0)
        !           511:                k = -k;
        !           512:        for(; k>0; k--)
        !           513:                x *= k3;
        !           514:        if(m & 1)
        !           515:                x *= k4;
        !           516: 
        !           517:        return(x);
        !           518: }
        !           519: 
        !           520: double
        !           521: cosx(coef,i,j,k,m,angle)
        !           522: double coef, angle;
        !           523: {
        !           524:        double x;
        !           525: 
        !           526:        x = i*mnom + j*msun + k*noded + m*dmoon + angle;
        !           527:        x = coef*cos(x*radian);
        !           528:        if(i < 0)
        !           529:                i = -i;
        !           530:        for(; i>0; i--)
        !           531:                x *= k1;
        !           532:        if(j < 0)
        !           533:                j = -j;
        !           534:        for(; j>0; j--)
        !           535:                x *= k2;
        !           536:        if(k < 0)
        !           537:                k = -k;
        !           538:        for(; k>0; k--)
        !           539:                x *= k3;
        !           540:        if(m & 1)
        !           541:                x *= k4;
        !           542: 
        !           543:        return(x);
        !           544: }

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