Introduction

Recently I posted about how to calculate the day of the week given the Julian Day Number and gave a reference to a Forth implementation, without giving the source code.

Why are Julian days (not “dates”, this refers to the Julian calendar, we’re talking about the astronomical time system here) useful?

  • It’s universal
  • It’s ignorant about leap second
  • However, conversion to the Gregorian date works just fine and even takes leap days into account
  • It’s proleptic, means, it even works before 1582 which might be kind of useful to other people than astronomers
  • You can easily do arithmetics since a Julian Day is nothing but the elapsed
  • time since January 1st, 4713 BC; Examples would be my mentioned modulo * thingy to calculate the day of the week or holiday calculations (“Rosenmontag”, the most important day of carnival here in Düsseldorf is Easter Sunday -48 days)
  • It can act as an intermediate to convert between different calendar systems

However, mind you, it does not correct for leap seconds; the Time Scale (e.g. UTC vs TT) is not defined, so you would need to correct for leap seconds yourself.

Here an example using Python and astropy to show what I mean with this:

In [23]: t = Time("2023-01-31", scale="utc")    # Set date using the UTC scale

In [24]: t, t.tt
Out[24]: 
(<Time object: scale='utc' format='iso' value=2023-01-31 00:00:00.000>,
 <Time object: scale='tt' format='iso' value=2023-01-31 00:01:09.184>)

So, mind your leap seconds and other secular effects.

Forth Usage Examples

Get Julian Day Number

Example: 2000-01-01

2000 1 1 jdn . 2451545 ok

Get Julian Day as fraction

The jd function will always return Midnight of this day, means the civil time when the day starts. The day must be given as a fraction, for example, when it’s day 25 at noon, one must pass 25.5 as the day.

Example: 2000-01-01 00:00:00:

2000 1 1e jd f.  2451544.5  ok

Example: 2000-01-01 at 13:00:00 (UTC); I use the supplied hms2dayfrac routine to convert HH:MM:SS to the fraction of a day.

2000 1 13 00 00 hms2dayfrac jd f.  2451544.04166667  ok

Source code

\ Julian Day conversion routines

\ if m < 2, then it's the preceeding year and
\ m=m+2
: jdn_s1 ( r r r -- r r r)
    { y m d }
    m 2 < if
        y 1 - ( y -- )
        m 12 + ( y m -- )
    else
        y
        m
    endif
    d
;

\ A = y / 100
\ B = 2 - A + ( A / 4 )
: jdn_s2 ( Y -- B )
  >l 2 @local0 100 / dup -rot - swap 4 / + lp+ ;

\ JD = INT(365.25 * ( Y + 4716 )) +
\      INT(30.6001 * ( M + 1 )) +
\      D + B - 1524.5
: jdn_s3
    { m d y b }
    y 4716 + s>f 365.25e f* f>s
    m 1 + s>f 30.6001e f* f>s
    +
    d +
    b +
    s>f
    1524.5e f-
;

\ Julian day number, Main routine, only accepts full Integer dates
: jdn
    { y m d }
    y m d jdn_s1
    rot
    dup
    jdn_s2
    jdn_s3
    .5E f+ f>s
;

\ jd Julian Day Main routine, accepts fractions of dates
: jd
  { y m F: d }
  y m d f>s jdn        \ get Julian Day Number (noon)
  s>f 0.5E f-          \ convert to midnight previous day
  d fdup floor f-      \ get fractional part of d
  fdup                 ( jdn frac frac)
  .5E f< IF
     f-                \ subtract fraction from noon if before noon
  ELSE
     f+                \ but add noon and fractional if not
  ENDIF
;

\ Convert H M S to fraction of a day.
\ Mind you, I noticed precision errors in gforth; or perhaps I have a bug
\ in my routine, I'm not quite sure and I need more coffee.
: hms2dayfrac
    { h m s }
    s s>f 216000e f/
    m s>f 1440e f/ f+
    h s>f f+
    24e f/
;