docs: function truephase (k, tphase) {
    // restrict tphase to (0, 1, 2, 3)
    tphase = tphase & PHASE_MASK
  
    // add phase to new moon time
    k = k + 0.25 * tphase
  
    // Time in Julian centuries from 1900 January 0.5
    const t = (1.0 / 1236.85) * k
  
    // Mean time of phase
    let pt = 2415020.75933 +
      SYNODIC_MONTH * k +
      (0.0001178 - 0.000000155 * t) * t * t +
      0.00033 * dsin(166.56 + (132.87 - 0.009173 * t) * t)
  
    // Sun's mean anomaly
    const m = 359.2242 + 29.10535608 * k - (0.0000333 - 0.00000347 * t) * t * t
  
    // Moon's mean anomaly
    const mprime = 306.0253 + 385.81691806 * k + (0.0107306 + 0.00001236 * t) * t * t
  
    // Moon's argument of latitude
    const f = 21.2964 + 390.67050646 * k - (0.0016528 - 0.00000239 * t) * t * t
  
    // use different correction equations depending on the phase being sought
    switch (tphase) {
      // new and full moon use one correction
      case NEW:
      case FULL:
        pt += (0.1734 - 0.000393 * t) * dsin(m) +
          0.0021 * dsin(2 * m) -
          0.4068 * dsin(mprime) +
          0.0161 * dsin(2 * mprime) -
          0.0004 * dsin(3 * mprime) +
          0.0104 * dsin(2 * f) -
          0.0051 * dsin(m + mprime) -
          0.0074 * dsin(m - mprime) +
          0.0004 * dsin(2 * f + m) -
          0.0004 * dsin(2 * f - m) -
          0.0006 * dsin(2 * f + mprime) +
          0.0010 * dsin(2 * f - mprime) +
          0.0005 * dsin(m + 2 * mprime)
        break
  
      // first and last quarter moon use a different correction
      case FIRST:
      case LAST:
        pt += (0.1721 - 0.0004 * t) * dsin(m) +
          0.0021 * dsin(2 * m) -
          0.6280 * dsin(mprime) +
          0.0089 * dsin(2 * mprime) -
          0.0004 * dsin(3 * mprime) +
          0.0079 * dsin(2 * f) -
          0.0119 * dsin(m + mprime) -
          0.0047 * dsin(m - mprime) +
          0.0003 * dsin(2 * f + m) -
          0.0004 * dsin(2 * f - m) -
          0.0006 * dsin(2 * f + mprime) +
          0.0021 * dsin(2 * f - mprime) +
          0.0003 * dsin(m + 2 * mprime) +
          0.0004 * dsin(m - 2 * mprime) -
          0.0003 * dsin(2 * m + mprime)
  
        // the sign of the last term depends on whether we're looking for a first
        // or last quarter moon!
        const sign = (tphase < FULL) ? +1 : -1
        pt += sign * (0.0028 - 0.0004 * dcos(m) + 0.0003 * dcos(mprime))
  
        break
    }
  
    return julian.toDate(pt)
  }