basebot.clj

;;; Project 1, 6.001, Spring 2005

(ns basebot
  (:use clojure.test))

(in-ns 'basebot)

;;; idea is to simulate a baseball robot

;; imagine hitting a ball with an initial velocity of v 
;; at an angle alpha from the horizontal, at a height h
;; we would like to know how far the ball travels.

;; as a first step, we can just model this with simple physics
;; so the equations of motion for the ball have a vertical and a 
;; horizontal component

;; the vertical component is governed by
;; y(t) = v sin alpha t + h - 1/2 g t^2 
;; where g is the gravitational constant of 9.8 m/s^2

;; the horizontal component is governed by
;; x(t) = v cos alpha t
;; assuming it starts at the origin

;; First, we want to know when the ball hits the ground
;; this is governed by the quadratic equation, so we just need to know when 
;; y(t)=0 (i.e. for what t_impact is y(t_impact)= 0).
;; note that there are two solutions, only one makes sense physically

(defn square [x] (* x x))

;; these are constants that will be useful to us
(def gravity 9.8)  ;; in m/s
(def pi 3.14159)

;; Problem 1

(defn position [a v u t]
  (+ (* 1/2 a (square t))
     (* v t)
     u))

;; you need to complete this procedure, then show some test cases

(deftest position-test
  (are [a v u t r] (= (position a v u t) r)
       0 0 0 0 0
       0 0 20 0 20
       0 5 10 10 60
       2 2 2 2 10
       5 5 5 5 185/2))

(position-test)

;; Problem 2

(defn root1 [a b c]
  (if (zero? a) (- (/ c b))
      (let [D (- (square b) (* 4 a c))]
        (when (>= D 0)
          (/ (+ (- b) (Math/sqrt D))
             (* 2 a))))))

(defn root2 [a b c]
  (if (zero? a) (- (/ c b))
      (let [D (- (square b) (* 4 a c))]
        (when (>= D 0)
          (/ (- (- b) (Math/sqrt D))
             (* 2 a))))))

;; complete these procedures and show some test cases

(deftest root-test
  (are [a b c r] (= (root1 a b c) r)
       0 2 1 -1/2
       5 3 6 nil
       1 6 3 -0.5505102572168221
       2 -5 0 5/2
       3 0 -27 3)
  (are [a b c r] (= (root2 a b c) r)
       2 -5 0 0
       3 0 -27 -3))

(root-test)

;; Problem 3

(defn time-to-impact [vertical-velocity elevation]
  (let [valid-time (fn [t] (and (number? t) (pos? t) t))
        time1 (valid-time (root1 (* -1/2 gravity)
                                 vertical-velocity
                                 elevation))
        time2 (valid-time (root2 (* -1/2 gravity)
                                 vertical-velocity
                                 elevation))]
    (or (when (and time1 time2) (min time1 time2))
        time1 
        time2
        0)))

;; Note that if we want to know when the ball drops to a particular height r 
;; (for receiver), we have

(defn time-to-height [vertical-velocity elevation target-elevation]
  (time-to-impact vertical-velocity (- elevation target-elevation)))

;; Problem 4

;; once we can solve for t_impact, we can use it to figure out how far the ball went

;; conversion procedure
(defn degree2radian [deg] (/ (*  deg pi) 180.))

(defn travel-distance-simple [elevation velocity angle]
  (let [vx (* velocity (Math/cos angle))
        vy (* velocity (Math/sin angle))]
    (* (time-to-impact vy elevation) vx)))

;; let's try this out for some example values.  Note that we are going to 
;; do everything in metric units, but for quaint reasons it is easier to think
;; about things in English units, so we will need some conversions.

(defn meters-to-feet [m] (/ (* m 39.6) 12))

(defn feet-to-meters [f] (/ (* f 12) 39.6))

(defn hours-to-seconds [h] (* h 3600))

(defn seconds-to-hours [s] (/ s 3600))

;; what is time to impact for a ball hit at a height of 1 meter
;; with a velocity of 45 m/s (which is about 100 miles/hour)
;; at an angle of 0 (straight horizontal)
;; at an angle of (/ pi 2) radians or 90 degrees (straight vertical)
;; at an angle of (/ pi 4) radians or 45 degrees

;; what is the distance traveled in each case?
;; record both in meters and in feet

(map (fn [angle] [angle
                  (time-to-impact (* 45 (Math/sin angle)) 1)
                  (let [d (travel-distance-simple 1 45 angle)]
                    [d (meters-to-feet d)])])
     [0 (/ pi 2) (/ pi 4)])

;; Problem 5

;; these sound pretty impressive, but we need to look at it more carefully

;; first, though, suppose we want to find the angle that gives the best
;; distance
;; assume that angle is between 0 and (/ pi 2) radians or between 0 and 90
;; degrees

(def alpha-increment 0.01)

(defn find-best-angle [velocity elevation]
  (first (reduce (fn [[a1 d1 :as td1] [a2 d2 :as td2]]
                    (if (> d1 d2) td1 td2))
                  (for [angle (range 0 (/ pi 2) alpha-increment)]
                    [angle (travel-distance-simple elevation velocity angle)]))))

;; find best angle
;; try for other velocities
;; try for other heights

(find-best-angle 100 0)   ; 0.79
(find-best-angle 1 1)     ; 0.22
(find-best-angle 1 10)    ; 0.07

;; Problem 6

;; problem is that we are not accounting for drag on the ball (or on spin 
;; or other effects, but let's just stick with drag)
;;
;; Newton's equations basically say that ma = F, and here F is really two 
;; forces.  One is the effect of gravity, which is captured by mg.  The
;; second is due to drag, so we really have
;;
;; a = drag/m + gravity
;;
;; drag is captured by 1/2 C rho A vel^2, where
;; C is the drag coefficient (which is about 0.5 for baseball sized spheres)
;; rho is the density of air (which is about 1.25 kg/m^3 at sea level 
;; with moderate humidity, but is about 1.06 in Denver)
;; A is the surface area of the cross section of object, which is pi D^2/4 
;; where D is the diameter of the ball (which is about 0.074m for a baseball)
;; thus drag varies by the square of the velocity, with a scaling factor 
;; that can be computed

;; We would like to again compute distance , but taking into account 
;; drag.
;; Basically we can rework the equations to get four coupled linear 
;; differential equations
;; let u be the x component of velocity, and v be the y component of velocity
;; let x and y denote the two components of position (we are ignoring the 
;; third dimension and are assuming no spin so that a ball travels in a plane)
;; the equations are
;;
;; dx/dt = u
;; dy/dt = v
;; du/dt = -(drag_x/m + g_x)
;; dv/dt = -(drag_y/m + g_y)
;; we have g_x = - and g_y = - gravity
;; to get the components of the drag force, we need some trig.
;; let speeed = (u^2+v^2)^(1/2), then
;; drag_x = - drag * u /speed
;; drag_y = - drag * v /speed
;; where drag = beta speed^2
;; and beta = 1/2 C rho pi D^2/4
;; note that we are taking direction into account here

;; we need the mass of a baseball -- which is about .15 kg.

;; so now we just need to write a procedure that performs a simple integration
;; of these equations -- there are more sophisticated methods but a simple one 
;; is just to step along by some step size in t and add up the values

;; dx = u dt
;; dy = v dt
;; du = - 1/m speed beta u dt
;; dv = - (1/m speed beta v + g) dt

;; initial conditions
;; u_0 = V cos alpha
;; v_0 = V sin alpha
;; y_0 = h
;; x_0 = 0

;; we want to start with these initial conditions, then take a step of size dt
;; (which could be say 0.1) and compute new values for each of these parameters
;; when y reaches the desired point (<= 0) we stop, and return the distance (x)

(def drag-coeff 0.5)
(def density 1.25)  ; kg/m^3
(def mass 0.145)  ; kg
(def diameter 0.074)  ; m

(def the-beta #(* 0.5 drag-coeff density (* 3.14159 0.25 (square diameter))))

(def beta (the-beta))

(def time-epsilon 0.1)
(def angle-epsilon 0.001)
(def distance-epsilon 5) ;; use so big epsilon cause of very hungry calculations for best value

(defn integrate [x0 y0 u0 v0 dt g m beta]
  (let [speed (fn [u v] (Math/sqrt (+ (square u) (square v))))
        dx (fn [u] (* u dt))
        dy (fn [v] (* v dt))
        du (fn [u v] (* -1 dt (/ 1 m) (speed u v) beta u))
        dv (fn [u v] (* -1 dt (+ g (* (/ 1 m) (speed u v) beta v))))]
    (loop [x x0, y y0, u u0, v v0, t 0]
      (if (neg? y) {:x x, :y y, :u u, :v v, :t t}
          (recur (+ x (dx u))
                 (+ y (dy v))
                 (+ u (du u v))
                 (+ v (dv u v))
                 (+ t dt))))))

(defn travel [elevation velocity angle]
  (integrate 0 elevation
             (* velocity (Math/cos angle))
             (* velocity (Math/sin angle))
             time-epsilon gravity mass beta))

(defn travel-distance [elevation velocity angle]
  (:x (travel elevation velocity angle)))

;; RUN SOME TEST CASES

(travel-distance 1 45 (/ pi 4)) ; 92.5

;; what about Denver?

(binding [density 1.06]
  (binding [beta (the-beta)]
    (travel-distance 1 45 (/ pi 4)))) ; 100.9

;; Problem 7
 
;; now let's turn this around.  Suppose we want to throw the ball.  The same
;; equations basically hold, except now we would like to know what angle to 
;; use, given a velocity, in order to reach a given height (receiver) at a 
;; given distance

(defn travel-time [elevation velocity angle]
  (:t (travel elevation velocity angle)))

(defn best-angle [elevation velocity distance]
  (let [angles-and-times 
        (for [angle (range (/ (- pi) 2) (/ pi 2) angle-epsilon)
              :let [travel (travel elevation velocity angle)
                    travel-distance (:x travel)
                    travel-time (:t travel)]
              :when (< (Math/abs (- distance travel-distance)) distance-epsilon)]
          [angle travel-time])]
    (when (seq angles-and-times)
      (first (reduce (fn [[a1 t1 :as at1] [a2 t2 :as at2]]
                       (if (< t1 t2) at1 at2))
                     angles-and-times)))))

;; a cather trying to throw someone out at second has to get it roughly 36 m
;; (or 120 ft) how quickly does the ball get there, if he throws at 55m/s,
;;  at 45m/s, at 35m/s?

(defn travel-time-on-distance [elevation velocity distance]
  (when-let [angle (best-angle elevation velocity distance)]
    (travel-time elevation velocity angle)))

(travel-time-on-distance 1 55 36) ; 0.7
(travel-time-on-distance 1 45 36) ; 0.8
(travel-time-on-distance 1 35 36) ; 1.01
      
;; try out some times for distances (30, 60, 90 m) or (100, 200, 300 ft) 
;; using 45m/s

(travel-time-on-distance 1 45 30) ; 0.7
(travel-time-on-distance 1 45 60) ; 1.6
(travel-time-on-distance 1 45 90) ; 3.2
(travel-time-on-distance 1 35 90) ; nil

;; Problem 8

(defn half [x] (/ x 2))

(defn travel-with-bounces [elevation velocity angle bounces]
  (loop [travel (travel elevation velocity angle)
         bounce bounces]
    (if (zero? bounce) travel
        (let [{:keys [x y u v]} travel]
          (recur (integrate x 0 (half u) (half (Math/abs v)) time-epsilon gravity mass beta)
                 (dec bounce))))))

(travel-with-bounces 1 45 (/ pi 4) 0) ; 92.5
(travel-with-bounces 1 45 (/ pi 4) 1) ; 103.7
(travel-with-bounces 1 45 (/ pi 4) 2) ; 106.5