orrery.clj

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; A simple orrery in clojure.  Models a 2-dimensional,
;; orthogonal solar system, with the planets most definitely
;; not to scale.  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(ns orrery
  (:import [java.awt Color Dimension]          
           [java.awt.event ActionListener]
           [javax.swing JFrame JPanel Timer]))

;Some constants you may want to play with
(def img-dims { :h 600, :w 600 })
(def frame-delay 20)
(def scale-factor 22)
(def home-planet "Earth") ;to view skymap from a different perspective

;------------------------------------------------------------
; Orbital calculations
;------------------------------------------------------------

(defn mean-anomaly 
  "t = time since perihelion; P = orbital period"
  [t P]
  (/ (* 2 Math/PI t) P))

(defn ecc-anomaly  
  "M = mean anomaly; e = eccentricity"
  [M e]
  ;approximate iteratively
  (loop [lastval M, newval (+ M (* e (Math/sin M)))]
    (if 
        (<= (Math/abs (- newval lastval)) 0.00001) newval
        (recur newval (+ M (* e (Math/sin newval)))))))

(defn true-anomaly 
  "Calculates the true anomaly as an angle in radians. 
   E = eccentric anomaly, e = eccentricity of orbit"
  [E e]
  (* 2 
     (Math/atan (* (Math/sqrt (/ (+ 1 e) (- 1 e)))
                   (Math/tan (/ E 2))))))

(defn radial-distance 
  "ta - true anomaly; a - semimajor axis; e - eccentricity
   Calculates the radial distance from the sun given the true anomaly in 
   radians"
  [ta a e]
  (* a (/ (- 1 (Math/pow e 2))
          (+ 1 (* e (Math/cos ta))))))

(defn angle-from-sides
  "law of cosines/determines angle C"
  [a b c] 
  (if (or (= a 0) (= b 0) (= c 0)) 0
    (Math/acos (/ ( - (+ (Math/pow a 2) 
                         (Math/pow b 2))
                      (Math/pow c 2))
                  (* 2 a b)))))

(defn get-coords-radial 
  "t = time since perihelion, a = semimajor axis, P = orbital period, 
   e = eccentricity"
  [t a e P]
  (let [theta (true-anomaly (ecc-anomaly (mean-anomaly t P) e) e)
        r (radial-distance theta a e)]
    [r theta]))
         
(defn radial-to-rect
  ([m a]    [ (* m (Math/cos a)) (* m (Math/sin a)) ])
  ([coords] (radial-to-rect (coords 0) (coords 1))))

;------------------------------------------------------------
; Planet definitions
;------------------------------------------------------------

;a = semi-major axis
;e = eccentricity
;P = orbital period; this is in days in the demo
;t = current time offset
;r = planet radius
;b = semi-minor axis; calculated.
;l = linear eccentricity; calculated.
;x = current-x position; calculated.
;y = current-y position; calculated.
(defstruct orbital-body :name :a :e :P :t :r :color :b :l :x :y)

(defn set-semiminor-axis! [planet] ;:a and :e must be set
  (reset! planet (assoc @planet
                        :b
                        (Math/sqrt (* (Math/pow (:a @planet) 2)
                                      (- 1 (Math/pow (:e @planet) 2)))))))

(defn set-linear-ecc! [planet] ;:b must be set
  (reset! planet (assoc @planet 
                        :l 
                        (Math/sqrt (- (Math/pow (:a @planet) 2)
                                      (Math/pow (:b @planet) 2))))))

(defn init-planets [planets]
  (doseq [planet planets]
    (set-semiminor-axis! planet)
    (set-linear-ecc! planet)))

;loosely from http://en.wikipedia.org/wiki/Attributes_of_the_largest_solar_system_bodies
(defn make-planets []
  (let [result 
        [(atom (struct orbital-body "Sol"     0    0          1      0 14 (Color. 255 255 0) 0 0 0 0))  
         (atom (struct orbital-body "Mercury" 0.77 0.20563069 86.704 0 3  (Color. 200 35 35) 0 0 0 0))
         (atom (struct orbital-body "Earth"   2    0.01671022 360    0 8  (Color. 0 200 255) 0 0 0 0))   
         (atom (struct orbital-body "Venus"   1.44 0.00677323 221.47 0 7  (Color. 100 150 200) 0 0 0 0))
         (atom (struct orbital-body "Mars"    3.04 0.09341233 677.11 0 5  (Color. 255 0 0) 0 0 0 0))
         (atom (struct orbital-body "Jupiter" 10.4 0.04839266 4270.5 0 12 (Color. 200 150 50) 0 0 0 0))
         ]]
    (do (init-planets result)
        result)))

