Moon illumination percent from moon phase age in Ruby

I have the following module to calculate some moon properties. Very simple according to wikipedia article. How do I calculate illumination percent though?

I can make it linear where 100% is the mid of LUNAR_MONTH. But I expect it might be incorrect to make it linear due to the oval shape of the body.

Or should it be liner?

module LunarPhase
  LUNAR_MONTH = 29.530588853
  KNOWN_NEW_MOON = DateTime.new(2000, 1, 6, 6, 14, 0).amjd.to_f

  module ClassMethods
    def percent(time = DateTime.now)
      age(time) / LUNAR_MONTH
    end

    def age(time = DateTime.now)
      (time.amjd - KNOWN_NEW_MOON).abs % LUNAR_MONTH
    end

    def illumination(time = DateTime.now)
      # this is linear but maybe we need something more complicated
      # because I see a difference with online calculators
      age_full = LUNAR_MONTH / 2
      diff = age_full - age(time)
      (age_full - diff.abs) / age_full
    end
  end

  extend ClassMethods
end

-- you can use code as WTFPL

Update: I see some calculations in Greg Miller's page but it is full of constants I don't understand. Some theory behind that would be very useful.

Solution

With a little clues from Greg Miller, I've adapted his routine and this is what I've got.

require "time"

module LunarPhase
  LUNAR_MONTH = 29.530588853
  JULIAN_YEAR = 365.25
  J2000 = DateTime.new(2000, 1, 1, 12, 0, 0).amjd.to_f
  KNOWN_NEW_MOON = DateTime.new(2000, 1, 6, 6, 14, 0).amjd.to_f

  module ClassMethods
    # credits to https://celestialprogramming.com/meeus-illuminated_fraction_of_the_moon.html
    # some more clues at https://en.wikipedia.org/wiki/Full_moon
    # rubocop:disable all
    def illumination(time = DateTime.now)
      t = centuries_since_j2000(time)

      d = rad(297.8501921 + 445267.1114034*t - 0.0018819*t*t + 1.0/545868.0*t*t*t - 1.0/113065000.0*t*t*t*t) #47.2
      m = rad(357.5291092 + 35999.0502909*t - 0.0001536*t*t + 1.0/24490000.0*t*t*t) #47.3
      mp = rad(134.9633964 + 477198.8675055*t + 0.0087414*t*t + 1.0/69699.0*t*t*t - 1.0/14712000.0*t*t*t*t) #47.4

      i = rad(180 - d*180/Math::PI - 6.289 * Math.sin(mp) + 2.1 * Math.sin(m) -1.274 * Math.sin(2*d - mp) -0.658 * Math.sin(2*d) -0.214 * Math.sin(2*mp) -0.11 * Math.sin(d)) #48.4
      (1+Math.cos(i))/2
      # rubocop:enable all
    end

    def percent(time = DateTime.now)
      age(time) / LUNAR_MONTH
    end

    def age(time = DateTime.now)
      age_from_known = time.to_datetime.amjd - KNOWN_NEW_MOON
      age_from_known += LUNAR_MONTH if age_from_known.negative?
      age_from_known % LUNAR_MONTH
    end

    def full?(time = DateTime.now)
      illumination(time) > 0.98
    end

    private

    def centuries_since_j2000(time = DateTime.now)
      (time.to_datetime.amjd.to_f - J2000) / (JULIAN_YEAR * 100)
    end

    def rad(num)
      angle = num % 360
      angle += 360 if angle.negative?
      angle * Math::PI / 180
    end
  end

  extend ClassMethods
end

As for accuracy, new moon is at

Jan 6 23:09 according to mooninfo.org
Jan 6 23:05 according to method above