resistor-cube 0.0.1 → 0.0.1.1
raw patch · 2 files changed
+37/−41 lines, 2 filesdep ~comfort-arraydep ~lapack
Dependency ranges changed: comfort-array, lapack
Files
- resistor-cube.cabal +9/−7
- src/Main.hs +28/−34
resistor-cube.cabal view
@@ -1,14 +1,16 @@ Name: resistor-cube-Version: 0.0.1+Version: 0.0.1.1 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.+ We obtain the voltages and currents+ by solving simultaneous linear equations according to this matrix. .- For a generalized version see the tests of the @linear-circuit@ package.+ * For an explanation see <http://code.henning-thielemann.de/bob2019/main.pdf>.+ .+ * For a generalized version see the tests of the @linear-circuit@ package. Homepage: http://hub.darcs.net/thielema/resistor-cube License: BSD3 License-File: LICENSE@@ -19,7 +21,7 @@ Cabal-Version: >=1.10 Source-Repository this- Tag: 0.0.1+ Tag: 0.0.1.1 Type: darcs Location: http://hub.darcs.net/thielema/resistor-cube @@ -30,8 +32,8 @@ Executable resistor-cube Main-is: Main.hs Build-Depends:- lapack >=0.2.2 && <0.3,- comfort-array >=0.3.1 && <0.4,+ lapack >=0.3 && <0.4,+ comfort-array >=0.4 && <0.5, base >=4.5 && <5 Hs-Source-Dirs: src Default-Language: Haskell2010
src/Main.hs view
@@ -9,8 +9,11 @@ import qualified Numeric.LAPACK.Matrix.Triangular as Triangular import qualified Numeric.LAPACK.Matrix as Matrix import qualified Numeric.LAPACK.Vector as Vector+import Numeric.LAPACK.Matrix.Triangular ((#%%%#))+import Numeric.LAPACK.Matrix ((#\|)) import Numeric.LAPACK.Format ((##)) +import qualified Data.Array.Comfort.Boxed as BoxedArray import qualified Data.Array.Comfort.Storable as Array import qualified Data.Array.Comfort.Shape as Shape import Data.Array.Comfort.Storable ((!))@@ -30,10 +33,10 @@ cornerShape :: CornerShape-cornerShape = (Shape.Enumeration, Shape.Enumeration, Shape.Enumeration)+cornerShape = Shape.static edgeShape :: EdgeShape-edgeShape = (Shape.Enumeration, Shape.Enumeration, Shape.Enumeration)+edgeShape = Shape.static {-@@ -43,55 +46,46 @@ -} voltageMatrix :: Matrix EdgeShape CornerShape Double voltageMatrix =- Matrix.fromRowMajor $- Array.sample- (edgeShape, cornerShape)- (\((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)-+ Matrix.fromRowArray cornerShape $+ fmap+ (\e ->+ Array.fromAssociations 0 cornerShape+ [(edgeCorner e C0, 1), (edgeCorner e C1, -1)]) $+ BoxedArray.indices edgeShape -selectCornerCoords :: Dim -> Corner -> ((Coord, Coord), Coord)-selectCornerCoords ed (cx,cy,cz) =+edgeCorner :: Edge -> Coord -> Corner+edgeCorner (ed,ex,ey) coord = case ed of- D0 -> ((cy, cz), cx)- D1 -> ((cx, cz), cy)- D2 -> ((cx, cy), cz)+ D0 -> (coord,ex,ey)+ D1 -> (ex,coord,ey)+ D2 -> (ex,ey,coord) sourceCorner, destCorner :: Corner sourceCorner = (C0,C0,C0) destCorner = (C1,C1,C1) resistances :: Vector.Vector EdgeShape Double-resistances = Vector.constant edgeShape 1+resistances = Vector.one edgeShape fullMatrix :: Triangular.Symmetric (():+:(EdgeShape:+:CornerShape)) Double fullMatrix =- Triangular.stackSymmetric- (Triangular.symmetricFromList MatrixShape.RowMajor () [0])- (Matrix.singleRow MatrixShape.RowMajor $- Vector.unit (edgeShape:+:cornerShape) (Right sourceCorner)) $- Triangular.stackSymmetric- (Triangular.diagonal MatrixShape.RowMajor resistances)- voltageMatrix- (Vector.constant- (MatrixShape.symmetric MatrixShape.RowMajor cornerShape) 0)+ let order = MatrixShape.RowMajor+ symmetricZero sh = Matrix.zero $ MatrixShape.symmetric order sh+ unit = Vector.unit (edgeShape:+:cornerShape) (Right sourceCorner)+ in (symmetricZero (), Matrix.singleRow order unit)+ #%%%#+ (Triangular.diagonal order resistances, voltageMatrix)+ #%%%#+ symmetricZero cornerShape main :: IO () main = do fullMatrix ## "%2.f"- let solutionVec =- Matrix.unliftColumn MatrixShape.ColumnMajor- (Triangular.solve fullMatrix) $- Vector.unit (():+:(edgeShape:+:cornerShape)) (Right (Right destCorner))- print $ - solutionVec ! Right (Right destCorner)+ let ix = Right (Right destCorner)+ let solutionVec = fullMatrix #\| Vector.unit (Triangular.size fullMatrix) ix+ print $ - solutionVec ! ix {- result: total resistance is 5/(2+2+2)