home *** CD-ROM | disk | FTP | other *** search
/ Chip 2001 Mobile / Chip_Mobile_2001.iso / palm / hobby / palmoon / palmoon.EXE / moontool.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  2000-07-21  |  12.2 KB  |  387 lines
  1. #include "moontool.h"
  2.  
  3.  
  4.  
  5. /*  KEPLER  --  Solve the equation of Kepler.  */
  6.  
  7. double kepler(double m, double ecc)
  8. {
  9.         double e, delta;
  10. #define EPSILON 1E-6
  11.  
  12.         e = m = torad(m);
  13.         do {
  14.            delta = e - ecc * sin(e) - m;
  15.            e -= delta / (1 - ecc * cos(e));
  16.         } while (abs(delta) > EPSILON);
  17.         return e;
  18. }
  19.  
  20.  
  21.  
  22. /*  PHASE  --  Calculate phase of moon as a fraction:
  23.  
  24.         The argument is the time for which the phase is requested,
  25.         expressed as a Julian date and fraction.  Returns the terminator
  26.         phase angle as a percentage of a full circle (i.e., 0 to 1),
  27.         and stores into pointer arguments the illuminated fraction of
  28.         the Moon's disc, the Moon's age in days and fraction, the
  29.         distance of the Moon from the centre of the Earth, and the
  30.         angular diameter subtended by the Moon as seen by an observer
  31.         at the centre of the Earth.
  32.  
  33. */
  34.  
  35. double phase(
  36.    double pdate,                      /* Date for which to calculate phase */
  37.    double *pphase,                    /* Illuminated fraction */
  38.    double *mage,                      /* Age of moon in days */
  39.    double *dist,                      /* Distance in kilometres */
  40.    double *angdia,                    /* Angular diameter in degrees */
  41.    double *sudist,                    /* Distance to Sun */
  42.    double *suangdia)                  /* Sun's angular diameter */
  43. {
  44.  
  45.         double Day, N, M, Ec, Lambdasun, ml, MM, MN, Ev, Ae, A3, MmP,
  46.                mEc, A4, lP, V, lPP, NP, y, x, Lambdamoon, BetaM,
  47.                MoonAge, MoonPhase,
  48.                MoonDist, MoonDFrac, MoonAng, MoonPar,
  49.                F, SunDist, SunAng;
  50.  
  51.         /* Calculation of the Sun's position */
  52.  
  53.         Day = pdate - epoch;        /* Date within epoch */
  54.         N = fixangle((360 / 365.2422) * Day); /* Mean anomaly of the Sun */
  55.         M = fixangle(N + elonge - elongp);    /* Convert from perigee
  56.                                        co-ordinates to epoch 1980.0 */
  57.         Ec = kepler(M, eccent);     /* Solve equation of Kepler */
  58.         Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2);
  59.         Ec = 2 * todeg(atan(Ec));   /* True anomaly */
  60.         Lambdasun = fixangle(Ec + elongp);  /* Sun's geocentric ecliptic
  61.                                                longitude */
  62.         /* Orbital distance factor */
  63.         F = ((1 + eccent * cos(torad(Ec))) / (1 - eccent * eccent));
  64.         SunDist = sunsmax / F;      /* Distance to Sun in km */
  65.         SunAng = F * sunangsiz;     /* Sun's angular size in degrees */
  66.  
  67.  
  68.         /* Calculation of the Moon's position */
  69.  
  70.         /* Moon's mean longitude */
  71.         ml = fixangle(13.1763966 * Day + mmlong);
  72.  
  73.         /* Moon's mean anomaly */
  74.         MM = fixangle(ml - 0.1114041 * Day - mmlongp);
  75.  
  76.         /* Moon's ascending node mean longitude */
  77.         MN = fixangle(mlnode - 0.0529539 * Day);
  78.  
  79.         /* Evection */
  80.         Ev = 1.2739 * sin(torad(2 * (ml - Lambdasun) - MM));
  81.  
  82.         /* Annual equation */
  83.         Ae = 0.1858 * sin(torad(M));
  84.  
  85.         /* Correction term */
  86.         A3 = 0.37 * sin(torad(M));
  87.  
  88.         /* Corrected anomaly */
  89.         MmP = MM + Ev - Ae - A3;
  90.  
  91.         /* Correction for the equation of the centre */
  92.         mEc = 6.2886 * sin(torad(MmP));
  93.  
  94.         /* Another correction term */
  95.         A4 = 0.214 * sin(torad(2 * MmP));
  96.  
  97.         /* Corrected longitude */
  98.         lP = ml + Ev + mEc - Ae + A4;
  99.  
  100.         /* Variation */
  101.         V = 0.6583 * sin(torad(2 * (lP - Lambdasun)));
  102.  
  103.         /* True longitude */
  104.         lPP = lP + V;
  105.  
  106.         /* Corrected longitude of the node */
  107.         NP = MN - 0.16 * sin(torad(M));
  108.  
  109.         /* Y inclination coordinate */
  110.         y = sin(torad(lPP - NP)) * cos(torad(minc));
  111.  
  112.         /* X inclination coordinate */
  113.         x = cos(torad(lPP - NP));
  114.  
  115.         /* Ecliptic longitude */
  116.         Lambdamoon = todeg(atan2(y, x));
  117.         Lambdamoon += NP;
  118.  
  119.         /* Ecliptic latitude */
  120.         BetaM = todeg(asin(sin(torad(lPP - NP)) * sin(torad(minc))));
  121.  
  122.         /* Calculation of the phase of the Moon */
  123.  
  124.         /* Age of the Moon in degrees */
  125.         MoonAge = lPP - Lambdasun;
  126.  
  127.         /* Phase of the Moon */
  128.         MoonPhase = (1 - cos(torad(MoonAge))) / 2;
  129.  
  130.         /* Calculate distance of moon from the centre of the Earth */
  131.  
  132.         MoonDist = (msmax * (1 - mecc * mecc)) /
  133.                    (1 + mecc * cos(torad(MmP + mEc)));
  134.  
  135.         /* Calculate Moon's angular diameter */
  136.  
  137.         MoonDFrac = MoonDist / msmax;
  138.         MoonAng = mangsiz / MoonDFrac;
  139.  
  140.         /* Calculate Moon's parallax */
  141.  
  142.         MoonPar = mparallax / MoonDFrac;
  143.  
  144.         *pphase = MoonPhase;
  145.         *mage = synmonth * (fixangle(MoonAge) / 360.0);
  146.         *dist = MoonDist;
  147.         *angdia = MoonAng;
  148.         *sudist = SunDist;
  149.         *suangdia = SunAng;
  150.         return fixangle(MoonAge) / 360.0;
  151. }
  152.  
  153.  
  154.  
  155. /* JTIME --    Convert internal GMT date and time to astronomical Julian
  156.                time (i.e. Julian date plus day fraction, expressed as
  157.                a double).  */
  158. /* Algorithm as given in Meeus, Astronomical Algorithms, Chapter 7, page 61 */
  159. FlpCompDouble *jtime(const DateTimeType *t, FlpCompDouble *j)
  160. {
  161.    Int32 year;
  162.    Int16 mon, mday, hour, min, sec;
  163.    Int16 a, b, m;
  164.    Int32 y;
  165.  
  166.    year = t->year;
  167.    mon  = t->month;
  168.    mday = t->day;
  169.    hour = t->hour;
  170.    min  = t->minute;
  171.    sec  = t->second;
  172.  
  173.     m = mon;
  174.     y = year;
  175.  
  176.    if (m <= 2) {
  177.       y--;
  178.       m += 12;
  179.    }
  180.  
  181.    /* Determine whether date is in Julian or Gregorian calendar based on
  182.       canonical date of calendar reform. */
  183.  
  184.    if ((year < 1582) || ((year == 1582) && ((mon < 9) ||
  185.        (mon == 9 && mday < 5)))) {
  186.       b = 0;
  187.    } else {
  188.       a = ((Int32) (y / 100));
  189.       b = 2 - a + (a / 4);
  190.    }
  191.    j->d = (((UInt32) (365.25 * (y + 4716))) + ((UInt16) (30.6001 * (m + 1))) +
  192.            mday + b - 1524.5) +
  193.            ((sec + 60L * (min + 60L * hour)) / 86400.0);
  194.    return j;
  195. }
  196.  
  197.  
  198.  
  199. /*  JYEAR  --  Convert Julian date to year, month, day, which are
  200.                returned via integer pointers to integers.  */
  201.  
  202. void jyear(double td, Int32 *yy, Int16 *mm, Int16 *dd)
  203. {
  204.    double z, f, a, alpha, b, c, d, e;
  205.  
  206.    td += 0.5;
  207.    z = floor(td);
  208.    f = td - z;
  209.  
  210.    if (z < 2299161.0) {
  211.       a = z;
  212.    } else {
  213.       alpha = floor((z - 1867216.25) / 36524.25);
  214.       a = z + 1 + alpha - floor(alpha / 4);
  215.    }
  216.  
  217.    b = a + 1524;
  218.    c = floor((b - 122.1) / 365.25);
  219.    d = floor(365.25 * c);
  220.    e = floor((b - d) / 30.6001);
  221.  
  222.    *dd = (Int16) (b - d - floor(30.6001 * e) + f);
  223.    *mm = (Int16) ((e < 14) ? (e - 1) : (e - 13));
  224.    *yy = (Int32) ((*mm > 2) ? (c - 4716) : (c - 4715));
  225. }
  226.  
  227.  
  228.  
  229. /*  JHMS  --  Convert Julian time to hour, minutes, and seconds.  */
  230.  
  231. void jhms(double j, Int16 *h, Int16 *m, Int16 *s)
  232. {
  233.    Int32 ij;
  234.  
  235.    j += 0.5;                  /* Astronomical to civil */
  236.    ij = (Int32) (((j - floor(j)) * 86400.0) + 0.5);  // Round to nearest second
  237.    *h = (Int16) (ij / 3600L);
  238.    *m = (Int16) ((ij / 60L) % 60L);
  239.    *s = (Int16) (ij % 60L);
  240. }
  241.  
  242.  
  243.  
  244. /*  MEANPHASE  --  Calculates  time  of  the mean new Moon for a given
  245.                    base date.  This argument K to this function is the
  246.                    precomputed synodic month index, given by:
  247.  
  248.                           K = (year - 1900) * 12.3685
  249.  
  250.                    where year is expressed as a year and fractional year.  */
  251.  
  252. double meanphase(double sdate, double k)
  253. {
  254.     double t, t2, t3, nt1;
  255.  
  256.     /* Time in Julian centuries from 1900 January 0.5 */
  257.     t = (sdate - 2415020.0) / 36525;
  258.     t2 = t * t;                       /* Square for frequent use */
  259.     t3 = t2 * t;                      /* Cube for frequent use */
  260.  
  261.     nt1 = 2415020.75933 + synmonth * k
  262.             + 0.0001178 * t2
  263.             - 0.000000155 * t3
  264.             + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);
  265.  
  266.     return nt1;
  267. }
  268.  
  269.  
  270. /*  TRUEPHASE  --  Given a K value used to determine the
  271.                    mean phase of the new moon, and a phase
  272.                    selector (0.0, 0.25, 0.5, 0.75), obtain
  273.                    the true, corrected phase time.  */
  274.  
  275. double truephase(double k, double phase)
  276. {
  277.         double t, t2, t3, pt, m, mprime, f;
  278.         Int16 apcor = 0;
  279.  
  280.         k += phase;                /* Add phase to new moon time */
  281.         t = k / 1236.85;           /* Time in Julian centuries from
  282.                                       1900 January 0.5 */
  283.         t2 = t * t;                /* Square for frequent use */
  284.         t3 = t2 * t;               /* Cube for frequent use */
  285.         pt = 2415020.75933         /* Mean time of phase */
  286.              + synmonth * k
  287.              + 0.0001178 * t2
  288.              - 0.000000155 * t3
  289.              + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2);
  290.  
  291.         m = 359.2242               /* Sun's mean anomaly */
  292.             + 29.10535608 * k
  293.             - 0.0000333 * t2
  294.             - 0.00000347 * t3;
  295.         mprime = 306.0253          /* Moon's mean anomaly */
  296.             + 385.81691806 * k
  297.             + 0.0107306 * t2
  298.             + 0.00001236 * t3;
  299.         f = 21.2964                /* Moon's argument of latitude */
  300.             + 390.67050646 * k
  301.             - 0.0016528 * t2
  302.             - 0.00000239 * t3;
  303.         if ((phase < 0.01) || (abs(phase - 0.5) < 0.01)) {
  304.  
  305.            /* Corrections for New and Full Moon */
  306.  
  307.            pt +=     (0.1734 - 0.000393 * t) * dsin(m)
  308.                     + 0.0021 * dsin(2 * m)
  309.                     - 0.4068 * dsin(mprime)
  310.                     + 0.0161 * dsin(2 * mprime)
  311.                     - 0.0004 * dsin(3 * mprime)
  312.                     + 0.0104 * dsin(2 * f)
  313.                     - 0.0051 * dsin(m + mprime)
  314.                     - 0.0074 * dsin(m - mprime)
  315.                     + 0.0004 * dsin(2 * f + m)
  316.                     - 0.0004 * dsin(2 * f - m)
  317.                     - 0.0006 * dsin(2 * f + mprime)
  318.                     + 0.0010 * dsin(2 * f - mprime)
  319.                     + 0.0005 * dsin(m + 2 * mprime);
  320.            apcor = 1;
  321.         } else if ((abs(phase - 0.25) < 0.01 || (abs(phase - 0.75) < 0.01))) {
  322.            pt +=     (0.1721 - 0.0004 * t) * dsin(m)
  323.                     + 0.0021 * dsin(2 * m)
  324.                     - 0.6280 * dsin(mprime)
  325.                     + 0.0089 * dsin(2 * mprime)
  326.                     - 0.0004 * dsin(3 * mprime)
  327.                     + 0.0079 * dsin(2 * f)
  328.                     - 0.0119 * dsin(m + mprime)
  329.                     - 0.0047 * dsin(m - mprime)
  330.                     + 0.0003 * dsin(2 * f + m)
  331.                     - 0.0004 * dsin(2 * f - m)
  332.                     - 0.0006 * dsin(2 * f + mprime)
  333.                     + 0.0021 * dsin(2 * f - mprime)
  334.                     + 0.0003 * dsin(m + 2 * mprime)
  335.                     + 0.0004 * dsin(m - 2 * mprime)
  336.                     - 0.0003 * dsin(2 * m + mprime);
  337.            if (phase < 0.5)
  338.               /* First quarter correction */
  339.               pt += 0.0028 - 0.0004 * dcos(m) + 0.0003 * dcos(mprime);
  340.            else
  341.               /* Last quarter correction */
  342.               pt += -0.0028 + 0.0004 * dcos(m) - 0.0003 * dcos(mprime);
  343.            apcor = 1;
  344.         }
  345.         ErrFatalDisplayIf((!apcor),
  346.                           "TRUEPHASE called with invalid phase selector");
  347.         return pt;
  348. }
  349.  
  350.  
  351.  
  352. /*  PHASEHUNT  --  Find time of phases of the moon which surround
  353.                    the current date.  Five phases are found, starting
  354.                    and ending with the new moons which bound the
  355.                    current lunation.  */
  356.  
  357. void phasehunt(double sdate, double phases[5])
  358. {
  359.     double adate, k1, k2, nt1, nt2;
  360.     Int32 yy;
  361.     Int16 mm, dd;
  362.  
  363.     adate = sdate - 45;
  364.  
  365.     jyear(adate, &yy, &mm, &dd);
  366.     k1 = floor((yy + ((mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685);
  367.  
  368.     adate = nt1 = meanphase(adate, k1);
  369.     while (1) {
  370.         adate += synmonth;
  371.         k2 = k1 + 1;
  372.         nt2 = meanphase(adate, k2);
  373.         if (nt1 <= sdate && nt2 > sdate)
  374.             break;
  375.         nt1 = nt2;
  376.         k1 = k2;
  377.     }
  378.     phases[0] = truephase(k1, 0.0);
  379.     phases[1] = truephase(k1, 0.25);
  380.     phases[2] = truephase(k1, 0.5);
  381.     phases[3] = truephase(k1, 0.75);
  382.     phases[4] = truephase(k2, 0.0);
  383. }
  384.  
  385.  
  386.  
  387.