diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Main.hs b/Main.hs new file mode 100644 --- /dev/null +++ b/Main.hs @@ -0,0 +1,8 @@ +{-# OPTIONS_GHC -Wall #-} + +module Main where + +import Numeric.HFoil.Repl + +main :: IO () +main = run diff --git a/Numeric/HFoil.hs b/Numeric/HFoil.hs new file mode 100644 --- /dev/null +++ b/Numeric/HFoil.hs @@ -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 diff --git a/Numeric/HFoil/Drawing.hs b/Numeric/HFoil/Drawing.hs new file mode 100644 --- /dev/null +++ b/Numeric/HFoil/Drawing.hs @@ -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 () diff --git a/Numeric/HFoil/Flow.hs b/Numeric/HFoil/Flow.hs new file mode 100644 --- /dev/null +++ b/Numeric/HFoil/Flow.hs @@ -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') diff --git a/Numeric/HFoil/Foil.hs b/Numeric/HFoil/Foil.hs new file mode 100644 --- /dev/null +++ b/Numeric/HFoil/Foil.hs @@ -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><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 diff --git a/Numeric/HFoil/Naca4.hs b/Numeric/HFoil/Naca4.hs new file mode 100644 --- /dev/null +++ b/Numeric/HFoil/Naca4.hs @@ -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 diff --git a/Numeric/HFoil/Repl.hs b/Numeric/HFoil/Repl.hs new file mode 100644 --- /dev/null +++ b/Numeric/HFoil/Repl.hs @@ -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 () diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/hfoil.cabal b/hfoil.cabal new file mode 100644 --- /dev/null +++ b/hfoil.cabal @@ -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