MECH2400 assignment code.

AssignmentE.hs

-- Main module so that this file can be compiled.
module Main where

-- Import module defined in StressAnalysis.hs (https://gist.github.com/3979818)
import StressAnalysis

-- Main function detemermines what the compiled executable does.
-- At the moment, I want to graph the stress at each element of a section.
main = graphElements

-- Generate data suitable for plotting the factor of safety against different
-- beam diameters. mapM_ just maps the 'displayPair' function over the
-- 'safeties' list in a way suitable for IO.
graphSafeties = mapM_ displayPair safeties
    where
        -- List of tuples (diameter, F.S.) created using a list comprehension
        safeties = [(d, f d) | d <- [17e-3, 17.01e-3 .. 18e-3]]
        -- Function that converts a diameter to a factor of safety, used above
        f = safe . (evaluateDiameter sections)

-- Generate the principle stress for all elements in the 'inside pinion'
-- Section. Output suitable for plotting with gnuplot.
graphElements = mapM_ display [(e, f e) | e <- elements] where
    elements = (getElements (defaults {
        dd = 0.0001,
        df = (2.22e-2)/2,
        t0 = pi / 50,
        dt = pi / 50
     }))
    f = (principalStress (Circle 2.22e-2) (sections !! 3))
    display r = putStrLn (
        (show $ distance (fst r)) ++ " " ++
        (show $ angle (fst r)) ++ " " ++
        (show (snd r))
     )

-- Find the diameter at which the factor of safety becomes greater than FS.
-- Recursively increases the given diameter until it meets the fs given in a
-- constant.
findDiameter :: [Section] -> Float -> Result
findDiameter ss d =
    if (safe result) > fs
        then result
        else findDiameter ss (d + dd')
    where
        result = evaluateDiameter ss d

-- Finds the maximum principal stress 
evaluateDiameter :: [Section] -> Float -> Result
evaluateDiameter ss d = result {
    -- Fill in additional fields in the Result
    diam = d,
    safe = sy / str result
 } where
    result = (maxPrincipalStress (Circle d) ss)

-- Convenience function so I can quickly get the result for a beam with a
-- particular diameter.
test = evaluateDiameter sections

-- Find the maximum principle stress given a Beam and several Sections.
maxPrincipalStress :: Beam -> [Section] -> Result
maxPrincipalStress b ss =
    -- Find the maximum over a constructed list
    maximum [
        -- Finds the maximum principle stress over a list of elements
        (maxPrincipalStress'
            -- Fundtion to convert an element to a shear stress
            (principalStress b s)
            -- Generated list of elements
            (getElements defaults { df = (diameter b) / 2 }))
         -- Append the current section to this result
         { sec = s }
        -- Try each section
        | s <- ss
     ]

-- Find the element that gives the maximum stress condition and return it with
-- the principal stress it sustains.
maxPrincipalStress' :: (Element -> Float) -> [Element] -> Result
maxPrincipalStress' f es
    | null es   = newResult
    | otherwise = newResult { str = fst max, el = snd max }
    where
        max = maximum [(f e, e) | e <- es]

-- Creates a list of elements with all variations of distance and angle given
-- in a Parameters object.
getElements :: Parameters -> [Element]
getElements p = [
    -- This is a list comprehension. Each element looks like:
    Element { distance = d, angle = a }
    -- Given all combinations of d and a governed by:
    | d <- [d0 p, d0 p + dd p .. df p],
      a <- [t0 p, t0 p + dt p .. tf p]
 ]

-- Parameters governing the way Elements are created on a Section.
data Parameters = Parameters {
    d0 :: Float, --Initial distance
    df :: Float, -- Final distance
    dd :: Float,
    t0 :: Float, -- Initial theta
    tf :: Float, -- Final theta
    dt :: Float
 } deriving (Show, Read)

defaults = Parameters {
    d0 = 0,
    df = rp,
    dd = 0.0001,
    t0 = pi / 50,
    tf = 2 * pi,
    dt = pi / 50
 }

-- A result that is constructed by several different function calls.
data Result = Result {
    str :: Float,
    diam :: Float,
    safe :: Float,
    el :: Element,
    sec :: Section
 } deriving (Eq, Ord)

instance Show Result where
    show r = "Result:" ++
        "\n    Principal stress: " ++ (show $ str r) ++
        "\n    Factor of safety: " ++ (show $ safe r) ++
        "\n    Beam diameter:    " ++ (show $ diam r) ++
        "\n    Element: " ++ (show $ el r) ++
        "\n    Section: " ++ (show $ sec r)

newResult = Result {
    str = 0/0,
    diam = 0/0,
    safe = 0/0,
    el = newElement,
    sec = newSection
 }

 -- Constants!
d0'  :: Float
d0' = 0.005
dd' :: Float
dd' = 0.0001

-- http://www.matweb.com/search/DataSheet.aspx?MatGUID=7573afc5e06c4a518c01efba4690182c
sy :: Float
sy = 180000000

fs :: Float
fs = 5

rp :: Float
rp = 0.054 - 1.25 * 0.002

-- The four critical points of interest along the shaft.
sections = [
    Section {
        name = "Outside gear",
        shearY = 43.14146621,
        shearZ = 32.62078375,
        momentX = 0,
        momentY = 0.9786235125,
        momentZ = 1.294243986
     },
    Section {
        name = "Inside gear",
        shearY = 43.14146621-107.6923077,
        shearZ = 32.62078375-26.96400044,
        momentX = 5.815384615,
        momentY = 0.9786235125,
        momentZ = 1.294243986
     },
    Section {
        name = "Inside pinion",
        shearY = 43.14146621-107.6923077,
        shearZ = 32.62078375-26.96400044,
        momentX = 5.815384615,
        momentY = 1.510361144,
        momentZ = -4.773535114
     },
    Section {
        name = "Outside pinion",
        shearY = 159.1178053,
        shearZ = -32.62078375,
        momentX = 0,
        momentY = 1.510361144,
        momentZ = -4.773535114
     }
 ]

-- Utility function to display a pair nicely.
displayPair :: (Show a, Show b) => (a, b) -> IO ()
displayPair p = putStrLn ((show (fst p)) ++ " " ++ (show (snd p)))