diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,4 +1,5 @@
 -*-change-log-*-
+	0.1.0.14 Added pendulum example
 	0.1.0.13 Adding ability to create linear operators from matrices; made linear operators a category and monoid
 	0.1.0.12 Adding type signatures to top level stuff
 	0.1.0.11 More inlining and specialisation
@@ -6,7 +7,7 @@
 	0.1.0.9 Inlined/specialised a bunch of function, hueg speed increase
         0.1.0.8 Implemented algebraic/transcendental typeclasses
 	0.1.0.7 Adding basic linear operators; made multivector a field
-	0.1.0.6 Memoising the blade index comparision function for a 20% speed increase 
+	0.1.0.6 Memoising the blade index comparision function for a 20% speed increase
 	0.1.0.5 Adding hspec tests, fixed blade comparison to order blades in terms of grade first
 
 	0.1.0.4 Made multivectors have a (p,q) metric signature at the type level
diff --git a/clifford.cabal b/clifford.cabal
--- a/clifford.cabal
+++ b/clifford.cabal
@@ -10,7 +10,7 @@
 -- PVP summary:      +-+------- breaking API changes
 --                   | | +----- non-breaking API additions
 --                   | | | +--- code changes with no API change
-version:             0.1.0.13
+version:             0.1.0.14
 
 -- A short (one-line) description of the package.
 synopsis:            A Clifford algebra library
@@ -57,11 +57,11 @@
   -- other-modules:       
   
   -- LANGUAGE extensions used by modules in this package.
-  other-extensions:    NoImplicitPrelude, ScopedTypeVariables, GADTs, DataKinds,KindSignatures,UnicodeSyntax,FlexibleContexts, RankNTypes,TemplateHaskell,NoMonomorphismRestriction,MultiParamTypeClasses,FlexibleInstances, TypeOperators
+  other-extensions:    NoImplicitPrelude, ScopedTypeVariables, GADTs, DataKinds,KindSignatures,UnicodeSyntax,FlexibleContexts, RankNTypes,TemplateHaskell,NoMonomorphismRestriction,MultiParamTypeClasses,FlexibleInstances, TypeOperators, UnicodeSyntax, ConstraintKinds
   -- Other library packages from which modules are imported.
   build-depends:       base >=4.6 && <4.9, numeric-prelude >= 0.4.0.1 && < 0.5.0, permutation >= 0.4.1 && < 0.5, 
                        data-ordlist >= 0.4.5 && < 0.5,  converge >= 0.1.0.1 && < 0.2, lens >= 4.0.3 && < 4.1, 
-                       deepseq >= 1.3.0.1 && < 1.4, vector >= 0.10.0.1 && < 0.11, stream-fusion >= 0.1 && < 0.2, criterion >= 0.8.0.0 && < 0.9, derive, QuickCheck, nats, tagged, cereal,hspec, MemoTrie >= 0.6 && < 0.7, semigroupoids >= 4.0 && < 4.1, monoid-extras >=0.3 && <0.4
+                       deepseq >= 1.3.0.1 && < 1.4, vector >= 0.10.0.1 && < 0.11, stream-fusion >= 0.1 && < 0.2, criterion >= 0.8.0.0 && < 0.9, derive, QuickCheck, nats, tagged, cereal,hspec, MemoTrie >= 0.6 && < 0.7, semigroupoids >= 4.0 && < 4.1, monoid-extras >=0.3 && <0.4, reflection >=1.4 && < 1.5
   
   ghc-options: -fllvm -fexcess-precision -optlc-O=3 -O3
   -- Directories containing source files.
@@ -70,6 +70,15 @@
   -- Base language which the package is written in.
   default-language:    Haskell2010
   
+
+executable pendulum
+  build-depends: base, clifford, numeric-prelude,gnuplot,stream-fusion, Chart, Chart-cairo, colour, lens, data-default-class
+  other-extensions: NoImplicitPrelude
+  main-is: Pendulum.hs
+  hs-source-dirs: examples
+  default-language: Haskell2010
+  ghc-options: -fllvm -fexcess-precision -optlc-O=3 -O3
+	
 test-suite spec
   type: exitcode-stdio-1.0
   default-extensions: DataKinds, ScopedTypeVariables
@@ -86,5 +95,5 @@
   main-is: benchmarks.hs
   build-depends: base, clifford, criterion, numeric-prelude, stream-fusion
   -- -optlo-O3
-  ghc-options: -fllvm -fexcess-precision -optlc-O=3 -O3
+  ghc-options: -fllvm -fexcess-precision -optlc-O=3 -O3 
   default-language: Haskell2010
