Volatile Yard 挥洒自如 Write anything I want to write...

I don't like to eat meats since I was a kid as I really don't like their taste, while my family members are ok with them. I vomited a few times trying to take meats. Wonder if this is fated? I recall an interesting incident many years ago that might probably tell.

One day, my family went to buy satay to eat together at home. We hardly had satay at home in my hometown, and this is in fact the only time I also ate together as I can remember. (I just enjoy chewing the satay but not really its taste. Odd huh?) Something interesting then happened. Well, you normally will not find any bone from satay, unless the hawker had missed it accidentally. Now the incident:

We bought 20 satay sticks altogether. I took 4 sticks to eat, and I found real 'genuine' chicken bones (later examined by my family members too) in 3 out of that 4 sticks! While they didn't find any bone in other sticks! I swear that I did not choose them purposely as I didn't see the bones before I eat them. I treat it as a sign that I was fated to be discouraged from taking meats, ha... So now the mathematical questions:

• There are 20 satay sticks, and bones are found in 3 of them. What is the probability that, when you choose 4 of them randomly, you picked up the 3 sticks with bone?
• Given the same scenario, what is the probability if you choose only 3 of them randomly but happened to picked up all the 3 sticks with bone? (I might have picked those 3 too if I only choose 3 of them. Unfortunately this cannot be verified now.)

Check below for the solution. Think first if you can!

Puzzle

I got a 'mathematical dream' last night, which involves an original mathematical puzzle. In the dream, I was in a book fair looking for some inexpensive second-hand books. Then, I saw a thick book with a title like "Adfferv...", well, which should not be a proper English word. The book covers a lot of things, and is written like "Sophie's World". The main character in the book somehow met a wise person, who taught him or her some good knowledge in literatures, physics, mathematics, etc. I browsed the book and saw a page with an eye-catching mathematical puzzle. It looked something like this:

The question is: "These people need to pay the money. But how much?" I somehow figured out the answer, and verified it with that listed at the back of the book. Are you able to figure it out? Read on for the solution.

Notes:

• From left to right: A half-burnt candle on a table; A fire hydrant [en.wikipedia.org] (also known as a fire plug or a johnny bump); Several people;
• Disclaimer: Don't try too hard! While the puzzle and the answer do make sense, they are not really rigorous and can be controversial. After all, all these things including the 'correct' answer are just from my dream!
• Hint: Trying to think in English may help here.

We normally need to learn some important trigonometric identities in high-school. For the identities among the trigonometric functions, some are 'manageable' and are not difficult to memorize. E.g.,

• Reciprocal Identities (E.g., $\sin \alpha = \frac{1}{\csc \alpha}$)
• Pythagorean Identities (E.g., $\sin^2 \alpha + \cos^2 \alpha = 1$)
• Quotient Identities (E.g., $\tan \alpha = \frac{\sin \alpha}{\cos \alpha}$)

But some are not easily to memorize. In particular, I used to get confused with the following identities:

• Sum-Difference Formulas
• Sum-to-Product Formulas
• Product-to-Sum Formulas

Why is it so? I analyzed them and made the following observations:

• The former formulas can be easily derived from the definitions of the trigonometric functions, normally involve only one angle, can obtain them by 'intuition' with the help of a right triangle, or some rather straightforward methods.
• The latter formulas involve two angles and there's a barrier to obtain the 'intuition' of manipulating them.

In particular, all the latter formulas can in fact be derived easily once you have $\sin (\alpha \pm \beta) = \sin \alpha \cos \beta \pm \cos \alpha \sin \beta$. Unfortunately, the identity makes some sense but not much. (Proving this is in fact tedious although not difficult.) I have not found any existing trick for helping us in memorizing them. Deriving all of them from this formula is feasible but not practical.

To solve the problem, I analyzed their patterns and 'invented' my own way. I came out with the following diagram:

    +   -
S  SC  CS
C  CC -SS
-> 2 /2 /2
<- /2

With this I have been able to memorize those formulas with ease. I can't guarantee that this method helps you, because the diagram can be as confusing as the formulas themselves! If you have a better way, please share! Read on if you are interested in my method.

1. Introductions

This is a rough guide on how to build or implement your own ephemeris program or library routines. There are some existing programs (and some of them are free) to do such a job. However, there are a few advantages if you build your own one.

