#lang racket
(require gigls/unsafe)

;;; File:
;;;   project/mandelbrot.rkt
;;; Authors:
;;;   Samuel A. Rebelsky
;;;   Logan Goldberg
;;; Summary:
;;;   A simulated CSC 151 project that illustates one approach to
;;;   making Mandelbrot sets.  Included to discourage students from
;;;   making Mandelbrot sets.
;;; References:
;;;   Complex numbers in Racket:
;;;     http://docs.racket-lang.org/reference/generic-numbers.html#%28part._.Complex_.Numbers%29
;;;   Wikipedia on the Mandelbrot set
;;;     http://en.wikipedia.org/wiki/Mandelbrot_set
;;;   The Julia Set (not implemented, but another similar approach)
;;;     http://mathworld.wolfram.com/JuliaSet.html

; +-------+----------------------------------------------------------
; | Notes |
; +-------+

; Converting points
;   Images have x coordinates that range from 0 to width-1 and y
;   coordinates that range from 0 to height-1.  The Mandelbrot set
;   normally works with x coordinates in the range -2 to 1 and y
;   coordinates in the range -1 to 1.  And the Mandelbrot computation
;   uses complex numbers.
;
;   Rather than converting directly, we do a multi-step conversion.
;
;     * We convert x coordinates to the range.  0..1 (yes, that's supposed 
;       to be a range of real numbers, not just integers).  This computation 
;       and the next allow us to do everything else assuming a consistent 
;       range that is independent of the width and height.
;     * We convert y coordinates to the range 0..1.
;     * We multiply the x coordinate by HSCALE and add HOFFSET.
;     * We multiply the y coordinate by VSCALE and add VOFFSET.
;     * We convert the new (x,y) pair to a complex number.
;
;   Constants (defined below) help with this conversion.  The constants
;   are used within make-complex, which is defined in the 
;   Utilities section.

; +-----------+------------------------------------------------------
; | Constants |
; +-----------+

;;; Constant:
;;;   MAX_ITERATIONS
;;; Type:
;;;   Integer
;;; Description:
;;;   The maximum number of iterations of the function we do before
;;;   giving up.
(define MAX_ITERATIONS 50)

; +---------+--------------------------------------------------------
; | Helpers |
; +---------+

;;; Procedure:
;;;   unit-coord-x
;;; Parameters:
;;;   x, an integer
;;;   width, an integer
;;; Purpose:
;;;   Convert x to the range 0..1.
;;; Produces:
;;;   ux, a real number
;;; Preconditions:
;;;   0 <= x < width
;;; Postconditions:
;;;   0 <= ux < 1
;;;   (* ux width) is approximately x.
(define unit-coord-x
  (lambda (x width)
    (/ x width)))

;;; Procedure:
;;;   unit-coord-y
;;; Parameters:
;;;   y, an integer
;;;   height, an integer
;;; Purpose:
;;;   Convert y to the range 0..1.
;;; Produces:
;;;   uy, a real number
;;; Preconditions:
;;;   0 <= y < height
;;; Postconditions:
;;;   0 <= uy < 1
;;;   (* uy height) is approximately y.
(define unit-coord-y
  (lambda (y height)
    (/ y height)))

;;; Procedure:
;;;   hsv2irgb
;;; Parameters:
;;;   hue, a real number
;;;   saturation, a real number
;;;   value, a real number
;;; Purpose:
;;;   Build an integer-encoded RGB color from a hue, saturation, 
;;;   and value
;;; Produces:
;;;   irgb, an integer-encoded RGB color
;;; Preconditions:
;;;   0 <= saturation <= 1
;;;   0 <= value <= 1
;;; Postconditions:
;;;   (irgb->hue irgb) is close to hue
;;;   (irgb->saturation irgb) is close to saturation
;;;   (irgb->value irgb) is close to value
;;; Philosophy:
;;;   The built-in hsv->irgb is somewhat broken and I'm too lazy
;;;   to fix it an propagate the changes everywhere.  This is a
;;;   temporary fix.
;;; Props:
;;;   Based on code developed by Janet Davis's research students.
;;;   Their names have been lost.
(define hsv2irgb
  (lambda (hue saturation value)
    (let* ([hi (mod (floor (/ hue 60)) 6)]
           [v value]
           [f (- (/ hue 60) hi)]
           [p (* value (- 1 saturation))]
           [q (* value (- 1 (* f saturation)))]
           [t (* value (- 1 (* saturation (- 1 f))))])
      (cond
        [(= hi 0) (irgb (* 255 v) (* 255 t) (* 255 p))]
        [(= hi 1) (irgb (* 255 q) (* 255 v) (* 255 p))]
        [(= hi 2) (irgb (* 255 p) (* 255 v) (* 255 t))]
        [(= hi 3) (irgb (* 255 p) (* 255 q) (* 255 v))]
        [(= hi 4) (irgb (* 255 t) (* 255 p) (* 255 v))]
        [(= hi 5) (irgb (* 255 v) (* 255 p) (* 255 q))]
        [else (irgb 0 0 0)]))))