diff --git a/examples/Pendulum.hs b/examples/Pendulum.hs
new file mode 100644
--- /dev/null
+++ b/examples/Pendulum.hs
@@ -0,0 +1,53 @@
+{-# LANGUAGE NoImplicitPrelude, NoMonomorphismRestriction, DataKinds #-}
+import NumericPrelude (Double, fst, snd, ($), (.), seq)
+import Prelude (getLine, putStrLn)
+import Algebra.Transcendental
+import Data.List.Stream
+import Numeric.Clifford.Multivector
+import Algebra.Ring 
+import Algebra.Additive 
+import Algebra.Field
+import Numeric.Clifford.NumericIntegration
+import Numeric.Clifford.NumericIntegration.DefaultIntegrators
+import Numeric.Clifford.Blade
+
+import Control.Lens
+import Graphics.Rendering.Chart
+import Data.Colour
+import Data.Colour.Names
+import Data.Default.Class
+import Graphics.Rendering.Chart.Backend.Cairo
+
+m = scalar 3 :: E3Vector
+l = scalar 20 :: E3Vector
+g = scalar 9.81 :: E3Vector
+
+hamil _ [p,theta] =   [ (-m*g*l)* sin theta, p / (m*l*l)] 
+
+integrator = gaussLegendreFourthOrder 0.1 hamil
+
+
+
+tenSeconds :: [(Double,Double,Double)]
+tenSeconds = map ((\ (t, ([BladeSum [Blade a []],BladeSum [Blade b []]])) -> (t,a,b) ) )  $ take 5001 $ iterate integrator (0,[zero,one/10])
+
+
+chart = toRenderable layout
+  where
+    momentum = plot_lines_values .~ [ ( map (\(t,p,_) -> (t,p)) tenSeconds )]
+                     $ plot_lines_style  . line_color .~ opaque blue
+                     $ plot_lines_title .~ "momentum"
+                     $ def
+
+    angle = plot_lines_style . line_color .~ (opaque red)
+                  $ plot_lines_values .~ [ ( map (\(t,_,theta) -> (t,theta)) tenSeconds )]
+                  $ plot_lines_title .~ "angle"
+                  $ def
+
+    layout = layout_title .~ "Pendulum"
+           $ layout_plots .~ [toPlot momentum,
+                              toPlot angle]
+           $ def
+
+main = renderableToFile def  chart "pendulum.png"
+
diff --git a/src/Numeric/Clifford/Blade.lhs b/src/Numeric/Clifford/Blade.lhs
--- a/src/Numeric/Clifford/Blade.lhs
+++ b/src/Numeric/Clifford/Blade.lhs
@@ -63,6 +63,7 @@
 \texttt{bScale} is the amplitude of the blade. \texttt{bIndices} are the indices for the basis. 
 \begin{code}
 
+
 data Blade (p :: Nat) (q :: Nat) f where
     Blade :: forall (p::Nat) (q::Nat) f . (Algebra.Field.C f, SingI p, SingI q) => {_scale :: f, _indices :: [Natural]} -> Blade p q f
 
@@ -138,7 +139,7 @@
 {-#INLINE bladeNormalForm#-}
 {-#SPECIALISE INLINE bladeNormalForm::E3Blade -> E3Blade #-}
 {-#SPECIALISE INLINE bladeNormalForm :: STBlade -> STBlade #-}
-bladeNormalForm :: forall (p::Nat) (q::Nat) f.  Blade p q f -> Blade p q f
+bladeNormalForm :: ∀ (p::Nat) (q::Nat) f.  Blade p q f -> Blade p q f
 bladeNormalForm (Blade scale indices)  = result 
         where
              result = if (any (\i -> (GHC.Real.toInteger i) >= d) indices) then zeroBlade else Blade scale' newIndices
diff --git a/src/Numeric/Clifford/ClassicalMechanics.lhs b/src/Numeric/Clifford/ClassicalMechanics.lhs
--- a/src/Numeric/Clifford/ClassicalMechanics.lhs
+++ b/src/Numeric/Clifford/ClassicalMechanics.lhs
@@ -19,12 +19,10 @@
 {-# LANGUAGE FlexibleInstances, TypeOperators #-}
 {-# LANGUAGE TemplateHaskell #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# OPTIONS_HADDOCK show-extensions #-}
 \end{code}
 %if False
 \begin{code}
-{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}
--- OPTIONS_GHC -Odph -fvectorise -package dph-lifted-vseg 
---  LANGUAGE ParallelArrays
 \end{code}
 %endif
 
@@ -34,7 +32,7 @@
 import Numeric.Clifford.Blade
 import GHC.TypeLits
 import Data.Proxy
-import NumericPrelude hiding (iterate, head, map, tail, reverse, scanl, zipWith, drop, (++), filter, null, length, foldr, foldl1, zip, foldl, concat, (!!), concatMap,any, repeat, replicate, elem, replicate, all)
+import NumericPrelude hiding (iterate, head, map, tail, reverse, scanl, zipWith, drop, (++), filter, null, length, foldr, foldl1, zip, foldl, concat, (!!), concatMap,any, repeat, replicate, elem, replicate, all, (.) )
 import Algebra.Absolute
 import Algebra.Algebraic
 import Algebra.Additive
@@ -56,9 +54,15 @@
 import Data.Maybe
 import Data.DeriveTH
 import Data.Word
-import Debug.Trace
---trace _ a = a
+import Numeric.Clifford.Internal
+import Numeric.Clifford.LinearOperators
+import Control.Applicative
+import Data.Monoid
+import Control.Category
 
+nonEqualFrames = "Non-equal reference frames! Insert code here to translate between them! :) Should really make reference frames as types and then have type families to convert between them :v :v :v"
+
+
 data EnergyMethod (p::Nat) (q::Nat) f = Hamiltonian{ _dqs :: [DynamicSystem p q f -> Multivector p q f], _dps :: [DynamicSystem p q f -> Multivector p q f]}
 
 data DynamicSystem (p::Nat) (q::Nat) f = DynamicSystem {_time :: f, coordinates :: [Multivector p q f], _momenta :: [Multivector p q f], _energyFunction :: EnergyMethod p q f, _projector :: DynamicSystem p q f -> DynamicSystem p q f}
@@ -80,23 +84,49 @@
 
 Now to make a physical object.
 \begin{code}
-data ReferenceFrame (p::Nat) (q::Nat) t = ReferenceFrame {basisFrame :: [Multivector p q t]}
-psuedoScalar' :: forall f (p::Nat) (q::Nat). (Ord f, Algebra.Field.C f, SingI p, SingI q) => ReferenceFrame p q f -> Multivector p q f
-psuedoScalar'  = multiplyList . basisFrame
+data ReferenceFrame (p::Nat) (q::Nat) t = RelativeFrame {frameName :: String, euclideanMove :: EuclideanMove p q t, velocityRelToParentFrame :: Multivector p q t, angularVelocityRelToParentFrame :: Multivector p q t, parent :: ReferenceFrame p q t}
+                                        |GlobalAbsoluteFrame deriving (Eq, Show)
 
 
+getRigidDisplacementRelToInertial :: (Algebra.Field.C t, Ord t, SingI p, SingI q) =>  ReferenceFrame p q t -> EuclideanMove   p q t
+getRigidDisplacementRelToInertial GlobalAbsoluteFrame = mempty
+--getRigidDisplacementRelToInertial (RelativeFrame _ displacement mother) = displacement <> (getRigidDisplacementRelToInertial mother)
 
-a `cross` b = (negate $ one)`e`[0,1,2] * (a ∧ b)
+getFrameTransformation :: forall (p::Nat) (q::Nat) t . (Algebra.Field.C t, Ord t, SingI p, SingI q) =>  ReferenceFrame p q t -> ReferenceFrame p q t -> EuclideanMove p q t
+getFrameTransformation r' r = undefined
+
+{-data InertialFrame (p::Nat) (q::Nat) f t = InertialFrame {objects :: t, changeFrame :: t -> EuclideanMove p q f -> t, frame :: ReferenceFrame p q f}
+
+
+instance Functor (InertialFrame p q f) where
+    fmap func (InertialFrame objs changeFrame frame) = InertialFrame (func objs) (changeFrame . func) frame
+
+instance (SingI p, SingI q) => Applicative (InertialFrame p q f) where
+    pure a = InertialFrame a GlobalAbsoluteFrame where 
+    (<*>) (InertialFrame func trans1 frame1) (InertialFrame objs trans2 frame2) = if (name frame1)==(name frame2)
+                                                                    then InertialFrame (func objs) frame1 
+                                                                    else InertialFrame (trans2 (func objs) (getFrameTransformation frame2 frame1)) trans2 frame1
+
+
+instance (SingI p, SingI q, Algebra.Field.C f, Ord f) => Monad (InertialFrame p q f) where
+    return = pure
+    (>>=) (InertialFrame objA changeFrameA frameA)  func = undefined where
+        (InertialFrame objB changeFrameB frameB) = func objA
+        
+-}
+
+a `cross` b = (negate psuedoScalar) * (a ∧ b)
+
+
+
 data PhysicalVector (p::Nat) (q::Nat) t = PhysicalVector {r :: Multivector p q t, referenceFrame :: ReferenceFrame p q t}
 
 
+
+
+
 data RigidBody (p::Nat) (q::Nat) f where
- RigidBody:: (Algebra.Field.C f, Algebra.Module.C f (Multivector p q f)) =>  {position :: PhysicalVector p q f,
-                              momentum :: PhysicalVector p q f,
-                              mass :: f,
-                              attitude :: PhysicalVector p q f,
-                              angularMomentum :: PhysicalVector p q f,
-                              inertia :: PhysicalVector p q f
+ RigidBody:: (Algebra.Field.C f, Algebra.Module.C f (Multivector p q f)) =>  {bodyName::String, frame::ReferenceFrame p q f, mass :: f, inertia :: Multivector p q f, position :: Multivector p q f, attitude :: Multivector p q f, velocity :: Multivector p q f, angularVelocity :: Multivector p q f
                              } -> RigidBody p q f
 
 --makeLenses ''RigidBody doesn't actually work
@@ -105,8 +135,8 @@
 5. figure a way to take exterior product of 1 forms at a type level so i can just go like: omega = df1 ^ df2 ^ df ; omega a b c
 -}
 
-type Vector3 f =  Multivector 3 0 f
-type STVector f = Multivector 3 1 f
+
+
 \end{code}
 \bibliographystyle{IEEEtran}
 \bibliography{biblio.bib}
diff --git a/src/Numeric/Clifford/Internal.hs b/src/Numeric/Clifford/Internal.hs
--- a/src/Numeric/Clifford/Internal.hs
+++ b/src/Numeric/Clifford/Internal.hs
@@ -1,6 +1,6 @@
 {-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}
-{-# LANGUAGE TypeOperators, TypeFamilies,CPP #-}
-module Numeric.Clifford.Internal (myTrace, trie, untrie, enumerate, dimension) where
+{-# LANGUAGE TypeOperators, TypeFamilies,CPP, ConstraintKinds, RankNTypes, DataKinds #-}
+module Numeric.Clifford.Internal (myTrace, trie, untrie, enumerate, dimension, DefaultField, AllowableCliffordType) where
 import Numeric.Natural
 import Prelude hiding (head,tail, null)
 import Data.MemoTrie
@@ -9,12 +9,19 @@
 import Data.Bits
 import Test.QuickCheck
 import Data.Word
+import GHC.TypeLits
+import Algebra.Field
 import qualified Debug.Trace as DebugTrace
 #ifdef DEBUG
 myTrace = DebugTrace.trace
 #else
 myTrace _ x = x
 #endif
+
+
+type AllowableCliffordType p q f = forall (p::Nat) (q::Nat) f. (Ord f, Algebra.Field.C f, SingI p, SingI q)
+type DefaultField = Double
+
 instance HasTrie Natural where
     newtype Natural :->: a = NaturalTrie ((Bool,[Bool]) :->: a)
     trie f = NaturalTrie (trie (f . unbitsZ)) 
diff --git a/src/Numeric/Clifford/LinearOperators.lhs b/src/Numeric/Clifford/LinearOperators.lhs
--- a/src/Numeric/Clifford/LinearOperators.lhs
+++ b/src/Numeric/Clifford/LinearOperators.lhs
@@ -1,8 +1,9 @@
 \begin{code}
-{-# LANGUAGE NoImplicitPrelude, RankNTypes, KindSignatures, DataKinds, GADTs, FlexibleInstances, UndecidableInstances, InstanceSigs, MultiParamTypeClasses #-}
+{-# LANGUAGE NoImplicitPrelude, RankNTypes, KindSignatures, DataKinds, GADTs, FlexibleInstances, UndecidableInstances, InstanceSigs, MultiParamTypeClasses, PolyKinds, ConstraintKinds, UnicodeSyntax, StandaloneDeriving #-}
+{-# OPTIONS_HADDOCK show-extensions #-}
 module Numeric.Clifford.LinearOperators where
 import qualified NumericPrelude as NP ((.), id)
-import NumericPrelude hiding ((.), id, (!!), zipWith, map, length)
+import NumericPrelude hiding ((.), id, (!!), zipWith, map, length, zipWith3, and)
 import Numeric.Clifford.Multivector
 import Algebra.Algebraic
 import Algebra.Field
@@ -20,6 +21,7 @@
 import Data.Semigroupoid
 import Numeric.Natural
 import Data.Word
+import Algebra.Additive
 import Numeric.Clifford.Internal
 import qualified Numeric.Clifford.Blade
 \end{code}
@@ -31,38 +33,87 @@
 data LinearOperator' p q f g where
     LinearOperator' :: {_operator' :: Multivector p q f -> Multivector p q g} -> LinearOperator' p q f g
     LinearOperator :: {_operator :: Multivector p q f -> Multivector p q f} -> LinearOperator' p q f f
+
+getFuncFromOperator :: LinearOperator' p q f g → (Multivector p q f → Multivector p q g)
+getFuncFromOperator (LinearOperator' op) = op
+getFuncFromOperator (LinearOperator op) = op
+
 type LinearOperator p q f = LinearOperator' p q f f
 type LinearOperatorCreator p q f = (Algebra.Algebraic.C f, Ord f, SingI p, SingI q) => Multivector p q f -> LinearOperator p q f
 
+instance (Show g, f ~ g) => Show (LinearOperator' p q f g) where
+    show = show . getMatrixElementsFromOperator
+
+instance (Algebra.Field.C f, Algebra.Field.C g, Ord f,  Ord g, SingI p, SingI q) => Eq (LinearOperator' p q f g) where
+    a == b = and (map (\ e → (f1 e) == (f2 e)) basisVectors) where
+           f1 = getFuncFromOperator a
+           f2 = getFuncFromOperator b
+
 instance Category (LinearOperator' p q) where
     id = LinearOperator' NP.id
     (.) (LinearOperator' a) (LinearOperator' b)  = LinearOperator' (a NP.. b)
 
+ 
 instance (Algebra.Field.C f, Ord f,Algebra.Field.C g, Ord g, SingI p, SingI q, f~g) => Monoid (LinearOperator' p q f g) where
     mempty = id
     mappend = (.)
 
+data EuclideanMove p q f where
+    EuclideanMove :: ∀ (p::Nat) (q::Nat) f.  (Algebra.Field.C f, Ord f, SingI p, SingI q) => { _rotation :: Multivector p q f, _translation :: Multivector p q f} -> EuclideanMove p q f
 
+deriving instance Eq(EuclideanMove p q f)
+deriving instance (Show f) => Show (EuclideanMove p q f)
 
+applyEuclideanMove (EuclideanMove r a) x = (rotate r x) + a
 
-class LinearOperatorClass' (p::Nat) (q::Nat) f g where
 
+
+
+instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Monoid (EuclideanMove p q f) where
+    mempty = EuclideanMove one zero
+    mappend (EuclideanMove s b) (EuclideanMove r a) = EuclideanMove rot trans where
+          rot = r*s
+          trans = (rotate s a) + b
+
+{-instance ∀ a b (p::Nat) (q::Nat).(Algebra.Field.C a, SingI p, SingI q, Ord a, Algebra.Field.C b, Ord b) => Category (AffineOperator' p q) where
+    id:: (Algebra.Field.C c) => AffineOperator' p q c c
+    id = AffineOperator id zero
+    (.) = undefined -}
+
 {-
+--GHC 7.8: The Control.Category module now has the PolyKinds extension enabled, meaning that instances of Category no longer need be of kind * -> * -> *.
+class Operator (p::Nat) (q::Nat) f g where
+    apply :: Operator p q f g ->  Multivector p q f -> Multivector p q g
+
+instance forall (p::Nat) (q::Nat) . Category (Operator p q) where
+    id = NP.id
+    (.) a b = a (NP..) b
+-}  
+
+{-
 [[f11, f12, f13],
  [f21, f22, f21],
  [f31, f32, f33]]
 -}
-createFunctionalFromElements :: forall (p::Nat) (q::Nat) f . (Algebra.Field.C f, Ord f, SingI p, SingI q) => [[f]] ->(Multivector p q f -> Multivector p q f)
+
+getMatrixElementsFromOperator :: LinearOperator' p q f' f'-> [[f']]
+getMatrixElementsFromOperator operator = error "Numeric.Clifford.LinearOperator.getMatrixElementsFromOperator not implemented yet!" where
+    f = getFuncFromOperator operator
+    
+   
+
+createFunctionalFromElements :: ∀  (p::Nat) (q::Nat) f . (Algebra.Field.C f, Ord f, SingI p, SingI q) => [[f]] ->(Multivector p q f -> Multivector p q f)
 createFunctionalFromElements elements = (\x -> f*x) where
     d = (length elements) - 1
     f = sumList $ map elementsForK [0..d]
     column k = let transposed = transpose elements in transposed !! k   
     elementsForK k =sumList $   zipWith (scaleRight) basisVectors (column k) 
     
-createLinearOperatorFromElements :: forall (p::Nat) (q::Nat) f . (Algebra.Field.C f, Ord f, SingI p, SingI q) => [[f]] -> LinearOperator p q f
+createLinearOperatorFromElements :: ∀ (p::Nat) (q::Nat) f . (Algebra.Field.C f, Ord f, SingI p, SingI q) => [[f]] -> LinearOperator p q f
 createLinearOperatorFromElements  = LinearOperator .  createFunctionalFromElements
 
 
+reflect :: (Algebra.Algebraic.C f, Ord f, SingI p, SingI q) => Multivector p q f -> Multivector p q f -> Multivector p q f
 reflect u x = (-u)*x*recip u
 
 makeReflectionOperator ::LinearOperatorCreator p q f
@@ -73,10 +124,10 @@
 
 makeRotationOperator :: LinearOperatorCreator p q f
 makeRotationOperator u = LinearOperator (rotate u)
-
 makeRotationOperatorFromPlaneAngle plane angle = LinearOperator (rotatePlaneAngle plane angle)
-project u x = inverse u * (u `dot` x)
 
+
+project u x = inverse u * (u `dot` x)
 makeProjectionOperator :: LinearOperatorCreator p q f
 makeProjectionOperator u = LinearOperator (project u)
 
diff --git a/src/Numeric/Clifford/Multivector.lhs b/src/Numeric/Clifford/Multivector.lhs
--- a/src/Numeric/Clifford/Multivector.lhs
+++ b/src/Numeric/Clifford/Multivector.lhs
@@ -21,16 +21,11 @@
 {-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-}
 {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs#-}
 {-# LANGUAGE FlexibleInstances, StandaloneDeriving, KindSignatures, DataKinds #-}
-{-# LANGUAGE TemplateHaskell, TypeOperators, DeriveFunctor #-}
+{-# LANGUAGE TemplateHaskell, TypeOperators, DeriveFunctor, DeriveFoldable, DeriveTraversable#-}
 {-# LANGUAGE MultiParamTypeClasses, UndecidableInstances #-}
-\end{code}
-%if False
-\begin{code}
-{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}
--- OPTIONS_GHC -Odph -fvectorise -package dph-lifted-vseg 
---  LANGUAGE ParallelArrays
+{-# OPTIONS_HADDOCK show-extensions #-}
 \end{code}
-%endif
+
 Clifford algebras are a module over a ring. They also support all the usual transcendental functions.
 \begin{code}
 module Numeric.Clifford.Multivector where
@@ -77,7 +72,8 @@
 import Data.Word
 import Control.Applicative ((<$>))
 import Numeric.Clifford.Internal
-
+import Data.Traversable
+import Data.Foldable (Foldable)
 
 \end{code}
 
@@ -86,8 +82,9 @@
 
 \begin{code}
 data Multivector (p::Nat) (q::Nat) f where
-    BladeSum :: forall (p::Nat) (q::Nat) f . (Ord f, Algebra.Field.C f,SingI p, SingI q) => { _terms :: [Blade p q f]} -> Multivector p q f
+    BladeSum :: ∀ (p::Nat) (q::Nat) f . (Ord f, Algebra.Field.C f,SingI p, SingI q) => { _terms :: [Blade p q f]} → Multivector p q f
 
+
 type STVector = Multivector 3 1 Double
 type E3Vector = Multivector 3 0 Double
 
@@ -98,12 +95,10 @@
        d = fromIntegral (p' + q')
 
 deriving instance Eq (Multivector p q f)
---instance  (SingI p, SingI q) => Functor (Multivector p q) where
---    fmap func x =  func x--((terms x) & scale %~ func)
 deriving instance Ord (Multivector p q f)
 deriving instance (Show f) => Show (Multivector p q f)
---deriving instance (Read f) => Read (Multivector p q f)
 
+
 signature :: forall (p::Nat) (q::Nat) f. (SingI p, SingI q) => Multivector p q f ->  (Natural,Natural)
 signature _ = (toNatural  ((fromIntegral $ fromSing (sing :: Sing p))::Word),toNatural  ((fromIntegral $ fromSing (sing :: Sing q))::Word))
 
@@ -118,10 +113,14 @@
 terms = lens _terms (\bladeSum v -> bladeSum {_terms = v})
 
 {-# INLINE mvNormalForm #-}
-mvNormalForm (BladeSum terms) = BladeSum $ if null resultant then [scalarBlade Algebra.Additive.zero] else resultant  where
+mvNormalForm (BladeSum terms) = BladeSum $ mvNormalForm' terms
+
+
+mvNormalForm' terms =  if null resultant then [scalarBlade Algebra.Additive.zero] else resultant  where
     resultant = filter bladeNonZero $ addLikeTerms' $ Data.List.Ordered.sortBy compare $  map bladeNormalForm $ terms
+
 {-#INLINE mvTerms #-}
-mvTerms m = m^.terms
+mvTerms m = _terms m
 
 {-# INLINE addLikeTerms' #-}
 addLikeTerms' = sumLikeTerms . groupLikeTerms
@@ -133,7 +132,8 @@
 compareTol :: (Algebra.Algebraic.C f, Algebra.Absolute.C f, Ord f, SingI p, SingI q) => Multivector p q f -> Multivector p q f -> f -> Bool
 compareTol x y tol = ((NPN.abs $ magnitude (x-y) ) <= tol)
 
-{-#INLINE compensatedSum' #-}
+{-#NOINLINE compensatedSum' #-}
+{-#SPECIALISE compensatedSum' :: [Double] -> Double #-}
 compensatedSum' :: (Algebra.Additive.C f) => [f] -> f
 compensatedSum' xs = kahan zero zero xs where
     kahan s _ [] = s
@@ -153,20 +153,6 @@
         t = s + y
         newc = (t - s) - y
             
---multiplyAdd a b c = a*b + c
---twoProduct a b = (x,y) where
---    x = a*b
---z    y = multiplyAdd a b (negate x)
-
---multiplyList [] = []
---multiplyList a@(x:[])=a
---multiplyList (a:b:xs) = loop a (b:xs) zero where
---  loop pm [] ei = pm+ei
---  loop pm1 (ai:remaining) eim1= loop pi remaining ei where
---      (pi, pii) = twoProduct pm1 ai
---      ei = multiplyAdd eim1 ai pii
-
-
 multiplyOutBlades :: (SingI p, SingI q, Algebra.Ring.C a) => [Blade p q a] -> [Blade p q a] -> [Blade p q a]
 multiplyOutBlades x y = [bladeMul l r | l <-x, r <- y]
 
@@ -179,11 +165,9 @@
 multiplyList1 l = mvNormalForm $ BladeSum listOfBlades where
     expandedBlades a = foldl1 multiplyOutBlades a
     listOfBlades = expandedBlades $ map mvTerms l
---things to test: is 1. adding blades into a map based on indices 2. adding errything together 3. sort results quicker than
---                   1. sorting by indices 2. groupBy-ing on indices 3. adding the lists of identical indices
 
 {-#INLINE sumList #-}
-sumList xs = mvNormalForm $ BladeSum $ concat $ map mvTerms xs
+sumList xs = BladeSum $ mvNormalForm' $ concatMap _terms xs 
 
 {-#INLINE sumLikeTerms #-}
 {-#SPECIALISE INLINE sumLikeTerms :: [[STBlade]] -> [STBlade] #-}
@@ -260,7 +244,6 @@
 
     a ^ 2 = a * a
     a ^ 0 = one
-    a ^ 1 = a
     --a ^ n  --n < 0 = Clifford.recip $ a ^ (negate n)
     a ^ n  =  multiplyList (replicate (NPN.fromInteger n) a)
 
@@ -293,7 +276,6 @@
     signum (BladeSum []) = scalar Algebra.Additive.zero
 
 instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Module.C f (Multivector p q f) where
---    (*>) zero v = Algebra.Additive.zero
     {-#INLINE (*>) #-}
     {-#SPECIALISE INLINE (*>) :: Double -> STVector -> STVector #-}
     {-#SPECIALISE INLINE (*>) :: Double -> E3Vector -> E3Vector #-}
@@ -315,8 +297,8 @@
 scaleRight v s = BladeSum $ map (bladeScaleRight s) $ mvTerms v
 {-#INLINE divideRight #-}
 divideRight v s = scaleRight v (recip s)
---integratePoly c x = c : zipWith (Numeric.Clifford.Multivector./) x progression
 
+
 {-# INLINE converge#-}
 converge [] = error "converge: empty list"
 converge xs = fromMaybe empty (convergeBy checkPeriodic Just xs) 
@@ -539,23 +521,6 @@
 normalised a = (a `scaleRight` ( recip $ mag),mag) where
     mag = magnitude a
 
-
-\end{code}
-
-Now let's do (slow as fuck probably) numerical integration! :D~! Since this is gonna be used for physical applications, it's we're gonna start off with a Hamiltonian structure and then a symplectic integrator.
-
-\begin{code}
-
-
-
-{- $(derive makeSerialize ''Blade)
-$(derive makeSerialize ''Multivector)
-$(derive makeData ''Blade)
-$(derive makeTypeable ''Blade)
-$(derive makeData ''Multivector)
-$(derive makeTypeable ''Multivector)-}
-
--- $(derive makeArbitrary ''Multivector)
 
 \end{code}
 \bibliographystyle{IEEEtran}
diff --git a/src/Numeric/Clifford/NumericIntegration.lhs b/src/Numeric/Clifford/NumericIntegration.lhs
--- a/src/Numeric/Clifford/NumericIntegration.lhs
+++ b/src/Numeric/Clifford/NumericIntegration.lhs
@@ -18,17 +18,7 @@
 {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs, DataKinds, KindSignatures #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE TemplateHaskell #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-\end{code}
-%if False
-\begin{code}
-{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}
--- OPTIONS_GHC -Odph -fvectorise -package dph-lifted-vseg 
---  LANGUAGE ParallelArrays
-\end{code}
-%endif
-
-\begin{code}
+{-# LANGUAGE MultiParamTypeClasses, BangPatterns #-}
 
 module Numeric.Clifford.NumericIntegration where
 import Numeric.Clifford.Multivector as MV
@@ -57,6 +47,7 @@
 import Numeric.Clifford.Internal
 import Data.DeriveTH
 
+import qualified Debug.Trace
 
 elementAdd = zipWith (+)
 elementScale = zipWith (*>) 
@@ -123,7 +114,7 @@
 type AdaptiveStepSizeFunction f state = f -> state -> f 
 
 data RKAttribute f state = Explicit
-                 | HamiltonianFunction
+                 | HamiltonianFunction {totalEnergy :: state -> f}
                  | AdaptiveStepSize {sigma :: AdaptiveStepSizeFunction f state}
                  | ConvergenceTolerance {epsilon :: f}
                  | ConvergenceFunction {converger :: ConvergerFunction f } 
@@ -172,7 +163,7 @@
                         Just (AdaptiveStepSize sigma) -> sigma
                         Nothing -> (\_ _ -> one)
 
-    {-#INLINE rkMethodImplicitFixedPoint#-}
+
 --    {-#SPECIALISE rkMethodImplicitFixedPoint :: RKStepper 3 0 Double stateType #-}
 --    {-#SPECIALISE rkMethodImplicitFixedPoint :: RKStepper 3 0 Double [E3Vector] #-}
 --    {-#SPECIALISE rkMethodImplicitFixedPoint :: RKStepper 3 1 Double stateType #-}
@@ -192,7 +183,8 @@
             guessTime :: t
             guessTime = time + h'
             zkp1 :: NPN.Int -> [Multivector p q t] -> [Multivector p q t]
-            zkp1 i zk = map (h*>) (sumOfJs i zk) where
+            zkp1 i !zk = map (h*>) (sumOfJs i zk) where
+                
                 {-#INLINE sumOfJs#-}
                 sumOfJs :: Int -> [Multivector p q t] -> [Multivector p q t]
                 sumOfJs i zk =  sumListOfLists $ map (scaledByAij zk) (a i) where 
diff --git a/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs b/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs
--- a/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs
+++ b/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs
@@ -36,9 +36,9 @@
 import           Test.QuickCheck
 --trace _ a = a
 
-gaussLegendreFourthOrderTableau = ButcherTableau [[0.25::NPN.Double, 0.25 - (1.0 NPN./6)* sqrt 3], [0.25 + (1.0 NPN./ 6) * sqrt 3, 0.25]] [0.5, 0.5] [0.5 - (1.0 NPN./6)* sqrt 3, 0.5 + (1.0 NPN./ 6)* sqrt 3]
+gaussLegendreFourthOrderTableau = ButcherTableau [[0.25::NPN.Double, 0.25 - ((sqrt 3.0) NPN./6.0)],[0.25 + ((sqrt 3.0) NPN./ 6.0) , 0.25]] [0.5, 0.5] [0.5 - ((sqrt 3.0) NPN./6.0), 0.5 + ((sqrt 3.0) NPN./ 6.0)]
 gaussLegendreFourthOrder h f (t, state) = impl h f id id (t,state) where
-    impl= genericRKMethod gaussLegendreFourthOrderTableau []
+    impl= genericRKMethod gaussLegendreFourthOrderTableau [ConvergenceTolerance 1.0e-8]
 
 
 rk4ClassicalFromTableau h f (t,state) = impl h f id id (t, state) where