• It's your thing. You have complete control of the licensing and the copyright of the program. You can also improve the programs anyhow you like.
• Good way to learn Spherical Astronomy. Of course, you need to like mathematics and astronomy. If you have no interest in these subjects at all, it's better to use an existing program or libraries routines for your purpose.
• Extend it for many other purposes. Since you are familiar with your program, you can easily extend it for many purposes such as (a) calculating the phase of the Moon; (b) predicting eclipses; (c) building your birth charts for many different traditions of astrology; (d) creating calendars (such as the Chinese calendar) that are based on astronomical events; (e) many other things.

While these sound interesting, you have to note that building your very own ephemeris is not for the faint of heart! If you are not familiar with spherical astronomy, be ready to learn a lot of concepts and technical terms related to the subject. As for implementation, expect to write

• a few hundred lines just for very basic set of features and low-precision calculations;
• a few thousand lines for more features and satisfactory accuracy;
• a few ten thousands lines for more complete features and high-precision calculations.

There are some easier ways to tackle this problem. For example, you can interpolate the positions of the planets from some existing ephemerides (such as DE405, DE406). Such an approach is definitely feasible and realistic. But here we focus more on performing the calculations "from scratch" which is useful in many cases. E.g., to calculate just the rough rising and setting time of the sun for a few centuries, it's better to stick to a simpler and less space-consuming low-precision calculation rather than to make use of the gigantic ephemerides.

Let's start.

2. Knowledge Required

Before starting your implementation, there are many concepts you need to understand and study. Here's a list of some important ones.

Mathematical Software Library Routines. In particular, you need to dial with a lot of trigonometric and logarithm calculations. Be very familiar with these to avoid making stupid mistakes. For example, a trigonometric function generally expect an angle expressed in radians, so don't use an angle expressed in degrees directly.

Concepts in Spherical Astronomy. You need to know a lot of thing specific to this subject. Knowing some popular concepts in astronomy like the big bang or the dark energy won't help you here, as we are dealing with locations here! Below is a list of some common concepts just to give you an idea what to expect. Note that the list is NOT exhaustive! The higher precision you want to calculate and implement, the more and the deeper of these concepts you need to understand.

• Time Tracking. TT/TDT (Terrestrial Time); ET (Ephemeris Time); ST (Sidereal Time); UTC (Coordinated Universal Time); UT (Universal Time); Delta-T; Julian day; Julian ephemeris day; Gregorian calendar; etc.
• Coordinate Systems. Horizontal coordinates (altitude, azimuth); Equatorial coordinates (right ascension, declination, equinox, solstice, hour-angle, transit, meridian); Ecliptic coordinates (longitude, latitude, obliquity); Galactic coordinates; Heliocentric; Geocentric; Topocentric; etc.
• Reference Frames. ICRF (International Celestial Reference Frame); FK5 (Fundamental Katalog 5); epoch; mean equator and equinox, etc.
• Orbital Elements. Perihelion; Aphelion; Perigee; Apogee; Eccentricity; Kepler's equation; Mean longitude; Mean anomaly; True anomaly; Equation of the center; etc.
• Miscellaneous. Astronomical Unit; Radius Vector; Aberration; Precession; Nutation; Parallax; Refraction; Twilight; Conjunction; Opposition; True/Geometric position; Mean position; Apparent position; etc.

Mathematics. To understand some basic calculations, you need only high-school mathematics. Be ready to deal with a lot of trigonometric functions and analytical geometry. But to deal with very high precision calculations, you may need to know more on calculus and astrophysics.

Where can you learn these concepts? Look for some good books or Internet resources related to this subject. I myself have learned a lot from [PDS] and various Internet resources. To learn spherical astronomy seriously, you can start with [WMS].

3. Guides for the Algorithms

3.1. Notes on Accuracy

Decide carefully the level of precision you need based on the purpose of your implementation. If you just want to use it to draw a simple sundial, then low-precision calculations are sufficient. But if you want to launch a rocket to the Moon or to visit your ET friends many light years away, then stick to the best algorithms known to human beings.

You may wonder, why don't we just go straight to the 'exact'/'perfect'-precision implementation? With such an implementation, we can then extend it anyhow we like. It's true, but there are some pragmatic issues. First, the state-of-the-art theories generally involve a lot of tedious calculations. E.g., ELP2000-82B is one of the best theories for calculating the position of the Moon. To calculate the position of the Moon in any given instant, you need to perform over 30 thousand trigonometric calculations. So... Second, storing such information takes space. You may not want to include such data especially for a portable device. Third, implementing the best solutions involves much more efforts in implementation and testing.

We will briefly mention the kind of applications you can do for calculations with different levels of precision. An interesting trade-off on the precision is to implement some calculations with different levels of precisions. Use only the one you need for your applications. E.g., to predict solar eclipse, you can first estimate the date and time a solar eclipse may occur by using low-precision (but faster) calculations. Then, use high-precision (but slower) algorithms to calculate the exact times related to the eclipse.

