hfoil 0.1.2 → 0.2.0
raw patch · 15 files changed
+960/−783 lines, 15 filesdep +hfoildep +lineardep +not-glossdep −MissingHdep −gloss
Dependencies added: hfoil, linear, not-gloss, parsec
Dependencies removed: MissingH, gloss
Files
- Main.hs +0/−8
- Numeric/HFoil.hs +0/−26
- Numeric/HFoil/Drawing.hs +0/−168
- Numeric/HFoil/Flow.hs +0/−206
- Numeric/HFoil/Foil.hs +0/−157
- Numeric/HFoil/Naca4.hs +0/−68
- Numeric/HFoil/Repl.hs +0/−125
- apps/Main.hs +8/−0
- hfoil.cabal +28/−25
- src/HFoil.hs +27/−0
- src/HFoil/Drawing.hs +149/−0
- src/HFoil/Flow.hs +205/−0
- src/HFoil/Foil.hs +246/−0
- src/HFoil/Naca4.hs +71/−0
- src/HFoil/Repl.hs +226/−0
− Main.hs
@@ -1,8 +0,0 @@-{-# OPTIONS_GHC -Wall #-}--module Main where--import Numeric.HFoil.Repl--main :: IO ()-main = run
− Numeric/HFoil.hs
@@ -1,26 +0,0 @@-{- |- Module : Numeric.Hfoil- Description : Top level module-- This is the top level module which exports the API- -}--{-# OPTIONS_GHC -Wall #-}--module Numeric.HFoil( -- * Airfoils- module Numeric.HFoil.Foil- -- * Flow solution- , module Numeric.HFoil.Flow- -- * Naca 4 utilities- , module Numeric.HFoil.Naca4- -- * Drawing utilities (gloss backend)- , module Numeric.HFoil.Drawing- -- * Interactive command line application- , module Numeric.HFoil.Repl- ) where--import Numeric.HFoil.Naca4-import Numeric.HFoil.Foil-import Numeric.HFoil.Drawing-import Numeric.HFoil.Flow-import Numeric.HFoil.Repl
− Numeric/HFoil/Drawing.hs
@@ -1,168 +0,0 @@-{-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE FlexibleContexts #-}--module Numeric.HFoil.Drawing( drawLine- , drawLineV- , drawSolution- , drawFoil- , drawOnce- , drawNormals- , drawForces- , drawKuttas- ) where--import Graphics.Gloss hiding(Vector,dim)-import Numeric.LinearAlgebra hiding(Element, scale,i)-import Foreign.Storable(Storable)-import qualified Numeric.LinearAlgebra as LA-import Text.Printf--import Numeric.HFoil.Flow-import Numeric.HFoil.Foil--xSize, ySize :: Int-xSize = 800-ySize = 500--cpScale :: Double-cpScale = -0.25--border :: Float-border = 0.7--normalLengths :: Double-normalLengths = 0.01--drawText :: Color -> (Float, Float) -> Float -> String -> Picture-drawText col (x,y) size str = translate (0.5*x*(fromIntegral xSize)) (0.5*y*(fromIntegral ySize))- $ scale size size- $ color col- $ Text str--drawLine :: Real a => Color -> [(a,a)] -> Picture-drawLine col coords = scale (border*(fromIntegral xSize)) (border*(fromIntegral xSize))- $ translate (-0.5) 0- $ color col- $ line $ map (\(x,y) -> (realToFrac x, realToFrac y)) coords--drawCircle :: Real a => Color -> (a, a) -> Float -> Picture-drawCircle col (x,y) size = scale (border*(fromIntegral xSize)) (border*(fromIntegral xSize))- $ translate (-0.5 + realToFrac x) (realToFrac y)- $ color col- $ Circle size--drawLineV :: (Real a, Storable a) => Color -> (Vector a, Vector a) -> Picture-drawLineV col (vx, vy) = drawLine col $ zip (toList vx) (toList vy)--drawFoil :: (Real a, Storable a, Show a) => Foil a -> Picture-drawFoil (Foil elements _) = pictures $ map drawElement elements--drawElement :: (Real a, Storable a) => Element a -> Picture-drawElement element = drawLineV white (fNodes element)--drawNormals :: Foil Double -> Picture-drawNormals (Foil elements _) = pictures $ map (\(xy0, xy1) -> drawLine green [xy0, xy1]) (zip xy0s xy1s)- where- xy0s = zip (toList xm) (toList ym)- xy1s = zip (toList (xm + (LA.scale normalLengths xUnitNormal))) (toList (ym + (LA.scale normalLengths yUnitNormal)))- - (xUnitNormal, yUnitNormal) = (\(x,y) -> (join x, join y)) $ unzip $ map fUnitNormals elements- (xm, ym) = (\(x,y) -> (join x, join y)) $ unzip $ map fMidpoints elements--colorFun :: (Fractional a, Real a) => a -> a -> a -> Color-colorFun min' max' x' = makeColor (1-x) (1-x) x 1- where- x = realToFrac $ (x' - min')/(max'-min')--drawKuttas :: (Real a, Storable a) => FlowSol a -> Picture-drawKuttas flow = pictures $ concatMap (\(k0,k1) -> [circ k0, circ k1]) kis- where- kis = solKuttaIndices flow- (xs',ys') = unzip $ map fMidpoints $ (\(Foil els _) -> els) (solFoil flow)- xs = join xs'- ys = join ys'- circ k = drawCircle yellow (xs @> k, ys @> k) 0.006--drawForces :: FlowSol Double -> Picture-drawForces flow = pictures $ map (\(xy0, xy1, cp) -> drawLine (colorFun minCp maxCp cp) [xy0, xy1])- $ zip3 xy0s xy1s (toList (solCps flow))- where- xy0s = zip (toList xm) (toList ym)- xy1s = zip (toList (xm + xPressures)) (toList (ym + yPressures))- (xPressures, yPressures) = (\(x,y) -> (LA.scale c x/lengths, LA.scale c y/lengths)) (solForces flow)- lengths = join $ map fLengths $ (\(Foil x _) -> x) $ solFoil flow- (xm, ym) = (\(x,y) -> (join x, join y)) $ unzip $ map fMidpoints $ (\(Foil x _) -> x) $ solFoil flow- - c = 0.1- - maxCp = maxElement (solCps flow)- minCp = minElement (solCps flow)--drawColoredFoil :: [Color] -> Foil Double -> Picture-drawColoredFoil colors foil@(Foil elements _) = pictures $ zipWith drawColoredElement colors' elements- where- colors' = groupSomethingByFoil foil colors--drawColoredElement :: [Color] -> Element Double -> Picture-drawColoredElement colors element = pictures $ map (\(xy0, xy1, col) -> drawLine col [xy0, xy1]) (zip3 xy0s xy1s colors)- where- xys = (\(x,y) -> zip (toList x) (toList y)) $ fNodes element- xy0s = tail xys- xy1s = init xys--groupSomethingByFoil :: Storable a => Foil a -> [b] -> [[b]]-groupSomethingByFoil (Foil elements _) somethings = f somethings (map (dim . fAngles) elements)- where- f xs (n:ns) = (take n xs):(f (drop n xs) ns)- f [] []= []- f _ _ = error "uh oh (groupSomethingByFoil)"--drawSolution :: FlowSol Double -> Picture-drawSolution flow = pictures $ onscreenText ++- [ drawColoredFoil colors foil- , drawCircle white (fst $ solCenterPressure flow, snd $ solCenterPressure flow) 0.006- , drawCircle white (fst $ solCenterPressure flow, 0) 0.006- , drawCircle green (0.25,0::Double) 0.006- ] ++ zipWith (\x y -> drawLineV red (x, y)) xs- (takesV (map dim xs) (LA.scale cpScale cps)) -- cp graph- where- foil@(Foil elements name) = solFoil flow- cps = solCps flow- - xs = map (fst . fMidpoints) elements- - onscreenText = zipWith (\m y -> drawText white (0.45,y) 0.15 m) msgs- $ take (length msgs) [0.8,0.65..]- msgs = [ name- , printf ("alpha: %.6f deg") ((solAlpha flow)*180/pi)- , printf ("Cl: %.6f") (solCl flow)- , printf ("Cd: %.6f") (solCd flow)- , printf ("Cm: %.6f (c/4, 0)") (solCm flow)- ]-- colors = map (colorFun (minElement cps) (maxElement cps)) (toList cps)----drawOnce :: [Picture] -> IO ()-drawOnce pics = do- display - (InWindow- "hfoil" -- window title- (xSize, ySize) -- window size- (10, 650)) -- window position- black -- background color- (pictures pics) -- picture to display---- let line = plot_lines_values ^= [[ (xc, yt (naca4 "0012") xc)--- | xc <- [0,0.01..0.99::Double]]]--- $ plot_lines_title ^= "naca 0012"--- $ defaultPlotLines--- --- chart = layout1_title ^= "naca yo"--- $ layout1_plots ^= [Left (toPlot line)]--- $ defaultLayout1--- --- renderableToWindow (toRenderable chart) 640 480--- _ <- renderableToPNGFile (toRenderable chart) 640 480 "mDiv_vs_tc.png"--- return ()
− Numeric/HFoil/Flow.hs
@@ -1,206 +0,0 @@-{-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE FlexibleContexts #-}--module Numeric.HFoil.Flow( FlowSol(..)- , solveFlow- ) where--import Data.Packed.ST(newMatrix, writeMatrix, freezeMatrix)-import Control.Monad(forM_)-import Control.Monad.ST(runST)-import Numeric.LinearAlgebra hiding(i)-import Foreign.Storable--import Numeric.HFoil.Foil--data FlowSol a = FlowSol { solFoil :: Foil a- , solVs :: Vector a- , solCps :: Vector a- , solVorticities :: [a]- , solAlpha :: a- , solForces :: (Vector a, Vector a)- , solForce :: (a,a)- , solCl :: a- , solCd :: a- , solCm :: a -- moment about (0,0.25)- , solCenterPressure :: (a,a)- , solKuttaIndices :: [(Int,Int)]- }--solveFlow :: (Num (Vector a), RealFloat a, Field a) => Foil a -> a -> FlowSol a-solveFlow foil@(Foil elements _) alpha =- FlowSol { solFoil = foil- , solVs = vs- , solCps = cps- , solVorticities = vorticities- , solAlpha = alpha- , solForces = (xForces, yForces)- , solForce = (xf,yf)- , solCl = cl- , solCd = cd- , solCm = cm- , solCenterPressure = (xCp, yCp)- , solKuttaIndices = kuttaIndices- }- where- -- surface velocities- (vs,vorticities) = ( (mapVector (\q -> cos(q - alpha)) angles) + (mV <> qsGamma)- , map (qsGamma @>) [n-(length elements),n-1]- )- where- (mA, mV, b) = getAVb geometries angles kuttaIndices alpha- qsGamma = flatten $ linearSolve mA b- n = dim angles- -- pressure coefficients- cps = 1 - vs*vs- - -- forces and force coefficients- xForces = -cps*(join $ map (fst . fNormals) elements)- yForces = -cps*(join $ map (snd . fNormals) elements)- xf = sumElements xForces- yf = sumElements yForces- - cd = xf*(cos alpha) + yf*(sin alpha)- cl = -xf*(sin alpha) + yf*(cos alpha)- - -- midpoints- (xs,ys) = (\(x,y) -> (join x, join y)) $ unzip $ map fMidpoints elements- - -- centers of pressure- (xCp, yCp) = ( (sumElements (xs*yForces)) / (sumElements yForces)- , (sumElements (ys*xForces)) / (sumElements xForces)- )-- -- moment about cp and quarter chord- cm = sumElements $ (mapVector (\z -> z - 0.25) xs)*yForces - ys*xForces-- kuttaIndices = ki 0 (map dim angles')- where- ki _ [] = []- ki n0 (n:ns) = (n0,n0+n-1):(ki (n0+n) ns)-- angles = join angles'- angles' = map fAngles elements-- -- (mSL,mCL,mSB,mCB)- geometries = setupGeometries angles (xms, yms) (xns, yns) (xnps, ynps)- where- joinTuples (x,y) = (join x, join y)- (xms,yms) = joinTuples $ unzip $ map fMidpoints elements- (xns,yns) = joinTuples $ unzip $ map fInits elements- (xnps,ynps) = joinTuples $ unzip $ map fTails elements----- make influence matrix A and also the matrix V where V*[sources; vortex] == tangential speeds---getAVb :: (Num (Vector a), RealFloat a, Container Vector a) =>--- Foil a -> a -> (Matrix a, Matrix a)-getAVb :: (Floating a, Num (Vector a), Container Vector a, Storable a) =>- (Matrix a, Matrix a, Matrix a, Matrix a)- -> Vector a- -> [(Int, Int)]- -> a- -> (Matrix a, Matrix a, Matrix a)-getAVb (mSL,mCL,mSB,mCB) angles kuttaIndices alpha =- ( scale (1/(2*pi)) $ fromBlocks [[ mAij, fromColumns vsAin]- ,[fromRows vsAnj, fromLists ann]]- , scale (1/(2*pi)) $ fromBlocks [[ mSB - mCL- , fromColumns vsTin]]- , asColumn $ join $ (mapVector (\q -> sin $ q - alpha) angles):- (map (\(n0,nf) -> fromList [-cos((angles @> n0) - alpha) - cos((angles @> nf) - alpha)]) kuttaIndices)- ) -- (for middle entry mV, mSL + mCB == mAij)- where- n = rows mSL- - -- sources influence matrix Aij- mAij = mSL + mCB- - -- vortices influence vectors Ains- vsAin = map (\(n0,nf) -> fromList $ map sumElements $ toRows $ subMatrix (0,n0) (n, nf-n0+1) m) kuttaIndices- where- m = (mCL - mSB)- - -- vortices tangential velocity outputs- vsTin = map (\(n0,nf) -> fromList $ map sumElements $ toRows $ subMatrix (0,n0) (n, nf-n0+1) m) kuttaIndices- where- m = (mSL + mCB)-- -- sources kutta condition influence vector Anj- vsAnj = map (\(n0,nf) -> (getRow n0 mSB) - (getRow n0 mCL) + (getRow nf mSB) - (getRow nf mCL)) kuttaIndices-- -- vortices kutta condition influence scalars ann- ann = (flip map) kuttaIndices $ \(ni0,nif) -> (flip map) kuttaIndices $ \(nj0,njf) ->- sumElements $ (getSubRow ni0 (nj0,njf) mSL) + (getSubRow ni0 (nj0,njf) mCB) + (getSubRow nif (nj0,njf) mSL) + (getSubRow nif (nj0,njf) mCB)-- getRow i mat = flatten $ subMatrix (i,0) (1,n) mat- getSubRow i (j0,jf) mat = flatten $ subMatrix (i,j0) (1,jf-j0+1) mat---{- -calcuate 4 matrices which will be useful--mSL(i,j) = sin(qi - ij) * ln(rij+1/rij)-mCL(i,j) = cos(qi - ij) * ln(rij+1/rij)-mSB(i,j) = sin(qi - ij) * beta(i,j)-mCB(i,j) = cos(qi - ij) * beta(i,j)--where for i==j: ln(r/r) = 0- beta = pi-are explicitly set--}-setupGeometries :: (RealFloat t, Storable t) =>- Vector t- -> (Vector t, Vector t)- -> (Vector t, Vector t)- -> (Vector t, Vector t)- -> (Matrix t, Matrix t, Matrix t, Matrix t)-setupGeometries angles (xms, yms) (xns, yns) (xnps, ynps) = runST $ do- let n = dim angles- mSL' <- newMatrix 0 n n- mCL' <- newMatrix 0 n n- mSB' <- newMatrix 0 n n- mCB' <- newMatrix 0 n n- - forM_ [0..n-1] $ \i -> forM_ [0..n-1] $ \j -> do- let qi = angles @> i- qj = angles @> j- - xmi = xms @> i- ymi = yms @> i- xnj = xns @> j- ynj = yns @> j- xnjp = xnps @> j- ynjp = ynps @> j- - s = sin (qi - qj)- c = cos (qi - qj)- - l - | i == j = 0- | otherwise = log(r1/r0)- where- r1 = distance (xmi, ymi) (xnjp, ynjp)- r0 = distance (xmi, ymi) (xnj, ynj )- distance (x1,y1) (x0,y0) = sqrt $ dx*dx + dy*dy- where- dx = x1 - x0- dy = y1 - y0-- b - | i == j = pi- | otherwise = atan2 (dyjp*dxj - dxjp*dyj)- (dxjp*dxj + dyjp*dyj)- where- dyj = ymi - ynj- dxj = xmi - xnj- dyjp = ymi - ynjp- dxjp = xmi - xnjp- - _ <- writeMatrix mSL' i j (s*l)- _ <- writeMatrix mCL' i j (c*l)- _ <- writeMatrix mSB' i j (s*b)- writeMatrix mCB' i j (c*b)- sl' <- freezeMatrix mSL'- cl' <- freezeMatrix mCL'- sb' <- freezeMatrix mSB'- cb' <- freezeMatrix mCB'- return (sl',cl',sb',cb')
− Numeric/HFoil/Foil.hs
@@ -1,157 +0,0 @@-{-# OPTIONS_GHC -Wall #-}-{-# LANGUAGE FlexibleContexts #-}--module Numeric.HFoil.Foil( Foil(..)- , Element(..)- , panelizeNaca4- , loadFoil- , getUIUCFoil- ) where--import System.Directory(doesFileExist)-import Numeric.LinearAlgebra hiding (Element)-import Data.Tuple.Utils(fst3)-import Network.HTTP(simpleHTTP, getRequest, getResponseBody)-import Foreign.Storable(Storable)--import qualified Numeric.HFoil.Naca4 as Naca4--data Foil a = Foil [Element a] String--instance (Storable a) => Show (Foil a) where- show (Foil [el] name) = "{"++name++": " ++ show (1 + dim (fLengths el)) ++ " nodes}"- show (Foil els name) = "{"++name++": " ++ show (length els) ++ " elements, " ++- show nodesPerEl ++ " nodes == "++show (sum nodesPerEl)++" total nodes}"- where- nodesPerEl = map (\x -> 1 + dim (fLengths x)) els--data Element a = Element { fNodes :: (Vector a, Vector a)- , fLengths :: Vector a- , fAngles :: Vector a- , fMidpoints :: (Vector a, Vector a)- , fTangents :: (Vector a, Vector a)- , fNormals :: (Vector a, Vector a)- , fUnitNormals :: (Vector a, Vector a)- , fInits :: (Vector a, Vector a)- , fTails :: (Vector a, Vector a)- }--toElement :: (Num (Vector a), RealFloat a, Container Vector a) =>- [(a, a)] -> Element a-toElement xynodes = Element { fNodes = (xNodes, yNodes)- , fLengths = lengths- , fAngles = zipVectorWith atan2 yTangents xTangents- , fMidpoints = (xMids, yMids)- , fTangents = (xTangents, yTangents)- , fNormals = (xNormals, yNormals)- , fUnitNormals = (xUnitNormals, yUnitNormals)- , fInits = (xInits, yInits)- , fTails = (xTails, yTails)- }- where- n = (dim xNodes) - 1- (xNodes, yNodes) = (\(xs,ys) -> (fromList xs, fromList ys)) $ unzip xynodes- (xInits, yInits) = (subVector 0 n xNodes, subVector 0 n yNodes)- (xTails, yTails) = (subVector 1 n xNodes, subVector 1 n yNodes)- (xTangents, yTangents) = (xTails - xInits, yTails - yInits)- (xMids, yMids) = (0.5*(xInits + xTails), 0.5*(yInits + yTails))- (xNormals, yNormals) = (-yTangents, xTangents)- lengths = mapVector sqrt $ xTangents*xTangents + yTangents*yTangents- (xUnitNormals, yUnitNormals) = (xNormals/lengths, yNormals/lengths)----- why isn't this standard???-poorMansStrip :: String -> String-poorMansStrip str = reverse $ dropWhile (== ' ') $ reverse $ dropWhile (== ' ') str--getUIUCFoil :: (Num (Vector a), Read a, RealFloat a, Container Vector a) =>- String -> IO (Either String (Foil a))-getUIUCFoil name' = do- let name = poorMansStrip name'- let file = "http://www.ae.illinois.edu/m-selig/ads/coord/" ++ name ++ ".dat"- dl <- simpleHTTP (getRequest file) >>= getResponseBody- return (parseRawFoil dl name)--parseRawFoil :: (Num (Vector a), Read a, RealFloat a, Container Vector a) =>- String -> String -> Either String (Foil a)-parseRawFoil raw name- -- bad data- | any (\x -> 2 /= length x) (concat groupsOfLines) = Left (raw ++ "\nError parsing the above data")- -- single element- | length elements > 0 = Right (Foil (map toElement elements) name)- | otherwise = Left (raw ++ "\nError parsing the above data")- where- elements = map (map (\(x:y:[]) -> (read x, read y))) groupsOfLines- - -- group lines by splitting at empty lines- groupsOfLines :: [[[String]]]- groupsOfLines = f (lines raw)- where- f [] = []- f ([]:xs) = f xs- f xs = (map words fst'):(f snd')- where- (fst', snd') = break (\x -> 0 == length x) xs- -loadFoil :: FilePath -> IO (Either String (Foil Double))-loadFoil filename' = do- let filename = poorMansStrip filename'- exists <- doesFileExist filename- if exists- then do rawData <- readFile filename- return (parseRawFoil rawData filename) -- use filename as name- else return (Left ("file \"" ++ filename ++ "\" couldn't be found"))--panelizeNaca4 :: (Enum a, Floating (Vector a), RealFloat a, Field a) =>- Naca4.Naca4 a -> Int -> Foil a-panelizeNaca4 foil nPanels = Foil [toElement $ [(1,0)]++reverse lower++[(0,0)]++upper++[(1,0)]]- (Naca4.naca4_name foil)- where- (upper, lower) = unzip $ map (Naca4.coords foil) xcs- xcs = toList $ fst3 $ bunchPanels (Naca4.yt foil) (Naca4.dyt foil) xcs0 0 0- xcs0 = fromList $ init $ tail $ toList $ linspace nXcs (0,1)- nXcs = (nPanels + (nPanels `mod` 2)) `div` 2 + 1- -bunchPanels :: (Enum a, Floating (Vector a), Floating a, Ord a, Field a) =>- (a -> a) -> (a -> a) -> Vector a -> Int -> Int -> (Vector a, Int, Int)-bunchPanels yt dyt xcs nIter nBadSteps - | nIter > 300 = error "panel buncher exceeded 300 iterations"- | nBadSteps > 1000 = error "panel buncher exceeded 1000 bad steps"- | sum (toList (abs deltaXcs)) < 1e-12 = (xcs, nIter, nBadSteps)- | otherwise = bunchPanels yt dyt goodXcs (nIter+1) (nBadSteps + length badOnes)- where- (badOnes, goodXcs:_) = break (\xs -> all (>0) (toList xs)) nextXcs- nextXcs = map (\alpha -> xcs + (scale alpha deltaXcs)) (map (2.0**) [0,-1..])- deltaXcs = xcsStep yt dyt xcs--xcsStep :: (Enum a, Floating (Vector a), Field a) => (a -> a) -> (a -> a) -> Vector a -> Vector a-xcsStep yt dyt xcs = flatten $ -(linearSolveLS mat2 (asColumn rs))- where- n = dim xcs- - xs = join [fromList [0], xcs, fromList[1]]- ys = join [fromList [0], mapVector yt xcs, fromList[0]]- dxs = (subVector 1 (n+1) xs) - (subVector 0 (n+1) xs)- dys = (subVector 1 (n+1) ys) - (subVector 0 (n+1) ys)- deltas = sqrt (dxs*dxs + dys*dys)- dydxs = mapVector dyt xcs- - mat = r1 - r0- where- r0 = fromBlocks [[(1><n)[0,0..]],[diagRect 0 diag0 n n]]- r1 = fromBlocks [ [diagRect 0 diag1 n n], [(1><n)[0,0..]]]- - (diag0, diag1) = (subVector 1 n d0, subVector 0 n d1)- where- d0 = (dxs + dys*dy0dxs)/deltas- d1 = (dxs + dys*dy1dxs)/deltas- dy0dxs = join [fromList [0], dydxs]- dy1dxs = join [dydxs, fromList [0]]- - zeros = (n><1)[0,0..]- eye = ident n- frontBunchingParam = 2.0- diff = (fromBlocks [[zeros, eye]]) - (fromBlocks [[eye*(1+frontBunchingParam/(fromIntegral n)), zeros]])-- mat2 = diff <> mat- rs = diff <> deltas
− Numeric/HFoil/Naca4.hs
@@ -1,68 +0,0 @@-{-# OPTIONS_GHC -Wall #-}--module Numeric.HFoil.Naca4( Naca4(..)- , coords- , yt- , dyt- , naca4- ) where--naca4 :: (Read a, Fractional a) => String -> Naca4 a-naca4 name@(m_:p_:t0:t1:[]) = Naca4 m p t ("NACA "++name)- where- m = 0.01 * read [m_]- p = 0.1 * read [p_]- t = 0.01 * read (t0:[t1])-naca4 _ = error "not a 4 digit airfoil"---- m: max camber in hundredths of chord--- p: position of max camber in tenths of chord--- t: max thickness in hundredths of chord-data Naca4 a = Naca4 { naca4_m :: a- , naca4_p :: a- , naca4_t :: a- , naca4_name :: String- } deriving Show---- xc: x/chord-yc :: (Ord a, Fractional a) => Naca4 a -> a -> a-yc (Naca4 {naca4_p = p, naca4_m = m}) xc- | xc < 0 = error "xc < 0"- | xc <= p = m/(p*p)*(2*p*xc - xc*xc)- | xc <= 1 = m/((1-p)*(1-p))*((1-2*p) + 2*p*xc - xc*xc)- | otherwise = error "xc > 1"--dyc :: (Ord a, Fractional a) => Naca4 a -> a -> a-dyc (Naca4 {naca4_p = p, naca4_m = m}) xc- | xc < 0 = error "xc < 0"- | xc <= p = m/(p*p)*(2*p - 2*xc)- | xc <= 1 = m/((1-p)*(1-p))*(2*p - 2*xc)- | otherwise = error "xc > 1"--yt :: (Ord a, Floating a) => Naca4 a -> a -> a-yt (Naca4 {naca4_t = t}) xc- | xc < 0 = error "xc < 0"- | xc > 1 = error $ "xc > 1"- | otherwise = 5*t*(0.2969*sqrt(xc) - 0.1260*(xc) - 0.3537*(xc)**2 + 0.2843*(xc)**3 - 0.1015*(xc)**4)--dyt :: (Ord a, Floating a) => Naca4 a -> a -> a-dyt (Naca4 {naca4_t = t}) xc- | xc < 0 = error "xc < 0"- | xc > 1 = error "xc > 1"- | otherwise = 5*t*(0.5*0.2969/sqrt(xc) - 0.1260 - 2*0.3537*(xc) + 3*0.2843*(xc)**2 - 4*0.1015*(xc)**3)--coords :: (Ord a, Floating a) => Naca4 a -> a -> ((a,a),(a,a))-coords foil xc - | naca4_m foil == 0 = ((xc,yt_), (xc,-yt_))- | otherwise = ((xu,yu ), (xl, yl ))- where- yt_ = yt foil xc- yc_ = yc foil xc- - yu = yc_ + yt_ * (cos theta)- yl = yc_ - yt_ * (cos theta)-- xu = xc - yt_ * (sin theta)- xl = xc + yt_ * (sin theta)-- theta = atan $ dyc foil xc
− Numeric/HFoil/Repl.hs
@@ -1,125 +0,0 @@-{-# OPTIONS_GHC -Wall #-}--module Numeric.HFoil.Repl( run- ) where--import System.Console.Haskeline hiding(display)-import Graphics.Gloss.Interface.IO.Animate hiding(scale, Vector)-import Control.Monad.IO.Class-import Control.Concurrent(forkIO)-import Control.Concurrent.MVar(newMVar, readMVar, swapMVar)--import Numeric.HFoil.Foil-import Numeric.HFoil.Naca4-import Numeric.HFoil.Drawing-import Numeric.HFoil.Flow------ configuration-nPanels :: Int-nPanels = 200--xSize, ySize :: Int-xSize = 800-ySize = 500--data Config = Config { confForces :: Bool- , confKuttas :: Bool- , confNormals :: Bool- }--defaultConfig :: Config-defaultConfig = Config { confForces = False- , confKuttas = False- , confNormals = False- }--run :: IO ()-run = do- let naca0 = "2412"- alfaDeg0 = 4- flow0 = solveFlow (panelizeNaca4 (naca4 naca0) nPanels) (pi/180*alfaDeg0)- mpics <- newMVar $ [drawSolution flow0]- - putStrLn "Welcome to hfoil\n"- - _ <- forkIO $ runInputT defaultSettings- $ topLoop (\pics -> swapMVar mpics pics >>= (\_ -> return ())) defaultConfig-- animateIO- (InWindow- "hfoil" -- window title- (xSize, ySize) -- window size- (10, 650)) -- window position- black -- background color- (\_ -> readMVar mpics >>= return . pictures) -- draw function--foilLoop :: ([Picture] -> IO ()) -> Config -> Foil Double -> InputT IO ()-foilLoop draw conf foil@(Foil _ name) = do- minput <- getInputLine $ "\ESC[1;32m\STXhfoil."++name++">> \ESC[0m\STX"- case minput of- Nothing -> return ()- Just "quit" -> do outputStrLn "gloss won't let you quit :(\ntry ctrl-c or hit ESC in drawing window"- foilLoop draw conf foil- Just ('a':'l':'f':'a':' ':[]) -> do outputStrLn $ "unrecognized command"- foilLoop draw conf foil- Just ('a':'l':'f':'a':' ':alphaDeg) -> do let flow = solveFlow foil (pi/180*(read alphaDeg))- forces = case (confForces conf) of- True -> [drawForces flow]- False -> []- kuttas = case (confKuttas conf) of- True -> [drawKuttas flow]- False -> []- normals = case (confNormals conf) of- True -> [drawNormals (solFoil flow)]- False -> []- liftIO $ draw $ forces++kuttas++normals++[drawSolution flow]- foilLoop draw conf foil- Just ('f':'o':'r':'c':'e':'s':[]) -> do- let newConf = conf {confForces = not (confForces conf)}- outputStrLn $ "force drawing set to "++ show (not (confForces conf))- foilLoop draw newConf foil- Just ('k':'u':'t':'t':'a':'s':[]) -> do - let newConf = conf {confKuttas = not (confKuttas conf)}- outputStrLn $ "kutta drawing set to "++ show (not (confKuttas conf))- foilLoop draw newConf foil- Just ('n':'o':'r':'m':'a':'l':'s':[]) -> do- let newConf = conf {confNormals = not (confNormals conf)}- outputStrLn $ "normals drawing set to "++ show (not (confNormals conf))- foilLoop draw newConf foil- Just "" -> return ()- Just input -> do outputStrLn $ "unrecognized command \"" ++ input ++ "\""- foilLoop draw conf foil--topLoop :: ([Picture] -> IO ()) -> Config -> InputT IO ()-topLoop draw conf = do- minput <- getInputLine "\ESC[1;32m\STXhfoil>> \ESC[0m\STX"- case minput of- Nothing -> return ()- Just "quit" -> do outputStrLn "gloss won't let you quit :(\ntry ctrl-c or hit ESC in drawing window"- topLoop draw conf- Just ('n':'a':'c':'a':' ':spec) -> do parseNaca draw conf spec- topLoop draw conf- Just ('l':'o':'a':'d':' ':name) -> do- foil <- liftIO (loadFoil name)- case foil of Left errMsg -> outputStrLn errMsg- Right foil' -> do liftIO $ draw [drawFoil foil', drawNormals foil']- foilLoop draw conf foil'- topLoop draw conf- Just ('u':'i':'u':'c':' ':name) -> do- efoil <- liftIO (getUIUCFoil name)- case efoil of Left errMsg -> outputStrLn errMsg- Right foil -> do liftIO $ draw [drawFoil foil, drawNormals foil]- foilLoop draw conf foil- topLoop draw conf-- Just "" -> topLoop draw conf- Just input -> do outputStrLn $ "unrecognized command \"" ++ input ++ "\""- topLoop draw conf--parseNaca :: ([Picture] -> IO ()) -> Config -> String -> InputT IO ()-parseNaca draw conf str - | length str == 4 = do let foil = panelizeNaca4 (naca4 str :: Naca4 Double) nPanels- liftIO $ draw [drawFoil foil, drawNormals foil]- foilLoop draw conf foil- | otherwise = do outputStrLn $ "Not 4 digits"- return ()
+ apps/Main.hs view
@@ -0,0 +1,8 @@+{-# OPTIONS_GHC -Wall #-}++module Main where++import HFoil.Repl++main :: IO ()+main = run
hfoil.cabal view
@@ -1,31 +1,31 @@ Name: hfoil-Version: 0.1.2+Version: 0.2.0 Stability: Experimental Synopsis: Hess-Smith panel code for inviscid 2-d airfoil analysis License: BSD3 License-file: LICENSE Author: Greg Horn Maintainer: gregmainland@gmail.com-Copyright: Greg Horn, 2012+Copyright: (c) Greg Horn 2012, 2015 Category: Numerical, Science Build-type: Simple-Cabal-version: >=1.6+Cabal-version: >= 1.8 Description: { Library and command line REPL with plotting to do simple inviscid hess-smith panel code. . Features include: .-* Single and multi-element airfoils+* Cheap and shameless xfoil ripoff for relp/plotting interface . * Naca 4-series support with Gauss-Newton paneling . * Broken UIUC database integration (type \"uiuc [foilname]\") .-* Shameless xfoil ripoff for relp/plotting interface+* Haskeline interface with tab-completion (oooh) .-* Only works with development version of Gloss that's not on Hackage yet+* Single and multi-element airfoils .-* Haskeline interface with tab-completion+* Inviscid, incompressible, 2-dimensional flow only . . To get started, do cabal install or whatever, then run the \"hfoil\" binary.@@ -38,22 +38,24 @@ Default: True Library- Exposed-modules: Numeric.HFoil- Numeric.HFoil.Naca4- Numeric.HFoil.Foil- Numeric.HFoil.Drawing- Numeric.HFoil.Flow- Numeric.HFoil.Repl+ hs-source-dirs: src+ Exposed-modules: HFoil+ HFoil.Naca4+ HFoil.Foil+ HFoil.Drawing+ HFoil.Flow+ HFoil.Repl -- Other-modules: - Build-depends: gloss >= 1.7 && <= 1.8,- base >= 4 && < 5,- MissingH >= 1.1,- hmatrix >= 0.12,- haskeline >= 0.6,- transformers >= 0.2,- directory >= 1.1,- HTTP >= 4000.2+ Build-depends: base >= 4 && < 5+ , not-gloss >= 0.7.2.0+ , linear+ , hmatrix >= 0.12+ , haskeline >= 0.6+ , transformers >= 0.2+ , directory >= 1.1+ , HTTP >= 4000.2+ , parsec >= 3 Ghc-Options: -Wall Ghc-Prof-Options: -prof -auto-all@@ -65,13 +67,14 @@ else Buildable: False + Build-depends: base >= 4 && < 5+ , hfoil >= 0.1.2++ hs-source-dirs: apps Main-Is: Main.hs- -- (gloss needs -O2 -threaded)+ -- (not-gloss needs -O2 -threaded) -- does it really? Ghc-Options: -O2 -threaded -Wall GHC-Prof-Options: -auto-all -caf-all -rtsopts-- if os(OSX)- extra-lib-dirs: /usr/lib source-repository head type: git
+ src/HFoil.hs view
@@ -0,0 +1,27 @@+{- |+ Module : Hfoil+ Description : Top level module++ This is the top level module which exports the API+ -}++{-# OPTIONS_GHC -Wall #-}++module HFoil+ ( -- * Airfoils+ module HFoil.Foil+ -- * Flow solution+ , module HFoil.Flow+ -- * Naca 4 utilities+ , module HFoil.Naca4+ -- * Drawing utilities (gloss backend)+ , module HFoil.Drawing+ -- * Interactive command line application+ , module HFoil.Repl+ ) where++import HFoil.Naca4+import HFoil.Foil+import HFoil.Drawing+import HFoil.Flow+import HFoil.Repl
+ src/HFoil/Drawing.hs view
@@ -0,0 +1,149 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE FlexibleContexts #-}++module HFoil.Drawing+ ( drawSolution+ , drawFoil+ , drawNormals+ , drawForces+ , drawKuttas+ , drawOnce+ ) where++import Numeric.LinearAlgebra hiding( Element, scale, i )+import Foreign.Storable ( Storable )+import qualified Numeric.LinearAlgebra as LA+import Text.Printf ( printf )+import Linear ( V3(..) )++import Vis++import HFoil.Flow+import HFoil.Foil++cpScale :: Fractional a => a+cpScale = -0.25++normalLengths :: Fractional a => a+normalLengths = 0.01++drawLine :: Num a => Color -> [(a,a)] -> VisObject a+drawLine col coords = Line (map (\(x,y) -> V3 x y 0) coords) col++drawCircle :: Num a => Color -> (a, a) -> a -> VisObject a+drawCircle col (x,y) size = Trans (V3 x y 0) $ Sphere size Solid col+--drawCircle col (x,y) size = scale (border*(fromIntegral xSize)) (border*(fromIntegral xSize))+-- $ translate (-0.5 + realToFrac x) (realToFrac y)+-- $ color col+-- $ Circle size++drawLineV :: (Num a, Storable a) => Color -> (Vector a, Vector a) -> VisObject a+drawLineV col (vx, vy) = Line (zipWith (\x y -> V3 x y 0) (toList vx) (toList vy)) col++drawFoil :: (Num a, Storable a) => Foil a -> VisObject a+drawFoil (Foil elements _) = VisObjects $ map drawElement elements++drawElement :: (Num a, Storable a) => Element a -> VisObject a+drawElement element = drawLineV white (fNodes element)++drawNormals :: Foil Double -> VisObject Double+drawNormals (Foil elements _) = VisObjects $ map (\(xy0, xy1) -> drawLine green [xy0, xy1]) (zip xy0s xy1s)+ where+ xy0s = zip (toList xm) (toList ym)+ xy1s = zip (toList (xm + (LA.scale normalLengths xUnitNormal))) (toList (ym + (LA.scale normalLengths yUnitNormal)))++ (xUnitNormal, yUnitNormal) = (\(x,y) -> (vjoin x, vjoin y)) $ unzip $ map fUnitNormals elements+ (xm, ym) = (\(x,y) -> (vjoin x, vjoin y)) $ unzip $ map fMidpoints elements++colorFun :: (Fractional a, Real a) => a -> a -> a -> Color+colorFun min' max' x' = makeColor (1-x) (1-x) x 1+ where+ x = realToFrac $ (x' - min') / (max'-min')++drawKuttas :: (Real a, Fractional a, Storable a) => FlowSol a -> VisObject a+drawKuttas flow = VisObjects $ concatMap (\(k0,k1) -> [circ k0, circ k1]) kis+ where+ kis = solKuttaIndices flow+ (xs',ys') = unzip $ map fMidpoints $ (\(Foil els _) -> els) (solFoil flow)+ xs = vjoin xs'+ ys = vjoin ys'+ circ k = drawCircle yellow (xs @> k, ys @> k) 0.006++drawForces :: FlowSol Double -> VisObject Double+drawForces flow = VisObjects $ map (\(xy0, xy1, cp) -> drawLine (colorFun minCp maxCp cp) [xy0, xy1])+ $ zip3 xy0s xy1s (toList (solCps flow))+ where+ xy0s = zip (toList xm) (toList ym)+ xy1s = zip (toList (xm + xPressures)) (toList (ym + yPressures))+ (xPressures, yPressures) = (\(x,y) -> (LA.scale c x/lengths, LA.scale c y/lengths)) (solForces flow)+ lengths = vjoin $ map fLengths $ (\(Foil x _) -> x) $ solFoil flow+ (xm, ym) = (\(x,y) -> (vjoin x, vjoin y)) $ unzip $ map fMidpoints $ (\(Foil x _) -> x) $ solFoil flow++ c = 0.1++ maxCp = maxElement (solCps flow)+ minCp = minElement (solCps flow)+++drawColoredFoil :: (Num a, Storable a) => [Color] -> Foil a -> VisObject a+drawColoredFoil colors foil@(Foil elements _) = VisObjects $ zipWith drawColoredElement colors' elements+ where+ colors' = groupSomethingByFoil foil colors++drawColoredElement :: (Num a, Storable a) => [Color] -> Element a -> VisObject a+drawColoredElement colors element = VisObjects $ map (\(xy0, xy1, col) -> drawLine col [xy0, xy1]) (zip3 xy0s xy1s colors)+ where+ xys = (\(x,y) -> zip (toList x) (toList y)) $ fNodes element+ xy0s = tail xys+ xy1s = init xys++groupSomethingByFoil :: Storable a => Foil a -> [b] -> [[b]]+groupSomethingByFoil (Foil elements _) somethings = f somethings (map (LA.dim . fAngles) elements)+ where+ f xs (n:ns) = (take n xs):(f (drop n xs) ns)+ f [] []= []+ f _ _ = error "uh oh (groupSomethingByFoil)"++drawSolution :: FlowSol Double -> VisObject Double+drawSolution flow = VisObjects $ onscreenText +++ [ drawColoredFoil colors foil+ , drawCircle white (fst $ solCenterPressure flow, snd $ solCenterPressure flow) 0.006+ , drawCircle white (fst $ solCenterPressure flow, 0) 0.006+ , drawCircle green (0.25,0) 0.006+ ] ++ zipWith (\x y -> drawLineV red (x, y)) xs+ (takesV (map LA.dim xs) (LA.scale cpScale cps)) -- cp graph+ where+ foil@(Foil elements name) = solFoil flow+ cps = solCps flow++ xs = map (fst . fMidpoints) elements++ onscreenText =+ zipWith (\s k -> Text2d s (30,fromIntegral $ 30*k) Fixed9By15 (makeColor 1 1 1 1))+ msgs (reverse [1..length msgs])++ msgs = [ name+ , printf ("alpha: %.6f deg") ((solAlpha flow)*180/pi)+ , printf ("Cl: %.6f") (solCl flow)+ , printf ("Cd: %.6f") (solCd flow)+ , printf ("Cm: %.6f (c/4, 0)") (solCm flow)+ ]++ colors = map (colorFun (minElement cps) (maxElement cps)) (toList cps)+++drawOnce :: Real a => [VisObject a] -> IO ()+drawOnce pics = display (defaultOpts {optWindowName = "hfoil"}) (VisObjects pics)++-- let line = plot_lines_values ^= [[ (xc, yt (naca4 "0012") xc)+-- | xc <- [0,0.01..0.99::Double]]]+-- $ plot_lines_title ^= "naca 0012"+-- $ defaultPlotLines+--+-- chart = layout1_title ^= "naca yo"+-- $ layout1_plots ^= [Left (toPlot line)]+-- $ defaultLayout1+--+-- renderableToWindow (toRenderable chart) 640 480+-- _ <- renderableToPNGFile (toRenderable chart) 640 480 "mDiv_vs_tc.png"+-- return ()
+ src/HFoil/Flow.hs view
@@ -0,0 +1,205 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE FlexibleContexts #-}++module HFoil.Flow+ ( FlowSol(..)+ , solveFlow+ ) where++import Data.Packed.ST (newMatrix, writeMatrix, freezeMatrix)+import Control.Monad ( forM_ )+import Control.Monad.ST ( runST )+import Numeric.LinearAlgebra hiding(i)+import Foreign.Storable ( Storable )++import HFoil.Foil++data FlowSol a = FlowSol { solFoil :: Foil a+ , solVs :: Vector a+ , solCps :: Vector a+ , solVorticities :: [a]+ , solAlpha :: a+ , solForces :: (Vector a, Vector a)+ , solForce :: (a,a)+ , solCl :: a+ , solCd :: a+ , solCm :: a -- moment about (0,0.25)+ , solCenterPressure :: (a,a)+ , solKuttaIndices :: [(Int,Int)]+ }++solveFlow :: (Num (Vector a), RealFloat a, Field a) => Foil a -> a -> FlowSol a+solveFlow foil@(Foil elements _) alpha =+ FlowSol { solFoil = foil+ , solVs = vs+ , solCps = cps+ , solVorticities = vorticities+ , solAlpha = alpha+ , solForces = (xForces, yForces)+ , solForce = (xf,yf)+ , solCl = cl+ , solCd = cd+ , solCm = cm+ , solCenterPressure = (xCp, yCp)+ , solKuttaIndices = kuttaIndices+ }+ where+ -- surface velocities+ (vs,vorticities) = ( (mapVector (\q -> cos(q - alpha)) angles) + (mV <> qsGamma)+ , map (qsGamma @>) [n-(length elements),n-1]+ )+ where+ (mA, mV, b) = getAVb geometries angles kuttaIndices alpha+ qsGamma = flatten $ linearSolve mA b+ n = dim angles+ -- pressure coefficients+ cps = 1 - vs*vs++ -- forces and force coefficients+ xForces = -cps*(vjoin $ map (fst . fNormals) elements)+ yForces = -cps*(vjoin $ map (snd . fNormals) elements)+ xf = sumElements xForces+ yf = sumElements yForces++ cd = xf*(cos alpha) + yf*(sin alpha)+ cl = -xf*(sin alpha) + yf*(cos alpha)++ -- midpoints+ (xs,ys) = (\(x,y) -> (vjoin x, vjoin y)) $ unzip $ map fMidpoints elements++ -- centers of pressure+ (xCp, yCp) = ( (sumElements (xs*yForces)) / (sumElements yForces)+ , (sumElements (ys*xForces)) / (sumElements xForces)+ )++ -- moment about cp and quarter chord+ cm = sumElements $ (mapVector (\z -> z - 0.25) xs)*yForces - ys*xForces++ kuttaIndices = ki 0 (map dim angles')+ where+ ki _ [] = []+ ki n0 (n:ns) = (n0,n0+n-1):(ki (n0+n) ns)++ angles = vjoin angles'+ angles' = map fAngles elements++ -- (mSL,mCL,mSB,mCB)+ geometries = setupGeometries angles (xms, yms) (xns, yns) (xnps, ynps)+ where+ vjoinTuples (x,y) = (vjoin x, vjoin y)+ (xms,yms) = vjoinTuples $ unzip $ map fMidpoints elements+ (xns,yns) = vjoinTuples $ unzip $ map fInits elements+ (xnps,ynps) = vjoinTuples $ unzip $ map fTails elements+++-- make influence matrix A and also the matrix V where V*[sources; vortex] == tangential speeds+getAVb :: (Floating a, Num (Vector a), Container Vector a, Storable a) =>+ (Matrix a, Matrix a, Matrix a, Matrix a)+ -> Vector a+ -> [(Int, Int)]+ -> a+ -> (Matrix a, Matrix a, Matrix a)+getAVb (mSL,mCL,mSB,mCB) angles kuttaIndices alpha =+ ( scale (1/(2*pi)) $ fromBlocks [[ mAij, fromColumns vsAin]+ ,[fromRows vsAnj, fromLists ann]]+ , scale (1/(2*pi)) $ fromBlocks [[ mSB - mCL+ , fromColumns vsTin]]+ , asColumn $ vjoin $ (mapVector (\q -> sin $ q - alpha) angles):+ (map (\(n0,nf) -> fromList [-cos((angles @> n0) - alpha) - cos((angles @> nf) - alpha)]) kuttaIndices)+ ) -- (for middle entry mV, mSL + mCB == mAij)+ where+ n = rows mSL++ -- sources influence matrix Aij+ mAij = mSL + mCB++ -- vortices influence vectors Ains+ vsAin = map (\(n0,nf) -> fromList $ map sumElements $ toRows $ subMatrix (0,n0) (n, nf-n0+1) m) kuttaIndices+ where+ m = (mCL - mSB)++ -- vortices tangential velocity outputs+ vsTin = map (\(n0,nf) -> fromList $ map sumElements $ toRows $ subMatrix (0,n0) (n, nf-n0+1) m) kuttaIndices+ where+ m = (mSL + mCB)++ -- sources kutta condition influence vector Anj+ vsAnj = map (\(n0,nf) -> (getRow n0 mSB) - (getRow n0 mCL) + (getRow nf mSB) - (getRow nf mCL)) kuttaIndices++ -- vortices kutta condition influence scalars ann+ ann = (flip map) kuttaIndices $ \(ni0,nif) -> (flip map) kuttaIndices $ \(nj0,njf) ->+ sumElements $ (getSubRow ni0 (nj0,njf) mSL) + (getSubRow ni0 (nj0,njf) mCB) + (getSubRow nif (nj0,njf) mSL) + (getSubRow nif (nj0,njf) mCB)++ getRow i mat = flatten $ subMatrix (i,0) (1,n) mat+ getSubRow i (j0,jf) mat = flatten $ subMatrix (i,j0) (1,jf-j0+1) mat+++{-+calcuate 4 matrices which will be useful++mSL(i,j) = sin(qi - ij) * ln(rij+1/rij)+mCL(i,j) = cos(qi - ij) * ln(rij+1/rij)+mSB(i,j) = sin(qi - ij) * beta(i,j)+mCB(i,j) = cos(qi - ij) * beta(i,j)++where for i==j: ln(r/r) = 0+ beta = pi+are explicitly set+-}+setupGeometries :: (RealFloat t, Storable t) =>+ Vector t+ -> (Vector t, Vector t)+ -> (Vector t, Vector t)+ -> (Vector t, Vector t)+ -> (Matrix t, Matrix t, Matrix t, Matrix t)+setupGeometries angles (xms, yms) (xns, yns) (xnps, ynps) = runST $ do+ let n = dim angles+ mSL' <- newMatrix 0 n n+ mCL' <- newMatrix 0 n n+ mSB' <- newMatrix 0 n n+ mCB' <- newMatrix 0 n n++ forM_ [0..n-1] $ \i -> forM_ [0..n-1] $ \j -> do+ let qi = angles @> i+ qj = angles @> j++ xmi = xms @> i+ ymi = yms @> i+ xnj = xns @> j+ ynj = yns @> j+ xnjp = xnps @> j+ ynjp = ynps @> j++ s = sin (qi - qj)+ c = cos (qi - qj)++ l+ | i == j = 0+ | otherwise = log(r1/r0)+ where+ r1 = distance (xmi, ymi) (xnjp, ynjp)+ r0 = distance (xmi, ymi) (xnj, ynj )+ distance (x1,y1) (x0,y0) = sqrt $ dx*dx + dy*dy+ where+ dx = x1 - x0+ dy = y1 - y0++ b+ | i == j = pi+ | otherwise = atan2 (dyjp*dxj - dxjp*dyj)+ (dxjp*dxj + dyjp*dyj)+ where+ dyj = ymi - ynj+ dxj = xmi - xnj+ dyjp = ymi - ynjp+ dxjp = xmi - xnjp++ _ <- writeMatrix mSL' i j (s*l)+ _ <- writeMatrix mCL' i j (c*l)+ _ <- writeMatrix mSB' i j (s*b)+ writeMatrix mCB' i j (c*b)+ sl' <- freezeMatrix mSL'+ cl' <- freezeMatrix mCL'+ sb' <- freezeMatrix mSB'+ cb' <- freezeMatrix mCB'+ return (sl',cl',sb',cb')
+ src/HFoil/Foil.hs view
@@ -0,0 +1,246 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE FlexibleContexts #-}++module HFoil.Foil+ ( Foil(..)+ , Element(..)+ , panelizeNaca4+ , loadFoil+ , getUIUCFoil+ ) where++import System.Directory ( doesFileExist )+import Numeric.LinearAlgebra hiding ( Element )+import Network.HTTP ( simpleHTTP, getRequest, getResponseBody )+import Foreign.Storable ( Storable )+import Text.Read ( readMaybe )+import Text.Parsec ( Parsec, ParsecT, Stream, (<?>), parse+ , option, char, try, lookAhead, manyTill, digit, endOfLine+ , many1, optional, anyChar, sepBy1, many, oneOf, skipMany )++import qualified HFoil.Naca4 as Naca4++data Foil a = Foil [Element a] String++instance (Storable a) => Show (Foil a) where+ show (Foil [el] name) = "{"++name++": " ++ show (1 + dim (fLengths el)) ++ " nodes}"+ show (Foil els name) = "{"++name++": " ++ show (length els) ++ " elements, " +++ show nodesPerEl ++ " nodes == "++show (sum nodesPerEl)++" total nodes}"+ where+ nodesPerEl = map (\x -> 1 + dim (fLengths x)) els++data Element a = Element { fNodes :: (Vector a, Vector a)+ , fLengths :: Vector a+ , fAngles :: Vector a+ , fMidpoints :: (Vector a, Vector a)+ , fTangents :: (Vector a, Vector a)+ , fNormals :: (Vector a, Vector a)+ , fUnitNormals :: (Vector a, Vector a)+ , fInits :: (Vector a, Vector a)+ , fTails :: (Vector a, Vector a)+ }+++-- make sure the nodes aren't reversed+toElement :: (Num (Vector a), RealFloat a, Container Vector a)+ => [(a, a)] -> Element a+toElement nodes+ | diff < 0 = el+ | otherwise = toElement' (reverse nodes)+ where+ el = toElement' nodes+ angles = toList (fAngles el)+ diff = sum $ zipWith (-) angles (drop 1 angles)++toElement' :: (Num (Vector a), RealFloat a, Container Vector a)+ => [(a, a)] -> Element a+toElement' xynodes =+ Element { fNodes = (xNodes, yNodes)+ , fLengths = lengths+ , fAngles = zipVectorWith atan2 yTangents xTangents+ , fMidpoints = (xMids, yMids)+ , fTangents = (xTangents, yTangents)+ , fNormals = (xNormals, yNormals)+ , fUnitNormals = (xUnitNormals, yUnitNormals)+ , fInits = (xInits, yInits)+ , fTails = (xTails, yTails)+ }+ where+ n = (dim xNodes) - 1+ (xNodes, yNodes) = (\(xs,ys) -> (fromList xs, fromList ys)) $ unzip xynodes+ (xInits, yInits) = (subVector 0 n xNodes, subVector 0 n yNodes)+ (xTails, yTails) = (subVector 1 n xNodes, subVector 1 n yNodes)+ (xTangents, yTangents) = (xTails - xInits, yTails - yInits)+ (xMids, yMids) = (0.5*(xInits + xTails), 0.5*(yInits + yTails))+ (xNormals, yNormals) = (-yTangents, xTangents)+ lengths = mapVector sqrt $ xTangents*xTangents + yTangents*yTangents+ (xUnitNormals, yUnitNormals) = (xNormals/lengths, yNormals/lengths)+++-- why isn't this standard???+poorMansStrip :: String -> String+poorMansStrip str = reverse $ dropWhile (== ' ') $ reverse $ dropWhile (== ' ') str++getUIUCFoil :: String -> IO (Either String (Foil Double))+getUIUCFoil name' = do+ let name = poorMansStrip name'+ file = "http://m-selig.ae.illinois.edu/ads/coord/" ++ name ++ ".dat"+ dl <- simpleHTTP (getRequest file) >>= getResponseBody+ return (parseRawFoil dl name)++loadFoil :: FilePath -> IO (Either String (Foil Double))+loadFoil filename' = do+ let filename = poorMansStrip filename'+ exists <- doesFileExist filename+ if exists+ then do rawData <- readFile filename+ return (parseRawFoil rawData filename) -- use filename as name+ else return (Left ("file \"" ++ filename ++ "\" couldn't be found"))++fst3 :: (a,b,c) -> a+fst3 (x,_,_) = x++panelizeNaca4 :: (Enum a, Floating (Vector a), RealFloat a, Field a) =>+ Naca4.Naca4 a -> Int -> Foil a+panelizeNaca4 foil nPanels = Foil [toElement $ [(1,0)]++reverse lower++[c0]++upper++[(1,0)]]+ (Naca4.naca4_name foil)+ where+ c0 = fst (Naca4.coords foil 0)+ (upper, lower) = unzip $ map (Naca4.coords foil) xcs+ xcs = toList $ fst3 $ bunchPanels (Naca4.yt foil) (Naca4.dyt foil) xcs0 0 0+ xcs0 = fromList $ init $ tail $ toList $ linspace nXcs (0,1)+ nXcs = (nPanels + (nPanels `mod` 2)) `div` 2 + 1++bunchPanels :: forall a . (Enum a, Floating (Vector a), Floating a, Ord a, Field a) =>+ (a -> a) -> (a -> a) -> Vector a -> Int -> Int -> (Vector a, Int, Int)+bunchPanels yt dyt xcs nIter nBadSteps+ | nIter > 300 = error "panel buncher exceeded 300 iterations"+ | nBadSteps > 1000 = error "panel buncher exceeded 1000 bad steps"+ | sum (toList (abs deltaXcs)) < 1e-12 = (xcs, nIter, nBadSteps)+ | otherwise = bunchPanels yt dyt goodXcs (nIter+1) (nBadSteps + length badOnes)+ where+ (badOnes, goodXcs:_) = break (\xs -> all (>0) (toList xs)) nextXcs+ nextXcs :: [Vector a]+ nextXcs = map (\alpha -> xcs + (scale alpha deltaXcs)) (map (2.0**) [0,-1..])+ deltaXcs :: Vector a+ deltaXcs = xcsStep yt dyt xcs++xcsStep :: (Enum a, Floating (Vector a), Field a) => (a -> a) -> (a -> a) -> Vector a -> Vector a+xcsStep yt dyt xcs = flatten $ -(linearSolveLS mat2 (asColumn rs))+ where+ n = dim xcs++ -- x coords+ xs = vjoin [fromList [0], xcs, fromList [1]]+ -- y coords+ ys = vjoin [fromList [0], mapVector yt xcs, fromList [0]]+ -- x deltas+ dxs = (subVector 1 (n+1) xs) - (subVector 0 (n+1) xs)+ -- y deltas+ dys = (subVector 1 (n+1) ys) - (subVector 0 (n+1) ys)+ -- delta magnitures+ deltas = sqrt (dxs*dxs + dys*dys)+ -- slopes+ dydxs = mapVector dyt xcs++ mat = r1 - r0+ where+ r0 = fromBlocks [[(1><n) [0,0..]], [diagRect 0 diag0 n n]]+ r1 = fromBlocks [ [diagRect 0 diag1 n n], [(1><n)[0,0..]]]++ (diag0, diag1) = (subVector 1 n d0, subVector 0 n d1)+ where+ d0 = (dxs + dys*dy0dxs)/deltas+ d1 = (dxs + dys*dy1dxs)/deltas+ dy0dxs = vjoin [fromList [0], dydxs]+ dy1dxs = vjoin [dydxs, fromList [0]]++ zeros = (n><1)[0,0..]+ eye = ident n+ frontBunchingParam = 2.0+ diff = (fromBlocks [[zeros, eye]]) - (fromBlocks [[eye*(1+frontBunchingParam/(fromIntegral n)), zeros]])++ mat2 = diff <> mat+ rs = diff <> deltas++discardIdenticalNeighbors :: Eq a => [a] -> [a]+discardIdenticalNeighbors (x0:x1:xs)+ | x0 == x1 = discardIdenticalNeighbors (x1:xs)+ | otherwise = x0 : discardIdenticalNeighbors (x1:xs)+discardIdenticalNeighbors xs = xs++stitchIncompleteElements :: [[(Double, Double)]] -> [[(Double,Double)]]+stitchIncompleteElements (as:bs:others)+ | and [ x0a < 0.5+ , x0b < 0.5+ , 0.5 < xFa+ , 0.5 < xFb+ ] = (discardIdenticalNeighbors ((reverse as) ++ bs)) : stitchIncompleteElements others+ | and [ 0.5 < x0a+ , 0.5 < x0b+ , xFa < 0.5+ , xFb < 0.5+ ] = (discardIdenticalNeighbors (as ++ reverse bs)) : stitchIncompleteElements others+ | otherwise = as : bs : stitchIncompleteElements others+ where+ x0a = fst $ head as+ xFa = fst $ last as++ x0b = fst $ head bs+ xFb = fst $ last bs+stitchIncompleteElements [] = []+stitchIncompleteElements xs = xs++parseRawFoil :: String -> String -> Either String (Foil Double)+parseRawFoil raw name = foil+ where+ foil = case parse foilP name raw of+ Left pe -> Left (raw ++ "\nError parsing the above data: " ++ (show pe))+ Right (_, els) -> Right (Foil (map toElement (stitchIncompleteElements els)) name)++headerP :: Parsec String () String+headerP = manyTill anyChar (try (lookAhead elementsP))++foilP :: Parsec String () (String, [[(Double, Double)]])+foilP = do+ s <- headerP <?> "header"+ els <- elementsP <?> "elements"+ return (s, els)++elementP :: Parsec String () [(Double, Double)]+elementP = many2 coordP+ where+ many2 :: (Stream s m t) => ParsecT s u m a -> ParsecT s u m [a]+ many2 p = do+ x0 <- p+ x1 <- p+ xs <- many p+ return (x0:x1:xs)+++elementsP :: Parsec String () [[(Double, Double)]]+elementsP = sepBy1 elementP endOfLine++doubleP :: Parsec String () Double+doubleP = do+ neg <- option ' ' (char '-')+ leading <- option "0" (many1 digit)+ dot' <- char '.'+ trailing <- option "0" (many1 digit)+ let doublish = neg : leading ++ dot' : trailing+ case readMaybe doublish of+ Just x -> return x+ Nothing -> error $ "failed to read this supposed double: " ++ show doublish++coordP :: Parsec String () (Double, Double)+coordP = do+ let space' = oneOf [' ', '\t']+ spaces' = skipMany space' <?> "space-like"+ spaces' <?> "leading whitespace"+ x <- doubleP <?> "first coordinate"+ spaces' <?> "middle whitespace"+ y <- doubleP <?> "second coordinate"+ spaces' <?> "trailing whitespace"+ _ <- optional endOfLine <?> "end of line"+ return (x,y)
+ src/HFoil/Naca4.hs view
@@ -0,0 +1,71 @@+{-# OPTIONS_GHC -Wall #-}++module HFoil.Naca4+ ( Naca4(..)+ , coords+ , yt+ , dyt+ , naca4+ ) where++import Text.Read ( readMaybe )++naca4 :: (Read a, Fractional a) => String -> Maybe (Naca4 a)+naca4 name@(m_:p_:t0:t1:[]) = do+ m <- fmap (0.01 *) $ readMaybe [m_]+ p <- fmap (0.1 *) $ readMaybe [p_]+ t <- fmap (0.01 *) $ readMaybe [t0,t1]+ return $ Naca4 m p t ("NACA " ++ name)+naca4 _ = Nothing++-- m: max camber in hundredths of chord+-- p: position of max camber in tenths of chord+-- t: max thickness in hundredths of chord+data Naca4 a = Naca4 { naca4_m :: a+ , naca4_p :: a+ , naca4_t :: a+ , naca4_name :: String+ } deriving Show++-- xc: x/chord+yc :: (Ord a, Fractional a) => Naca4 a -> a -> a+yc (Naca4 {naca4_p = p, naca4_m = m}) xc+ | xc < 0 = error "xc < 0"+ | xc < p = m/(p*p)*(2*p*xc - xc*xc)+ | xc <= 1 = m/((1-p)*(1-p))*((1-2*p) + 2*p*xc - xc*xc)+ | otherwise = error "xc > 1"++dyc :: (Ord a, Fractional a) => Naca4 a -> a -> a+dyc (Naca4 {naca4_p = p, naca4_m = m}) xc+ | xc < 0 = error "xc < 0"+ | xc < p = m/(p*p)*(2*p - 2*xc)+ | xc <= 1 = m/((1-p)*(1-p))*(2*p - 2*xc)+ | otherwise = error "xc > 1"++yt :: (Ord a, Floating a) => Naca4 a -> a -> a+yt (Naca4 {naca4_t = t}) xc+ | xc < 0 = error "xc < 0"+ | xc > 1 = error $ "xc > 1"+ | otherwise = 5*t*(0.2969*sqrt(xc) - 0.1260*(xc) - 0.3537*(xc)**2 + 0.2843*(xc)**3 - 0.1015*(xc)**4)++dyt :: (Ord a, Floating a) => Naca4 a -> a -> a+dyt (Naca4 {naca4_t = t}) xc+ | xc < 0 = error "xc < 0"+ | xc > 1 = error "xc > 1"+ | otherwise = 5*t*(0.5*0.2969/sqrt(xc) - 0.1260 - 2*0.3537*(xc) + 3*0.2843*(xc)**2 - 4*0.1015*(xc)**3)++coords :: (Ord a, Floating a) => Naca4 a -> a -> ((a,a),(a,a))+coords foil xc+ | naca4_m foil == 0 = ((xc,yt_), (xc,-yt_))+ | otherwise = ((xu,yu ), (xl, yl ))+ where+ yt_ = yt foil xc+ yc_ = yc foil xc++ yu = yc_ + yt_ * (cos theta)+ yl = yc_ - yt_ * (cos theta)++ xu = xc - yt_ * (sin theta)+ xl = xc + yt_ * (sin theta)++ theta = atan $ dyc foil xc
+ src/HFoil/Repl.hs view
@@ -0,0 +1,226 @@+{-# OPTIONS_GHC -Wall #-}+-- {-# Language FlexibleContexts #-}++module HFoil.Repl+ ( run+ ) where++import Control.Concurrent ( forkIO )+import Control.Concurrent.MVar ( newMVar, readMVar, swapMVar )+import Control.Monad.IO.Class ( MonadIO, liftIO )+import Control.Monad.Trans.Class ( lift )+import Control.Monad.Trans.State.Strict ( StateT, evalStateT, get, modify )+import Data.List ( isPrefixOf )+import Linear ( Quaternion(..), V3(..) )+import System.Console.Haskeline ( InputT, runInputT, defaultSettings, getInputLine, outputStrLn, setComplete )+import System.Console.Haskeline.Completion ( CompletionFunc, Completion, completeWord, simpleCompletion )+import Text.Read ( readMaybe )++import Vis++import HFoil.Foil+import HFoil.Naca4+import HFoil.Drawing+import HFoil.Flow++nPanels :: Int+nPanels = 200++-- configuration+data Config = Config { confForces :: Bool+ , confKuttas :: Bool+ , confNormals :: Bool+ }++defaultConfig :: Config+defaultConfig = Config { confForces = False+ , confKuttas = False+ , confNormals = False+ }++data Mode = TopMode | FoilMode++foilCommands :: [(String, String)]+foilCommands =+ [ ("alfa", "alfa [#]")+ , ("forces", "forces")+ , ("kuttas", "kuttas")+ , ("normals", "normals")+ , ("help", "help")+ ]++topCommands :: [(String, String)]+topCommands =+ [ ("naca", "naca xxxx")+ , ("load", "load [filename]")+ , ("uiuc", "uiuc [foil name]")+ ]++topHelp :: InputT (StateT Mode IO) ()+topHelp = mapM_ (outputStrLn . snd) topCommands++foilHelp :: InputT (StateT Mode IO) ()+foilHelp = mapM_ (outputStrLn . snd) foilCommands++comp :: CompletionFunc (StateT Mode IO)+comp = completeWord Nothing " \t" searchFunc+ where+ searchFunc :: String -> (StateT Mode IO) [Completion]+ searchFunc str = do+ mode <- get+ let wordList = case mode of+ TopMode -> map fst topCommands+ FoilMode -> map fst foilCommands+ return $ map simpleCompletion $ filter (str `isPrefixOf`) wordList++run :: IO ()+run = do+ mpics <- newMVar $ []++ putStrLn "Welcome to hfoil\n"++ let go :: InputT (StateT Mode IO) ()+ go = topLoop (\pics -> swapMVar mpics pics >>= (\_ -> return ())) defaultConfig++ settings = setComplete comp defaultSettings+ _ <- forkIO $ flip evalStateT TopMode $ runInputT settings go++ let toScreen xs =+ RotQuat (Quaternion 0 (V3 1 0 0))+ $ Trans (V3 (-0.5) 0 0)+ $ VisObjects xs+ -- $ VisObjects (Axes (0.3,15) : xs)+ cam0 =+ Camera0+ { phi0 = 90+ , theta0 = 90+ , rho0 = 2+ }++ animateIO+ (defaultOpts {optWindowName = "hfoil", optInitialCamera = Just cam0})+ (\_ -> fmap toScreen (readMVar mpics))++data FoilState =+ FoilState+ { fsFlowSol :: Maybe (FlowSol Double)+ , fsConf :: Config+ , fsFoil :: Foil Double+ }++drawPicture :: MonadIO m => ([VisObject Double] -> IO ()) -> StateT FoilState m ()+drawPicture draw = do+ fs <- get+ let conf = fsConf fs+ foil = fsFoil fs+ normals = case confNormals conf of+ True -> [drawNormals foil]+ False -> []++ case fsFlowSol fs of+ Nothing -> liftIO $ draw (drawFoil foil: normals)+ Just flow -> do+ let forces = case (confForces conf) of+ True -> [drawForces flow]+ False -> []+ kuttas = case (confKuttas conf) of+ True -> [drawKuttas flow]+ False -> []+ liftIO $ draw $ forces++kuttas++normals++[drawSolution flow]++strip :: String -> String+strip = rstrip . lstrip+ where+ lstrip (' ':xs) = lstrip xs+ lstrip x = x++ rstrip = reverse . lstrip . reverse++foilLoop :: ([VisObject Double] -> IO ()) -> StateT FoilState (InputT (StateT Mode IO)) ()+foilLoop draw = do+ lift (lift (modify (const FoilMode)))+ fs <- get+ let foil@(Foil _ name) = fsFoil fs+ conf = fsConf fs+ drawPicture draw+ minput <- lift $ getInputLine $ "\ESC[1;32m\STXhfoil."++name++">> \ESC[0m\STX"++ case fmap strip minput of+ Nothing -> return ()+ Just "quit" -> do lift $ outputStrLn "not-gloss won't let you quit :(\ntry ctrl-c or hit ESC in drawing window"+ foilLoop draw+ Just ('a':'l':'f':'a':' ':alphaDeg') -> do+ case readMaybe alphaDeg' of+ Nothing -> do lift $ outputStrLn $ "parse fail on " ++ show alphaDeg'+ foilLoop draw+ Just alphaDeg -> do let flow :: FlowSol Double+ flow = solveFlow foil (pi/180*alphaDeg)+ modify (\fs' -> fs' {fsFlowSol = Just flow})+ foilLoop draw+ Just "forces" -> do+ let newConf = conf {confForces = not (confForces conf)}+ lift $ outputStrLn $ "force drawing set to "++ show (not (confForces conf))+ modify (\fs' -> fs' {fsConf = newConf})+ foilLoop draw+ Just "kuttas" -> do+ let newConf = conf {confKuttas = not (confKuttas conf)}+ lift $ outputStrLn $ "kutta drawing set to "++ show (not (confKuttas conf))+ modify (\fs' -> fs' {fsConf = newConf})+ foilLoop draw+ Just "normals" -> do+ let newConf = conf {confNormals = not (confNormals conf)}+ lift $ outputStrLn $ "normals drawing set to "++ show (not (confNormals conf))+ modify (\fs' -> fs' {fsConf = newConf})+ foilLoop draw+ Just "help" -> lift foilHelp >> foilLoop draw+ Just "h" -> lift foilHelp >> foilLoop draw+ Just "?" -> lift foilHelp >> foilLoop draw+ Just "" -> return ()+ Just input -> do lift $ outputStrLn $ "unrecognized command \"" ++ input ++ "\""+ foilLoop draw+++topLoop :: ([VisObject Double] -> IO ()) -> Config -> InputT (StateT Mode IO) ()+topLoop draw conf = do+ lift (modify (const TopMode))+ minput <- getInputLine "\ESC[1;32m\STXhfoil>> \ESC[0m\STX"+ case minput of+ Nothing -> return ()+ Just msg -> do runTop draw conf msg+ topLoop draw conf++runTop :: ([VisObject Double] -> IO ()) -> Config -> String -> InputT (StateT Mode IO) ()+runTop draw conf msg = case strip msg of+ "quit" -> outputStrLn "not-gloss won't let you quit :(\ntry ctrl-c or hit ESC in drawing window"+ ('n':'a':'c':'a':' ':spec) -> do+ case naca4 spec :: Maybe (Naca4 Double) of+ Nothing -> outputStrLn "not a valid naca4"+ Just n4 -> runFoil draw conf (panelizeNaca4 n4 nPanels)+ ('l':'o':'a':'d':' ':name) -> do+ mfoil <- liftIO (loadFoil name)+ case mfoil of Left errMsg -> outputStrLn errMsg+ Right foil -> do runFoil draw conf foil+ ('u':'i':'u':'c':' ':name) -> do+ efoil <- liftIO (getUIUCFoil name)+ case efoil of Left errMsg -> outputStrLn errMsg+ Right foil -> do let Foil els _ = foil+ outputStrLn $ "got " ++ show (length els) ++ " elements"+ runFoil draw conf foil+ "help" -> topHelp+ "h" -> topHelp+ "?" -> topHelp+ "" -> return ()+ other -> outputStrLn $ "unrecognized command \"" ++ other ++ "\""+++runFoil :: ([VisObject Double] -> IO ()) -> Config -> Foil Double -> InputT (StateT Mode IO) ()+runFoil draw conf foil = do+ let state0 =+ FoilState+ { fsFlowSol = Nothing+ , fsConf = conf+ , fsFoil = foil+ }+ let go :: StateT FoilState (InputT (StateT Mode IO)) ()+ go = foilLoop draw+ flip evalStateT state0 go