diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2015, Henning Thielemann + +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 Henning Thielemann 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/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/resistor-cube.cabal b/resistor-cube.cabal new file mode 100644 --- /dev/null +++ b/resistor-cube.cabal @@ -0,0 +1,36 @@ +Name: resistor-cube +Version: 0.0 +Synopsis: Compute total resistance of a cube of resistors +Description: + This is an example of how to compute the total resistance + of a non-trivial circuit of resistors. + It demonstrates how to build the necessary matrix. + The computed voltages and currents + are elements of the null vector of that matrix. +Homepage: http://hub.darcs.net/thielema/resistor-cube +License: BSD3 +License-File: LICENSE +Author: Henning Thielemann +Maintainer: haskell@henning-thielemann.de +Category: Math +Build-Type: Simple +Cabal-Version: >=1.10 + +Source-Repository this + Tag: 0.0 + Type: darcs + Location: http://hub.darcs.net/thielema/resistor-cube + +Source-Repository head + Type: darcs + Location: http://hub.darcs.net/thielema/resistor-cube + +Executable resistor-cube + Main-is: Main.hs + Build-Depends: + hmatrix >=0.16 && <0.17, + transformers >=0.3 && <0.4, + utility-ht >=0.0 && <0.1, + base >=4.5 && <4.6 + Hs-Source-Dirs: src + Default-Language: Haskell2010 diff --git a/src/Main.hs b/src/Main.hs new file mode 100644 --- /dev/null +++ b/src/Main.hs @@ -0,0 +1,110 @@ +{- | +Consider a cube of resistors of equal resistance. +What is the overall resistance from one corner to the opposite one? +-} +module Main where + +import qualified Numeric.Container as NC +import qualified Data.Packed.Matrix as Matrix +import qualified Data.Packed.Vector as Vector +import qualified Numeric.LinearAlgebra.HMatrix as HMatrix +import Data.Packed.Matrix (Matrix) + +import Control.Applicative (liftA3) +import Control.Functor.HT (outerProduct) + + +{- +Set resistance of a primitive resistor to 1. +This way, voltage equals current. +-} +data Coord = C0 | C1 deriving (Eq, Ord, Show, Enum, Bounded) +data Dim = D0 | D1 | D2 deriving (Eq, Ord, Show, Enum, Bounded) +data Corner = Corner Coord Coord Coord deriving (Eq, Ord, Show) +data Edge = Edge Dim Coord Coord deriving (Eq, Ord, Show) + +allCoords :: [Coord] +allCoords = [minBound .. maxBound] + +allDims :: [Dim] +allDims = [minBound .. maxBound] + +flattenCornerIndex :: Corner -> Int +flattenCornerIndex (Corner x y z) = + (fromEnum x * 2 + fromEnum y) * 2 + fromEnum z + +flattenEdgeIndex :: Edge -> Int +flattenEdgeIndex (Edge d x y) = + (fromEnum d * 2 + fromEnum x) * 2 + fromEnum y + + +voltageMatrix :: Matrix Double +voltageMatrix = + Matrix.fromLists $ + outerProduct + (\(Edge ed ex ey) c -> + let ((cx, cy), cz) = selectCornerCoords ed c + in if ex==cx && ey==cy + then + case cz of + C0 -> 1 + C1 -> -1 + else 0) + (liftA3 Edge allDims allCoords allCoords) + (liftA3 Corner allCoords allCoords allCoords) + + +-- ToDo: How about cyclic arrangement of dimensions? +selectCornerCoords :: Dim -> Corner -> ((Coord, Coord), Coord) +selectCornerCoords ed (Corner cx cy cz) = + case ed of + D0 -> ((cy, cz), cx) + D1 -> ((cx, cz), cy) + D2 -> ((cx, cy), cz) + +sourceCorner, destCorner :: Corner +sourceCorner = Corner C0 C0 C0 +destCorner = Corner C1 C1 C1 + +currentMatrix :: Matrix Double +currentMatrix = + Matrix.fromLists $ + outerProduct + (\c@(Corner _ _ _) (Edge ed ex ey) -> + let ((cx, cy), cz) = selectCornerCoords ed c + in if ex==cx && ey==cy + then + case cz of + C0 -> 1 + C1 -> -1 + else 0) + (filter (/= sourceCorner) $ + filter (/= destCorner) $ + liftA3 Corner allCoords allCoords allCoords) + (liftA3 Edge allDims allCoords allCoords) + + +fullMatrix :: Matrix Double +fullMatrix = + Matrix.fromBlocks + [[NC.konst 0 (1,12), Matrix.asRow $ Vector.fromList $ 1 : replicate 7 0], + [HMatrix.ident 12, voltageMatrix], + [currentMatrix, NC.konst 0 (6,8)]] + +main :: IO () +main = do + print fullMatrix + let [currentVec, potentialVec] = + Vector.takesV [12,8] $ HMatrix.null1 fullMatrix + let totalCurrent = + sum $ + map + (\d -> NC.atIndex currentVec $ flattenEdgeIndex $ Edge d C0 C0) + [minBound .. maxBound] + let cornerPot c = NC.atIndex potentialVec (flattenCornerIndex c) + let totalVoltage = cornerPot destCorner - cornerPot sourceCorner + print $ totalVoltage / totalCurrent + +{- +result: total resistance is 5/(2+2+2) +-}