//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 2.0
 *  @date    Fri Jul 29 14:33:21 EDT 2011
 *  @see     LICENSE (MIT style license file).
 *
 *  @note    Golf Ball Flight Dynamics
 */

package scalation
package dynamics

import scala.math.{cos, sin, Pi}
import scala.util.control.Breaks.{breakable, break}

import scalation.mathstat._

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ballFlight` main function is used to illustrate the `RungeKutta` 'RK' and
 *  `DormandPrince` 'DP' ODE solvers by applying them to Newton's Second Law of Motion,
 *  'f = ma = -gm'.
 *  The flight of a golf ball is simulated from impact until the ball hits the
 *  ground.  Note, a more realistic simulation would take additional forces into
 *  account: drag, lift and spin.
 *  @see home2.fvcc.edu/~dhicketh/DiffEqns/Spring11projects/Brett_Burglund_Ryan_Street/
 *  @see Diff%20Q/pdfscreen/projectoutline.pdf
 *  @see claymore.engineer.gvsu.edu/~lait/312/golfball.pdf
 *  The accuracies of 'RK' and 'DP' versus the exact solution (EX) are compared.
 *  > runMain scalation.dynamics.ballFlight
 */
@main def ballFlight (): Unit =

    val n  = 200                                          // maximum number of time points
    val tm =   5.0                                        // simulate for a maximum of tm seconds
    val g  =   9.80665                                    // gravitational force (meters/second^2)
//  val m  =  45.93                                       // mass of a golf ball in grams
    val aa =  15.00                                       // launch angle in degrees
    val ss = 100.00                                       // swing speed in miles/hour
    val sf =   1.49                                       // smash factor
    val s  = ss * sf * 1609.344 / 3600                    // initial ball speed in meters/second
    val a  = aa * Pi / 180.0                              // launch angle in radians
    val p0 = VectorD (0.0, 0.0)                           // initial position (x, y) at time t0=0
    val v0 = VectorD (s * cos(a), s * sin(a))             // initial velocity (v_x, v_y) at t0

    println ("ball speed    s  = " + s)
    println ("launch angle  a  = " + a)
    println ("ball velocity v0 = " + v0)

    // define the system of Ordinary Differential Equations (ODEs)
    def dx_dt (t: Double, x: Double) = v0(0)              // ODE 1
    def dy_dt (t: Double, y: Double) = v0(1) - g * t      // ODE 2
    val odes: Array [Derivative] = Array (dx_dt, dy_dt)

    def exactSolution (t: Double) = VectorD (v0(0) * t, v0(1) * t - .5 * g * t * t)

    val p_r = new MatrixD (n, 2); p_r(0) = p0
    val p_d = new MatrixD (n, 2); p_d(0) = p0
    val p_e = new MatrixD (n, 2); p_e(0) = p0
    val tt  = new VectorD (n);    tt(0)  = 0.0

    val dt = tm / n                                              // time step
    var t  = dt                                                  // next time point to examine

    breakable {
        for i <- 1 to n do
            p_e(i) = exactSolution (t)                           // compute new position using EX
            if p_e(i, 1) < 0.0 then { p_e(i, 1) = 0; break () }  // quit after hitting the ground

            p_r(i) = RungeKutta.integrateV (odes, p0, t)         // compute new position using RK
            p_d(i) = DormandPrince.integrateV (odes, p0, t)      // compute new position using DP

            println (s"> at t = ${"%4.1f".format(t)}, p_r = ${p_r(i)}, p_d = ${p_d(i)}, p_e = ${p_e(i)}")
            tt(i) = t
            t += dt
        end for
    } // breakable

    new Plot (tt, p_r(?, 0), p_d(?, 0), "Plot x vs. t (black-RK, red-DP)")
    new Plot (tt, p_r(?, 1), p_d(?, 1), "Plot y vs. t (black-RK, red-DP)")
    new Plot (tt, p_r(?, 1), p_e(?, 1), "Plot y vs. t (black-RK, red-EX)")
    new Plot (tt, p_d(?, 1), p_e(?, 1), "Plot y vs. t (black-DP, red-EX)")
    new Plot (p_d(?, 0), p_d(?, 1))

end ballFlight