In this post we’ll continue the previous one about vectors and take a look at calculating the path of a projectile and rendering that path to the screen.
Imagine a single particle in three dimensional space, we can characterise a state for it as its position and velocity at some instant in time.
type State = (Time, Displacement, Velocity)
From basic mechanics and Newtons laws we know that if no forces act on it then not a lot happens really! Time will increase and depending on your point of view, not much else will change. However, if some forces are acting then things become more interesting. Knowing the forces acting is equivalent to knowing the acceleration of the particle and we can suggest that acceleration is a function of State.
|
1 |
type Accnf = State -> Vector |
We know from simple dynamics the relationships between velocity, acceleration and displacement
Converting this to our Vector notation and taking small time intervals we can derive a function that takes an Accnf function, a small time interval and a State and returns a new State.
|
1 2 3 4 5 |
step :: Accnf -> Float -> State -> State step f dt st@(t, r, v) = (t', r', v') where t' = t + dt -- advance t by a small ammount r' = r ^+^ v ^* dt -- from equation 1 above v' = v ^+^ f st ^* dt -- from equation 2 above |
So for a given acceleration function, a (small) time step and a start State we can approximate a numerical solution with
|
1 2 |
solution :: Accnf -> Float -> State -> [State] solution a dt = iterate (step a dt) |
this gives a lazily infinite list of State and so at some point we use some form of ‘take‘ function to collect a list of State that can be further processed and plotted.
As a simple example of applying this consider a projectile fired from a gun with velocity v at an angle theta to the horizontal. A detailed analytic treatment can be found here.
For our purposes we will ignore air resistance and assume a mass of 1 unit so that the only force acting on the particle is the acceleration due to gravity – this gives
|
1 2 3 |
-- Only accn due to gravity acting down projectile :: Accnf -- i.e. State -> Vector projectile (_, _, _) = V (0, -9.8, 0) |
To initialise the projectile with velocity v at angle theta we need to consider the horizontal and vertical components of the velocity as shown here
|
1 2 3 4 |
-- Projectile fire at angle theta with velocity v has x, y componenst v cos and v sin projectileInit :: Float -> Float -> V.Vector -- 0.0174533 is conversion factor degrees to radians projectileInit v theta = V (v * cos (theta * 0.0174533), v * sin (theta * 0.0174533), 0) |
We can now create a solution giving a list of State
|
1 2 3 4 5 |
-- Projectile at velocity v, angle theta. Keep evaluating while the projectile is above ground projectileAtVandTheta :: Float -> Float -> [State] projectileAtVandTheta v theta = takeWhile heightPositive $ solution projectile 0.01 (0, p0, projInit) where projInit = projectileInit v theta p0 = V (0, 0 ,0) |
which computes a solution and just takes values while the projectile is above ground.
For display purposes I want to ‘fire’ several projectiles and vary the angle to the horizontal keeping the initial velocity the same. This function will do that for us
|
1 2 3 4 5 |
-- calculate projectile at v theta and keep reducing theta until < 0 severalProjectiles :: Float -> Float -> Float -> [[State]] severalProjectiles v theta dtheta | theta > 0 = projectileAtVandTheta v theta : severalProjectiles v (theta - dtheta) dtheta | otherwise = [] |
Rather than go through the code function by function I’ll put the ‘finished’ code below.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 |
{-# LANGUAGE OverloadedStrings #-} {-# OPTIONS -Wall -fno-warn-type-defaults #-} module Dynamics where import Text.Printf import Vectors as V import Graphics.Gloss import Views type Time = Float type Displacement = V.Vector type Velocity = V.Vector type State = (Time, Displacement, Velocity) type Accnf = State -> V.Vector -- get the x, y component from the Displacement vector displacement :: State -> Point displacement (_, V (x, y, _) , V (_, _, _)) = (x, y) visualPathFromStates :: [State] -> Path visualPathFromStates = map displacement -- get the x, y component from the Velocity vector velocity :: State -> Point velocity (_, V (_, _, _) , V (velx, vely, _)) = (velx, vely) visualVelFromStates :: [State] -> Path visualVelFromStates = map velocity step :: Accnf -> Float -> State -> State step f dt st@(t, r, v) = (t', r', v') where t' = t + dt r' = r ^+^ v ^* dt v' = v ^+^ f st ^* dt solution :: Accnf -> Float -> State -> [State] solution a dt = iterate (step a dt) compareYVal :: V.Vector -> Float -> (Float -> Float -> Bool) -> Bool compareYVal (V (_, y, _)) d p = p y d yNZero :: V.Vector -> Bool yNZero v = compareYVal v 0 (>=) heightPositive :: State -> Bool heightPositive (_, d, _) = yNZero d displayStates :: [State] -> String displayStates = foldr f "" where f (t, p, v) str = str ++ "T:" ++ (printf "%.4f" t :: String) ++ " P:" ++ show p ++ " V:" ++ show v ++ "\n" -- Projectile fire at angle theta with velocity v has x, y componenst v cos and v sin projectileInit :: Float -> Float -> V.Vector projectileInit v theta = V (v * cos (theta * 0.0174533), v * sin (theta * 0.0174533), 0) -- Only accn due to gravity acting down projectile :: Accnf projectile (_, _, _) = V (0, -9.8, 0) -- Projectile at velocity v, angle theta. Keep evaluating while the projectile is above ground projectileAtVandTheta :: Float -> Float -> [State] projectileAtVandTheta v theta = takeWhile heightPositive $ solution projectile 0.01 (0, p0, projInit) where projInit = projectileInit v theta p0 = V (0, 0 ,0) -- calculate projectile at v theta and keep reducing theta until < 0 severalProjectiles :: Float -> Float -> Float -> [[State]] severalProjectiles v theta dtheta | theta > 0 = projectileAtVandTheta v theta : severalProjectiles v (theta - dtheta) dtheta | otherwise = [] paths :: [[State]] -> [Path] paths = map visualPathFromStates -- make a picture from each path pathPlots :: [Path] -> [Picture] pathPlots = map (color blue . Line ) plotSeveralProjectiles :: Float -> Float -> Float -> IO () plotSeveralProjectiles v theta dtheta = drawPics . pathPlots . paths $ severalProjectiles v theta dtheta plot :: IO () plot = plotSeveralProjectiles 50 90 5 |
The ‘plot’ function calculates and shows the path of a particle fired at angles from 90 to 5 degrees varying by 5 degrees. and here are the results!
which agrees quite well with an analytical solution that would show an angle of 45 degrees having the maximum range and complementary angles giving the same range.
And just for the pretty picture doing plotSeveralProjectiles 50 180 5 gives
And finally… In the last post I showed a function
|
1 2 3 4 5 6 |
vecsAtOrigin :: Int -> V.Vector -> IO () vecsAtOrigin n = drawPics . take n . map lineVectorO . iterate (V.rotateXY (2*pi / fromIntegral n)) |
that would render n equally spaced vectors about the origin.
Because vecsAtOrigin is a specific case of a more general one of rendering vectors about a point and then translating them all to the origin it became an ‘itch’ that I really needed to ‘scratch’! Here’s the version after scratching the itch.
|
1 2 3 4 5 6 7 8 9 |
vecsAtOrigin :: Int -> V.Vector -> IO () vecsAtOrigin n = vecsAtPos n origin vecsAtPos :: Int -> V.Vector -> V.Vector -> IO () vecsAtPos n p = drawPics . take n . map (\x -> lineVector p (x ^+^ p)) . iterate (V.rotateXY (2*pi / fromIntegral n) ) |
and here is an example using vecsAtPos
I think that’s enough for now. Next time we’ll maybe look at other, slightly more complex, forces and see if the naive step function is capable of handling them. As usual all the code is available of GitHub.
Thanks for reading!