DEV Community

KunStudio
KunStudio

Posted on

Lichun Is Not Always February 4: Computing Solar Terms With Newton's Method

Most calendar code that touches East Asian astrology has the same bug hiding in it: someone wrote down that Lichun (立春, the start of spring) falls on February 4, and moved on.

It does not. In 2021 and 2025 it fell on February 3. In 2025 the exact moment was 23:10 KST — fifty minutes before midnight. Anyone born in that fifty-minute window gets a different year pillar and a different month pillar than a hardcoded table would give them. For a birth-chart engine, that is not a rounding error. That is the wrong chart.

Here is what a solar term actually is, and how to compute one properly.

A solar term is an angle, not a date

The 24 solar terms (절기 / 節氣) are not calendar conventions. Each one is the instant at which the sun's apparent ecliptic longitude reaches a multiple of 15 degrees. Lichun is 315 degrees. The vernal equinox is 0. There are 24 of them because 360 / 15 = 24.

That definition has two consequences that matter for implementation:

  1. The instant is global. It is one moment in physical time, the same for every observer on Earth. Only its local clock representation differs.
  2. The instant does not land on a tidy calendar boundary. It drifts year to year because the tropical year is not 365 days, and leap years shove it around.

So the question "what date is Lichun in year Y" is really "solve for the time t where apparentSolarLongitude(t) == 315, then render t in the target timezone."

That is a root-finding problem.

Step 1: apparent solar longitude

I use the low-precision solar series from Jean Meeus, Astronomical Algorithms (2nd ed.). Given a Julian Ephemeris Day:

function sunApparentLongitude(jde) {
  const T = (jde - J2000) / 36525.0;
  const L0 = rev360(280.46646 + 36000.76983 * T + 0.0003032 * T * T);
  const M  = rev360(357.52911 + 35999.05029 * T - 0.0001537 * T * T);
  const Mr = M * D2R;
  const C =
    (1.914602 - 0.004817 * T - 0.000014 * T * T) * Math.sin(Mr) +
    (0.019993 - 0.000101 * T)                    * Math.sin(2 * Mr) +
     0.000289                                    * Math.sin(3 * Mr);
  const trueLong = L0 + C;
  const omega = 125.04 - 1934.136 * T;
  return rev360(trueLong - 0.00569 - 0.00478 * Math.sin(omega * D2R));
}
Enter fullscreen mode Exit fullscreen mode

L0 is the geometric mean longitude, M the mean anomaly, C the equation of the center. The final line is what makes it apparent rather than true: -0.00569 is aberration, and the omega term is nutation in longitude. Skip those two and you are off by roughly 20 arcseconds, which is about seven minutes of time. Usually harmless, occasionally the difference between February 3 and February 4.

Step 2: Newton's method on a circle

Now solve sunApparentLongitude(jde) == targetLon. The derivative is well-behaved (the sun moves about 1 degree per day, always forward), so plain Newton converges in a handful of iterations:

function solveLongitude(targetLon, jdeGuess) {
  let jde = jdeGuess;
  for (let i = 0; i < 12; i++) {
    const lon = sunApparentLongitude(jde);
    let diff = lon - targetLon;
    while (diff >  180) diff -= 360;
    while (diff < -180) diff += 360;

    const h = 0.02;
    let d1 = sunApparentLongitude(jde + h) - sunApparentLongitude(jde - h);
    while (d1 >  180) d1 -= 360;
    while (d1 < -180) d1 += 360;

    const rate = d1 / (2 * h);
    const corr = diff / rate;
    jde -= corr;
    if (Math.abs(corr) < 1e-7) break;
  }
  return jde;
}
Enter fullscreen mode Exit fullscreen mode

The while loops wrapping diff and d1 into [-180, 180] are the entire trick, and they are the thing people get wrong. Longitude is modular. Near 0 degrees the naive difference between 359.9 and 0.1 is -359.8, which sends Newton flying off to the wrong year. Normalizing the difference rather than the values keeps the step honest across the wrap point. The same normalization is needed on the numerical derivative for the same reason.

I use a central difference (jde + h and jde - h) instead of deriving the series analytically. It costs one extra evaluation per iteration and removes a whole class of transcription errors.

Step 3: the part everyone forgets, delta-T

Newton gives an answer in Terrestrial Time. Civil clocks run on Universal Time. The gap between them, ΔT = TT − UT, is not a constant — it tracks the irregular rotation of the Earth, and it is currently around 69 seconds and growing.

Sixty-nine seconds sounds ignorable until the term lands at 23:59:30 local. Then it moves the date by a day, and the date is the thing the chart is keyed on. So it gets applied as a piecewise polynomial in the year:

if (year >= 2005 && year < 2050) {
  const t = y - 2000;
  dt = 62.92 + 0.32217 * t + 0.005589 * t * t;
}
Enter fullscreen mode Exit fullscreen mode

with separate branches for earlier epochs, then:

const dt = deltaTSeconds(approx.year, approx.month);
const jdKST = jde - dt / 86400 + 9 / 24;   // TT -> UT -> KST
Enter fullscreen mode Exit fullscreen mode

Subtract ΔT to get UT, add 9 hours for KST, and only then convert to a calendar date.

Note the ordering. Converting to a date first and adjusting afterwards is the bug that produces off-by-one days at year boundaries.

Step 4: guess the right year

One more trap. Seeding Newton from a day-of-year approximation can converge into the neighbouring year's occurrence of the same longitude, because 315 degrees happens every year. The fix is to probe and verify:

for (const probeYear of [year - 1, year, year + 1]) {
  let approxDoy = 79 + (lon / 360) * 365.2422;
  if (approxDoy > 365.2422) approxDoy -= 365.2422;
  const jde = solveLongitude(lon, gregorianToJD(probeYear, 1, 0.5) + approxDoy);
  // ...convert to KST...
  if (g.year !== year) continue;   // converged into a neighbour, try again
  return g;
}
Enter fullscreen mode Exit fullscreen mode

The 79 offset is the approximate day-of-year of the vernal equinox, which is where ecliptic longitude 0 sits.

Validation

Claiming astronomical accuracy without checking it is how these bugs survive. I ran all 12 major terms across 1920–2050 against Korea's public calendar API and the dates match exactly across the whole range.

On precision: the date is exact throughout. The minute can differ by up to roughly 15 minutes against a full VSOP87 solution, because the Meeus low-precision series carries about 0.01 degrees of error. For pillar boundaries that is irrelevant — the boundary is evaluated at day granularity. If you are building something that needs the minute, use the full series.

Why bother

The honest answer is that the hardcoded version is right most of the time. Lichun is February 4 in the majority of years, so the table passes casual testing and fails silently for the people born on the edge — who are exactly the people most likely to notice their chart is wrong.

Roughly 130 lines of astronomy replaces the table and removes the failure mode. It is self-contained, has no dependencies, and runs fine at the edge; ours executes inside a Cloudflare Function.

I build Cheonmyeongdang, a Korean Saju (four pillars) reading service, and this calculation sits under every chart it produces. If you want to see the boundary behaviour without installing anything, the English chart calculator is free and needs no signup — try a birth date on February 3 or 4 of 2021 or 2025 and watch the month pillar change.

The general lesson transfers beyond astrology: when a date in your domain is defined by a physical event rather than a calendar rule, storing it as a constant is a bug with a delayed fuse.

Top comments (0)