;------------------------------------------------------------
; Display some planets in orbit
;------------------------------------------------------------

(defn draw-skymap-planet [g home target]
  ;determine the sky angle between the home and target planets
  (let [c (- (:y home) (:y target))
        b (- (:x home) (:x target))
        a (Math/sqrt (+ (Math/pow c 2) (Math/pow b 2)))
        phi (angle-from-sides a b c)
        H (if (> c 0) phi ;0-PI degrees
              (- (* Math/PI 2) phi))
        ]
    (doto g
      (.setColor (:color target))
      (.fillOval (/ (* H (:w img-dims)) (* Math/PI 2))
                 ( - 25 (/ (:r target) 2))
                 (:r target)
                 (:r target)))))

(defn draw-skymap [g planets]
  (let [signs 12
        unit-width (/ (:w img-dims) signs)
        home (first (filter #(= (:name @%) home-planet) planets))]
    (do          
      ;draw the constellation lines
      (. g setColor (Color. 200 200 200))
      (doseq [r (range 0 signs)] 
        (. g drawRect 
           (* r unit-width)
           0 
           unit-width 
           50))

      ;plot each planet in the map based on the angle relative to home
      (doseq [planet (filter #(not (identical? home %)) planets)]
        (draw-skymap-planet g @home @planet)))))

(defn draw-orbit-planet [g {:keys [a b color r l x y] :as planet}]  
  (doto g
    ;draw the orbital path
    (.setColor (Color. 50 50 50))
    (.drawOval (- (/ (:w img-dims) 2) 
                  (* a scale-factor)
                  (* l scale-factor))    ;x
               (- (/ (:h img-dims) 2)
                  (* b scale-factor))    ;y
               (* 2 a scale-factor)      ;w
               (* 2 b scale-factor))     ;h

    ;plot the planet
    (.setColor (:color planet))
    (.fillOval (- x (/ r 2)) (- y (/ r 2)) r r)
    ))

(defn draw-orbits [g planets]
  (doseq [planet planets]
    (draw-orbit-planet g @planet)))

(defn update-planets!
"planets - a collection of planet atoms
Increments the planets' time counter and x, y positions in the orbital 
plane"
  [planets]
  (doseq [planet planets]
    (let [p @planet
          new-time (if (>= (:t p) (:P p)) 
                     (- (:t p) (:P p)) ;0 
                     (inc (:t p)))
          new-pos (map #(+ (/ (:h img-dims) 2) (* scale-factor %))
                       (radial-to-rect (get-coords-radial (:t p) 
                                                          (:a p)
                                                          (:e p)
                                                          (:P p))))
          ]
      (reset! planet (assoc p 
                       :t new-time
                       :x (nth new-pos 0)
                       :y (nth new-pos 1))))))

(defn make-canvas [frame planets]
  (proxy [JPanel ActionListener] []
    (paintComponent [g]
      (proxy-super paintComponent g)
      
      (doto g
        ;background
        (.setColor (Color. 0 20 20)) 
        (.fillRect 0 0 (:w img-dims) (:h img-dims)))
  
      ;paint the planets
      (draw-orbits g planets)
      (draw-skymap g planets))

    (actionPerformed [e]
      (update-planets! planets)
      (.repaint this))
    (getPreferredSize []
      (Dimension. (:w img-dims) (:h img-dims)))))

(defn main [& args]
  ;load and display the canvas
  (let [frame (JFrame. "Orrery")
        planets (make-planets)
        canvas (make-canvas frame planets)
        timer (Timer. frame-delay canvas)]
    (. canvas setFocusable true)
    (doto frame
      (.setSize (:w img-dims) (:h img-dims))
      (.add canvas)
      (.pack)
      (.setVisible true))
    (.start timer)
  planets))

(main)