;;; Procedure:
;;;   indexed-color
;;; Parameters:
;;;   i, an integer
;;; Purpose:
;;;   Produce a color based on an index
;;; Produces:
;;;   irgb, an integer-encoded RGB color
;;; Preconditions:
;;;   [No additional]
;;; Postconditions:
;;;   If j != i and (abs (- j i)) < MAX_ITERATIONS, then
;;;     (indexed-color i) and (index-color j) are likely to
;;;     be different.
(define index-color
  (lambda (i)
    (hsv2irgb (* i (/ 360.0 MAX_ITERATIONS))
              (+ .5 (* .1 (mod i 6)))
              (+ .5 (* .05 (mod i 11))))))

;;; Procedure:
;;;   complex-distance-squared
;;; Parameters:
;;;   c1, a complex number
;;;   c2, a complex number
;;; Purpose:
;;;   Computes the square of the distance from the point represented
;;;   by c1 to the point represented by c2
;;; Produces:
;;;   distance-squared, a real number
;;; Preconditions:
;;;   [No additional]
;;; Postconditions:
;;;   distance-squared is the distance using the standard formula.
;;; Ponderings:
;;;   We compute the squared distance, rather than the distance, because
;;;   computing distance normally requires computing a square root, and
;;;   that's likely to be expensive.
(define complex-distance-squared
  (lambda (c1 c2)
    (+ (square (- (real-part c1) (real-part c2)))
       (square (- (imag-part c1) (imag-part c2))))))

;;; Procedure:
;;;   steps-to-escape
;;; Parameters:
;;;   c, a complex number.
;;; Purpose:
;;;   Counts the number of times we can repeatedly compute z[i+1] =
;;;   z[i]*z[i]+c before hitting MAX_ITERATIONS or getting a number 
;;;   that is more than 2 away from the complex origin.  If we hit
;;;   MAX_ITERATIONS, we return a special value (-1).
;;; Produces:
;;;   s, an integer
;;; Preconditions:
;;;   Let f(z) be z*z+c.  (We use this in the postconditions.)
;;; Postconditions:
;;;   If s >= 0, then c "escapes" after exactly s iterations.  That is,
;;;     distance((f^s)(c), 0) >= 2 and for all i, 0 <= i < s, 
;;;     distance((f^i)(c), 0) < 2.
;;;   If s = -1, then c does not "escape" within MAX_ITERATIONS. That is,
;;;     for all i, 0 <= i <= MAX_ITERATIONS, distance((f^i)(c), 0) < 2.
(define steps-to-escape
  (lambda (c)
    (let kernel ([z c]
                 [s 0])
      (cond
        [(>= (complex-distance-squared z 0) 4)
         s]
        [(>= s MAX_ITERATIONS)
         -1]
        [else
         (kernel (+ (* z z) c) (+ s 1))]))))

; +--------------------+---------------------------------------------
; | Primary Procedures |
; +--------------------+

; An version of image-series that generates various portions of the
; mandelbrot set.  Not documented with the six P's because that's
; part of the project, and this is sample code for the project.
(define mandelbrot-image-series
  (lambda (n width height)
    (let* (; Our default color.  Used when we hit MAX_ITERATIONS.
           [DEFAULT (irgb 0 0 0)]
           ; The horizontal offset of the center from top-left.
           ; Use -2.0 for normal.  Might be based on n.
           [HOFFSET (+ -2.0 (* 0.5 (mod n 3)))]
           ; The vertical offset of the center from top-left.
           ; Use -1.0 for normal.  Might be based on n.
           [VOFFSET (+ -1.0 (* 0.25 (mod n 5)))]
           ; The horizontal scale.  Might be based on n
           [HSCALE (- 1.0 HOFFSET)]
           ; The Vertical scale.  Might be based on n.
           [VSCALE (- 1.0 VOFFSET)])
      (let (; Convert a point in the unit plane to a complex number within
            ; some more interesting range.  That range is determined by the
            ; parameters above.
            [complicate (lambda (x y)
                          (+ (+ HOFFSET (* HSCALE x))
                             (* 0+i (+ VOFFSET (* VSCALE y)))))])
        (image-compute
         (lambda (col row)
           (let* ([c (complicate (unit-coord-x col width) 
                                 (unit-coord-y row height))]
                  [s (steps-to-escape c)]
                  [color (if (= s -1)
                             DEFAULT
                             (index-color (+ n s)))])
             color))
         width height)))))

; The image-series procedure, which is described elsewhere.
(define image-series mandelbrot-image-series)