packages feed

rsagl-math (empty) → 0.6.0.0

raw patch · 18 files changed

+2082/−0 lines, 18 filesdep +OpenGLdep +OpenGLRawdep +Vecsetup-changed

Dependencies added: OpenGL, OpenGLRaw, Vec, Vec-OpenGLRaw, array, base, containers, parallel, parsec, random

Files

+ LICENSE view
@@ -0,0 +1,26 @@+Copyright (c) 2006, 2007, Christopher Lane Hinson+ 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 Christopher Lane Hinson nor the names of its 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.
+ RSAGL/Math.hs view
@@ -0,0 +1,25 @@+module RSAGL.Math (+    module RSAGL.Math.AbstractVector,+    module RSAGL.Math.Affine,+    module RSAGL.Math.Angle,+    module RSAGL.Math.BoundingBox,+    module RSAGL.Math.Curve,+    module RSAGL.Math.Interpolation,+    module RSAGL.Math.Matrix,+    module RSAGL.Math.Orthogonal,+    module RSAGL.Math.Ray,+    module RSAGL.Math.Types,+    module RSAGL.Math.Vector)+  where++import RSAGL.Math.AbstractVector+import RSAGL.Math.Affine+import RSAGL.Math.Angle+import RSAGL.Math.BoundingBox+import RSAGL.Math.Curve+import RSAGL.Math.Interpolation+import RSAGL.Math.Matrix+import RSAGL.Math.Orthogonal+import RSAGL.Math.Ray+import RSAGL.Math.Types+import RSAGL.Math.Vector
+ RSAGL/Math/AbstractVector.hs view
@@ -0,0 +1,230 @@+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, UndecidableInstances #-}++-- | Provides generic typeclasses for common operations among many types: addition, subtraction, scalar multiplication, magnitude, and zero.+module RSAGL.Math.AbstractVector+    (AbstractVector,+     AbstractZero(..),+     AbstractAdd(..),+     AbstractSubtract(..),+     AbstractScale(..),+     AbstractMagnitude(..),+     abstractScaleTo,+     abstractSum,+     abstractAverage,+     abstractDistance)+    where++import Data.Fixed+import Control.Applicative+import RSAGL.Math.Types++-- | A data type that has an additive identity.+class AbstractZero a where+    zero :: a++-- | A data type that supports addition.+--+-- * @a `add` zero = a@+class AbstractAdd p v | p -> v where+    add :: p -> v -> p++-- | A data type that supports subtraction.+--+-- * @a `sub` a = zero@+class AbstractSubtract p v | p -> v where+    sub :: p -> p -> v++-- | A data type that supports scalar multiplication.+--+-- * @scalarMultiply 0 a = zero@+class AbstractScale v where+    scalarMultiply :: RSdouble -> v -> v++-- | A data type that supports scalar magnitude.+--+-- * @magnitude (scalarMultiply (recip $ magnitude a) a) = 1@+class AbstractMagnitude v where+    magnitude :: v -> RSdouble++-- | A convenience class for many vector types.+class (AbstractZero v,AbstractAdd v v,AbstractSubtract v v,AbstractScale v) => AbstractVector v where++-- Integer++instance AbstractZero Integer where+    zero = 0++instance AbstractAdd Integer Integer where+    add = (+)++instance AbstractSubtract Integer Integer where+    sub = (-)++instance AbstractMagnitude Integer where+    magnitude = abs . fromInteger++-- Float++instance AbstractZero Float where+    zero = 0++instance AbstractAdd Float Float where+    add = (+)++instance AbstractSubtract Float Float where+    sub = (-)++instance AbstractScale Float where+    scalarMultiply d = (f2f d *)++instance AbstractMagnitude Float where+    magnitude = abs . f2f++instance AbstractVector Float++-- Double++instance AbstractAdd Double Double where+    add = (+)++instance AbstractSubtract Double Double where+    sub = (-)++instance AbstractScale Double where+    scalarMultiply d = (f2f d *)++instance AbstractMagnitude Double where+    magnitude = abs . f2f++instance AbstractVector Double++instance AbstractZero Double where+    zero = 0++-- Fixed++instance (HasResolution a) => AbstractZero (Fixed a) where+    zero = 0++instance (HasResolution a) => AbstractAdd (Fixed a) (Fixed a) where+    add = (+)++instance (HasResolution a) => AbstractSubtract (Fixed a) (Fixed a) where+    sub = (-)++instance (HasResolution a) => AbstractScale (Fixed a) where+    scalarMultiply d = (realToFrac d *)++instance (HasResolution a) => AbstractMagnitude (Fixed a) where+    magnitude = abs . realToFrac++instance (HasResolution a) => AbstractVector (Fixed a)++-- Tuples++instance (AbstractZero a,AbstractZero b) => AbstractZero (a,b) where+    zero = (zero,zero)++instance (AbstractAdd a a',AbstractAdd b b') => AbstractAdd (a,b) (a',b') where+    add (a,b) (c,d) = (add a c,add b d)++instance (AbstractSubtract a a',AbstractSubtract b b') => AbstractSubtract (a,b) (a',b') where+    sub (a,b) (c,d) = (sub a c,sub b d)++instance (AbstractScale a,AbstractScale b) => AbstractScale (a,b) where+    scalarMultiply d (a,b) = (scalarMultiply d a,scalarMultiply d b)++instance (AbstractMagnitude a,AbstractMagnitude b) => AbstractMagnitude (a,b) where+    magnitude (a,b) = sqrt $ magnitude a ^ 2 + magnitude b ^ 2++instance (AbstractVector a,AbstractVector b) => AbstractVector (a,b)++-- Functions++instance (AbstractAdd a a') => AbstractAdd ((->) x a) ((->) x a') where+    add a b = \x -> a x `add` b x++instance (AbstractSubtract a a') => AbstractSubtract ((->) x a) ((->) x a') where+    sub a b = \x -> a x `sub` b x++instance (AbstractScale a) => AbstractScale ((->) x a) where+    scalarMultiply d f = scalarMultiply d . f++-- RSfloat++instance AbstractAdd RSfloat RSfloat where+    add = (+)++instance AbstractSubtract RSfloat RSfloat where+    sub = (-)++instance AbstractScale RSfloat where+    scalarMultiply d = (f2f d *)++instance AbstractMagnitude RSfloat where+    magnitude = abs . f2f++instance AbstractVector RSfloat++instance AbstractZero RSfloat where+    zero = 0++-- RSdouble++instance AbstractAdd RSdouble RSdouble where+    add = (+)++instance AbstractSubtract RSdouble RSdouble where+    sub = (-)++instance AbstractScale RSdouble where+    scalarMultiply = (*)++instance AbstractMagnitude RSdouble where+    magnitude = abs . f2f++instance AbstractVector RSdouble++instance AbstractZero RSdouble where+    zero = 0++-- Lists++instance (AbstractZero a) => AbstractZero [a] where+    zero = repeat zero++instance (AbstractAdd a b) => AbstractAdd [a] [b] where+    add = zipWith add++instance (AbstractSubtract a b) => AbstractSubtract [a] [b] where+    sub = zipWith sub++instance (AbstractScale a) => AbstractScale [a] where+    scalarMultiply d = map (scalarMultiply d)++instance (AbstractMagnitude a) => AbstractMagnitude [a] where+    magnitude = sqrt . sum . map ((^2) . magnitude)++instance (AbstractVector a) => AbstractVector [a] where++-- Generic functions.++-- | Force a vector to the specified magnitude.+abstractScaleTo :: (AbstractScale v,AbstractMagnitude v) => RSdouble -> v -> v+abstractScaleTo _ v | magnitude v == 0 = v+abstractScaleTo x v = scalarMultiply (x / magnitude v) v++-- | Sum of a list.+abstractSum :: (AbstractAdd p v,AbstractZero p) => [v] -> p+abstractSum = foldr (flip add) zero++-- | Average of a list.+abstractAverage :: (AbstractAdd p v,AbstractSubtract p v,AbstractVector v,AbstractZero p) => [p] -> p+abstractAverage vs = zero `add` scalarMultiply (recip $ fromInteger total_count) total_sum+    where f y (i,x) = i `seq` x `seq` (i+1,y `add` x)+          (total_count,total_sum) = foldr f (0,zero) $ map (`sub` zero) vs++-- | Distance between two points, based on the 'magnitude' of the difference.+abstractDistance :: (AbstractMagnitude v,AbstractSubtract p v) => p -> p -> RSdouble+abstractDistance x y = magnitude $ x `sub` y+
+ RSAGL/Math/Affine.hs view
@@ -0,0 +1,145 @@+-- | Affine Transformations of Arbitrary Geometric Objects+module RSAGL.Math.Affine+    (AffineTransformable(..),+     scale',+     inverseTransform,+     withTransformation,+     transformAbout,+     translateToFrom,+     rotateToFrom,+     scaleAlong)+    where++import Graphics.Rendering.OpenGL.GL as GL hiding (R)+import RSAGL.Math.Vector+import RSAGL.Math.Matrix+import RSAGL.Math.Angle+import RSAGL.Math.Types+import Data.Maybe++-- | 'AffineTransformable' objects are subject to affine transformations using matrix multiplication.+class AffineTransformable a where+    -- | Apply an affine transformation, defined by a 4x4 matrix.  (This is the only required method.)+    transform :: RSAGL.Math.Matrix.Matrix -> a -> a+    -- | Scale an entity along the @x@ @y@ and @z@ axes.  For example, @scale (Vector3D 2 3 4)@ will make an object twice as wide, three times as tall,+    -- and four times as deep.  It may be helpful to think of the vector as a control point on the vertex of a unit cube.+    scale :: Vector3D -> a -> a+    scale vector = transform $ scaleMatrix vector+    -- | Translate an entity along the specified vector.+    translate :: Vector3D -> a -> a+    translate vector = transform $ translationMatrix vector+    -- | Rotate an entity about the origin, using the specified vector as the axis of rotation.+    -- See also 'transformAbout' to rotate around an arbitrary point.+    rotate :: Vector3D -> Angle -> a -> a+    rotate vector angle = transform $ rotationMatrix vector angle+    -- | Specific rotation around the x-axis.+    rotateX :: Angle -> a -> a+    rotateX = RSAGL.Math.Affine.rotate (Vector3D 1 0 0)+    -- | Specific rotation around the y-axis.+    rotateY :: Angle -> a -> a+    rotateY = RSAGL.Math.Affine.rotate (Vector3D 0 1 0)+    -- | Specific rotation around the z-axis.+    rotateZ :: Angle -> a -> a+    rotateZ = RSAGL.Math.Affine.rotate (Vector3D 0 0 1)++-- | Apply the inverse of an affine transformation, defined by a 4x4 matrix.+{-# INLINE inverseTransform #-}+inverseTransform :: (AffineTransformable a) => RSAGL.Math.Matrix.Matrix -> a -> a+inverseTransform m = transform (matrixInverse m)++-- | Specific scale preserving proportions.+{-# INLINE scale' #-}+scale' :: (AffineTransformable a) => RSdouble -> a -> a+scale' x = RSAGL.Math.Affine.scale (Vector3D x x x)++-- | Apply a function under an affine transformation.  @withTransformation m id@ is an identity if @m@ is invertable.+{-# INLINE withTransformation #-}+withTransformation :: (AffineTransformable a) => RSAGL.Math.Matrix.Matrix -> (a -> a) -> a -> a+withTransformation m f = inverseTransform m . f . transform m++-- | Apply a function treating a particular point as the origin.  For example, combining 'transformAbout' with 'RSAGL.Math.Affine.rotate'+-- performs a rotation about an arbitrary point rather than the origin.+{-# INLINE transformAbout #-}+transformAbout :: (AffineTransformable a) => Point3D -> (a -> a) -> a -> a+transformAbout center f = withTransformation (translateToFrom origin_point_3d center identity_matrix) f++-- | Specific translation along the vector between two points.+-- This ordinary use is to set the second point as the center of a model (typically origin_point_3d)+-- and the first point as the desired position of the model.+{-# INLINE translateToFrom #-}+translateToFrom :: (AffineTransformable a) => Point3D -> Point3D -> a -> a+translateToFrom a b = RSAGL.Math.Affine.translate (vectorToFrom a b)++-- | Specific rotation along the shortest path that brings the second vector in line with the first.+{-# INLINE rotateToFrom #-}+rotateToFrom :: (AffineTransformable a) => Vector3D -> Vector3D -> a -> a+rotateToFrom u v = RSAGL.Math.Affine.rotate c a+    where c = vectorNormalize $ vectorScale (-1) $ fromMaybe (fst $ orthos u) $ aNonZeroVector $ crossProduct u v+          a = angleBetween u v++-- | Specific scale along an arbitary axis.+{-# INLINE scaleAlong #-}+scaleAlong :: (AffineTransformable a) => Vector3D -> RSdouble -> a -> a+scaleAlong v u = withTransformation (rotateToFrom (Vector3D 0 1 0) v identity_matrix) (RSAGL.Math.Affine.scale (Vector3D 1 u 1))++instance AffineTransformable a => AffineTransformable (Maybe a) where+    transform m = fmap (transform m)++instance AffineTransformable a => AffineTransformable [a] where+    transform m = map (transform m)++instance (AffineTransformable a,AffineTransformable b) => AffineTransformable (a,b) where+    transform m (a,b) = (transform m a,transform m b)++instance (AffineTransformable a,AffineTransformable b,AffineTransformable c) => AffineTransformable (a,b,c) where+    transform m (a,b,c) = (transform m a,transform m b,transform m c)++instance AffineTransformable RSAGL.Math.Matrix.Matrix where+    transform mat = matrixMultiply mat++instance AffineTransformable Vector3D where+    transform m (Vector3D x y z) = transformHomogenous x y z 0 Vector3D m+    scale (Vector3D x1 y1 z1) (Vector3D x2 y2 z2) = Vector3D (x1*x2) (y1*y2) (z1*z2)+    translate _ = id+    rotateX a (Vector3D x y z) = Vector3D x (c*y-s*z) (c*z+s*y)+        where s = sine a+              c = cosine a+    rotateY a (Vector3D x y z) = Vector3D (c*x+s*z) y (c*z-s*x)+        where s = sine a+              c = cosine a+    rotateZ a (Vector3D x y z) = Vector3D (c*x-s*y) (c*y+s*x) z+        where s = sine a+              c = cosine a++instance AffineTransformable Point3D where+    transform m (Point3D x y z) = transformHomogenous x y z 1 Point3D m+    scale (Vector3D x1 y1 z1) (Point3D x2 y2 z2) = Point3D (x1*x2) (y1*y2) (z1*z2)+    translate (Vector3D x1 y1 z1) (Point3D x2 y2 z2) = Point3D (x1+x2) (y1+y2) (z1+z2)+    rotateX a (Point3D x y z) = Point3D x (c*y-s*z) (c*z+s*y)+        where s = sine a+              c = cosine a+    rotateY a (Point3D x y z) = Point3D (c*x+s*z) y (c*z-s*x)+        where s = sine a+              c = cosine a+    rotateZ a (Point3D x y z) = Point3D (c*x-s*y) (c*y+s*x) z+        where s = sine a+              c = cosine a++instance AffineTransformable SurfaceVertex3D where+    transform m (SurfaceVertex3D p v) = SurfaceVertex3D (RSAGL.Math.Affine.transform m p) (RSAGL.Math.Affine.transform (matrixTranspose $ matrixInverse m) v)+    translate vector (SurfaceVertex3D p v) = SurfaceVertex3D (RSAGL.Math.Affine.translate vector p) v++-- | The IO monad itself is AffineTransformable.  This is done by wrapping the IO action in an OpenGL transformation.+instance AffineTransformable (IO a) where+    transform mat iofn = preservingMatrix $ do mat' <- newMatrix RowMajor $ map f2f $ concat $ rowMajorForm mat+                                               multMatrix (mat' :: GLmatrix GLdouble)+                                               iofn+    translate (Vector3D x y z) iofn = preservingMatrix $ +        do GL.translate $ Vector3 (f2f x) (f2f y) (f2f z :: GLdouble)+           iofn+    scale (Vector3D x y z) iofn = preservingMatrix $ +        do GL.scale (f2f x) (f2f y) (f2f z :: GLdouble)+           iofn+    rotate (Vector3D x y z) angle iofn = preservingMatrix $ +        do GL.rotate (f2f $ toDegrees_ angle) (Vector3 (f2f x) (f2f y) (f2f z :: GLdouble))+           iofn
+ RSAGL/Math/Angle.hs view
@@ -0,0 +1,193 @@+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}+module RSAGL.Math.Angle+    (Angle,+     BoundAngle(..),+     fromDegrees,+     fromRadians,+     fromRotations,+     fromTimeOfDayHMS,+     fromArcMinutes,+     fromArcSeconds,+     sine,+     arcSine,+     cosine,+     arcCosine,+     tangent,+     arcTangent,+     cartesianToPolar,+     polarToCartesian,+     toRadians,+     toRadians_,+     toDegrees,+     toDegrees_,+     toRotations,+     toRotations_,+     scaleAngle,+     supplementaryAngle,+     zero_angle,+     angularIncrements,+     angleAdd,+     angleSubtract,+     angleNegate,+     absoluteAngle,+     unboundAngle)+    where++import RSAGL.Math.FMod+import RSAGL.Math.AbstractVector+import RSAGL.Math.Types++-- | An angular value.+newtype Angle = Radians RSdouble deriving (Read,Show)++-- | An angular value.  'BoundAngle's are always in the range between -180 and 180 degrees, inclusive.+newtype BoundAngle = BoundAngle Angle deriving (Read,Show)++zero_angle :: Angle+zero_angle = Radians 0++instance Eq Angle where+    (==) a b = case (toRadians_ a,toRadians_ b) of+                   (x,y) | abs x == pi && abs y == pi -> True+                   (x,y) | x == y -> True+                   _ -> False++instance Ord Angle where+    compare x y = case () of+                      _ | x == y -> EQ+                      _ -> compare (toRadians_ x) (toRadians_ y)++instance AbstractZero Angle where+    zero = zero_angle++instance AbstractZero BoundAngle where+    zero = BoundAngle zero_angle++instance AbstractAdd Angle Angle where+    add = angleAdd++instance AbstractAdd BoundAngle Angle where+    add (BoundAngle a) x = BoundAngle $ boundAngle $ a `add` x++instance AbstractSubtract Angle Angle where+    sub = angleSubtract++instance AbstractSubtract BoundAngle Angle where+    sub (BoundAngle a) (BoundAngle b) = boundAngle $ a `sub` b++instance AbstractScale Angle where+    scalarMultiply = scaleAngle++instance AbstractVector Angle++instance AbstractMagnitude Angle where+    magnitude = toRotations_ . absoluteAngle++-- | angularIncrements answers n evenly distributed angles from 0 to 2*pi.+angularIncrements :: Integer -> [Angle]+angularIncrements subdivisions = map (fromRadians . (2*pi*) . (/ fromInteger subdivisions) . fromInteger) [0 .. subdivisions - 1]++-- | There are 2*pi radians in a circle.+fromRadians :: RSdouble -> Angle+fromRadians = Radians++-- | There are 260 degrees in a circle.+fromDegrees :: RSdouble -> Angle+fromDegrees = Radians . ((*) (pi/180))++-- | There is 1 rotation in a circle.+fromRotations :: RSdouble -> Angle+fromRotations = Radians . ((*) (2*pi))++-- | Get an angle based on time of day, hours, minutes, seconds, where noon is considered a zero angle.+fromTimeOfDayHMS :: RSdouble -> RSdouble -> RSdouble -> Angle+fromTimeOfDayHMS h m s = fromRotations (((s/60+m)/60+h)/24)++-- | There are 21600 arc minutes in a circle, 60 arc minutes in a degree.+fromArcMinutes :: RSdouble -> Angle+fromArcMinutes = fromDegrees . (/60)++-- | There are 1296000 arc seconds in a circle, 60 arc seconds in an arc minutes.+fromArcSeconds :: RSdouble -> Angle+fromArcSeconds = fromArcMinutes . (/60)++-- | Answers the angle in the range of -180 to 180, inclusive.+toDegrees :: Angle -> RSdouble+toDegrees x = let x' = toRadians x+                  in x' * 180 / pi+++-- | 'toDegrees_' answers the angle in degrees with no range limitation.+toDegrees_ :: Angle -> RSdouble+toDegrees_ (Radians x) = x * 180 / pi++-- | 'toRadians' answers the angle in the range of -pi .. pi, inclusive.+toRadians :: Angle -> RSdouble+toRadians x = let (Radians x') = boundAngle x+                  in x'+++-- | toRadians answers the angle in radians with no range limitation.+toRadians_ :: Angle -> RSdouble+toRadians_ (Radians x) = x++-- | 'toRotations' answers the angle in the range of -0.5 to 0.5, inclusive.+toRotations :: Angle -> RSdouble+toRotations x= let x' = toRadians x+                   in x' / pi / 2+++-- | 'toRotations' answers the angle in rotations with no range limitation.+toRotations_ :: Angle -> RSdouble+toRotations_ (Radians x) = x / pi / 2++scaleAngle :: RSdouble -> Angle -> Angle+scaleAngle x = Radians . (*x) . toRadians_++supplementaryAngle :: Angle -> Angle+supplementaryAngle (Radians x) = Radians $ pi - x++angleAdd :: Angle -> Angle -> Angle+angleAdd (Radians x) (Radians y) = Radians $ x + y++angleSubtract :: Angle -> Angle -> Angle+angleSubtract (Radians x) (Radians y) = Radians $ x - y++angleNegate :: Angle -> Angle+angleNegate (Radians x) = Radians $ negate x++-- | Absolute value of an angle.+absoluteAngle :: Angle -> Angle+absoluteAngle (Radians x) = Radians $ abs x++sine :: Angle -> RSdouble+sine (Radians x) = sin x++arcSine :: RSdouble -> Angle+arcSine = fromRadians . asin++cosine :: Angle -> RSdouble+cosine (Radians x) = cos x++arcCosine :: RSdouble -> Angle+arcCosine = fromRadians . acos++tangent :: Angle -> RSdouble+tangent (Radians x) = tan x++arcTangent :: RSdouble -> Angle+arcTangent = fromRadians . atan++cartesianToPolar :: (RSdouble,RSdouble) -> (Angle,RSdouble)+cartesianToPolar (u,v) = (fromRadians $ atan2 v u,sqrt $ u*u + v*v)++polarToCartesian :: (Angle,RSdouble) -> (RSdouble,RSdouble)+polarToCartesian (a,d) = (cosine a*d,sine a*d)++-- | 'boundAngle' forces the angle into the range (-pi..pi).+boundAngle :: Angle -> Angle+boundAngle (Radians x) = Radians $ if bounded > pi then bounded - 2*pi else bounded+    where bounded = x `fmod` (2*pi)++unboundAngle :: BoundAngle -> Angle+unboundAngle (BoundAngle a) = a
+ RSAGL/Math/BoundingBox.hs view
@@ -0,0 +1,87 @@+module RSAGL.Math.BoundingBox+    (BoundingBox,+     Bound3D(..),+     boundingCenterRadius,+     minimalDistanceToBoundingBox)+   where++import RSAGL.Math.Vector+import RSAGL.Math.Interpolation+import RSAGL.Math.Affine+import RSAGL.Math.Types++-- | A simple bounding box.  Operations on bounding boxes+-- are designed to be fast, not accurate.  The only+-- guarantee on any bounding box operation is that+-- objects reported to be outside a bounding box,+-- are.+data BoundingBox = BoundingBox {+    bbox_bottom, bbox_top, bbox_left,+    bbox_right, bbox_far, bbox_near :: !RSdouble }+        deriving (Show)++-- | A convenience class for any finite geometry.+-- In particular, it's easy to concatenate the bounding+-- box of multiple geometries by placing them in a list+-- and taking the bounding box of the entire list.+class Bound3D a where+    boundingBox :: a -> BoundingBox++instance Bound3D Point3D where+    boundingBox (Point3D x y z) = BoundingBox {+                                      bbox_bottom = y,+                                      bbox_top = y,+                                      bbox_right = x,+                                      bbox_left = x,+                                      bbox_near = z,+                                      bbox_far = z }++instance Bound3D SurfaceVertex3D where+    boundingBox = boundingBox . sv3d_position++instance (Bound3D a) => Bound3D [a] where+    boundingBox [] = error "instance Bound3D [a], boundingBox: can't construct boundingBox []"+    boundingBox (x:xs) = foldr combineBoundingBoxes (boundingBox x) $ xs++instance Bound3D BoundingBox where+    boundingBox = id++instance AffineTransformable BoundingBox where+    transform m = boundingBox . transform m . boundingBoxToPointCloud++combineBoundingBoxes :: (Bound3D a,Bound3D b) => a -> b -> BoundingBox+combineBoundingBoxes x y =+    BoundingBox {+        bbox_left = min (bbox_left b1) (bbox_left b2),+        bbox_right = max (bbox_right b1) (bbox_right b2),+        bbox_bottom = min (bbox_bottom b1) (bbox_bottom b2),+        bbox_top = max (bbox_top b1) (bbox_top b2),+        bbox_near = min (bbox_near b1) (bbox_near b2),+        bbox_far = max (bbox_far b1) (bbox_far b2) }+    where b1 = boundingBox x+          b2 = boundingBox y++boundingBoxToPointCloud :: BoundingBox -> [Point3D]+boundingBoxToPointCloud bbox =+    [Point3D (bbox_left bbox)  (bbox_bottom bbox) (bbox_near bbox),+     Point3D (bbox_right bbox) (bbox_bottom bbox) (bbox_near bbox),+     Point3D (bbox_left bbox)  (bbox_top bbox)    (bbox_near bbox),+     Point3D (bbox_right bbox) (bbox_top bbox)    (bbox_near bbox),+     Point3D (bbox_left bbox)  (bbox_bottom bbox) (bbox_far bbox),+     Point3D (bbox_right bbox) (bbox_bottom bbox) (bbox_far bbox),+     Point3D (bbox_left bbox)  (bbox_top bbox)    (bbox_far bbox),+     Point3D (bbox_right bbox) (bbox_top bbox)    (bbox_far bbox)]++-- | View of a bounding box in the form of a bounding spehre.+boundingCenterRadius :: BoundingBox -> (Point3D,RSdouble)+boundingCenterRadius bbox = (lerp 0.5 (nlb,frt),distanceBetween nlb frt / 2)+    where nlb = Point3D (bbox_near bbox) (bbox_left bbox) (bbox_bottom bbox)+          frt = Point3D (bbox_far bbox) (bbox_right bbox) (bbox_top bbox)++-- | Estimates distance between a point and the outside surface of a bounding+-- box.  If the value is negative, then the point lies inside the bound+-- region.+minimalDistanceToBoundingBox :: Point3D -> BoundingBox -> RSdouble+minimalDistanceToBoundingBox p bbox = distanceBetween p c - r+    where (c,r) = boundingCenterRadius bbox+
+ RSAGL/Math/Curve.hs view
@@ -0,0 +1,302 @@+{-# LANGUAGE GeneralizedNewtypeDeriving, ExistentialQuantification, Rank2Types #-}+module RSAGL.Math.Curve+    (Curve,+     zipCurve,+     iterateCurve,+     transposeCurve,+     curve,+     Surface,+     surface,+     wrapSurface,+     unwrapSurface,+     transposeSurface,+     zipSurface,+     iterateSurface,+     halfIterateSurface,+     flipTransposeSurface,+     translateCurve,+     scaleCurve,+     clampCurve,+     loopCurve,+     controlCurve,+     transformCurve2,+     uv_identity,+     translateSurface,+     scaleSurface,+     transformSurface,+     transformSurface2,+     surfaceDerivative,+     curveDerivative,+     orientationLoops,+     newellCurve,+     surfaceNormals3D,+     SamplingAlgorithm,+     IntervalSample,+     intervalRange,+     intervalSize,+     intervalSample,+     intervalSingleIntegral,+     linearSamples,+     adaptiveMagnitudeSamples,+     integrateCurve)+    where++import RSAGL.Math.ListUtils+import RSAGL.Math.Vector+import RSAGL.Math.Angle+import RSAGL.Math.Affine+import Data.List+import Data.Maybe+import Control.Parallel.Strategies+import Control.Applicative+import RSAGL.Math.AbstractVector+import Debug.Trace+import RSAGL.Math.Interpolation+import RSAGL.Math.FMod+import RSAGL.Math.Types+import RSAGL.Math.BoundingBox++-- | A parametric function that is aware of it's own sampling interval.  The first parameter is the sampling interval, while the second is the curve input parameter.+type CurveF a = (RSdouble,RSdouble) -> a+-- | A surface is a curve of curves.+type SurfaceF a = CurveF (CurveF a)+-- | A 'Curve' is a parametric function from a one-dimensional space into a space of an arbitrary datatype.  The key feature of a 'Curve' is that it is aware of it's own+-- sampling interval.  Using this information and appropriate arithmetic and scalar multiplication functions provided by RSAGL.AbstractVector, a 'Curve' can be differentiated or integrated.+newtype Curve a = Curve { fromCurve :: CurveF a }++instance Functor Curve where+   fmap g (Curve f) = Curve $ g . f++instance Applicative Curve where+    pure a = Curve $ const a+    f <*> a = zipCurve ($) f a++instance (AffineTransformable a) => AffineTransformable (Curve a) where+    scale v = fmap (scale v)+    translate v = fmap (translate v)+    rotate vector angle = fmap (rotate vector angle)+    transform m = fmap (transform m)+    rotateX angle = fmap (rotateX angle)+    rotateY angle = fmap (rotateY angle)+    rotateZ angle = fmap (rotateZ angle)++instance NFData (Curve a) where+    rnf (Curve f) = seq f ()++-- | Sample a specific point on a curve given the sampling interval (first parameter) and the point on the curve (second parameter).+sampleCurve :: Curve a -> RSdouble -> RSdouble -> a+sampleCurve (Curve f) = curry f++-- | Sample a curve at regular intervals in the range 0..1 inclusive.+iterateCurve :: Integer -> Curve x -> [x]+iterateCurve n c = map f $ zeroToOne n+    where f = sampleCurve c (0.25/fromInteger n)++-- | Combine two curves using an arbitrary function.+zipCurve :: (x -> y -> z) -> Curve x -> Curve y -> Curve z+zipCurve f (Curve x) (Curve y) = Curve $ \hu -> f (x hu) (y hu)++-- | Arbitrarily transform a 'Curve'.+mapCurve :: (CurveF a -> CurveF a) -> Curve a -> Curve a+mapCurve f = Curve . f . fromCurve++-- | Arbitrarily transform a surface 'Curve'.+mapCurve2 :: (SurfaceF a -> SurfaceF a) -> Curve (Curve a) -> Curve (Curve a)+mapCurve2 f = Curve . (Curve .) . f . (fromCurve .) . fromCurve++-- | Translate a curve along the axis of the input parameter.+translateCurve :: RSdouble -> Curve a -> Curve a+translateCurve x = mapCurve $ \f (h,u) -> f (h,u-x)++-- | Scale a curve along the axis of the input parameter.  Factors greater than one have a "zoom in" effect, while+-- factors less than one have a "zoom out" effect.+scaleCurve :: RSdouble -> Curve a -> Curve a+scaleCurve s = mapCurve (\f (h,u) -> f (h/s,u/s))++-- | Clamp lower and upper bounds of a curve along the axis of the input parameter.+clampCurve :: (RSdouble,RSdouble) -> Curve a -> Curve a+clampCurve (a,b) | b < a = clampCurve (b,a)+clampCurve (a,b) = mapCurve $ \f (h,u) -> f (h,min b $ max a u)++-- | Loop a curve onto itself at the specified bounds.+loopCurve :: (RSdouble,RSdouble) -> Curve a -> Curve a+loopCurve (a,b) | b < a = loopCurve (b,a)+loopCurve (a,b) = mapCurve $ \f (h,u) -> f (h,(u-a) `fmod` (b-a) + a)++-- | Transform a curve by manipulating control points.+controlCurve :: (RSdouble,RSdouble) -> (RSdouble,RSdouble) -> Curve a -> Curve a+controlCurve (u,v) (u',v') = translateCurve u' . scaleCurve ((u'-v') / (u-v)) . translateCurve (negate u)++-- | Transpose the inner and outer components of a curve.+transposeCurve :: Curve (Curve a) -> Curve (Curve a)+transposeCurve = mapCurve2 flip++-- | Lift two curve transformations onto each axis of a second order curve.+transformCurve2 :: (forall u. Curve u -> Curve u) -> (forall v. Curve v -> Curve v) -> Curve (Curve a) -> Curve (Curve a)+transformCurve2 fu fv = transposeCurve . fu . transposeCurve . fv++-- | Define a simple curve.+curve :: (RSdouble -> a) -> Curve a+curve = Curve . uncurry . const++-- | A 'Surface' is a based on a 'Curve' with an output of another 'Curve'.+newtype Surface a = Surface (Curve (Curve a)) deriving (NFData,AffineTransformable)++-- | Define a simple surface.+surface :: (RSdouble -> RSdouble -> a) -> Surface a+surface f = Surface $ curve (\u -> curve $ flip f u)++wrapSurface :: Curve (Curve a) -> Surface a+wrapSurface = Surface++unwrapSurface :: Surface a -> Curve (Curve a)+unwrapSurface (Surface s) = s++-- | Transpose the axes of a 'Surface'.+transposeSurface :: Surface a -> Surface a+transposeSurface (Surface s) = Surface $ transposeCurve s++-- | Sample a surface at regularly spaced lattice points in the range 0..1 inclusive.+iterateSurface :: (Integer,Integer) -> Surface a -> [[a]]+iterateSurface (u,v) (Surface s) = map (iterateCurve u) $ iterateCurve v s++-- | Sample the outer 'Curve' of a 'Surface' at regularly spaced intervals.+halfIterateSurface :: Integer -> Surface a -> [Curve a]+halfIterateSurface u = iterateCurve u . unwrapSurface++instance Functor Surface where+    fmap f (Surface x) = Surface $ fmap (fmap f) x++instance Applicative Surface where+    pure a = surface (const $ const a)+    f <*> a = zipSurface ($) f a++-- | Combine two surfaces using an arbitrary function.+zipSurface :: (x -> y -> z) -> Surface x -> Surface y -> Surface z+zipSurface f (Surface x) (Surface y) = Surface $ zipCurve (zipCurve f) x y++-- | Lift a transformation on a second order 'Curve' onto a Surface.+transformSurface :: (Curve (Curve a) -> Curve (Curve a)) -> Surface a -> Surface a+transformSurface f = Surface . f . unwrapSurface++-- | Lift two curve transformations onto each axis of a Surface.+transformSurface2 :: (forall u. Curve u -> Curve u) -> (forall v. Curve v -> Curve v) -> Surface a -> Surface a+transformSurface2 fu fv = transformSurface (transformCurve2 fu fv)++-- | Translate a surface over each of its input parameter axes, as translateCurve.+translateSurface :: (RSdouble,RSdouble) -> Surface a -> Surface a+translateSurface (u,v) = transformSurface2 (translateCurve u) (translateCurve v)++-- | Scale a surface along each of its input parameter axes, as scaleCurve.+scaleSurface :: (RSdouble,RSdouble) -> Surface a -> Surface a+scaleSurface (u,v) = transformSurface2 (scaleCurve u) (scaleCurve v)++-- | Transpose a surface while flipping the inner curve, so that that orientable surfaces retain their original orientation.+flipTransposeSurface :: Surface a -> Surface a+flipTransposeSurface = transformSurface (mapCurve2 $ \f (hu,u) (hv,v) -> f (hu,u) (hv,1-v)) . transposeSurface++-- | An identity 'Surface'.+uv_identity :: Surface (RSdouble,RSdouble)+uv_identity = surface (curry id)++-- | Take the derivative of a 'Curve'.+curveDerivative :: (AbstractSubtract p v,AbstractScale v) => Curve p -> Curve v+curveDerivative (Curve f) = Curve $ \(h,u) -> scalarMultiply (recip $ 2 * h) $ f (h/2,u+h) `sub` f (h/2,u-h)++-- | Take the piecewise derivative of a 'Surface' along the inner and outer curves.+surfaceDerivative :: (AbstractSubtract p v,AbstractScale v) => Surface p -> Surface (v,v)+surfaceDerivative s = zipSurface (,) (curvewiseDerivative s) (transposeSurface $ curvewiseDerivative $ transposeSurface s)+    where curvewiseDerivative (Surface t) = Surface $ fmap curveDerivative t++-- | Determine the orientation of a 'Surface' by passing very small circles centered on each sampled point as the parametric input.+--+-- A gotchya with this operation is that as much as 3/4ths of the orientation loop may lie outside of the 0..1 range that is normally+-- sampled.  Depending on how the surface is constructed, this may produce unexpected results.  The solution is to clamp the+-- the problematic parametric inputs at 0 and 1 using 'clampSurface'.+--+-- As a rule, do clamp longitudinal axes that come to a singularity at each end.+-- Do not clamp latitudinal axes that are connected at each end.+--+orientationLoops :: Surface p -> Surface (Curve p)+orientationLoops (Surface s) = Surface $ Curve $ \(uh,u) -> Curve $ \(vh,v) ->+                                     curve $ \t -> f (uh/2,u + uh*(sine $ fromRotations t)) +				                     (vh/2,v + vh*(cosine $ fromRotations t))+   where f = fromCurve . fromCurve s++-- | Try to determine the normal vector to a curve.+newellCurve :: Curve Point3D -> Maybe Vector3D+newellCurve c = newell $ iterateCurve 16 c++-- | Try to determine the normal vectors of a surface using orientation loops.  This is usually slower but more successful than 'surfaceNormals3DByPartialDerivatives'.+-- This generate a warning message if it cannot determine the normal vector at a sampled point.+-- See also 'orientationLoops'.+surfaceNormals3DByOrientationLoops :: Surface Point3D -> Surface SurfaceVertex3D+surfaceNormals3DByOrientationLoops s = SurfaceVertex3D <$> s <*> ((\c -> errmsg c (newellCurve c)) <$> orientationLoops s)+    where errmsg c = fromMaybe (trace ("surfaceNormals3DByOrientationLoops: zero normal gave up: " ++ show (iterateCurve 16 c)) (Vector3D 0 0 0))++-- | Try to determine the normal vectors of a surface using partial derivatives.+surfaceNormals3DByPartialDerivatives :: Surface Point3D -> Surface (Maybe Vector3D)+surfaceNormals3DByPartialDerivatives s = safeCrossProduct <$> surfaceDerivative s+    where x = snd $ boundingCenterRadius $ boundingBox $ concat $ iterateSurface (8,8) s+          safeCrossProduct (u_,v_) =+              do u <- aLargeVector (x/100) u_+	         v <- aLargeVector (x/100) v_+		 return $ vectorNormalize $ crossProduct u v++-- | Try to determine the normal vectors of a surface using multiple techniques.+surfaceNormals3D :: Surface Point3D -> Surface SurfaceVertex3D+surfaceNormals3D s = (\p by_pd by_newell -> case by_pd of+                           Just v -> SurfaceVertex3D p v+			   Nothing -> by_newell) <$>+                     s <*> (surfaceNormals3DByPartialDerivatives s) <*> (surfaceNormals3DByOrientationLoops s)++-- | An interval of a curve, including the curve, lower and upper bounds of the interval, and an instantaneous sample value for that interval.+data IntervalSample a = IntervalSample (Curve a) RSdouble RSdouble a++intervalSample :: Curve a -> RSdouble -> RSdouble -> IntervalSample a+intervalSample c l h = IntervalSample c l h $ sampleCurve c ((abs $ l - h) / 2) ((l+h) / 2) ++-- | Lower and upper bounds of an 'IntervalSample'.+intervalRange :: IntervalSample a -> (RSdouble,RSdouble)+intervalRange (IntervalSample _ l h _) = (l,h)++-- | Size of the range of an 'IntervalSample'.+intervalSize :: IntervalSample a -> RSdouble+intervalSize (IntervalSample _ l h _) = abs $ h - l++-- | Instantaneous sample value of an 'Interval'.+intervalValue :: IntervalSample a -> a+intervalValue (IntervalSample _ _ _ a) = a++-- | Integral of the sample value over the range of the 'IntervalSample'.+intervalSingleIntegral :: (AbstractScale a) => IntervalSample a -> a+intervalSingleIntegral x = scalarMultiply (intervalSize x) $ intervalValue x++-- | Split an interval into three equal parts.+splitInterval :: IntervalSample a -> [IntervalSample a]+splitInterval (IntervalSample c l h a) = [intervalSample c l l',IntervalSample c l' h' a,intervalSample c h' h]+    where l' = lerp (1/3) (l,h)+          h' = lerp (2/3) (l,h)++type SamplingAlgorithm a = Curve a -> [IntervalSample a]++-- | Definite integral of a curve.+integrateCurve :: (AbstractAdd p v,AbstractScale v,AbstractZero p) => SamplingAlgorithm v -> Curve v -> p -> p+integrateCurve samplingAlgorithm c initial_value = foldl' add initial_value $ map intervalSingleIntegral $ samplingAlgorithm c++-- | Sampling algorithm that takes a fixed count of samples.+linearSamples :: Integer -> SamplingAlgorithm a+linearSamples n c = map (\(l,h) -> intervalSample c l h) $ doubles $ zeroToOne (n+1)++-- | Sampling algorithm that takes increasing numbers of samples over intervals where the magnitude of the sample is large.+adaptiveMagnitudeSamples :: (AbstractMagnitude a) => Integer -> SamplingAlgorithm a+adaptiveMagnitudeSamples n c = resampleLoop (\xs -> if genericLength xs > n then Nothing else Just $ newSamples xs) $ linearSamples (max 1 $ n `div` 10) c+    where newSamples xs = let a = abstractAverage $ map intervalMagnitude xs+                              in flip concatMap xs $ \x -> if intervalMagnitude x >= a then splitInterval x else [x]+          intervalMagnitude :: (AbstractMagnitude a) => IntervalSample a -> RSdouble+          intervalMagnitude (IntervalSample _ l h a) = magnitude a * (abs $ h-l)++-- | Loop to keep generating samples until finished.+resampleLoop :: (b -> Maybe b) -> b -> b+resampleLoop nextPass initial_value = f $ initial_value+    where f x = maybe x f $ nextPass x
+ RSAGL/Math/CurveExtras.lhs view
@@ -0,0 +1,110 @@+\section{Specific and Interpolated Curves}++\begin{code}+module RSAGL.Math.CurveExtras+    (sphericalCoordinates,+     cylindricalCoordinates,+     toroidalCoordinates,+     circularCoordinates,+     polarCoordinates,+     transformUnitSquareToUnitCircle,+     transformUnitCubeToUnitSphere,+     circleXY,+     regularPolygon,+     linearInterpolation,+     loopedLinearInterpolation,+     smoothCurve,+     loopCurve)+    where++import RSAGL.Math.Curve+import RSAGL.Math.Interpolation+import RSAGL.Math.Vector+import RSAGL.Math.Angle+import RSAGL.Math.AbstractVector+import RSAGL.Math.Affine+import RSAGL.Math.ListUtils+import Control.Arrow+import RSAGL.Math.Types+\end{code}++\subsection{Alternate Coordinate Systems for Models}++\begin{code}+sphericalCoordinates :: ((Angle,Angle) -> a) -> Surface a+sphericalCoordinates f = transformSurface2 id (clampCurve (0,1)) $ surface $ curry (f . (\(u,v) -> (fromRadians $ u*2*pi,fromRadians $ ((pi/2) - v*pi))))++cylindricalCoordinates :: ((Angle,RSdouble) -> a) -> Surface a+cylindricalCoordinates f = transformSurface2 id (clampCurve (0,1)) $ surface $ curry (f . (\(u,v) -> (fromRadians $ u*2*pi,v)))++toroidalCoordinates :: ((Angle,Angle) -> a) -> Surface a+toroidalCoordinates f = surface $ curry (f . (\(u,v) -> (fromRadians $ u*2*pi,fromRadians $ negate $ v*2*pi)))++circularCoordinates :: ((RSdouble,RSdouble) -> a) -> Surface a+circularCoordinates f = surface $ curry $ (f . second negate . transformUnitSquareToUnitCircle)++polarCoordinates :: ((Angle,RSdouble) -> a) -> Surface a+polarCoordinates f = circularCoordinates (f . cartesianToPolar)+\end{code}++\subsection{Transformations Between Unit Volumes}++\begin{code}+transformUnitSquareToUnitCircle :: (RSdouble,RSdouble) -> (RSdouble,RSdouble)+transformUnitSquareToUnitCircle (u,v) = (x,z)+    where (Point3D x _ z) = transformUnitCubeToUnitSphere (Point3D u 0.5 v)++transformUnitCubeToUnitSphere :: Point3D -> Point3D+transformUnitCubeToUnitSphere p =+    let p_centered@(Point3D x y z) = scale' 2.0 $ translate (Vector3D (-0.5) (-0.5) (-0.5)) p+        p_projected = scale' (minimum [recip $ abs x,recip $ abs y,recip $ abs z]) p_centered+        k = recip $ distanceBetween origin_point_3d p_projected+        w = maximum $ [abs x, abs y, abs z] -- 'w' could be 1, but this gives a smoother tesselation+        in if p_centered == origin_point_3d then origin_point_3d else lerp w (p_centered,scale' k p_centered)+\end{code}++\subsection{Circles}++\begin{code}+circleXY :: Curve Point3D+circleXY = curve $ \u_ -> let u = fromRotations u_ in Point3D (cosine u) (sine u) 0+\end{code}++\subsection{Regular Polygons}++A regular polygon, centered at the origin, in the XY plane.++\begin{code}+regularPolygon :: Integer -> Curve Point3D+regularPolygon n = loopedLinearInterpolation $ map (flip rotateZ (Point3D 0 1 0) . fromRotations) $ zeroToOne n+\end{code}++\subsection{Piecewise Length Normalization}++\begin{code}+normalizePolyline :: (AbstractSubtract p v,AbstractMagnitude v) => [p] -> [(RSdouble,p)]+normalizePolyline pts = zip (map (/ total_dist) accumulated_dists) pts+    where dists = map (uncurry abstractDistance) $ doubles pts+          total_dist = last accumulated_dists+          accumulated_dists = scanl (+) 0 dists+\end{code}++\subsection{Interpolated Curves}++\begin{code}+linearInterpolation :: (AbstractSubtract p v,AbstractAdd p v,AbstractMagnitude v,AbstractScale v) => [p] -> Curve p+linearInterpolation = curve . lerpMap . normalizePolyline++loopedLinearInterpolation :: (AbstractSubtract p v,AbstractAdd p v,AbstractMagnitude v,AbstractScale v) => [p] -> Curve p+loopedLinearInterpolation = loopCurve (0,1) . linearInterpolation . (\a -> last a:a)+\end{code}++\subsection{Smoothing Curves}++\texttt{smoothCurve i h} takes i samples of an h-long piece of a 'Curve' at each point to smooth it.  This is not an interpolation function and will tend to shrink shapes+toward their center of gravity.++\begin{code}+smoothCurve :: (AbstractAdd p v,AbstractSubtract p v,AbstractVector v,AbstractZero p) => Integer -> RSdouble -> Curve p -> Curve p+smoothCurve i h c = curve $ \u -> abstractAverage $ iterateCurve i $ controlCurve (u-h/2,u+h/2) (0,1) c+\end{code}
+ RSAGL/Math/FMod.hs view
@@ -0,0 +1,24 @@+module RSAGL.Math.FMod+    (FMod(..))+    where++import RSAGL.Math.Types++class FMod f where+    fmod :: f -> f -> f++instance FMod Float where+    fmod n d = n - fromInteger f * d+        where f = floor $ n / d++instance FMod Double where+    fmod n d = n - fromInteger f * d+        where f = floor $ n / d++instance FMod RSfloat where+    fmod n d = n - fromInteger f * d+        where f = floor $ n / d++instance FMod RSdouble where+    fmod n d = n - fromInteger f * d+        where f = floor $ n / d
+ RSAGL/Math/Interpolation.lhs view
@@ -0,0 +1,109 @@+\section{RSAGL.Interpolation}++\begin{code}+module RSAGL.Math.Interpolation+    (lerp,+     lerpClamped,+     lerpBetween,+     lerpBetweenMutated,+     lerpBetweenClamped,+     lerpBetweenClampedMutated,+     lerp_mutator_continuous_1st,+     lerpMap)+    where++import RSAGL.Math.AbstractVector+import Data.Map as Map+import Data.Maybe+import RSAGL.Math.Types+\end{code}++\subsection{The Lerpable typeclass}++Implements linear interpolation.++\begin{code}+{-# INLINE lerp #-}+lerp :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r) => r -> (p,p) -> p+lerp u (a,b) = a `add` scalarMultiply (f2f u) (b `sub` a)+\end{code}++\subsection{Non-linear interpolations}++The mutated versions of the lerp functions takes an arbitrary ``mutator function'' to define non-linear interpolation curves.+The ``between'' versions of the lerp functions allow the u-value to lie between any two numbers, as opposed to between 0 and 1.+The ``clamped'' versions of the lerp functions clamp the u-value to lie between its boundaries.  Otherwise, with non-clamped+interpolations, the u-value may lie outside of its boundaries.++\begin{code}+{-# INLINE lerpClamped #-}++lerpClamped :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r) => r -> (p,p) -> p+lerpClamped u = lerpBetweenClamped (0,u,1)++{-# INLINE lerpBetween #-}+lerpBetween :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r) => (r,r,r) -> (p,p) -> p+lerpBetween = lerpBetweenMutated id++{-# INLINE lerpBetweenMutated #-}+lerpBetweenMutated :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r) => (r -> r) -> (r,r,r) -> (p,p) -> p+lerpBetweenMutated _ (l,_,r) | l == r = lerp 0.5+lerpBetweenMutated mutator (l,u,r) = lerp $ mutator $ (u-l) / (r-l)++{-# INLINE lerpBetweenClamped #-}+lerpBetweenClamped :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r,Ord r) => (r,r,r) -> (p,p) -> p+lerpBetweenClamped = lerpBetweenClampedMutated id++{-# INLINE lerpBetweenClampedMutated #-}+lerpBetweenClampedMutated :: (AbstractScale v,AbstractSubtract p v,AbstractAdd p v,RealFloat r,Ord r) => (r -> r) -> (r,r,r) -> (p,p) -> p+lerpBetweenClampedMutated mutator (l,u,r) = lerpBetweenMutated (lerp_mutator_clamp . mutator) (l,u,r)+\end{code}++\subsection{Lerp mutators}++\texttt{lerp\_mutator\_clamp} implements clamping between 0 and 1.++\begin{code}+{-# INLINE lerp_mutator_clamp #-}++lerp_mutator_clamp :: (RealFloat r) => r -> r+lerp_mutator_clamp = min 1 . max 0+\end{code}++\texttt{lerp\_mutator\_continuous\_1st} implements clamping between 0 and 1, but such that the 1st derivative of the result is continuous.++\begin{code}+{-# INLINE lerp_mutator_continuous_1st #-}++lerp_mutator_continuous_1st :: (RealFloat r) => r -> r+lerp_mutator_continuous_1st x | x < 0 = 0+lerp_mutator_continuous_1st x | x > 1 = 1+lerp_mutator_continuous_1st x | x <= 0.5 = 2*x^2+lerp_mutator_continuous_1st x = 4*x - 2*x^2 - 1+\end{code}++\subsection{lerpMap}++Given many entities, lerp between the two entities closest to the given point+on either side.  For example, if we wanted to lerp between colors on a rainbow,+we might use the map [(0,red),(1,orange),(2,yellow),(3,green),(4,blue),(5,indigo),(6,violet)].+lerpMap 3.5 would result in a blue-green color.++\begin{code}+{-# INLINE lerpMap #-}++lerpMap :: (RealFloat r,AbstractScale v,AbstractSubtract p v,AbstractAdd p v) => [(r,p)] -> r -> p+lerpMap pts = lerpMapSorted $ Map.fromList pts++{-# INLINE lerpMapSorted #-}+lerpMapSorted :: (RealFloat r,AbstractScale v,AbstractSubtract p v,AbstractAdd p v) => Map r p -> r -> p+lerpMapSorted m _ | Map.null m = error "lerpMapSorted: empty map"+lerpMapSorted m u = case () of+        () | isJust exactly -> fromJust exactly+        () | Map.null lowers -> higher_a+        () | Map.null highers -> lower_a+        () -> lerpBetween (lower_r,u,higher_r) (lower_a,higher_a)+    where (lowers, exactly, highers) = splitLookup u m+          (lower_r,lower_a) = findMax lowers+          (higher_r,higher_a) = findMin highers+\end{code}
+ RSAGL/Math/ListUtils.hs view
@@ -0,0 +1,53 @@+module RSAGL.Math.ListUtils+    (doubles,+     loopedDoubles,+     consecutives,+     loopedConsecutives,+     zeroToOne)+    where++import RSAGL.Math.Types+import Debug.Trace++-- | Transforms a list to a list of adjacent elements.+--+-- @doubles [1,2,3,4,5] = [(1,2),(2,3),(3,4),(4,5)]@+doubles :: [a] -> [(a,a)]+doubles [] = []+doubles [_] = []+doubles (x:y:zs) = (x,y) : doubles (y:zs)++-- | loopedDoubles transforms a list to a list of adjacent elements, looping+-- back to the beginning of the list.+--+-- @loopedRSdoubles [1,2,3,4,5] = [(1,2),(2,3),(3,4),(4,5),(5,1)]@+loopedDoubles :: [a] -> [(a,a)]+loopedDoubles as = loopedDoubles_ (head as) as+    where loopedDoubles_ _ [] = []+          loopedDoubles_ a [x] = [(x,a)]+          loopedDoubles_ a (x:y:zs) = (x,y) : loopedDoubles_ a (y:zs)++-- | Answers a list containing every sequence of n consecutive+-- elements in the parameter.+--+-- @consecutives 3 [1,2,3,4] = [[1,2,3],[2,3,4]]@+consecutives :: Int -> [a] -> [[a]]+consecutives n xs = let taken = take n xs+                        in if (length taken == n)+                           then (taken : (consecutives n $ tail xs))+                           else []++-- | Answers a list containing every sequence of n consecutive+-- elements in the parameter, looping back to the beginning of the list.+--+-- @consecutives 3 [1,2,3,4] = [[1,2,3],[2,3,4],[3,4,1],[4,1,2]]@+loopedConsecutives :: Int -> [a] -> [[a]]+loopedConsecutives n xs = consecutives n $ take (n + length xs - 1) $ cycle xs++-- | Creates a list of numbers from 0.0 to 1.0, using n steps.+-- This can't be done with the enum-from-to method, due to roundoff errors.+zeroToOne :: Integer -> [RSdouble]+zeroToOne n | n > 100000 = trace ("Warning: zeroToOne was asked for " ++ show n ++ " subdivisions, which seems high.  Using 100,000 instead.") zeroToOne 100000+zeroToOne n = map (*x) [0..(fromInteger $ n-1)]+    where x = recip (fromInteger $ n - 1)+
+ RSAGL/Math/Matrix.lhs view
@@ -0,0 +1,190 @@+\section{RSAGL.Matrix}++\begin{code}++module RSAGL.Math.Matrix+    (Matrix,+     matrix,+     transformHomogenous,+     rowMajorForm,+     colMajorForm,+     rowAt,+     matrixAt,+     identity_matrix,+     translationMatrix,+     rotationMatrix,+     scaleMatrix,+     xyzMatrix,+     matrixAdd,+     matrixMultiply,+     matrixTranspose,+     matrixInverse,+     determinant,+     matrixInversePrim,+     matrixTransposePrim,+     matrixInverseTransposePrim,+     determinantPrim)+    where++import Data.List as List+import RSAGL.Math.Angle+import RSAGL.Math.Vector+import Data.Vec as Vec+import RSAGL.Math.Types+\end{code}++A 4-by-4 matrix with cached inverse, transpose, inverse transpose, and determinant.++\begin{code}+data Matrix = Matrix { matrix_data :: !(Mat44 RSdouble),+                       matrix_inverse :: Matrix,+                       matrix_transpose :: Matrix,+                       matrix_determinant :: RSdouble }++instance Eq Matrix where+    x == y = matrix_data x == matrix_data y++instance Show Matrix where+    show m = show $ rowMajorForm m++rowMajorForm :: Matrix -> [[RSdouble]]+rowMajorForm = matToLists . matrix_data++colMajorForm :: Matrix -> [[RSdouble]]+colMajorForm = List.transpose . rowMajorForm+\end{code}++rowAt answers the nth row of a matrix.++\begin{code}+rowAt :: Matrix -> Int -> [RSdouble]+rowAt m n = (rowMajorForm m) !! n+\end{code}++matrixAt answers the (i'th,j'th) element of a matrix.++\begin{code}+matrixAt :: Matrix -> (Int,Int) -> RSdouble+matrixAt m (i,j) = rowAt m i !! j+\end{code}++\subsection{Constructing matrices}++matrix constructs a matrix from row major list form.  (Such a list form can be formatted correctly in +monospaced font and haskell syntax, so that it looks like a matrix as it would be normally written.)++\begin{code}+matrix :: [[RSdouble]] -> Matrix+matrix = uncheckedMatrix . matFromLists++-- | Generate a column matrix of length 4, perform an affine transformation on it, and produce the resulting value.+{-# INLINE transformHomogenous #-}+transformHomogenous :: RSdouble -> RSdouble -> RSdouble -> RSdouble -> (RSdouble -> RSdouble -> RSdouble -> a) -> Matrix -> a+transformHomogenous x y z w f m = f x' y' z'+    where (x':.y':.z':._) = multmv (matrix_data m) (x:.y:.z:.w:.())++uncheckedMatrix :: Mat44 RSdouble -> Matrix+uncheckedMatrix dats = m+    where m_inverse = (matrixInversePrim m) { matrix_inverse = m, matrix_transpose = m_inverse_transpose, matrix_determinant = recip m_det }+          m_transpose = (matrixTransposePrim m) { matrix_inverse = m_inverse_transpose, matrix_transpose = m, matrix_determinant = m_det }+          m_inverse_transpose = (matrixTransposePrim m_inverse) { matrix_inverse = m_transpose, matrix_transpose = m_inverse, matrix_determinant = recip m_det }+          m_det = determinantPrim m+          m = Matrix { matrix_data=dats,+                       matrix_inverse = m_inverse,+                       matrix_transpose = m_transpose,+                       matrix_determinant = m_det }+\end{code}++identityMatrix constructs the n by n identity matrix++\begin{code}+identity_matrix :: Matrix+identity_matrix = Matrix {+    matrix_data = identity,+    matrix_inverse = identity_matrix,+    matrix_transpose = identity_matrix,+    matrix_determinant = 1 }+\end{code}++\begin{code}+translationMatrix :: Vector3D -> Matrix+translationMatrix (Vector3D x y z) = uncheckedMatrix $ (1:.0:.0:.x:.()) :.+                                                       (0:.1:.0:.y:.()) :.+                                                       (0:.0:.1:.z:.()) :.+                                                       (0:.0:.0:.1:.()) :. ()+\end{code}++\begin{code}+rotationMatrix :: Vector3D -> Angle -> Matrix+rotationMatrix vector angle = let s = sine angle+				  c = cosine angle+				  c' = 1 - c+				  (Vector3D x y z) = vectorNormalize vector+				  in uncheckedMatrix $ ((c+c'*x*x)   :. (c'*y*x - s*z) :. (c'*z*x + s*y) :. 0 :. ()) :.+					               ((c'*x*y+s*z) :. (c+c'*y*y)     :. (c'*z*y - s*x) :. 0 :. ()) :.+					               ((c'*x*z-s*y) :. (c'*y*z+s*x)   :. (c+c'*z*z)     :. 0 :. ()) :.+					               (0            :. 0              :. 0              :. 1 :. ()) :. ()+\end{code}++\begin{code}+scaleMatrix :: Vector3D -> Matrix+scaleMatrix (Vector3D x y z) = uncheckedMatrix $ (x :. 0 :. 0 :. 0) :.+				                 (0 :. y :. 0 :. 0) :.+				                 (0 :. 0 :. z :. 0) :.+				                 (0 :. 0 :. 0 :. 1) :. ()+\end{code}++\texttt{xyzMatrix} constructs the matrix in which the x y and z axis are transformed to point in the direction of the specified+vectors.++\begin{code}+xyzMatrix :: Vector3D -> Vector3D -> Vector3D -> Matrix+xyzMatrix (Vector3D x1 y1 z1) (Vector3D x2 y2 z2) (Vector3D x3 y3 z3) =+    uncheckedMatrix $ (x1 :. x2 :. x3 :. 0) :.+                      (y1 :. y2 :. y3 :. 0) :.+                      (z1 :. z2 :. z3 :. 0) :.+                      (0  :. 0  :. 0  :. 1.0) :. ()+\end{code}++\subsection{Matrix arithmetic}++\begin{code}+matrixAdd :: Matrix -> Matrix -> Matrix+matrixAdd m n = uncheckedMatrix $ Vec.zipWith (+) (matrix_data m) (matrix_data n)+\end{code}++\label{howtoUseMatrixMultiplyImpl}++\begin{code}+matrixMultiply :: Matrix -> Matrix -> Matrix+matrixMultiply m n = uncheckedMatrix $ multmm (matrix_data m) (matrix_data n)++matrixTranspose :: Matrix -> Matrix+matrixTranspose = matrix_transpose+\end{code}++\subsection{Inverse and determinant of a matrix}++\begin{code}+matrixInverse :: Matrix -> Matrix+matrixInverse = matrix_inverse+\end{code}++\begin{code}+determinant :: Matrix -> RSdouble+determinant = matrix_determinant+\end{code}++\begin{code}+matrixInverseTransposePrim :: Matrix -> Matrix+matrixInverseTransposePrim = uncheckedMatrix . Vec.transpose . fst . invertAndDet . matrix_data++matrixInversePrim :: Matrix -> Matrix+matrixInversePrim = uncheckedMatrix . fst . invertAndDet . matrix_data++determinantPrim :: Matrix -> RSdouble+determinantPrim = det . matrix_data++matrixTransposePrim :: Matrix -> Matrix+matrixTransposePrim = uncheckedMatrix . Vec.transpose . matrix_data+\end{code}
+ RSAGL/Math/Orthogonal.hs view
@@ -0,0 +1,116 @@++-- | It's useful to work with the set of coordinate systems restricted to those+-- that use orthogonal unit-scaled axes, that is, that are subject only to+-- rotation and translation.  This is because these coordinate systems are the+-- describe rigid objects.+module RSAGL.Math.Orthogonal+    (up,down,left,right,forward,backward,+     orthogonalFrame,+     modelLookAt,+     FUR)+    where++import RSAGL.Math.Affine+import RSAGL.Math.Vector+import RSAGL.Math.Matrix++-- | FUR stands for Forward Up Right.  It's used to specify arbitrary+-- orthogonal coordinate systems given any combination of forward up and right+-- vectors.  It also accepts down, left, and backward vectors.++-- | This modules uses a left-handed coordinate system.  Right is positive X,+-- up is positive Y, forward is positive Z.++-- When specifying FUR coordinate systems, the first vector is fixed+-- while the second vector will be adjusted as little as possible to guarantee+-- that it is orthogonal to the first.  The third vector never needs to be+-- specified, it can be deduced.+data FURAxis = ForwardAxis | UpAxis | RightAxis | DownAxis | LeftAxis | BackwardAxis++data FUR a = FUR FURAxis a++instance Functor FUR where+    fmap f (FUR a x) = FUR a $ f x++-- | A reference to the +Y axis.+up :: a -> FUR a+up = FUR UpAxis++-- | A reference to the -Y axis.+down :: a -> FUR a+down = FUR DownAxis++-- | A reference to the -X axis.+left :: a -> FUR a+left = FUR LeftAxis++-- | A reference to the +X axis.+right :: a -> FUR a+right = FUR RightAxis++-- | A reference to the +Z axis.+forward :: a -> FUR a+forward = FUR ForwardAxis++-- | A reference to the -Z axis.+backward :: a -> FUR a+backward = FUR BackwardAxis++-- | Combine two axial references to describe a rigid affine transformation.+-- Accepts any combination of non-coaxial references.+-- In the affine transformation, the old axes will be mapped onto the specified+-- freeform axes.+--+-- The first parameter is absolute, meaning that the source axis will always map+-- perfectly onto the destination axis.  The second parameter will be obeyed+-- on a "best effort" basis.+--+orthogonalFrame :: (AffineTransformable a) => FUR Vector3D -> FUR Vector3D -> a -> a+orthogonalFrame (FUR ForwardAxis f) (FUR RightAxis r) =+    let (r',u') = fixOrtho2 f r in transform (xyzMatrix r' u' (vectorNormalize f))+orthogonalFrame (FUR UpAxis u) (FUR ForwardAxis f) =+    let (f',r') = fixOrtho2 u f in transform (xyzMatrix r' (vectorNormalize u) f')+orthogonalFrame (FUR RightAxis r) (FUR UpAxis u) =+    let (u',f') = fixOrtho2 r u in transform (xyzMatrix (vectorNormalize r) u' f')+orthogonalFrame (FUR RightAxis r) (FUR ForwardAxis f) =+    let (f',u') = fixOrtho2Left r f in transform (xyzMatrix (vectorNormalize r) u' f')+orthogonalFrame (FUR ForwardAxis f) (FUR UpAxis u) =+    let (u',r') = fixOrtho2Left f u in transform (xyzMatrix r' u' (vectorNormalize f))+orthogonalFrame (FUR UpAxis u) (FUR RightAxis r) =+    let (r',f') = fixOrtho2Left u r in transform (xyzMatrix r' (vectorNormalize u) f')+orthogonalFrame (FUR ForwardAxis _) (FUR ForwardAxis _) =+    error "orthogonalFrame: two forward vectors"+orthogonalFrame (FUR UpAxis _) (FUR UpAxis _) =+    error "orthogonalFrame: two up vectors"+orthogonalFrame (FUR RightAxis _) (FUR RightAxis _) =+    error "orthogonalFrame: two right vectors"+orthogonalFrame x y = orthogonalFrame (furCorrect x) (furCorrect y)++-- | Function to transform down, left, or backward axes into+-- up, right, or forward axes.+furCorrect :: FUR Vector3D -> FUR Vector3D+furCorrect (FUR ForwardAxis f) = FUR ForwardAxis f+furCorrect (FUR UpAxis u) = FUR UpAxis u+furCorrect (FUR RightAxis r) = FUR RightAxis r+furCorrect (FUR DownAxis d) = FUR UpAxis $ vectorScale (-1) d+furCorrect (FUR LeftAxis l) = FUR RightAxis $ vectorScale (-1) l+furCorrect (FUR BackwardAxis b) = FUR ForwardAxis $ vectorScale (-1) b++-- | Translates and rotates a model to aim at a given position or in a+-- given direction from a given vantage point.  This is analogous+-- to camera look-at functions, and could be used, for example, to+-- cause a model of an eyeball to track a particular target.+-- The first parameter is the position of the model.  Typically the second+-- parameter will be the position of the target, and the third parameter will+-- @(up \$ Vector3D 0 1 0)@.+modelLookAt :: (AffineTransformable a) =>+    Point3D ->+    FUR (Either Point3D Vector3D) ->+    FUR (Either Point3D Vector3D) ->+    a -> a+modelLookAt pos primaryish secondaryish =+          RSAGL.Math.Affine.translate (vectorToFrom pos origin_point_3d) .+              orthogonalFrame primary secondary+    where primary = fmap (either (`vectorToFrom` pos) id) primaryish+          secondary = fmap (either (`vectorToFrom` pos) id) secondaryish+
+ RSAGL/Math/Ray.hs view
@@ -0,0 +1,47 @@+module RSAGL.Math.Ray+    (Ray3D(..),+     projectRay,+     distanceAlong,+     angleFrom,+     normalizeRay)+    where++import RSAGL.Math.Vector+import RSAGL.Math.Angle+import RSAGL.Math.Affine+import RSAGL.Math.Types++-- | Rays with endpoints and vectors.+--+-- Although a ray is isomorphic to a 'SurfaceVertex3D', it does not have the+-- same behavior.+--+data Ray3D = Ray3D { ray_endpoint :: Point3D,+                     ray_vector :: Vector3D }+    deriving (Read,Show)++instance AffineTransformable Ray3D where+    transform m (Ray3D p v) = Ray3D (transform m p) (transform m v)++-- | The parametric function of a ray.  The parameter is measured as a+-- proportion of the length of the vector.  @projectRay 0@ is the+-- endpoint of the ray.  @projectRay 1@ is the endpoint offset+-- by the ray's vector.+projectRay :: RSdouble -> Ray3D -> Point3D+projectRay t (Ray3D (Point3D x y z) (Vector3D u v w)) = Point3D (x+u*t) (y+v*t) (z+w*t)++-- | The inverse operation to 'projectRay'.  This could also be+-- understood as the height of the point above the plane defined+-- by the ray.+distanceAlong :: Ray3D -> Point3D -> RSdouble+distanceAlong (Ray3D p v) p' = dotProduct (vectorToFrom p' p) v / vectorLengthSquared v++-- | The angle between vector of the ray and the vector from the+-- endpoint of the ray to the specified point.+angleFrom :: Ray3D -> Point3D -> Angle+angleFrom (Ray3D p v) p' = angleBetween (vectorToFrom p' p) v++-- | A ray normalize to a length of 1.+normalizeRay :: Ray3D -> Ray3D+normalizeRay (Ray3D p v) = Ray3D p $ vectorNormalize v+
+ RSAGL/Math/Types.hs view
@@ -0,0 +1,70 @@+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+module RSAGL.Math.Types+    (RSfloat,+     RSdouble,+     toGLfloat,+     toGLdouble,+     fromGLfloat,+     fromGLdouble,+     f2f)+    where++import Graphics.Rendering.OpenGL.Raw.Core31+import Data.Vec as Vec+import Data.Vec.OpenGLRaw+import System.Random+import Control.Arrow++newtype RSfloat = RSfloat GLfloat+    deriving (Enum,+              Eq,+              Floating,+              Fractional,+              Num,+              Ord,+              Read,+              Real,+              RealFloat,+              RealFrac,+              Show,+              NearZero)+newtype RSdouble = RSdouble GLdouble+    deriving (Enum,+              Eq,+              Floating,+              Fractional,+              Num,+              Ord,+              Read,+              Real,+              RealFloat,+              RealFrac,+              Show,+              NearZero)++toGLfloat :: RSfloat -> GLfloat+toGLfloat (RSfloat x) = x++toGLdouble :: RSdouble -> GLdouble+toGLdouble (RSdouble x) = x++fromGLfloat :: GLfloat -> RSfloat+fromGLfloat = RSfloat++fromGLdouble :: GLdouble -> RSdouble+fromGLdouble = RSdouble++{-# RULES+"f2f/id"    f2f = id #-}+f2f :: (RealFloat a,RealFloat b) => a -> b+f2f = uncurry encodeFloat . decodeFloat+{-# INLINE f2f #-}++instance Random RSfloat where+    randomR (lo,hi) g = first f2f $ randomR (f2f lo, f2f hi :: Float) g+    random g = first (f2f :: Float -> RSfloat) $ random g++instance Random RSdouble where+    randomR (lo,hi) g = first f2f $ randomR (f2f lo, f2f hi :: Double) g+    random g = first (f2f :: Double -> RSdouble) $ random g+
+ RSAGL/Math/Vector.lhs view
@@ -0,0 +1,298 @@+\section{Points and Vectors: RSAGL.Vector}++\begin{code}+{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, FunctionalDependencies, PatternGuards #-}+module RSAGL.Math.Vector+    (Point3D(..),+     origin_point_3d,+     Vector3D(..),+     SurfaceVertex3D(..),+     zero_vector,+     point3d,+     points3d,+     point2d,+     points2d,+     vector3d,+     dotProduct,+     angleBetween,+     crossProduct,+     distanceBetween,+     distanceBetweenSquared,+     aNonZeroVector,+     aLargeVector,+     vectorAdd,+     vectorSum,+     vectorScale,+     vectorScaleTo,+     vectorToFrom,+     vectorNormalize,+     vectorAverage,+     vectorLength,+     vectorLengthSquared,+     newell,+     Xyz(..),+     XYZ,+     vectorString,+     randomXYZ,+     fixOrtho,+     fixOrtho2,+     fixOrtho2Left,+     orthos)+    where++import Control.Parallel.Strategies+import RSAGL.Math.Angle+import System.Random+import RSAGL.Math.AbstractVector+import RSAGL.Math.Types+import RSAGL.Math.ListUtils+\end{code}++\subsection{Generic 3-dimensional types and operations}++\begin{code}+type XYZ = (RSdouble,RSdouble,RSdouble)++class Xyz a where+    toXYZ :: a -> XYZ+    fromXYZ :: XYZ -> a++instance Xyz (RSdouble,RSdouble,RSdouble) where+    toXYZ = id+    fromXYZ = id++vectorString :: Xyz a => a -> String+vectorString xyz = let (x,y,z) = toXYZ xyz+		       in (show x) ++ "," ++ (show y) ++ "," ++ (show z)++uncurry3d :: (RSdouble -> RSdouble -> RSdouble -> a) -> XYZ -> a+uncurry3d fn (x,y,z) = fn x y z+\end{code}++\subsection{Points in 3-space}++\begin{code}+data Point3D = Point3D {-# UNPACK #-} !RSdouble {-# UNPACK #-} !RSdouble {-# UNPACK #-} !RSdouble+	     deriving (Read,Show,Eq)++origin_point_3d :: Point3D+origin_point_3d = Point3D 0 0 0++point3d :: (RSdouble,RSdouble,RSdouble) -> Point3D+point3d = uncurry3d Point3D++point2d :: (RSdouble,RSdouble) -> Point3D+point2d (x,y) = point3d (x,y,0)++points3d :: [(RSdouble,RSdouble,RSdouble)] -> [Point3D]+points3d = map point3d++points2d :: [(RSdouble,RSdouble)] -> [Point3D]+points2d = map point2d++instance Xyz Point3D where+    toXYZ (Point3D x y z) = (x,y,z)+    fromXYZ (x,y,z) = Point3D x y z++instance AbstractZero Point3D where+    zero = origin_point_3d++instance AbstractAdd Point3D Vector3D where+    add = displace++instance AbstractSubtract Point3D Vector3D where+    sub = vectorToFrom++instance NFData Point3D+\end{code}++\subsection{Vectors in 3-space}++\begin{code}+data Vector3D = Vector3D {-# UNPACK #-} !RSdouble {-# UNPACK #-} !RSdouble {-# UNPACK #-} !RSdouble+	      deriving (Read,Show,Eq)++zero_vector :: Vector3D+zero_vector = Vector3D 0 0 0++vector3d :: (RSdouble,RSdouble,RSdouble) -> Vector3D+vector3d = uncurry3d Vector3D++instance Xyz Vector3D where+    toXYZ (Vector3D x y z) = (x,y,z)+    fromXYZ (x,y,z) = Vector3D x y z++instance NFData Vector3D++instance AbstractZero Vector3D where+    zero = zero_vector++instance AbstractAdd Vector3D Vector3D where+    add = vectorAdd++instance AbstractSubtract Vector3D Vector3D where+    sub = vectorToFrom++instance AbstractScale Vector3D where+    scalarMultiply = vectorScale++instance AbstractMagnitude Vector3D where+    magnitude = vectorLength++instance AbstractVector Vector3D+\end{code}++A \texttt{SurfaceVertex3D} is a a point on an orientable surface, having a position and a normal vector.++\subsection{Surface Vertices in 3-space}++\begin{code}+data SurfaceVertex3D = SurfaceVertex3D { sv3d_position :: Point3D, +                                         sv3d_normal :: Vector3D }+    deriving (Read,Show)++instance NFData SurfaceVertex3D where+    rnf (SurfaceVertex3D p v) = rnf (p,v)+\end{code}++\subsection{Vector Arithmetic}++\begin{code}+aNonZeroVector :: Vector3D -> Maybe Vector3D+aNonZeroVector v = case vectorLength v of+    x | x <= 0 -> Nothing+    x | isDenormalized x -> Nothing+    x | isNaN x -> Nothing+    x | isInfinite x -> Nothing+    _ | otherwise -> Just v++aLargeVector :: RSdouble -> Vector3D -> Maybe Vector3D+aLargeVector x v_ =+    case aNonZeroVector v_ of+        Just v | vectorLength v > x -> Just v+        _ | otherwise -> Nothing++dotProduct :: Vector3D -> Vector3D -> RSdouble+dotProduct (Vector3D ax ay az) (Vector3D bx by bz) = +    (ax*bx) + (ay*by) + (az*bz)++angleBetween :: Vector3D -> Vector3D -> Angle+angleBetween a b = arcCosine $ dotProduct (vectorNormalize a) (vectorNormalize b)++crossProduct :: Vector3D -> Vector3D -> Vector3D+crossProduct (Vector3D ax ay az) (Vector3D bx by bz) = +    Vector3D (ay*bz - az*by) (az*bx - ax*bz) (ax*by - ay*bx)++distanceBetween :: (Xyz xyz) => xyz -> xyz -> RSdouble+distanceBetween a b = vectorLength $ vectorToFrom a b++distanceBetweenSquared :: (Xyz xyz) => xyz -> xyz -> RSdouble+distanceBetweenSquared a b = vectorLengthSquared $ vectorToFrom a b++displace :: (Xyz xyz) => xyz -> Vector3D -> xyz+displace xyz (Vector3D x2 y2 z2) =+    let (x1,y1,z1) = toXYZ xyz+        in fromXYZ (x1+x2,y1+y2,z1+z2)++vectorAdd :: Vector3D -> Vector3D -> Vector3D+vectorAdd = displace++vectorSum :: [Vector3D] -> Vector3D+vectorSum vectors = foldr vectorAdd zero_vector vectors++vectorToFrom :: (Xyz xyz) => xyz -> xyz -> Vector3D+vectorToFrom a b = +    let (ax,ay,az) = toXYZ a+        (bx,by,bz) = toXYZ b+        in Vector3D (ax - bx) (ay - by) (az - bz)++vectorLength :: Vector3D -> RSdouble+vectorLength = sqrt . vectorLengthSquared++vectorLengthSquared :: Vector3D -> RSdouble+vectorLengthSquared (Vector3D x y z) = (x*x + y*y + z*z)++vectorScale :: RSdouble -> Vector3D -> Vector3D+vectorScale s (Vector3D x y z) = Vector3D (x*s) (y*s) (z*s)+\end{code}++vectorScaleTo forces the length of a vector to a certain value, without changing the vector's direction.++\begin{code}+vectorScaleTo :: RSdouble -> Vector3D -> Vector3D+vectorScaleTo new_length vector = vectorScale new_length $ vectorNormalize vector+\end{code}++vectorNormalize forces the length of a vector to 1, without changing the vector's direction.++\begin{code}+vectorNormalize :: Vector3D -> Vector3D+vectorNormalize v = +    let l = vectorLength v+	in maybe (Vector3D 0 0 0) (vectorScale (1/l)) $ aNonZeroVector v+\end{code}++vectorAverage takes the average of a list of vectors.  The result is a normalized vector+where the length of the element vectors is not reflected in the result.++\begin{code}+vectorAverage :: [Vector3D] -> Vector3D+vectorAverage vects = vectorNormalize $ vectorSum $ map vectorNormalize vects+\end{code}++\subsection{Generating normal vectors}++\texttt{newell} calculates the normal vector of an arbitrary polygon.  If the points specified are non-coplanar,+\texttt{newell} often calculates a reasonable result; if they are colinear or singular \texttt{newell} will fail.++The result is a normalized vector.++\begin{code}+newell :: [Point3D] -> Maybe Vector3D+newell points = fmap vectorNormalize $ aNonZeroVector $ vectorSum $ map newell_ $ loopedDoubles points+    where newell_ (Point3D x0 y0 z0,Point3D x1 y1 z1) =+              (Vector3D+               ((y0 - y1)*(z0 + z1))+               ((z0 - z1)*(x0 + x1))+               ((x0 - x1)*(y0 + y1)))+\end{code}++\subsection{Randomly Generated Coordinates}++\texttt{randomXYZ} can generate random coordinates within the cube where x, y, and z are each in the range (lo,hi).++\begin{code}+randomXYZ :: (RandomGen g,Xyz p) => (RSdouble,RSdouble) -> g -> (p,g)+randomXYZ (lo,hi) g = (fromXYZ (f2f x,f2f y,f2f z),g')+    where (g_,g') = split g+          (x:y:z:_) = randomRs (f2f lo,f2f hi) g_ :: [Double]+\end{code}++\subsection{Orthagonal Vectors}++\texttt{fixOrtho a v} finds the vector, orthogonal to a, that has the least angle to v.++\texttt{fixOrtho2 right up} yields \texttt{(up,forward)}.++\texttt{fixOrtho2Left right up} yields \texttt{(up,backward)}.++\texttt{orthos} finds two arbitrary vectors orthogonal to the parameter.++\begin{code}+fixOrtho :: Vector3D -> Vector3D -> Vector3D+fixOrtho a = fst . fixOrtho2 a++fixOrtho2 :: Vector3D -> Vector3D -> (Vector3D,Vector3D)+fixOrtho2 a v = (vectorNormalize $ crossProduct a $ vectorScale (-1) b,vectorNormalize b)+    where b = crossProduct a v++fixOrtho2Left :: Vector3D -> Vector3D -> (Vector3D,Vector3D)+fixOrtho2Left a v = (vectorNormalize $ crossProduct a b,vectorNormalize b)+    where b = vectorScale (-1) $ crossProduct a v++orthos :: Vector3D -> (Vector3D,Vector3D)+orthos v@(Vector3D x y z) | abs y >= abs x && abs z >= abs x = fixOrtho2 v (Vector3D (abs x + abs y + abs z) y z)+orthos v@(Vector3D x y z) | abs x >= abs y && abs z >= abs y = fixOrtho2 v (Vector3D x (abs x + abs y + abs z) z)+orthos v@(Vector3D x y z) | abs x >= abs z && abs y >= abs z = fixOrtho2 v (Vector3D x y (abs x + abs y + abs z))+orthos v = error $ "orthos: (" ++ show v ++ ")"+\end{code}
+ Setup.hs view
@@ -0,0 +1,5 @@+#!/usr/bin/runhaskell++import Distribution.Simple++main = defaultMainWithHooks simpleUserHooks
+ rsagl-math.cabal view
@@ -0,0 +1,52 @@+name:                rsagl-math+version:             0.6.0.0+license:             BSD3+license-file:        LICENSE+author:              Christopher Lane Hinson+maintainer:          Christopher Lane Hinson <lane@downstairspeople.org>++category:            Math+synopsis:            The RogueStar Animation and Graphics Library: Mathematics+description:         RSAGL, the RogueStar Animation and Graphics Library,+                     was specifically designed for a computer game called+                     roguestar, but effort has been made to make it accessable+                     to other projects that might benefit from it.  This+                     package includes mathematical algorithms to support+                     the project.+cabal-version:       >= 1.2+homepage:            http://roguestar.downstairspeople.org/++build-type:          Simple+tested-with:         GHC==6.12.1++Library+    exposed-modules:     RSAGL.Math+                         RSAGL.Math.AbstractVector,+                         RSAGL.Math.Affine,+                         RSAGL.Math.Angle,+                         RSAGL.Math.BoundingBox,+                         RSAGL.Math.Curve,+                         RSAGL.Math.CurveExtras,+                         RSAGL.Math.FMod,+                         RSAGL.Math.Interpolation,+                         RSAGL.Math.ListUtils,+                         RSAGL.Math.Matrix,+                         RSAGL.Math.Orthogonal,+                         RSAGL.Math.Ray,+                         RSAGL.Math.Types+                         RSAGL.Math.Vector++    ghc-options:         -fno-warn-type-defaults -fexcess-precision+    ghc-prof-options:    -prof -auto-all++    build-depends:       base>=4 && <5,+                         random>= 1.0.0.2 && < 1.1,+                         array>= 0.3.0.0 && < 0.4,+                         containers>= 0.3.0.0 && < 0.4,+                         OpenGL>= 2.4.0.1 && < 2.5,+                         OpenGLRaw>= 1.1.0.1 && < 1.2,+                         parsec>=3.1.0 && < 3.2,+                         parallel>=2.2.0.1 && < 2.3,+                         Vec>=0.9.8 && < 0.10,+                         Vec-OpenGLRaw>=0.2.0.0 && < 0.3+