Last but not least, you need to know the maximum error introduced by each algorithm. With this, you can tell when the algorithm fails. For example, if an algorithm tells you that a new Moon starts from 11.59pm tonight, but its precision is low and may incur up to 10 minutes of maximum error in time, then you can't tell whether a new Moon starts today or tomorrow. Such a case may not affect you if you want to observe other stars on a new Moon night, but it can be extremely important if you want to calculate an astronomical lunar calendar. To resolve such a case, you need to resort to high-precision calculations.

3.2. Low-Precision Algorithms

To find out the rough locations of the planets at night, to design a normal sundial, to include a cute moon phase indicator in your PDA, to design a compass based on the sun and moon positions, or to know the rough rising and setting time of the sun and the moon, low-precision calculations are sufficient.

There are a few good places to learn such calculations. The book [PDS] covers a lot of useful concepts and low-precision algorithms. The following websites contain good resources too: Approximate astronomical positions, How to compute planetary positions.

With this level of precession, you normally don't need to distinguish between the TT (Terrestrial Time) and UTC (Coordinated Universal Time), among the various reference frames, and you need only simple linear interpolation.

3.3. High-Precision Algorithms

To calculate an accurate birth chart for the modern western astrology or chinese astrology, to calculate lunar and solar eclipses with satisfactory results, or to calculate a calendar that is based on astronomical events, you may need high-precision algorithms. In particular, knowing the maximum possible error that can incur could be critical for some applications.

A good book that includes such algorithms is [JM]. The book introduces various important concepts. The calculations of the planet positions are based on an bridged version of the VSOP87 theory, which is one of the best planetary theories we have. While the calculations of the moon position is based on an abridged version of the famous ELP 2000/82 lunar solution.

You need to distinguish among the time systems here. Don't get confused between TT and UTC. You need to implement the routines to convert between them. A proper interpolation may be necessary to obtain good conversion results. You also need to distinguish among the true, mean, and apparent positions. Be careful to deal with nutation, precession, aberration, parallax, refraction, etc.; and convert among the different reference frames properly.

3.4. "Platinum" Algorithms

If the calculations above don't suit your purpose because you want to design a bullet-proof astrological birth chart, to calculate astronomical calendars, to predict eclipses to an accuracy in the order of milliseconds, or to navigate your spacecraft to another galaxy, then you need to study more and implement very high-precision algorithms.

You can find a lot of great resources from the wonderful Internet. Some useful web sites are listed here.

• JPL Solar System Dynamics. A lot of great information including very high-precision ephemerides (DE405/DE406, probably the most precise ephemerides on Earth!), an online HORIZONS ephemeride system, etc.
• Institut de mecanique celeste et de calcul des ephemerides. Ephemeris generator, resources on the complete VSOP87 and ELP 2000-82B theories.
• International Earth Rotation & Reference Systems Service. Stuffs related to the Earth rotation such as the Precession-Nutation Model, Delta-T, announcements of leap seconds, etc.
• International Terrestrial Reference Frame.

There are more things you need to take care of. A cubic spline interpolation for Delta-T (= TT - UT1) would be great. You may also want to learn more about astrodynamics and numerical integration.

4. References

[JM] Jean Meeus. Astronomical Algorithms. 2nd Edition. December 1998. ISBN 0-943-39661-1.
[PDS] Peter Duffett-Smith. Practical Astronomy with Your Calculator. Third Edition. March 1989. ISBN 0-521-33699-7.
[WMS] W. M. Smart. Textbook on Spherical Astronomy. 6th Edition (Revised by G.M. Green). 1977. ISBN 0-521-21516-1/0-521-29180-1.
Also take a look at the links I've bookmarked under the Astronomy category on this site.

1. What is This?

calcn includes some handy C functions for converting a date in Chinese calendar to a date in Western calendar, and vice versa. It uses a small lookup table to achieve space as well as time efficient conversion. No floating point operation is used in the conversion. Thus the code can be used by many types of platforms.

2. Source Code

calcn-0.5.tar.bz2 (Beta Version)

There are only 3 files in the package.

• calcn.h: The header file. Include this to use the conversion functions.
• calcn.c: Implementation of the conversion functions. Link to this during the compilation.
• calcn_test.c: A simple test program that makes use of all functions from calcn.c. A good place to learn all functions within 3 minutes if you are lazy to read the documentation. Just do a gcc calcn_test.c calcn.c -o calcn_test to get start.

