Projectile.hs

{-# OPTIONS_GHC -Wall #-}

module Main where

import Graphics.Gnuplot.Simple
import Physics.Learn

--type StateTuple = (Double,Vec,Vec)
--type AccelerationFunction = StateTuple -> Vec

--eulerCromerSolution :: Double -> AccelerationFunction
--                    -> StateTuple -> StateTuple
--eulerCromerSolution

-- vertical direction is y direction
earthSurfaceGravity :: OneParticleAccelerationFunction
earthSurfaceGravity _state = vec 0 (-g) 0

g :: Double
g = 9.81

projectileTuples :: Double -> Double
                 -> OneParticleAccelerationFunction
                 -> [OneParticleSystemState]
projectileTuples v0 theta af
    = oneParticleRungeKuttaSolution af 0.01
      (0, St (cart 0 0 0) (vec vx0 vy0 0))
      where
        vx0 = v0 * cos theta
        vy0 = v0 * sin theta

yCoord :: Position -> Double
yCoord r = y
    where
      (_,y,_) = cartesianCoordinates r

inAir :: [OneParticleSystemState] -> [OneParticleSystemState]
inAir = takeWhile (\(_,St r _) -> yCoord r >= 0)

initialProjState :: Double -> Double -> OneParticleSystemState
initialProjState v0 theta
    = (0, St (cart 0 0 0) (vec vx0 vy0 0))
      where
        vx0 = v0 * cos theta
        vy0 = v0 * sin theta

-- air resistance quadratic in the speed
surfaceGravityAirResistance :: Double -> Double
                            -> OneParticleAccelerationFunction
surfaceGravityAirResistance m b (_t,St _r v)
    = netForce ^/ m
      where
        netForce      = gravity ^+^ airResistance
        gravity       = vec 0 (-m * g) 0
        airResistance = ((-b) * magnitude v) *^ v

trajectory :: [OneParticleSystemState] -> [(Double,Double)]
trajectory sts = [(x,y) | (_,St r _) <- sts, let (x,y,_) = cartesianCoordinates r]

traj :: Double -> [(Double,Double)]
traj b = trajectory $ inAir
         $ oneParticleRungeKuttaSolution
               (surfaceGravityAirResistance 2 b)
               0.01
               (initialProjState 30 (pi/6))

main :: IO ()
main = plotPathsStyle
       [Title "Trajectories of 2-kg object, initial speed 30 m/s, angle 30 degrees"
       ,XLabel "Range (m)"
       ,YLabel "Height (m)"
       ,PNG "learn-physics-Projectile.png"
       ] [(defaultStyle {lineSpec = CustomStyle [LineTitle "No air resistance"]}, traj 0)
         ,(defaultStyle {lineSpec = CustomStyle [LineTitle "Drag 0.01 kg/m"]}, traj 0.01)
         ,(defaultStyle {lineSpec = CustomStyle [LineTitle "Drag 0.02 kg/m"]}, traj 0.02)
         ] >> putStrLn "output sent to file learn-physics-Projectile.png"