#lang racket/base

(require racket/math racket/vector racket/match
         "./calendar.rkt")

(provide lunar-next-phase)

;; meeus, ch 47

(define (lunar-mean-phase+ k T)
  (+ 2451550.09766
     (* 29.530588861 k)
     (* 0.00015437 T T)
     (- (* 0.000000150 T T T))
     (* 0.00000000073 T T T T)))

;; phase must correspond to the decimal part of k
;; phase is one of 'new 'first 'full 'last
(define (lunar-precise-phase+ k phase)
  (define T (/ k 1236.85))
  (define lmm (lunar-mean-phase+ k T))

  ;; corrections
  (define E (- 1 (* 0.002516 T) (* 0.0000074 T T)))
  (define (d->r val)
    (degrees->radians val))

  (define M (d->r (+ 2.5534 (* 29.10535670 k)
                     (* -0.0000014 T T)
                     (* -0.00000011 T T T))))
  (define M* (d->r (+ 201.5643 (* 385.81693528 k)
                      (* 0.0107582 T T)
                      (* 0.00001238 T T T)
                      (* -0.000000058 T T T T))))
  (define F (d->r (+ 160.7108 (* 390.67050284 k)
                     (* -0.0016118 T T)
                     (* -0.00000227 T T T)
                     (* 0.000000011 T T T T))))
  (define Ω (d->r (+ 124.7746 (* -1.56375588 k)
                     (* 0.0020672 T T)
                     (* 0.00000215 T T T))))

  (define A
    (vector-map
      d->r
      (vector (+ 299.77 (*  0.107408 k) (* -0.009173 T T))
              (+ 251.88 (*  0.016321 k))
              (+ 251.83 (* 26.651886 k))
              (+ 349.42 (* 36.412478 k))
              (+  84.66 (* 18.206239 k))
              (+ 141.74 (* 53.303771 k))
              (+ 207.14 (*  2.453732 k))
              (+ 154.84 (*  7.306860 k))
              (+  34.52 (* 27.261239 k))
              (+ 207.19 (*  0.121824 k))
              (+ 291.34 (*  1.844379 k))
              (+ 161.72 (* 24.198154 k))
              (+ 239.56 (* 25.513099 k))
              (+ 331.55 (*  3.592518 k)))))

  (define periodic-multipliers
    (match phase
      ['new
       (vector    -0.40720
               (* +0.17241 E)
                  +0.01608
                  +0.01039
               (* +0.00739 E)
               (* -0.00514 E)
               (* +0.00208 E E)
                  -0.00111
                  -0.00057
               (* +0.00056 E)
                  -0.00042
               (* +0.00042 E)
               (* +0.00038 E)
               (* -0.00024 E)
                  -0.00017
                  -0.00007
                  +0.00004
                  +0.00004
                  +0.00003
                  +0.00003
                  -0.00003
                  +0.00003
                  -0.00002
                  -0.00002
                  +0.00002)]
      ['full
       (vector    -0.40614
               (* +0.17302 E)
                  +0.01614
                  +0.01043
               (* +0.00734 E)
               (* -0.00515 E)
               (* +0.00209 E E)
                  -0.00111
                  -0.00057
               (* +0.00056 E)
                  -0.00042
               (* +0.00042 E)
               (* +0.00038 E)
               (* -0.00024 E)
                  -0.00017
                  -0.00007
                  +0.00004
                  +0.00004
                  +0.00003
                  +0.00003
                  -0.00003
                  +0.00003
                  -0.00002
                  -0.00002
                  +0.00002)]
      [(or 'first 'last)
       (vector    -0.62801
               (* +0.17172 E)
               (* -0.01182 E)
                  +0.00862
                  +0.00804
               (* +0.00454 E)
               (* +0.00204 E E)
                  -0.00180
                  -0.00070
                  -0.00040
               (* -0.00034 E)
               (* +0.00032 E)
               (* +0.00032 E)
               (* -0.00028 E E)
               (* +0.00027 E)
                  -0.00017
                  -0.00005
                  +0.00004
                  -0.00004
                  +0.00004
                  +0.00004
                  +0.00004
                  +0.00002
                  +0.00002
                  -0.00002)]))
  (define periodic-sine-terms
    (match phase
      [(or 'new 'full)
       (vector    M*
                  M
               (* M* 2)
               (* F 2)
               (- M* M)
               (+ M* M)
               (* M 2)
               (- M* (* 2 F))
               (+ M* (* 2 F))
               (+ (* 2 M*) M)
               (* 3 M*)
               (+ M (* 2 F))
               (- M (* 2 F))
               (- (* 2 M*) M)
                  Ω
               (+ M* (* 2 M))
               (- (* 2 M*) (* 2 F))
               (* 3 M)
               (+ M* M (* -2 F))
               (+ (* 2 M*) (* 2 F))
               (+ M* M (* 2 F))
               (+ M* (- M) (* 2 F))
               (- M* M (* 2 F))
               (+ (* 3 M*) M)
               (* 4 M*))]
      [(or 'first 'last)
       (vector    M*
                  M
               (+ M* M)
               (* 2 M*)
               (* 2 F)
               (- M* M)
               (* 2 M)
               (- M* (* 2 F))
               (+ M* (* 2 F))
               (* 3 M*)
               (- (* 2 M*) M)
               (+ M (* 2 F))
               (- M (* 2 F))
               (+ M* (* 2 M))
               (+ (* 2 M*) M)
               Ω
               (- M* M (* 2 F))
               (+ (* 2 M*) (* 2 F))
               (+ M* M (* 2 F))
               (- M* (* 2 M))
               (+ M* M (* -2 F))
               (* 3 M)
               (- (* 2 M*) (* 2 F))
               (+ M* (- M) (* 2 F))
               (+ (* 3 M*) M))]))

  (define periodic-correction
    (for/sum ([multiplier (in-vector periodic-multipliers)]
              [sine-term (in-vector periodic-sine-terms)])
      (* multiplier (sin sine-term))))

  (define W (+ 0.00306
               (* -0.00038 E (cos M))
               (*  0.00026 (cos M*))
               (* -0.00002 (cos (- M* M)))
               (*  0.00002 (cos (+ M* M)))
               (*  0.00002 (cos (* 2 F)))))
  (define phase-correction
    (match phase
      ['first W]
      ['last (- W)]
      [(or 'new 'full) 0]))

  (define additional-multipliers
    (vector 0.000325
            0.000165
            0.000164
            0.000126
            0.000110
            0.000062
            0.000060
            0.000056
            0.000047
            0.000042
            0.000040
            0.000037
            0.000035
            0.000023))
  (define additional-correction
    (for/sum ([Ai (in-vector A)]
              [multipler (in-vector additional-multipliers)])
      (* multipler (sin Ai))))

  ; (define r->d radians->degrees)
  ; (printf "lmm ~a E ~a M ~a M* ~a F ~a Ω ~a\n" lmm E (r->d M) (r->d M*) (r->d F) (r->d Ω))
  ; (printf "C1 ~a C2 ~a C3 ~a\n" periodic-correction phase-correction additional-correction)

  (+ lmm periodic-correction phase-correction additional-correction))

;; 1-indexed again
(define (lunar-next-phase jd phase)
  (define fraction
    (match phase
      ['new   0.00]
      ['first 0.25]
      ['full  0.50]
      ['last  0.75]))

  (define seekb -3)
  (define seekf 10)

  (define-values [year month day] (jd->date-values jd))
  (define k~ (let ([k~ (* (- (+ year (/ (sub1 month) 12)) 2000) 12.3685)])
               (+ (floor k~) seekb fraction)))

  (define trials (for/vector #:length 10 ([offset (in-range seekb seekf)])
                   (lunar-precise-phase+ (+ k~ offset) phase)))

  (for/first ([x (in-vector trials)]
              #:when (>= x jd))
             x))