hfoil (empty) → 0.1.1
raw patch · 10 files changed
+809/−0 lines, 10 filesdep +HTTPdep +MissingHdep +basesetup-changed
Dependencies added: HTTP, MissingH, base, directory, gloss, haskeline, hmatrix, transformers
Files
- LICENSE +30/−0
- Main.hs +8/−0
- Numeric/HFoil.hs +26/−0
- Numeric/HFoil/Drawing.hs +161/−0
- Numeric/HFoil/Flow.hs +184/−0
- Numeric/HFoil/Foil.hs +157/−0
- Numeric/HFoil/Naca4.hs +68/−0
- Numeric/HFoil/Repl.hs +92/−0
- Setup.hs +2/−0
- hfoil.cabal +81/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c)2012, Greg Horn++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Greg Horn nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Main.hs view
@@ -0,0 +1,8 @@+{-# OPTIONS_GHC -Wall #-}++module Main where++import Numeric.HFoil.Repl++main :: IO ()+main = run
+ Numeric/HFoil.hs view
@@ -0,0 +1,26 @@+{- |+ 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 view
@@ -0,0 +1,161 @@+{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE FlexibleContexts #-}++module Numeric.HFoil.Drawing( drawLine+ , drawLineV+ , drawSolution+ , drawFoil+ , drawOnce+ , drawNormals+ , drawForces+ ) 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')++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 $ [ drawText white (0.45, 0.8) 0.15 m0+ , drawText white (0.45, 0.65) 0.15 m1+ , drawText white (0.45, 0.5) 0.15 m2+ , drawText white (0.45, 0.35) 0.15 m3+ , drawForces flow+ , drawColoredFoil colors foil+ , drawCircle white (fst $ solCenterPressure flow, snd $ solCenterPressure flow) 0.006+ , drawCircle white (fst $ solCenterPressure flow, 0) 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+ +-- ] ++ zipWith (\x y -> drawLine red (toList x, y)) xs (groupSomethingByFoil foil (toList (LA.scale cpScale cps))) -- cp graph+-- ++++ [m0,m1,m2,m3] = [ name+ , printf ("alpha: %.6f") ((solAlpha flow)*180/pi)+ , printf ("Cl: %.6f") (solCl flow)+ , printf ("Cd: %.6f") (solCd 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 view
@@ -0,0 +1,184 @@+{-# 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)+ , solCl :: a+ , solCd :: a+ , solCenterPressure :: (a,a)+ }++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)+ , solCl = cl+ , solCd = cd+ , solCenterPressure = (xCp, yCp)+ }+ 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)++ -- centers of pressure+ (xCp, yCp) = ( (sumElements (xs*yForces)) / (sumElements yForces)+ , (sumElements (ys*xForces)) / (sumElements xForces)+ )+ where+ (xs,ys) = (\(x,y) -> (join x, join y)) $ unzip $ map fMidpoints elements++ 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+++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 view
@@ -0,0 +1,157 @@+{-# 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 view
@@ -0,0 +1,68 @@+{-# 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 view
@@ -0,0 +1,92 @@+{-# 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++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 ())))++ 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 ()) -> Foil Double -> InputT IO ()+foilLoop draw 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 foil+ Just ('a':'l':'f':'a':' ':[]) -> do outputStrLn $ "unrecognized command"+ foilLoop draw foil+ Just ('a':'l':'f':'a':' ':alphaDeg) -> do let flow = solveFlow foil (pi/180*(read alphaDeg))+ liftIO $ draw $ [drawSolution flow]+ foilLoop draw foil+ Just "" -> return ()+ Just input -> do outputStrLn $ "unrecognized command \"" ++ input ++ "\""+ foilLoop draw foil++topLoop :: ([Picture] -> IO ()) -> InputT IO ()+topLoop draw = 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+ Just ('n':'a':'c':'a':' ':spec) -> do parseNaca draw spec+ topLoop draw+ 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 foil'+ topLoop draw+ 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 foil+ topLoop draw++ Just "" -> topLoop draw+ Just input -> do outputStrLn $ "unrecognized command \"" ++ input ++ "\""+ topLoop draw++parseNaca :: ([Picture] -> IO ()) -> String -> InputT IO ()+parseNaca draw str + | length str == 4 = do let foil = panelizeNaca4 (naca4 str :: Naca4 Double) nPanels+ liftIO $ draw [drawFoil foil, drawNormals foil]+ foilLoop draw foil+ | otherwise = do outputStrLn $ "Not 4 digits"+ return ()
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ hfoil.cabal view
@@ -0,0 +1,81 @@+Name: hfoil+Version: 0.1.1+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+Category: Numerical, Science+Build-type: Simple+Cabal-version: >=1.6+Description: {+Library and command line REPL with plotting to do simple inviscid hess-smith panel code.+.+Features include:+.+* Single and multi-element airfoils+.+* Naca 4-series support with Gauss-Newton paneling+.+* Broken UIUC database integration (type \"uiuc [foilname]\")+.+* Shameless xfoil ripoff for the relp/plotting+.+* Only works with development version of Gloss that's not yet on Hackage+.+* Haskeline interface with tab-completion+.+* And nothing else.+.+.+To get started, do cabal install or whatever, then run the \"hfoil\" binary.+.+Things to try: \"naca 2412\", \"alfa 4\", (hit enter before entering another airfoil), \"load [filename]\", \"uiuc e330\"+}++Flag repl+ Description: Build the command line interface+ Default: True++Library+ Exposed-modules: Numeric.HFoil+ Numeric.HFoil.Naca4+ Numeric.HFoil.Foil+ Numeric.HFoil.Drawing+ Numeric.HFoil.Flow+ Numeric.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++ Ghc-Options: -Wall+ Ghc-Prof-Options: -prof -auto-all+++Executable hfoil+ if flag(repl)+ Buildable: True+ else+ Buildable: False++ Main-Is: Main.hs+ -- (gloss needs -O2 -threaded)+ 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+ location: git://github.com/ghorn/hfoil.git+ tag: 0.1.1