docs: function phase (phase_date) {
    if (!phase_date) {
      phase_date = new Date()
    }
    phase_date = julian.fromDate(phase_date)
  
    const day = phase_date - EPOCH
  
    // calculate sun position
    const sun_mean_anomaly =
      (360.0 / 365.2422) * day +
      (ECLIPTIC_LONGITUDE_EPOCH - ECLIPTIC_LONGITUDE_PERIGEE)
    const sun_true_anomaly =
      2 * todeg(Math.atan(
        Math.sqrt((1.0 + ECCENTRICITY) / (1.0 - ECCENTRICITY)) *
        Math.tan(0.5 * kepler(sun_mean_anomaly, ECCENTRICITY))
      ))
    const sun_ecliptic_longitude =
      ECLIPTIC_LONGITUDE_PERIGEE + sun_true_anomaly
    const sun_orbital_distance_factor =
      (1 + ECCENTRICITY * dcos(sun_true_anomaly)) /
      (1 - ECCENTRICITY * ECCENTRICITY)
  
    // calculate moon position
    const moon_mean_longitude =
      MOON_MEAN_LONGITUDE_EPOCH + 13.1763966 * day
    const moon_mean_anomaly =
      moon_mean_longitude - 0.1114041 * day - MOON_MEAN_PERIGEE_EPOCH
    const moon_evection =
      1.2739 * dsin(
        2 * (moon_mean_longitude - sun_ecliptic_longitude) - moon_mean_anomaly
      )
    const moon_annual_equation =
      0.1858 * dsin(sun_mean_anomaly)
    // XXX: what is the proper name for this value?
    const moon_mp =
      moon_mean_anomaly +
      moon_evection -
      moon_annual_equation -
      0.37 * dsin(sun_mean_anomaly)
    const moon_equation_center_correction =
      6.2886 * dsin(moon_mp)
    const moon_corrected_longitude =
      moon_mean_longitude +
      moon_evection +
      moon_equation_center_correction -
      moon_annual_equation +
      0.214 * dsin(2.0 * moon_mp)
    const moon_age =
      fixangle(
        moon_corrected_longitude -
        sun_ecliptic_longitude +
        0.6583 * dsin(
          2 * (moon_corrected_longitude - sun_ecliptic_longitude)
        )
      )
    const moon_distance =
      (MOON_SMAXIS * (1.0 - MOON_ECCENTRICITY * MOON_ECCENTRICITY)) /
      (1.0 + MOON_ECCENTRICITY * dcos(moon_mp + moon_equation_center_correction))
  
    return {
      phase: (1.0 / 360.0) * moon_age,
      illuminated: 0.5 * (1.0 - dcos(moon_age)),
      age: (SYNODIC_MONTH / 360.0) * moon_age,
      distance: moon_distance,
      angular_diameter: MOON_ANGULAR_SIZE_SMAXIS / moon_distance,
      sun_distance: SUN_SMAXIS / sun_orbital_distance_factor,
      sun_angular_diameter: SUN_ANGULAR_SIZE_SMAXIS * sun_orbital_distance_factor
    }
  }