The code should be compilable by any C compiler that is compliant with the C89 or the C99 standard. (It has been tested to be compiled with gcc -c -O3 -Wall -W -ansi -pedantic calcn.c and gcc -c -O3 -Wall -W -std=c99 -pedantic calcn.c without any warning and error.)

The code is distributed under the terms of LGPL.

3. Sample Test Session

This is a sample session from the test program. Note that the test program uses UTF-8 encoding to display some Chinese characters.

[1] Chinese -> Gregorian
[2] Gregorian -> Chinese
[3] Chinese year -> Gregorian year
[4] Gregorian year -> Chinese year
[q] Quit
? 1
Cycle (1-)? 78
Year in the cycle (1-60)? 23
Month (1-12)? 10
A leap month (0=no; 1=leap)? 0
Day (1-30)? 18
Western: 2006/12/8

[1] Chinese -> Gregorian
[2] Gregorian -> Chinese
[3] Chinese year -> Gregorian year
[4] Gregorian year -> Chinese year
[q] Quit
? 2
Year (1929-2043)? 2006
Month (1-12)? 12
Day (1-31)? 8
Chinese: 第78元运丙戌年十月十八日

[1] Chinese -> Gregorian
[2] Gregorian -> Chinese
[3] Chinese year -> Gregorian year
[4] Gregorian year -> Chinese year
[q] Quit
? 3
Cycle (1-)? 78
Year in the cycle (1-60)? 23
Gregorian year: 2006

[1] Chinese -> Gregorian
[2] Gregorian -> Chinese
[3] Chinese year -> Gregorian year
[4] Gregorian year -> Chinese year
[q] Quit
? 4
Year? 2006
Cycle             : 78
Year in the cycle : 23

4. How do use it?

4.1. Date Representations

For a date in Western (Gregorian) calendar:

struct calgregorian {
  int year;    // The year  [CALCN_YEAR_FROM..CALCN_YEAR_TO]
  int month;   // The month [1..12]
  int day;     // The day of the month [1..31]
};

For a date in the Chinese calendar:

struct calchinese {
  int scycle;      /* The sexagenary cycle [1..] (Year 2637 BCE as epoch) */
  int syear;       /* The year in the sexagenary cycle [1..60] */
  int month;       /* The month [1..12] */
  int leap_month;  /* 0 for a normal month; 1 for a leap/intercalary month */
  int day;         /* The day of the month [1..30] */
};

4.2. Conversion Functions

int calcn_chinese_from_gregorian(struct calchinese *p_calchinese,
                                 const struct calgregorian *p_calgregorian);
int calcn_gregorian_from_chinese(struct calgregorian *p_calgregorian,
                                 const struct calchinese *p_calchinese);
void calcn_year_chinese_from_astro(int *p_scycle, int *p_syear, const int year);
int calcn_year_astro_from_chinese(const int scycle, const int syear);

calcn_chinese_from_gregorian() converts a date in the Gregorian calendar to a date in the Chinese calendar, while calcn_gregorian_from_chinese() does the reverse. The return value has the following meanings:

• 0: The conversion is successful.
• CALCN_ERROR_OUT_OF_RANGE: The date is out of the range that can be handled by the program. Only the Gregorian year between CALCN_YEAR_FROM and CALCN_YEAR_TO (inclusive) can be handled by the function.
• CALCN_ERROR_INVALID_DATE: The date is invalid. It can be that the day of the month or the month is out of range, or a Chinese month has been specified as a leap month while there is no such a leap month in that year.

calcn_year_chinese_from_astro() converts an astronomical year to the Chinese year in the sexagenary cycle, while calcn_year_astro_from_chinese() does the reverse. The astronomical year -2636 is used as the epoch (the 1st year in the 1st cycle) for the conversion.

The numbering of the astronomical year is based on the familiar AD/BCE scheme, but it follows the normal integer numbering strictly. An astronomical year with a positive value is simply equal to the Western AD year. However, there is an astronomical year 0, which corresponds to 1 BCE. So, year n BCE corresponds to the -(n-1) astronomical year.

5. Known Issues

5.1. Bugs

No known outstanding bugs. Please contact me if you find any.

5.2. Limitations

The date conversion routines currently work for Gregorian years 1929-2043 only. I use my own program to generate the table and verify it against the calendar calculated by the Purple Mountain Observatory (紫金山天文台). The situation is a bit 'chaotic' before 1929 when the modern meridian of 120 degree East was not used yet. I have yet to decide some issues on this. I will add the data for future years when I have time to verify them against a reliable source.

The functions are not C++ friendly. They should be C++ friendly in the next version. (From Version 0.6 onwards.)