diff --git a/Data/Manifold/Atlas.hs b/Data/Manifold/Atlas.hs
--- a/Data/Manifold/Atlas.hs
+++ b/Data/Manifold/Atlas.hs
@@ -73,20 +73,20 @@
   lookupAtlas = id
 instance Atlas S¹ where
   type ChartIndex S¹ = S⁰
-  chartReferencePoint NegativeHalfSphere = S¹ $ -pi/2
-  chartReferencePoint PositiveHalfSphere = S¹ $ pi/2
-  interiorChartReferencePoint _ NegativeHalfSphere = S¹ $ -pi/2
-  interiorChartReferencePoint _ PositiveHalfSphere = S¹ $ pi/2
-  lookupAtlas (S¹ φ) | φ<0        = NegativeHalfSphere
+  chartReferencePoint NegativeHalfSphere = S¹Polar $ -pi/2
+  chartReferencePoint PositiveHalfSphere = S¹Polar $ pi/2
+  interiorChartReferencePoint _ NegativeHalfSphere = S¹Polar $ -pi/2
+  interiorChartReferencePoint _ PositiveHalfSphere = S¹Polar $ pi/2
+  lookupAtlas (S¹Polar φ) | φ<0        = NegativeHalfSphere
                      | otherwise  = PositiveHalfSphere
 instance Atlas S² where
   type ChartIndex S² = S⁰
-  chartReferencePoint PositiveHalfSphere = S² 0 0
-  chartReferencePoint NegativeHalfSphere = S² pi 0
-  interiorChartReferencePoint _ PositiveHalfSphere = S² 0 0
-  interiorChartReferencePoint _ NegativeHalfSphere = S² pi 0
-  lookupAtlas (S² ϑ _) | ϑ<pi/2     = PositiveHalfSphere
-                       | otherwise  = NegativeHalfSphere
+  chartReferencePoint PositiveHalfSphere = S²Polar 0 0
+  chartReferencePoint NegativeHalfSphere = S²Polar pi 0
+  interiorChartReferencePoint _ PositiveHalfSphere = S²Polar 0 0
+  interiorChartReferencePoint _ NegativeHalfSphere = S²Polar pi 0
+  lookupAtlas (S²Polar ϑ _) | ϑ<pi/2     = PositiveHalfSphere
+                            | otherwise  = NegativeHalfSphere
 
 instance (LinearSpace (a n), Needle (a n) ~ a n, Interior (a n) ~ a n)
               => Atlas (LinAff.Point a n) where
diff --git a/Data/Manifold/DifferentialEquation.hs b/Data/Manifold/DifferentialEquation.hs
--- a/Data/Manifold/DifferentialEquation.hs
+++ b/Data/Manifold/DifferentialEquation.hs
@@ -92,7 +92,7 @@
 --   be an arbitrary one-dimensional space (i.e. basically real intervals or 'S¹').
 --   In these cases, there is always only one partial derivative: that which we
 --   integrate over, in the only possible direction for propagation.
-type ODE x y = DifferentialEqn AffineModel x y
+type ODE x y = DifferentialEqn QuadraticModel x y
 
 constLinearDEqn :: ∀ x y . ( SimpleSpace x
                            , SimpleSpace y, AffineManifold y
@@ -132,7 +132,7 @@
     ,LinearManifoldWitness BoundarylessWitness, DualSpaceWitness ) -> \bwt' ->
     let bwt'inv = pseudoInverse bwt'
     in \(Shade (_x,y) δxy) -> LocalDifferentialEqn
-            (\(AffineModel shy' _) ->
+            (\(QuadraticModel shy' _ _) ->
                     let shy = dualShade shy'
                     in ( return $ shy & shadeNarrowness %~ scaleNorm 0.01
                        , return $ projectShade (Embedding (arr bwt')
diff --git a/Data/Manifold/FibreBundle.hs b/Data/Manifold/FibreBundle.hs
new file mode 100644
--- /dev/null
+++ b/Data/Manifold/FibreBundle.hs
@@ -0,0 +1,360 @@
+-- |
+-- Module      : Data.Manifold.FibreBundle
+-- Copyright   : (c) Justus Sagemüller 2018
+-- License     : GPL v3
+-- 
+-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de
+-- Stability   : experimental
+-- Portability : portable
+-- 
+
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE MultiParamTypeClasses      #-}
+{-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE UndecidableInstances       #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE UnicodeSyntax              #-}
+{-# LANGUAGE GADTs                      #-}
+{-# LANGUAGE DefaultSignatures          #-}
+{-# LANGUAGE CPP                        #-}
+#if __GLASGOW_HASKELL__ >= 800
+{-# LANGUAGE UndecidableSuperClasses    #-}
+#endif
+
+
+module Data.Manifold.FibreBundle where
+
+
+import Data.AdditiveGroup
+import Data.VectorSpace
+import Math.LinearMap.Category
+
+import Data.Manifold.Types.Primitive
+import Data.Manifold.PseudoAffine
+    
+import qualified Prelude as Hask
+
+import Control.Category.Constrained.Prelude hiding ((^))
+import Control.Category.Discrete
+import Control.Arrow.Constrained
+
+import Linear.V2 (V2(V2))
+import Linear.V3 (V3(V3))
+
+import Data.Tagged
+
+
+data TransportOnNeedleWitness k m f where
+  TransportOnNeedle :: (ParallelTransporting (LinearFunction (Scalar (Needle m)))
+                                             (Needle m) (Needle f))
+                     => TransportOnNeedleWitness k m f
+
+data ForgetTransportProperties k m f where
+  ForgetTransportProperties :: ParallelTransporting (->) m f
+                     => ForgetTransportProperties k m f
+
+class (PseudoAffine m, m ~ Interior m, Category k, Object k f)
+           => ParallelTransporting k m f where
+  transportOnNeedleWitness :: TransportOnNeedleWitness k m f
+  default transportOnNeedleWitness
+      :: ParallelTransporting (LinearFunction (Scalar (Needle m))) (Needle m) (Needle f)
+           => TransportOnNeedleWitness k m f
+  transportOnNeedleWitness = TransportOnNeedle
+  forgetTransportProperties :: ForgetTransportProperties k m f
+  default forgetTransportProperties :: ParallelTransporting (->) m f
+           => ForgetTransportProperties k m f
+  forgetTransportProperties = ForgetTransportProperties
+  
+  parallelTransport :: m -> Needle m -> k f f
+  translateAndInvblyParTransport
+        :: m -> Needle m -> (m, (k f f, k f f))
+  translateAndInvblyParTransport p v
+              = (q, ( parallelTransport p v
+                    , parallelTransport q $ p.-~!q ))
+   where q = p.+~^v
+
+instance ∀ m s . (PseudoAffine m, m ~ Interior m, s ~ (Scalar (Needle m)), Num' s)
+      => ParallelTransporting Discrete m (ZeroDim s) where
+  transportOnNeedleWitness = case (pseudoAffineWitness :: PseudoAffineWitness m) of
+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) -> TransportOnNeedle
+  forgetTransportProperties = case (pseudoAffineWitness :: PseudoAffineWitness m) of
+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))
+        -> ForgetTransportProperties
+  parallelTransport _ _ = id
+instance ∀ m s . (PseudoAffine m, m ~ Interior m, s ~ (Scalar (Needle m)), Num' s)
+      => ParallelTransporting (LinearFunction s) m (ZeroDim s) where
+  transportOnNeedleWitness = case (pseudoAffineWitness :: PseudoAffineWitness m) of
+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) -> TransportOnNeedle
+  forgetTransportProperties = case (pseudoAffineWitness :: PseudoAffineWitness m) of
+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))
+        -> ForgetTransportProperties
+  parallelTransport _ _ = id
+instance ∀ m s . (PseudoAffine m, m ~ Interior m, s ~ (Scalar (Needle m)), Num' s)
+      => ParallelTransporting (->) m (ZeroDim s) where
+  transportOnNeedleWitness = case (pseudoAffineWitness :: PseudoAffineWitness m) of
+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)) -> TransportOnNeedle
+  forgetTransportProperties = case (pseudoAffineWitness :: PseudoAffineWitness m) of
+    (PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness))
+        -> ForgetTransportProperties
+  parallelTransport _ _ = id
+
+instance (Category k, Object k ℝ) => ParallelTransporting k ℝ ℝ where
+  parallelTransport _ _ = id
+instance (Category k, Object k ℝ²) => ParallelTransporting k ℝ² ℝ² where
+  parallelTransport _ _ = id
+instance (Category k, Object k ℝ³) => ParallelTransporting k ℝ³ ℝ³ where
+  parallelTransport _ _ = id
+instance (Category k, Object k ℝ⁴) => ParallelTransporting k ℝ⁴ ℝ⁴ where
+  parallelTransport _ _ = id
+
+instance (Category k, Object k ℝ) => ParallelTransporting k S¹ ℝ where
+  parallelTransport _ _ = id
+
+instance (EnhancedCat k (LinearMap ℝ), Object k ℝ²)
+             => ParallelTransporting k S² ℝ² where
+  parallelTransport p v = (fst . snd) (translateAndInvblyParTransport p v)
+  translateAndInvblyParTransport (S²Polar θ₀ φ₀) 𝐯
+     | d < pi     = (S²Polar θ₁ φ₁, (arr fwd, arr bwd))
+     | d < 2*pi   = translateAndInvblyParTransport (S²Polar θ₀ φ₀)
+                      $ 𝐯^*(-(2*pi-d)/d)
+     | otherwise  = translateAndInvblyParTransport (S²Polar θ₀ φ₀)
+                      $ let revolutions = floor $ d/(2*pi)
+                        in 𝐯^*((d - 2*pi*fromIntegral revolutions)/d)
+   where -- See images/constructions/sphericoords-needles.svg. Translation as in
+         -- "Data.Manifold.PseudoAffine" instance.
+         S¹Polar γc₀ = coEmbed 𝐯
+         γ₀ | θ₀ < pi/2   = γc₀ - φ₀
+            | otherwise   = γc₀ + φ₀
+         d = magnitude 𝐯
+         S¹Polar φ₁ = S¹Polar φ₀ .+~^ δφ
+         
+         -- Cartesian coordinates of p₁ in the system whose north pole is p₀
+         -- with φ₀ as the zero meridian
+         V3 bx by bz = embed $ S²Polar d γ₀
+         
+         sθ₀ = sin θ₀; cθ₀ = cos θ₀
+         -- Cartesian coordinates of p₁ in the system with the standard north pole,
+         -- but still φ₀ as the zero meridian
+         (qx,qz) = ( cθ₀ * bx + sθ₀ * bz
+                   ,-sθ₀ * bx + cθ₀ * bz )
+         qy      = by
+         
+         S²Polar θ₁ δφ = coEmbed $ V3 qx qy qz
+         
+         sθ₁ = sin θ₁; cθ₁ = cos θ₁
+         
+         γ₁
+          | sθ₀<=sθ₁  = let
+              -- Cartesian coordinates of the standard north pole in the system whose north
+              -- pole is p₀ with 𝐯 along the zero meridian
+              V3 nbx nby nbz = embed $ S²Polar θ₀ (pi-γ₀)
+              
+              sd = sin d; cd = cos d
+              -- Cartesian coordinates of the standard north pole in the system whose north
+              -- pole is p₁ with 𝐯 along the zero meridian
+              (ox,oz) = ( cd * nbx - sd * nbz
+                        , sd * nbx + cd * nbz )
+              oy      = nby
+
+           in atan2 oy (-ox)
+
+          | otherwise = let
+              -- Cartesian coordinates of p₀ in the system with the standard north pole,
+              -- with p₁ on the zero meridian
+              V3 gx gy gz = embed $ S²Polar θ₀ (-δφ)
+              
+              -- Cartesian coordinates of p₀ in the system whose north
+              -- pole is p₁ and the standard north pole on the zero meridian
+              (ux,uz) = ( cθ₁ * gx - sθ₁ * gz
+                        , sθ₁ * gx + cθ₁ * gz )
+              uy      = gy
+
+           in atan2 (-uy) (-ux)
+
+         γc₁ | θ₁ < pi/2  = γ₁ + φ₁
+             | otherwise  = γ₁ - φ₁
+
+         (sδγc, cδγc) = sin &&& cos $ γc₁ - γc₀
+
+         fwd = LinearMap (V2 (V2   cδγc  sδγc)
+                             (V2 (-sδγc) cδγc)) :: LinearMap ℝ ℝ² ℝ²
+         bwd = LinearMap (V2 (V2 cδγc (-sδγc))
+                             (V2 sδγc   cδγc )) :: LinearMap ℝ ℝ² ℝ²
+
+
+instance {-# OVERLAPS #-} ∀ k a b fa fb s .
+         ( ParallelTransporting k a fa, ParallelTransporting k b fb
+         , PseudoAffine fa, PseudoAffine fb
+         , Scalar (Needle a) ~ s, Scalar (Needle b) ~ s
+         , Scalar (Needle fa) ~ s, Scalar (Needle fb) ~ s
+         , Num' s
+         , Morphism k, ObjectPair k fa fb )
+              => ParallelTransporting k (a,b) (fa,fb) where
+  transportOnNeedleWitness = case
+         ( pseudoAffineWitness :: PseudoAffineWitness a
+         , pseudoAffineWitness :: PseudoAffineWitness b
+         , pseudoAffineWitness :: PseudoAffineWitness fa
+         , pseudoAffineWitness :: PseudoAffineWitness fb
+         , transportOnNeedleWitness :: TransportOnNeedleWitness k a fa
+         , transportOnNeedleWitness :: TransportOnNeedleWitness k b fb ) of
+     ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)
+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)
+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)
+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)
+      ,TransportOnNeedle, TransportOnNeedle)
+         -> TransportOnNeedle
+  forgetTransportProperties = case
+    ( forgetTransportProperties :: ForgetTransportProperties k a fa
+    , forgetTransportProperties :: ForgetTransportProperties k b fb ) of
+     (ForgetTransportProperties, ForgetTransportProperties) -> ForgetTransportProperties
+  parallelTransport (pa,pb) (va,vb)
+       = parallelTransport pa va  *** parallelTransport pb vb
+
+instance ∀ k a f g s .
+         ( ParallelTransporting k a f, ParallelTransporting k a g
+         , ParallelTransporting (LinearFunction s) (Needle a) (Needle f, Needle g)
+         , PseudoAffine f, PseudoAffine g
+         , Morphism k, ObjectPair k f g )
+              => ParallelTransporting k a (f,g) where
+  transportOnNeedleWitness = case
+         ( pseudoAffineWitness :: PseudoAffineWitness a
+         , pseudoAffineWitness :: PseudoAffineWitness f
+         , pseudoAffineWitness :: PseudoAffineWitness g
+         , transportOnNeedleWitness :: TransportOnNeedleWitness k a f
+         , transportOnNeedleWitness :: TransportOnNeedleWitness k a g ) of
+     ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)
+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)
+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)
+      ,TransportOnNeedle, TransportOnNeedle)
+         -> TransportOnNeedle
+  forgetTransportProperties = case
+    ( forgetTransportProperties :: ForgetTransportProperties k a f
+    , forgetTransportProperties :: ForgetTransportProperties k a g ) of
+     (ForgetTransportProperties, ForgetTransportProperties) -> ForgetTransportProperties
+  parallelTransport p v
+       = parallelTransport p v *** parallelTransport p v
+
+
+instance ( ParallelTransporting (LinearFunction (Scalar f)) m f, AdditiveGroup m
+         , VectorSpace f )
+                => AdditiveGroup (FibreBundle m f) where
+  zeroV = FibreBundle zeroV zeroV
+  FibreBundle p v ^+^ FibreBundle q w = FibreBundle (p^+^q) (v^+^w)
+  negateV (FibreBundle p v) = FibreBundle (negateV p) (negateV v)
+
+instance ∀ m f s .
+         ( ParallelTransporting (->) m (Interior f), Semimanifold f
+         , ParallelTransporting (LinearFunction s) (Needle m) (Needle f)
+         , s ~ Scalar (Needle m) )
+                => Semimanifold (FibreBundle m f) where
+  type Interior (FibreBundle m f) = FibreBundle m (Interior f)
+  type Needle (FibreBundle m f) = FibreBundle (Needle m) (Needle f)
+  toInterior (FibreBundle p f) = FibreBundle p <$> toInterior f
+  translateP = Tagged $ case ( translateP :: Tagged m (Interior m -> Needle m -> Interior m)
+                             , semimanifoldWitness :: SemimanifoldWitness f) of
+      (Tagged tpm, SemimanifoldWitness BoundarylessWitness)
+           -> \(FibreBundle p f) (FibreBundle v δf)
+                   -> FibreBundle (tpm p v) (parallelTransport p v f.+~^δf)
+  semimanifoldWitness = case ( semimanifoldWitness :: SemimanifoldWitness m
+                             , semimanifoldWitness :: SemimanifoldWitness f
+                             , forgetTransportProperties
+                               :: ForgetTransportProperties (LinearFunction s) (Needle m) (Needle f)
+                             ) of
+         (SemimanifoldWitness BoundarylessWitness, SemimanifoldWitness BoundarylessWitness
+          ,ForgetTransportProperties)
+           -> SemimanifoldWitness BoundarylessWitness
+  FibreBundle p f .+~^ FibreBundle v δf
+      = FibreBundle (p.+~^v) (parallelTransport p v f.+~^δf)
+
+instance ∀ m f s .
+         ( ParallelTransporting (->) m f, ParallelTransporting (->) m (Interior f)
+         , PseudoAffine f
+         , ParallelTransporting (LinearFunction s) (Needle m) (Needle f)
+         , s ~ Scalar (Needle m) )
+                => PseudoAffine (FibreBundle m f) where
+  pseudoAffineWitness = case ( pseudoAffineWitness :: PseudoAffineWitness m
+                             , pseudoAffineWitness :: PseudoAffineWitness f
+                             , forgetTransportProperties
+                               :: ForgetTransportProperties (LinearFunction s) (Needle m) (Needle f)
+                             ) of
+     ( PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)
+      ,PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)
+      ,ForgetTransportProperties)
+         -> PseudoAffineWitness (SemimanifoldWitness BoundarylessWitness)
+  FibreBundle p f .-~. FibreBundle q g = case p.-~.q of
+      Nothing -> Nothing
+      Just v  -> FibreBundle v <$> f .-~. parallelTransport p v g
+
+
+instance (AdditiveGroup f, x ~ Interior x) => NaturallyEmbedded x (FibreBundle x f) where
+  embed x = FibreBundle x zeroV
+  coEmbed (FibreBundle x _) = x
+
+instance (NaturallyEmbedded (Interior m) (Interior v), VectorSpace f)
+    => NaturallyEmbedded (FibreBundle m ℝ⁰) (FibreBundle v f) where
+  embed (FibreBundle x Origin) = FibreBundle (embed x) zeroV
+  coEmbed (FibreBundle u _) = FibreBundle (coEmbed u) Origin
+
+instance (AdditiveGroup (Interior y), AdditiveGroup g)
+           => NaturallyEmbedded (FibreBundle x f) (FibreBundle (x,y) (f,g)) where
+  embed (FibreBundle x δx) = FibreBundle (x,zeroV) (δx,zeroV)
+  coEmbed (FibreBundle (x,_) (δx,_)) = FibreBundle x δx
+
+instance NaturallyEmbedded v w
+      => NaturallyEmbedded (FibreBundle ℝ v) (FibreBundle ℝ w) where
+  embed (FibreBundle p v) = FibreBundle p $ embed v
+  coEmbed (FibreBundle p w) = FibreBundle p $ coEmbed w
+instance NaturallyEmbedded v w
+      => NaturallyEmbedded (FibreBundle ℝ² v) (FibreBundle ℝ² w) where
+  embed (FibreBundle p v) = FibreBundle p $ embed v
+  coEmbed (FibreBundle p w) = FibreBundle p $ coEmbed w
+instance NaturallyEmbedded v w
+      => NaturallyEmbedded (FibreBundle ℝ³ v) (FibreBundle ℝ³ w) where
+  embed (FibreBundle p v) = FibreBundle p $ embed v
+  coEmbed (FibreBundle p w) = FibreBundle p $ coEmbed w
+instance NaturallyEmbedded v w
+      => NaturallyEmbedded (FibreBundle ℝ⁴ v) (FibreBundle ℝ⁴ w) where
+  embed (FibreBundle p v) = FibreBundle p $ embed v
+  coEmbed (FibreBundle p w) = FibreBundle p $ coEmbed w
+
+instance NaturallyEmbedded (FibreBundle S¹ ℝ) (FibreBundle ℝ² ℝ²) where
+  embed (FibreBundle (S¹Polar φ) l) = FibreBundle (V2 cφ sφ) $ l*^(V2 (-sφ) cφ)
+   where (cφ, sφ) = (cos &&& sin) φ
+  coEmbed (FibreBundle (V2 0 0) (V2 _ δy)) = FibreBundle (S¹Polar 0) δy
+  coEmbed (FibreBundle p (V2 δx δy)) = FibreBundle (S¹Polar $ atan2 sφ cφ) $ cφ*δy - sφ*δx
+   where V2 cφ sφ = p^/r
+         r = magnitude p
+
+instance NaturallyEmbedded (FibreBundle S² ℝ²) (FibreBundle ℝ³ ℝ³) where
+  embed (FibreBundle (S²Polar θ φ) 𝐯@(V2 δξ δυ))
+       = FibreBundle (V3 (sθ*cφ) (sθ*sφ) cθ) 𝐯r
+   where [V2 cθ sθ, V2 cφ sφ] = embed . S¹Polar <$> [θ,φ]
+         S¹Polar γc = coEmbed 𝐯
+         γ | θ < pi/2   = γc - φ
+           | otherwise  = γc + φ
+         d = magnitude 𝐯
+
+         V2 δθ δφ = d *^ embed (S¹Polar γ)
+         
+         𝐞φ = V3 (-sφ) cφ 0
+         𝐞θ = V3 (cθ*cφ) (cθ*sφ) (-sθ)
+         𝐯r = δθ*^𝐞θ ^+^ δφ*^𝐞φ
+  
+  coEmbed (FibreBundle (V3 x y z) 𝐯r)
+           = FibreBundle (S²Polar θ φ) (magnitude (δθ,δφ) *^ embed (S¹Polar γc))
+   where r = sqrt $ x^2 + y^2 + z^2
+         rxy = sqrt $ x^2 + y^2
+         θ = atan2 rxy z
+         φ = atan2 y x
+         cθ = z / r
+         sθ = rxy / r
+         (cφ,sφ) | rxy>0      = (x,y)^/rxy
+                 | otherwise  = (1,0)
+         𝐞φ = V3 (-sφ) cφ 0
+         𝐞θ = V3 (cθ*cφ) (cθ*sφ) (-sθ)
+         δθ = 𝐞θ <.> 𝐯r
+         δφ = 𝐞φ <.> 𝐯r
+         γ = atan2 δφ δθ
+         γc | θ < pi/2   = γ + φ
+            | otherwise  = γ - φ
diff --git a/Data/Manifold/Function/Interpolation.hs b/Data/Manifold/Function/Interpolation.hs
new file mode 100644
--- /dev/null
+++ b/Data/Manifold/Function/Interpolation.hs
@@ -0,0 +1,106 @@
+-- |
+-- Module      : Data.Manifold.Function.Interpolation
+-- Copyright   : (c) Justus Sagemüller 2017
+-- License     : GPL v3
+-- 
+-- Maintainer  : (@) jsagemue $ uni-koeln.de
+-- Stability   : experimental
+-- Portability : portable
+-- 
+
+{-# LANGUAGE ScopedTypeVariables      #-}
+{-# LANGUAGE UnicodeSyntax            #-}
+{-# LANGUAGE TypeOperators            #-}
+{-# LANGUAGE TupleSections            #-}
+{-# LANGUAGE TypeFamilies             #-}
+{-# LANGUAGE UndecidableInstances     #-}
+{-# LANGUAGE FlexibleContexts         #-}
+{-# LANGUAGE StandaloneDeriving       #-}
+{-# LANGUAGE TemplateHaskell          #-}
+{-# LANGUAGE ConstraintKinds          #-}
+
+module Data.Manifold.Function.Interpolation (
+      InterpolationFunction
+    ) where
+
+
+import Data.Manifold.Types
+import Data.Manifold.Types.Primitive ((^))
+import Data.Manifold.PseudoAffine
+import Data.Manifold.Shade
+import Data.Manifold.Web
+import Data.Manifold.Web.Internal
+import Data.Manifold.Function.LocalModel
+
+import Data.VectorSpace
+import Math.LinearMap.Category
+
+import Data.Foldable (toList)
+import Data.List.NonEmpty (NonEmpty (..))
+import qualified Data.List.NonEmpty as NE
+
+import qualified Prelude as Hask
+
+import Control.Category.Constrained.Prelude
+import Control.Arrow.Constrained
+import Control.Monad.Constrained
+
+import Control.Lens
+import Control.Lens.TH
+
+
+newtype InterpolationFunction ㄇ x y = InterpolationFunction {
+      _interpWeb :: PointsWeb x (ㄇ x y)
+    }
+makeLenses ''InterpolationFunction
+
+
+fromPointsWeb :: (ModellableRelation x y, LocalModel ㄇ)
+                 => PointsWeb x (Shade' y) -> InterpolationFunction ㄇ x y
+fromPointsWeb = InterpolationFunction . localFmapWeb (
+                 \locInfo -> case fitLocally $
+                                    (zeroV, locInfo^.thisNodeData)
+                                  : [ (ngbx, ngb^.thisNodeData)
+                                    | (ngbx,ngb) <- concat $ localOnion locInfo []] of
+                                 Just locModl -> locModl )
+
+
+adjustMetricToModel :: ∀ x y ㄇ . (ModellableRelation x y, LocalModel ㄇ)
+                 => InterpolationFunction ㄇ x y -> InterpolationFunction ㄇ x y
+adjustMetricToModel = _interpWeb >>> webLocalInfo
+    >>> \(PointsWeb w) -> InterpolationFunction . PointsWeb $ fmap remetricise w
+ where remetricise :: Neighbourhood x (WebLocally x (ㄇ x y))
+             -> Neighbourhood x (ㄇ x y)
+       remetricise nd = nd & dataAtNode .~ localModel
+                           & localScalarProduct .~ newNorm
+        where localModel = nd^.dataAtNode.thisNodeData
+              newNorm = spanNorm
+                  [ dx ^/ ((0.1 + occlusion (ngb^.thisNodeData.tweakLocalOffset)
+                                            (fromInterior ySynth))
+                           * (dx<.>^δx))
+                  | (δx,ngb) <- concat . take 2 $ localOnion (nd^.dataAtNode) []
+                  , let dx = nd^.localScalarProduct<$|δx
+                        Shade' ySynth _ = evalLocalModel localModel δx ]
+                      :: Metric x
+
+
+upsampleAtLargeDist :: (ModellableRelation x y, LocalModel ㄇ)
+        => ℝ -> (x -> ㄇ x y -> Needle x -> Shade' y)
+            -> InterpolationFunction ㄇ x y -> PointsWeb x (Shade' y)
+upsampleAtLargeDist dmax gapFillFunc (InterpolationFunction web)
+     = fromWebNodes (\(Shade x _) -> case nearestNeighbour webI (fromInterior x) of
+                         Just (_,nearest) -> nearest ^. nodeLocalScalarProduct) $ do
+          local <- toList webI
+          (local^.thisNodeCoord, evalLocalModel (local^.thisNodeData) zeroV) : do 
+             (ngId, (δx, ngb)) <- local^.nodeNeighbours
+             guard (ngId > local^.thisNodeId
+                   && (local^.nodeLocalScalarProduct|$|δx) > dmax)
+             return ( local^.thisNodeCoord !+~^ δx^/2
+                    , gapFillFunc (local^.thisNodeCoord)
+                                  (local^.thisNodeData)
+                                  (δx^/2) )
+ where webI = webLocalInfo web
+
+autoUpsampleAtLargeDist :: (ModellableRelation x y, LocalModel ㄇ)
+        => ℝ -> InterpolationFunction ㄇ x y -> PointsWeb x (Shade' y)
+autoUpsampleAtLargeDist dmax = upsampleAtLargeDist dmax $ const evalLocalModel
diff --git a/Data/Manifold/Function/LocalModel.hs b/Data/Manifold/Function/LocalModel.hs
--- a/Data/Manifold/Function/LocalModel.hs
+++ b/Data/Manifold/Function/LocalModel.hs
@@ -260,6 +260,7 @@
                   => [(Needle x, Shade' y)] -> Maybe (ㄇ x y)
   tweakLocalOffset :: ModellableRelation x y
                   => Lens' (ㄇ x y) (Shade y)
+  evalLocalModel :: ModellableRelation x y => ㄇ x y -> Needle x -> Shade' y
 
 modelParametersOverdetMargin :: Int -> Int
 modelParametersOverdetMargin n = n + round (sqrt $ fromIntegral n) - 1
@@ -296,6 +297,12 @@
                      $ (p₀:|ps++[pω])
           | otherwise  = Nothing
   tweakLocalOffset = affineModelOffset
+  evalLocalModel = aEvL pseudoAffineWitness
+   where aEvL :: ∀ x y . ModellableRelation x y
+                => PseudoAffineWitness y -> AffineModel x y -> Needle x -> Shade' y
+         aEvL (PseudoAffineWitness (SemimanifoldWitness _)) (AffineModel shy₀ shj) δx
+          = convolveShade' (dualShade shy₀)
+                           (dualShade . linearProjectShade (lfun ($ δx)) $ shj)
 
 instance LocalModel QuadraticModel where
   fitLocally = qFitL
@@ -308,3 +315,13 @@
                      $ (p₀:|ps++[pω])
           | otherwise  = Nothing
   tweakLocalOffset = quadraticModelOffset
+  evalLocalModel = aEvL pseudoAffineWitness
+   where aEvL :: ∀ x y . ModellableRelation x y
+                => PseudoAffineWitness y -> QuadraticModel x y -> Needle x -> Shade' y
+         aEvL (PseudoAffineWitness (SemimanifoldWitness _))
+              (QuadraticModel shy₀ shj shjj) δx
+          = (dualShade shy₀)
+           `convolveShade'`
+            (dualShade . linearProjectShade (lfun ($ δx)) $ shj)
+           `convolveShade'`
+            (dualShade . linearProjectShade (lfun ($ squareV δx)) $ shjj)
diff --git a/Data/Manifold/PseudoAffine.hs b/Data/Manifold/PseudoAffine.hs
--- a/Data/Manifold/PseudoAffine.hs
+++ b/Data/Manifold/PseudoAffine.hs
@@ -55,7 +55,7 @@
             , PseudoAffine(..)
             -- * Type definitions
             -- ** Needles
-            , Local(..), (⊙+^)
+            , Local(..), (⊙+^), (!+~^)
             -- ** Metrics
             , Metric, Metric'
             , RieMetric, RieMetric'
@@ -107,6 +107,7 @@
 
 import Control.Lens (Lens', lens, (^.), (&), (%~), (.~))
 
+import Data.CallStack (HasCallStack)
 import GHC.Exts (Constraint)
 
 
@@ -385,46 +386,62 @@
   fromInterior = id
   toInterior = pure
   translateP = Tagged (.+~^)
-  S² ϑ₀ φ₀ .+~^ δv
-     | ϑ₀ < pi/2  = sphereFold PositiveHalfSphere $ ϑ₀*^embed(S¹ φ₀) ^+^ δv
-     | otherwise  = sphereFold NegativeHalfSphere $ (pi-ϑ₀)*^embed(S¹ φ₀) ^+^ δv
+  S²Polar θ₀ φ₀ .+~^ 𝐯 = S²Polar θ₁ φ₁
+   where -- See images/constructions/sphericoords-needles.svg.
+         S¹Polar γc = coEmbed 𝐯
+         γ | θ₀ < pi/2   = γc - φ₀
+           | otherwise   = γc + φ₀
+         d = magnitude 𝐯
+         S¹Polar φ₁ = S¹Polar φ₀ .+~^ δφ
+         
+         -- Cartesian coordinates of p₁ in the system whose north pole is p₀
+         -- with φ₀ as the zero meridian
+         V3 bx by bz = embed $ S²Polar d γ
+         
+         sθ₀ = sin θ₀; cθ₀ = cos θ₀
+         -- Cartesian coordinates of p₁ in the system with the standard north pole,
+         -- but still φ₀ as the zero meridian
+         (qx,qz) = ( cθ₀ * bx + sθ₀ * bz
+                   ,-sθ₀ * bx + cθ₀ * bz )
+         qy      = by
+         
+         S²Polar θ₁ δφ = coEmbed $ V3 qx qy qz
+
 instance PseudoAffine S² where
-  S² ϑ₁ φ₁ .-~. S² ϑ₀ φ₀
-     | ϑ₀ < pi/2  = pure ( ϑ₁*^embed(S¹ φ₁) ^-^ ϑ₀*^embed(S¹ φ₀) )
-     | otherwise  = pure ( (pi-ϑ₁)*^embed(S¹ φ₁) ^-^ (pi-ϑ₀)*^embed(S¹ φ₀) )
+  S²Polar θ₁ φ₁ .-~! S²Polar θ₀ φ₀ = d *^ embed(S¹Polar γc)
+   where -- See images/constructions/sphericoords-needles.svg.
+         V3 qx qy qz = embed $ S²Polar θ₁ (φ₁-φ₀)
 
-sphereFold :: S⁰ -> ℝ² -> S²
-sphereFold hfSphere v
-   | ϑ₀ > pi     = S² (inv $ tau - ϑ₀) (toS¹range $ φ₀+pi)
-   | otherwise  = S² (inv ϑ₀) φ₀
- where S¹ φ₀ = coEmbed v
-       ϑ₀ = magnitude v `mod'` tau
-       inv ϑ = case hfSphere of PositiveHalfSphere -> ϑ
-                                NegativeHalfSphere -> pi - ϑ
+         sθ₀ = sin θ₀; cθ₀ = cos θ₀
+         (bx,bz) = ( cθ₀ * qx - sθ₀ * qz
+                   , sθ₀ * qx + cθ₀ * qz )
+         by      = qy
 
+         S²Polar d γ = coEmbed $ V3 bx by bz
+         
+         γc | θ₀ < pi/2   = γ + φ₀
+            | otherwise   = γ - φ₀
 
+
+
+
 instance Semimanifold ℝP² where
   type Needle ℝP² = ℝ²
   fromInterior = id
   toInterior = pure
   translateP = Tagged (.+~^)
-  ℝP² r₀ φ₀ .+~^ V2 δr δφ
-   | r₀ > 1/2   = case r₀ + δr of
-                   r₁ | r₁ > 1     -> ℝP² (2-r₁) (toS¹range $ φ₀+δφ+pi)
-                      | otherwise  -> ℝP²    r₁  (toS¹range $ φ₀+δφ)
-  ℝP² r₀ φ₀ .+~^ δxy = let v = r₀*^embed(S¹ φ₀) ^+^ δxy
-                           S¹ φ₁ = coEmbed v
-                           r₁ = magnitude v `mod'` 1
-                       in ℝP² r₁ φ₁  
+  HemisphereℝP²Polar θ₀ φ₀ .+~^ v
+      = case S²Polar θ₀ φ₀ .+~^ v of
+          S²Polar θ₁ φ₁
+           | θ₁ > pi/2   -> HemisphereℝP²Polar (pi-θ₁) (-φ₁)
+           | otherwise   -> HemisphereℝP²Polar θ₁        φ₁
 instance PseudoAffine ℝP² where
-  ℝP² r₁ φ₁ .-~. ℝP² r₀ φ₀
-   | r₀ > 1/2   = pure `id` case φ₁-φ₀ of
-                          δφ | δφ > 3*pi/2  -> V2 (  r₁ - r₀) (δφ - 2*pi)
-                             | δφ < -3*pi/2 -> V2 (  r₁ - r₀) (δφ + 2*pi)
-                             | δφ > pi/2    -> V2 (2-r₁ - r₀) (δφ - pi  )
-                             | δφ < -pi/2   -> V2 (2-r₁ - r₀) (δφ + pi  )
-                             | otherwise    -> V2 (  r₁ - r₀) (δφ       )
-   | otherwise  = pure ( r₁*^embed(S¹ φ₁) ^-^ r₀*^embed(S¹ φ₀) )
+  HemisphereℝP²Polar θ₁ φ₁ .-~! HemisphereℝP²Polar θ₀ φ₀
+      = case S²Polar θ₁ φ₁ .-~! S²Polar θ₀ φ₀ of
+          v -> let r² = magnitudeSq v
+               in if r²>pi^2/4
+                   then S²Polar (pi-θ₁) (-φ₁) .-~! S²Polar θ₀ φ₀
+                   else v
 
 
 -- instance (PseudoAffine m, VectorSpace (Needle m), Scalar (Needle m) ~ ℝ)
@@ -460,6 +477,13 @@
 
 type DualNeedleWitness x = DualSpaceWitness (Needle x)
 
+
+
+infixl 6 !+~^
+-- | Boundary-unsafe version of `.+~^`.
+(!+~^) :: ∀ x . (Semimanifold x, HasCallStack) => x -> Needle x -> x
+p!+~^v = case toInterior p of
+           Just p' -> p'.+~^v
 
 
 infix 6 ⊙+^
diff --git a/Data/Manifold/Riemannian.hs b/Data/Manifold/Riemannian.hs
--- a/Data/Manifold/Riemannian.hs
+++ b/Data/Manifold/Riemannian.hs
@@ -164,16 +164,16 @@
   geodesicBetween _ _ = empty
 
 instance Geodesic S¹ where
-  geodesicBetween (S¹ φ) (S¹ ϕ)
-    | abs (φ-ϕ) < pi  = (>>> S¹) <$> geodesicBetween φ ϕ
-    | φ > 0           = (>>> S¹ . \ψ -> signum ψ*pi - ψ)
+  geodesicBetween (S¹Polar φ) (S¹Polar ϕ)
+    | abs (φ-ϕ) < pi  = (>>> S¹Polar) <$> geodesicBetween φ ϕ
+    | φ > 0           = (>>> S¹Polar . \ψ -> signum ψ*pi - ψ)
                         <$> geodesicBetween (pi-φ) (-ϕ-pi)
-    | otherwise       = (>>> S¹ . \ψ -> signum ψ*pi - ψ)
+    | otherwise       = (>>> S¹Polar . \ψ -> signum ψ*pi - ψ)
                         <$> geodesicBetween (-pi-φ) (pi-ϕ)
-  middleBetween (S¹ φ) (S¹ ϕ)
-    | abs (φ-ϕ) < pi  = S¹ <$> middleBetween φ ϕ
-    | φ > 0           = S¹ <$> middleBetween (pi-φ) (-ϕ-pi)
-    | otherwise       = S¹ <$> middleBetween (-pi-φ) (pi-ϕ)
+  middleBetween (S¹Polar φ) (S¹Polar ϕ)
+    | abs (φ-ϕ) < pi  = S¹Polar <$> middleBetween φ ϕ
+    | φ > 0           = S¹Polar <$> middleBetween (pi-φ) (-ϕ-pi)
+    | otherwise       = S¹Polar <$> middleBetween (-pi-φ) (pi-ϕ)
 
 
 -- instance Geodesic (Cℝay S⁰) where
diff --git a/Data/Manifold/TreeCover.hs b/Data/Manifold/TreeCover.hs
--- a/Data/Manifold/TreeCover.hs
+++ b/Data/Manifold/TreeCover.hs
@@ -365,26 +365,39 @@
     , (pre, (x,node):post) <- splitAt i lvs
               -> Right . fmap (PlainLeaves . (pre++) . (:post) . (x,)) $ f node
     | otherwise -> Left $ i-n
-treeLeaf i f (DisjointBranches n brs)
-    | i<n        = foldl (\case 
-                             Left i' -> (treeLeaf i' f)
-                             result  -> return result
-                         ) (Left i) brs
-    | otherwise  = Left $ i-n
-treeLeaf i f sh@(OverlappingBranches n _ brs)
-    | i<n        = foldl (\case 
-                             Left i' -> (treeLeaf i' f)
-                             result  -> return result
-                         ) (Left i) (toList brs>>=toList)
-    | otherwise  = Left $ i-n
+treeLeaf i f (DisjointBranches n _)
+    | i>=n   = Left $ i-n
+treeLeaf i f (DisjointBranches n (br:|[]))
+        = fmap (DisjointBranches n . pure) <$> treeLeaf i f br
+treeLeaf i f (DisjointBranches n (br:|br':brs))
+        = case treeLeaf i f br of
+            Left overshoot -> fmap (\(DisjointBranches _ (br'':|brs'))
+                                   -> DisjointBranches n (br:|br'':brs'))
+                  <$> treeLeaf overshoot f
+                     (DisjointBranches (n-nLeaves br) $ br':|brs)
+            Right done -> Right $ DisjointBranches n . (:|br':brs) <$> done
+treeLeaf i f (OverlappingBranches n extend (br@(DBranch dir (Hourglass t b)):|brs))
+    | i<nt       = fmap (OverlappingBranches n extend
+                         . (:|brs) . DBranch dir . (`Hourglass`b))
+                    <$> treeLeaf i f t
+    | i<nt+nb    = fmap (OverlappingBranches n extend
+                         . (:|brs) . DBranch dir . ( Hourglass t))
+                    <$> treeLeaf (i-nt) f b
+    | br':brs' <- brs
+                 = fmap (\(OverlappingBranches _ _ (br'':|brs''))
+                         -> OverlappingBranches n extend $ br:|br'':brs'')
+                    <$> treeLeaf (i-nt-nb) f (OverlappingBranches n extend $ br':|brs')
+    | otherwise  = Left $ i - nt - nb
+ where [nt,nb] = nLeaves<$>[t,b]
 
+
 -- | “Inverse indexing” of a tree. This is roughly a nearest-neighbour search,
 --   but not guaranteed to give the correct result unless evaluated at the
 --   precise position of a tree leaf.
 positionIndex :: ∀ x y . (WithField ℝ Manifold x, SimpleSpace (Needle x))
        => Maybe (Metric x)   -- ^ For deciding (at the lowest level) what “close” means;
                              --   this is optional for any tree of depth >1.
-        -> x`Shaded`y        -- ^ The tree to index into
+        -> (x`Shaded`y)      -- ^ The tree to index into
         -> x                 -- ^ Position to look up
         -> Maybe (Int, ([x`Shaded`y], (x,y)))
                    -- ^ Index of the leaf near to the query point, the “path” of
diff --git a/Data/Manifold/Types.hs b/Data/Manifold/Types.hs
--- a/Data/Manifold/Types.hs
+++ b/Data/Manifold/Types.hs
@@ -28,14 +28,16 @@
 {-# LANGUAGE TypeOperators            #-}
 {-# LANGUAGE ScopedTypeVariables      #-}
 {-# LANGUAGE UnicodeSyntax            #-}
+{-# LANGUAGE PatternSynonyms          #-}
 
 
 module Data.Manifold.Types (
         -- * Index / ASCII names
           Real0, Real1, RealPlus, Real2, Real3
         , Sphere0, Sphere1, Sphere2
-        , Projective1, Projective2
+        , Projective0, Projective1, Projective2
         , Disk1, Disk2, Cone, OpenCone
+        , FibreBundle(..), TangentBundle
         -- * Linear manifolds
         , ZeroDim(..)
         , ℝ, ℝ⁰, ℝ¹, ℝ², ℝ³, ℝ⁴
@@ -44,11 +46,11 @@
         , Stiefel1(..), stiefel1Project, stiefel1Embed
         -- ** Specific examples
         , HasUnitSphere(..)
-        , S⁰(..), S¹(..), S²(..)
+        , S⁰(..), S¹(..), pattern S¹, S²(..), pattern S²
         -- * Projective spaces
-        , ℝP¹,  ℝP²(..)
+        , ℝP⁰(..), ℝP¹(..), pattern ℝP¹,  ℝP²(..), pattern ℝP²
         -- * Intervals\/disks\/cones
-        , D¹(..), D²(..)
+        , D¹(..), D²(..), pattern D²
         , ℝay
         , CD¹(..), Cℝay(..)
         -- * Affine subspaces
diff --git a/Data/Manifold/Types/Primitive.hs b/Data/Manifold/Types/Primitive.hs
--- a/Data/Manifold/Types/Primitive.hs
+++ b/Data/Manifold/Types/Primitive.hs
@@ -28,23 +28,25 @@
 {-# LANGUAGE TypeOperators            #-}
 {-# LANGUAGE ScopedTypeVariables      #-}
 {-# LANGUAGE RecordWildCards          #-}
+{-# LANGUAGE PatternSynonyms          #-}
 
 
 module Data.Manifold.Types.Primitive (
         -- * Index / ASCII names
           Real0, Real1, RealPlus, Real2, Real3
         , Sphere0, Sphere1, Sphere2
-        , Projective1, Projective2
+        , Projective0, Projective1, Projective2
         , Disk1, Disk2, Cone, OpenCone
+        , FibreBundle(..), TangentBundle
         -- * Linear manifolds
         , ZeroDim(..)
         , ℝ, ℝ⁰, ℝ¹, ℝ², ℝ³, ℝ⁴
         -- * Hyperspheres
-        , S⁰(..), otherHalfSphere, S¹(..), S²(..)
+        , S⁰(..), otherHalfSphere, S¹(..), pattern S¹, S²(..), pattern S²
         -- * Projective spaces
-        , ℝP¹,  ℝP²(..)
+        , ℝP⁰(..), ℝP¹(..), pattern ℝP¹,  ℝP²(..), pattern ℝP²
         -- * Intervals\/disks\/cones
-        , D¹(..), fromIntv0to1, D²(..)
+        , D¹(..), fromIntv0to1, D²(..), pattern D²
         , ℝay
         , CD¹(..), Cℝay(..)
         -- * Tensor products
@@ -57,6 +59,7 @@
 
 
 import Math.Manifold.Core.Types
+import Math.Manifold.Core.PseudoAffine (FibreBundle(..), TangentBundle, Interior)
 
 import Data.VectorSpace
 import Data.VectorSpace.Free
@@ -67,6 +70,7 @@
 import Data.Basis
 import Data.Void
 import Data.Monoid
+import Data.Fixed (mod')
 import Math.LinearMap.Category (type (⊗)())
 
 import Control.Applicative (Const(..), Alternative(..))
@@ -80,7 +84,8 @@
 
 import Data.Embedding
 
-
+import qualified Test.QuickCheck as QC
+import qualified Text.Show.Pragmatic as SP
 
 
 
@@ -98,52 +103,9 @@
 
 
 
--- | The ordinary unit sphere.
-data S² = S² { ϑParamS² :: !Double -- ^ Range @[0, π[@.
-             , φParamS² :: !Double -- ^ Range @[-π, π[@.
-             } deriving (Show)
 
 
 
--- | The two-dimensional real projective space, implemented as a unit disk with
---   opposing points on the rim glued together.
-data ℝP² = ℝP² { rParamℝP² :: !Double -- ^ Range @[0, 1]@.
-               , φParamℝP² :: !Double -- ^ Range @[-π, π[@.
-               } deriving (Show)
-
-
-
--- | The standard, closed unit disk. Homeomorphic to the cone over 'S¹', but not in the
---   the obvious, &#x201c;flat&#x201d; way. (And not at all, despite
---   the identical ADT definition, to the projective space 'ℝP²'!)
-data D² = D² { rParamD² :: !Double -- ^ Range @[0, 1]@.
-             , φParamD² :: !Double -- ^ Range @[-π, π[@.
-             } deriving (Show)
-
--- | A (closed) cone over a space @x@ is the product of @x@ with the closed interval 'D¹'
---   of &#x201c;heights&#x201d;,
---   except on its &#x201c;tip&#x201d;: here, @x@ is smashed to a single point.
---   
---   This construct becomes (homeomorphic-to-) an actual geometric cone (and to 'D²') in the
---   special case @x = 'S¹'@.
-data CD¹ x = CD¹ { hParamCD¹ :: !Double -- ^ Range @[0, 1]@
-                 , pParamCD¹ :: !x      -- ^ Irrelevant at @h = 0@.
-                 } deriving (Show)
-
-
--- | An open cone is homeomorphic to a closed cone without the &#x201c;lid&#x201d;,
---   i.e. without the &#x201c;last copy&#x201d; of @x@, at the far end of the height
---   interval. Since that means the height does not include its supremum, it is actually
---   more natural to express it as the entire real ray, hence the name.
-data Cℝay x = Cℝay { hParamCℝay :: !Double -- ^ Range @[0, &#x221e;[@
-                   , pParamCℝay :: !x      -- ^ Irrelevant at @h = 0@.
-                   } deriving (Show)
-
-
-
-
-
-
 class NaturallyEmbedded m v where
   embed :: m -> v
   coEmbed :: v -> m
@@ -156,23 +118,33 @@
   embed x = (embed x, zeroV)
   coEmbed (x,_) = coEmbed x
 
+instance NaturallyEmbedded ℝ⁰ ℝ⁰ where embed = id; coEmbed = id
+instance NaturallyEmbedded ℝ  ℝ  where embed = id; coEmbed = id
+instance NaturallyEmbedded ℝ² ℝ² where embed = id; coEmbed = id
+instance NaturallyEmbedded ℝ³ ℝ³ where embed = id; coEmbed = id
+instance NaturallyEmbedded ℝ⁴ ℝ⁴ where embed = id; coEmbed = id
+
 instance NaturallyEmbedded S⁰ ℝ where
   embed PositiveHalfSphere = 1
   embed NegativeHalfSphere = -1
   coEmbed x | x>=0       = PositiveHalfSphere
             | otherwise  = NegativeHalfSphere
 instance NaturallyEmbedded S¹ ℝ² where
-  embed (S¹ φ) = V2 (cos φ) (sin φ)
-  coEmbed (V2 x y) = S¹ $ atan2 y x
+  embed (S¹Polar φ) = V2 (cos φ) (sin φ)
+  coEmbed (V2 x y) = S¹Polar $ atan2 y x
 instance NaturallyEmbedded S² ℝ³ where
-  embed (S² ϑ φ) = V3 (cos φ * sin ϑ) (sin φ * sin ϑ) (cos ϑ)
-  coEmbed (V3 x y z) = S² (acos $ z/r) (atan2 y x)
-   where r = sqrt $ x^2 + y^2 + z^2
+  embed (S²Polar ϑ φ) = V3 (cos φ * sϑ) (sin φ * sϑ) (cos ϑ)
+   where sϑ = sin ϑ
+  {-# INLINE embed #-}
+  coEmbed (V3 x y z) = S²Polar (atan2 rxy z) (atan2 y x)
+   where rxy = sqrt $ x^2 + y^2
+  {-# INLINE coEmbed #-}
  
 instance NaturallyEmbedded ℝP² ℝ³ where
-  embed (ℝP² r φ) = V3 (r * cos φ) (r * sin φ) (sqrt $ 1-r^2)
-  coEmbed (V3 x y z) = ℝP² (sqrt $ 1-(z/r)^2) (atan2 (y/r) (x/r))
-   where r = sqrt $ x^2 + y^2 + z^2
+  embed (HemisphereℝP²Polar θ φ) = V3 (cθ * cos φ) (cθ * sin φ) (sin θ)
+   where cθ = cos θ
+  coEmbed (V3 x y z) = HemisphereℝP²Polar (atan2 rxy z) (atan2 y x)
+   where rxy = sqrt $ x^2 + y^2
 
 instance NaturallyEmbedded D¹ ℝ where
   embed = xParamD¹
@@ -210,6 +182,7 @@
 type Sphere1 = S¹
 type Sphere2 = S²
 
+type Projective0 = ℝP⁰
 type Projective1 = ℝP¹
 type Projective2 = ℝP²
 
@@ -229,3 +202,48 @@
 
 
 
+instance QC.Arbitrary S⁰ where
+  arbitrary = (\hsph -> if hsph then PositiveHalfSphere else NegativeHalfSphere)
+               <$> QC.arbitrary
+instance SP.Show S⁰ where
+  showsPrec = showsPrec
+
+instance QC.Arbitrary S¹ where
+  arbitrary = S¹Polar . (pi-) . (`mod'`(2*pi))
+               <$> QC.arbitrary
+  shrink (S¹Polar φ) = S¹Polar . (pi/12*) <$> QC.shrink (φ*12/pi)
+instance SP.Show S¹ where
+  showsPrec p (S¹Polar φ) = showParen (p>9) $ ("S¹Polar "++) . SP.showsPrec 10 φ
+
+instance QC.Arbitrary S² where
+  arbitrary = ( \θ φ -> S²Polar (θ`mod'`pi) (pi - (φ`mod'`(2*pi))) )
+               <$> QC.arbitrary<*>QC.arbitrary
+  shrink (S²Polar θ φ) = uncurry S²Polar . (pi/12*^) <$> QC.shrink (θ*12/pi, φ*12/pi)
+instance SP.Show S² where
+  showsPrec p (S²Polar θ φ) = showParen (p>9) $ ("S²Polar "++)
+                           . SP.showsPrec 10 θ . (' ':) . SP.showsPrec 10 φ
+
+instance QC.Arbitrary ℝP⁰ where
+  arbitrary = pure ℝPZero
+
+instance QC.Arbitrary ℝP¹ where
+  arbitrary = ( \θ -> HemisphereℝP¹Polar (pi/2 - (θ`mod'`pi)) ) <$> QC.arbitrary
+  shrink (HemisphereℝP¹Polar θ) = HemisphereℝP¹Polar . (pi/6*) <$> QC.shrink (θ*6/pi)
+
+instance QC.Arbitrary ℝP² where
+  arbitrary = ( \θ φ -> HemisphereℝP²Polar (θ`mod'`pi/2) (pi - (φ`mod'`(2*pi))) )
+               <$> QC.arbitrary<*>QC.arbitrary
+  shrink (HemisphereℝP²Polar θ φ) = [ HemisphereℝP²Polar (θ'*pi/6) (φ'*pi/12)
+                                    | θ' <- QC.shrink (θ*6/pi)
+                                    , φ' <- QC.shrink (φ*12/pi) ]
+
+
+instance (SP.Show (Interior m), SP.Show f) => SP.Show (FibreBundle m f) where
+  showsPrec p (FibreBundle m v) = showParen (p>9)
+                $ ("FibreBundle "++) . SP.showsPrec 10 m
+                            . (' ':) . SP.showsPrec 10 v
+instance (QC.Arbitrary (Interior m), QC.Arbitrary f) => QC.Arbitrary (FibreBundle m f) where
+  arbitrary = FibreBundle <$> QC.arbitrary <*> QC.arbitrary
+  shrink (FibreBundle m v) = [ FibreBundle m' v'
+                             | m' <- QC.shrink m
+                             , v' <- QC.shrink v ]
diff --git a/Data/Manifold/Web.hs b/Data/Manifold/Web.hs
--- a/Data/Manifold/Web.hs
+++ b/Data/Manifold/Web.hs
@@ -49,7 +49,8 @@
             , localModels_CGrid
               -- * Differential equations
               -- ** Fixed resolution
-            , iterateFilterDEqn_static, iterateFilterDEqn_static_selective
+            , iterateFilterDEqn_static, iterateFilterDEqn_pathwise
+            , iterateFilterDEqn_static_selective
               -- ** Automatic resolution
             , filterDEqnSolutions_adaptive, iterateFilterDEqn_adaptive
               -- ** Configuration
@@ -59,7 +60,7 @@
             , PropagationInconsistency(..)
               -- * Misc
             , ConvexSet(..), ellipsoid, ellipsoidSet, coerceWebDomain
-            , rescanPDELocally, webOnions, knitShortcuts
+            , rescanPDELocally, localOnion, webOnions, knitShortcuts
             ) where
 
 
@@ -818,6 +819,28 @@
                            . fmap (shading $->)
 
 
+iterateFilterDEqn_pathwise
+     :: ( ModellableRelation x y, Hask.MonadPlus m, Hask.Traversable m, LocalModel ㄇ )
+       => InformationMergeStrategy [] m (x,Shade' y) iy
+           -> Embedding (->) (Shade' y) iy
+           -> DifferentialEqn ㄇ x y
+                 -> PointsWeb x (Shade' y) -> Cofree m (PointsWeb x (Shade' y))
+iterateFilterDEqn_pathwise strategy shading f
+            = fmap (fmap (shading >-$))
+            . (`evalState`7438)
+            . unfoldM (\oldWeb -> do
+               r <- get
+               let i = r `mod` nLeaves (webNodeRsc oldWeb)
+                   m = 2^31 - 1
+                   a = 963345    :: Int  -- revised Park-Miller
+               put $ (a*r)`mod`m
+               return ( oldWeb
+                      , filterDEqnSolutions_static strategy shading f
+                       =<<filterDEqnSolutions_pathsTowards i strategy shading f oldWeb
+                      ))
+            . fmap (shading $->)
+
+
 iterateFilterDEqn_static_selective :: ( ModellableRelation x y
                                       , Hask.MonadPlus m, badness ~ ℝ
                                       , LocalModel ㄇ )
@@ -869,6 +892,42 @@
               _ -> mergeInformation strategy oldValue empty
         )
 
+
+filterDEqnSolutions_pathsTowards :: ∀ x y ㄇ iy m .
+                     ( ModellableRelation x y, Hask.MonadPlus m, LocalModel ㄇ )
+       => WebNodeId
+          -> InformationMergeStrategy [] m  (x,Shade' y) iy
+          -> Embedding (->) (Shade' y) iy
+          -> DifferentialEqn ㄇ x y -> PointsWeb x iy -> m (PointsWeb x iy)
+filterDEqnSolutions_pathsTowards = case ( geodesicWitness :: GeodesicWitness y
+                                        , boundarylessWitness :: BoundarylessWitness x ) of
+   (GeodesicWitness _, BoundarylessWitness) -> \targetNode strategy shading f
+       -> traversePathsTowards targetNode
+            (\(PathStep stepStart stepEnd) -> StateT $
+              \odeState ->
+                let apriori = shading >-$ stepEnd^.thisNodeData
+                in case propagateDEqnSolution_loc
+                                f
+                                (LocalDataPropPlan{
+                                   _sourcePosition = stepStart^.thisNodeCoord
+                                 , _targetPosOffset = (stepEnd^.thisNodeCoord)
+                                                        .-~! (stepStart^.thisNodeCoord)
+                                 , _sourceData = odeState
+                                 , _targetAPrioriData = apriori
+                                 , _relatedData
+                                     = (fmap (second ((shading>-$) . _thisNodeData))
+                                               . concat . tail $ localOnion stepStart
+                                                                  [stepEnd^.thisNodeId])
+                                 }) of
+                          Nothing -> undefined
+                              <$> mergeInformation strategy (stepEnd^.thisNodeData) []
+                          Just propd -> (, propd)
+                                  <$> mergeInformation strategy
+                                        (stepEnd^.thisNodeData)
+                                        [ ( stepEnd^.thisNodeCoord, apriori )
+                                        , ( stepStart^.thisNodeCoord, propd ) ] )
+            (\startPoint pathTrav
+               -> evalStateT pathTrav $ shading >-$ startPoint^.thisNodeData)
 
 
 data Average a = Average { weight :: Int
diff --git a/Data/Manifold/Web/Internal.hs b/Data/Manifold/Web/Internal.hs
--- a/Data/Manifold/Web/Internal.hs
+++ b/Data/Manifold/Web/Internal.hs
@@ -63,6 +63,7 @@
 import Control.Monad (guard, forM_)
 import Control.Comonad
 import Control.Monad.Trans.State
+import Control.Monad.Trans.Writer
 
 import Control.DeepSeq
 
@@ -511,6 +512,16 @@
 
 type WNIPath = [WebNodeId]
 type NodeSet = ℤSet.IntSet
+
+
+pathsTowards :: ∀ x y . (WithField ℝ Manifold x, HasCallStack)
+     => WebNodeId -> PointsWeb x y -> [[y]]
+pathsTowards target web = execWriter $ traversePathsTowards
+       target
+       (\(PathStep _ y) -> tell [y^.thisNodeData] >> return (y^.thisNodeData))
+       (\startNode (WriterT (Identity (ν, pathTrav)))
+            -> tell [startNode^.thisNodeData : pathTrav] >> return ν)
+       web
 
 traversePathInIWeb :: ∀ φ x y . (WithField ℝ Manifold x, Monad φ, HasCallStack)
      => [WebNodeId] -> (PathStep x y -> φ y)
diff --git a/Math/Manifold/Embedding/Simple/Class.hs b/Math/Manifold/Embedding/Simple/Class.hs
new file mode 100644
--- /dev/null
+++ b/Math/Manifold/Embedding/Simple/Class.hs
@@ -0,0 +1,26 @@
+-- |
+-- Module      : Math.Manifold.Embedding.Simple.Class
+-- Copyright   : (c) Justus Sagemüller 2018
+-- License     : GPL v3
+-- 
+-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de
+-- Stability   : experimental
+-- Portability : portable
+-- 
+-- Some manifolds are “naturally” embedded within some bigger space. For instance,
+-- the topological spheres are readily identified with the geometric unit spheres in
+-- real vector spaces.
+--
+-- An embedding is a pretty strong relationship, but often all that's needed is being
+-- able to map single points from the manifold to the enclosing space. This module offers
+-- a class which does just that.
+
+
+
+
+module Math.Manifold.Embedding.Simple.Class (
+          NaturallyEmbedded(..)
+   ) where
+
+
+import Data.Manifold.Types.Primitive
diff --git a/manifolds.cabal b/manifolds.cabal
--- a/manifolds.cabal
+++ b/manifolds.cabal
@@ -1,5 +1,5 @@
 Name:                manifolds
-Version:             0.4.4.0
+Version:             0.4.5.0
 Category:            Math
 Synopsis:            Coordinate-free hypersurfaces
 Description:         Manifolds, a generalisation of the notion of &#x201c;smooth curves&#x201d; or surfaces,
@@ -40,10 +40,10 @@
 
 Library
   Build-Depends:     base>=4.5 && < 6
-                     , manifolds-core == 0.4.4.0
+                     , manifolds-core == 0.4.5.0
                      , transformers
                      , vector-space>=0.8
-                     , free-vector-spaces>=0.1.1
+                     , free-vector-spaces>=0.1.5
                      , linear
                      , MemoTrie
                      , vector
@@ -59,8 +59,9 @@
                      , placeholders
                      , lens
                      , call-stack
-                     , constrained-categories >= 0.2.3 && < 0.3.1
+                     , constrained-categories >= 0.3.1 && < 0.4
                      , pragmatic-show
+                     , QuickCheck >=2.8 && <3
   other-extensions:  FlexibleInstances
                      , TypeFamilies
                      , FlexibleContexts
@@ -80,6 +81,7 @@
                      Data.Manifold.Web.Internal
                      Data.Manifold.DifferentialEquation
                      Data.Manifold.Function.LocalModel
+                     Data.Manifold.Function.Interpolation
                      Data.SimplicialComplex
                      Data.Function.Differentiable
                      Data.Function.Affine
@@ -87,7 +89,9 @@
                      Data.Manifold.Types.Stiefel
                      Data.Manifold.Griddable
                      Data.Manifold.Atlas
+                     Data.Manifold.FibreBundle
                      Data.Manifold.Riemannian
+                     Math.Manifold.Embedding.Simple.Class
   Other-modules:   Data.List.FastNub
                    Data.Manifold.Types.Primitive
                    Data.SetLike.Intersection
@@ -119,6 +123,7 @@
     , pragmatic-show
     , containers
     , vector-space
+    , linear
     , constrained-categories
     , linearmap-category
     , lens
diff --git a/test/tasty/test.hs b/test/tasty/test.hs
--- a/test/tasty/test.hs
+++ b/test/tasty/test.hs
@@ -9,17 +9,23 @@
 -- 
 
 {-# LANGUAGE OverloadedLists, TypeFamilies, FlexibleContexts, UndecidableInstances #-}
-{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE FlexibleInstances, AllowAmbiguousTypes  #-}
+{-# LANGUAGE TypeOperators, TypeApplications, ScopedTypeVariables, UnicodeSyntax #-}
 
 module Main where
 
 import Data.Manifold.Types
 import Data.Manifold.PseudoAffine
+import Data.Manifold.FibreBundle
 import Data.Manifold.TreeCover
 import Data.Manifold.Web
 import Data.Manifold.Web.Internal
 import Data.Manifold.Function.LocalModel
+import Math.Manifold.Embedding.Simple.Class
 import Data.VectorSpace
+import Data.Cross (cross3)
+import Linear.V2 (V2(V2))
+import Linear.V3 (V3(V3))
 import Math.LinearMap.Category
 import Prelude hiding (id, fst, snd)
 import Control.Category.Constrained (id)
@@ -35,7 +41,7 @@
 import qualified Data.Graph as Graph
 import qualified Data.Set as Set
 import Control.Arrow
-import Control.Lens
+import Control.Lens hiding ((<.>))
 
 import qualified Text.Show.Pragmatic as SP
 
@@ -44,7 +50,262 @@
 
 tests :: TestTree
 tests = testGroup "Tests"
- [ testGroup "Graph structure of webs"
+ [ testGroup "Semimanifold laws"
+  [ testGroup "Asymptotic associativity"
+   [ QC.testProperty "Real vector space" (nearlyAssociative @(ℝ,ℝ))
+   , QC.testProperty "1-sphere" (nearlyAssociative @S¹)
+   , QC.testProperty "Projective line" (nearlyAssociative @ℝP¹)
+   , QC.testProperty "2-sphere" (QC.expectFailure $ nearlyAssociative @S²)
+   , QC.testProperty "Projective plane" (QC.expectFailure $ nearlyAssociative @ℝP²)
+   ]
+  ]
+ , testGroup "Pseudoaffine laws"
+  [ testGroup "Displacement cancellation"
+   [ QC.testProperty "Real vector space" (originCancellation @(ℝ,ℝ))
+   , QC.testProperty "1-sphere" (originCancellation @S¹)
+   , QC.testProperty "Projective line" (originCancellation @ℝP¹)
+   , QC.testProperty "2-sphere" (originCancellation @S²)
+   , testGroup "2-sphere corner cases"
+    [ QC.testProperty "To north pole"
+        $ \(S¹Polar φ) p -> originCancellation (S²Polar 0 φ) p
+    , QC.testProperty "From north pole"
+        $ \(S¹Polar φ) p -> originCancellation p (S²Polar 0 φ)
+    , QC.testProperty "To south pole"
+        $ \(S¹Polar φ) p -> originCancellation (S²Polar pi φ) p
+    , QC.testProperty "From south pole"
+        $ \(S¹Polar φ) p -> originCancellation p (S²Polar pi φ)
+    , QC.testProperty "South- to north pole"
+        $ \(S¹Polar φ) (S¹Polar ψ) -> originCancellation (S²Polar 0 φ) (S²Polar pi ψ)
+    , QC.testProperty "North- to south pole"
+        $ \(S¹Polar φ) (S¹Polar ψ) -> originCancellation (S²Polar pi ψ) (S²Polar 0 φ)
+    , QC.testProperty "Along equator"
+        $ \(S¹Polar φ) (S¹Polar ψ) -> originCancellation (S²Polar (pi/2) ψ) (S²Polar (pi/2) φ)
+    , QC.testProperty "Just south of equator"
+        $ \(S¹Polar φ) (S¹Polar ψ) -> originCancellation (S²Polar (pi/2 + 1e-10) ψ) (S²Polar (pi/2 + 1e-10) φ)
+    , QC.testProperty "Just across the equator"
+        $ \(S¹Polar φ) (S¹Polar ψ) -> originCancellation (S²Polar (pi/2) ψ) (S²Polar (pi/2 + 1e-10) φ)
+    , QC.testProperty "To equator"
+        $ \(S¹Polar φ) p -> originCancellation (S²Polar (pi/2) φ) p
+    , QC.testProperty "From equator"
+        $ \(S¹Polar φ) p -> originCancellation p (S²Polar (pi/2) φ)
+    ]
+   , QC.testProperty "Projective plane" (originCancellation @ℝP²)
+   ]
+  ]
+ , testGroup "Natural embeddings"
+  [ testGroup "1-sphere"
+     [ testCase "North pole" $ embed (S¹Polar $ pi/2) @?≈ (V2 0 1 :: ℝ²)
+     , testCase "South pole" $ embed (S¹Polar $ -pi/2) @?≈ (V2 0 (-1) :: ℝ²)
+     ]
+  , testGroup "2-sphere"
+     [ testCase "North pole" $ embed (S²Polar 0 0) @?≈ (V3 0 0 1 :: ℝ³)
+     , testCase "South pole" $ embed (S²Polar pi 0) @?≈ (V3 0 0 (-1) :: ℝ³)
+     ]
+  , testGroup "1-sphere tangent bundle"
+     [ testCase "North pole"
+           $ embed (FibreBundle (S¹Polar $  pi/2) 1 :: TangentBundle S¹)
+               @?≈ (FibreBundle (V2 0 1) (V2 (-1) 0) :: TangentBundle ℝ²)
+     , testCase "South pole"
+           $ embed (FibreBundle (S¹Polar $ -pi/2) 1 :: TangentBundle S¹)
+               @?≈ (FibreBundle (V2 0 (-1)) (V2 1 0) :: TangentBundle ℝ²)
+     , testCase "45°"
+           $ embed (FibreBundle (S¹Polar $ pi/4) 1 :: TangentBundle S¹)
+               @?≈ (FibreBundle (V2 1 1^/sqrt 2) (V2 (-1) 1^/sqrt 2) :: TangentBundle ℝ²)
+     ]
+  , testGroup "2-sphere tangent bundle"
+     [ testCase "North pole, x-dir"
+           $ embed (FibreBundle (S²Polar 0 0) (V2 1 0) :: TangentBundle S²)
+               @?≈ (FibreBundle (V3 0 0 1) (V3 1 0 0) :: TangentBundle ℝ³)
+     , testCase "North pole (alternative φ), x-dir"
+           $ embed (FibreBundle (S²Polar 0 1.524) (V2 1 0) :: TangentBundle S²)
+               @?≈ (FibreBundle (V3 0 0 1) (V3 1 0 0) :: TangentBundle ℝ³)
+     , testCase "North pole, y-dir"
+           $ embed (FibreBundle (S²Polar 0 0) (V2 0 1) :: TangentBundle S²)
+               @?≈ (FibreBundle (V3 0 0 1) (V3 0 1 0) :: TangentBundle ℝ³)
+     , testCase "Close to north pole"
+           $ embed (FibreBundle (S²Polar 1e-11 0.602) (V2 3.7 1.1) :: TangentBundle S²)
+               @?≈ (FibreBundle (V3 0 0 1) (V3 3.7 1.1 0) :: TangentBundle ℝ³)
+     , testCase "South pole, x-dir"
+           $ embed (FibreBundle (S²Polar pi 0) (V2 1 0) :: TangentBundle S²)
+               @?≈ (FibreBundle (V3 0 0 (-1)) (V3 (-1) 0 0) :: TangentBundle ℝ³)
+     , testCase "South pole, y-dir"
+           $ embed (FibreBundle (S²Polar pi 0) (V2 0 1) :: TangentBundle S²)
+               @?≈ (FibreBundle (V3 0 0 (-1)) (V3 0 1 0) :: TangentBundle ℝ³)
+     , testCase "Close to south pole"
+           $ embed (FibreBundle (S²Polar (pi-1e-11) 0.602) (V2 3.7 1.1) :: TangentBundle S²)
+               @?≈ (FibreBundle (V3 0 0 (-1)) (V3 (-3.7) 1.1 0) :: TangentBundle ℝ³)
+     , testCase "Equator, y-dir"
+           $ embed (FibreBundle (S²Polar (pi/2) 0) (V2 0 1) :: TangentBundle S²)
+               @?≈ (FibreBundle (V3 1 0 0) (V3 0 1 0) :: TangentBundle ℝ³)
+     , testCase "Equator, x-dir"
+           $ embed (FibreBundle (S²Polar (pi/2) (pi/2)) (V2 1 0) :: TangentBundle S²)
+               @?≈ (FibreBundle (V3 0 1 0) (V3 (-1) 0 0) :: TangentBundle ℝ³)
+     , testCase "Equator, z-dir"
+           $ embed (FibreBundle (S²Polar (pi/2) 0) (V2 1 0) :: TangentBundle S²)
+               @?≈ (FibreBundle (V3 1 0 0) (V3 0 0 (-1)) :: TangentBundle ℝ³)
+     ]
+  ]
+ , testGroup "Embedding tangent bundles"
+  [ QC.testProperty "Real vector space" (embeddingTangentiality @ℝ² @ℝ² 1)
+  , QC.testProperty "1-sphere (unlimited)" (QC.expectFailure
+                                       $ embeddingTangentiality @ℝ² @S¹ 1)
+  , QC.testProperty "1-sphere" (embeddingTangentiality @ℝ² @S¹ 1e-5)
+  , QC.testProperty "2-sphere" (embeddingTangentiality @ℝ³ @S² 1e-5)
+  ]
+ , testGroup "Embedding back-projection"
+  [ QC.testProperty "Real vector space" (embeddingBackProject @(ℝ,ℝ) @ℝ)
+  , QC.testProperty "1-sphere" (embeddingBackProject @ℝ² @S¹)
+  , QC.testProperty "2-sphere" (embeddingBackProject @ℝ³ @S²)
+  , QC.testProperty "Vector space tangent bundle"
+       (embeddingBackProject @(TangentBundle (ℝ,ℝ)) @(TangentBundle ℝ) )
+  , QC.testProperty "S¹ tangent bundle"
+       (embeddingBackProject @(TangentBundle ℝ²) @(TangentBundle S¹) )
+  , QC.testProperty "S² tangent bundle"
+       (embeddingBackProject @(TangentBundle ℝ³) @(TangentBundle S²) )
+  ]
+ , testGroup "Special properties of translations"
+  [ testGroup "2-sphere"
+   [ QC.testProperty "S²-movement as rotation in ℝ³"
+      $ \p v -> let FibreBundle pCart vCart :: TangentBundle ℝ³
+                         = embed (FibreBundle p v :: TangentBundle S²)
+                    q = p .+~^ v :: S²
+                    qCart = embed q :: ℝ³
+                    axis = pCart `cross3` qCart
+                    FibreBundle _ axisProj :: TangentBundle S²
+                        = coEmbed (FibreBundle pCart axis :: TangentBundle ℝ³)
+                in vCart <.> axis + 1 ≈ 1    -- i.e. the movement vector is always
+                  && v <.> axisProj + 1 ≈ 1  -- orthogonal to the rotation axis.
+   ]
+  ]
+ , testGroup "Parallel transport"
+  [ testGroup "Displacement cancellation"
+   [ QC.testProperty "Real vector space" (parTransportAssociativity @(ℝ,ℝ))
+   , QC.testProperty "1-sphere" (parTransportAssociativity @S¹)
+   ]
+  , testGroup "Nearby tangent spaces of embedding"
+   [ QC.testProperty "Real vector space" (nearbyTangentSpaceEmbedding @(ℝ,ℝ) @ℝ 1)
+   , QC.testProperty "1-sphere (unlimited)"
+         $ QC.expectFailure (nearbyTangentSpaceEmbedding @ℝ² @S¹ 1)
+   , QC.testProperty "1-sphere" (nearbyTangentSpaceEmbedding @ℝ² @S¹ 1e-5)
+   , QC.testProperty "2-sphere" (nearbyTangentSpaceEmbedding @ℝ³ @S² 1e-5)
+   ]
+  , testGroup "2-sphere"
+   [ testCase "Non-movement on the equator"
+        $ sphereParallelTransportTest
+            (S²Polar (pi/2) 0) (S²Polar (pi/2) 0) [V3 0 0 1] [V3 0 0 1]
+   , testCase "Micro-movement on the equator"
+        $ sphereParallelTransportTest
+            (S²Polar (pi/2) 0) (S²Polar (pi/2) 1e-3) [V3 0 0 1] [V3 0 0 1]
+   , testCase "Small movement on the equator (ez)"
+        $ sphereParallelTransportTest
+            (S²Polar (pi/2) 0) (S²Polar (pi/2) (pi/2)) [V3 0 0 1, V3   0  1 0]
+                                                       [V3 0 0 1, V3 (-1) 0 0]
+   , testCase "Big movement on the equator"
+        $ sphereParallelTransportTest
+            (S²Polar (pi/2) 0) (S²Polar (pi/2) 3) [V3 0 0 1] [V3 0 0 1]
+   , testCase "Big negative movement on the equator"
+        $ sphereParallelTransportTest
+            (S²Polar (pi/2) 0) (S²Polar (pi/2) (-3)) [V3 0 0 1] [V3 0 0 1]
+   , testCase "Movement on the zero meridian from north pole"
+        $ sphereParallelTransportTest
+            (S²Polar 0 0) (S²Polar (pi/2) 0) [V3 0 1 0] [V3 0 1 0]
+   , testCase "Movement on the zero meridian to north pole"
+        $ sphereParallelTransportTest
+            (S²Polar (pi/2) 0) (S²Polar 0 0) [V3 0 1 0, V3   0  0 1]
+                                             [V3 0 1 0, V3 (-1) 0 0]
+   , testCase "Crossing the equator on the zero meridian"
+        $ sphereParallelTransportTest
+            (S²Polar (pi/4) 0) (S²Polar (3*pi/4) 0) [V3 0 1 0, V3 (-1) 0 1] 
+                                                    [V3 0 1 0, V3   1  0 1]
+   , testCase "Crossing the equator on the 90° meridian"
+        $ sphereParallelTransportTest
+            (S²Polar (pi/4) (pi/2)) (S²Polar (3*pi/4) (pi/2)) [V3 1 0 0, V3 0 (-1) 1]
+                                                              [V3 1 0 0, V3 0   1  1]
+   , testCase "Crossing the equator on the 180° meridian"
+        $ sphereParallelTransportTest
+            (S²Polar (pi/4) pi) (S²Polar (3*pi/4) pi) [V3 0 1 0, V3   1  0 1]
+                                                      [V3 0 1 0, V3 (-1) 0 1]
+   , testCase "Crossing the equator on the -90° meridian"
+        $ sphereParallelTransportTest
+            (S²Polar (pi/4) (-pi/2)) (S²Polar (3*pi/4) (-pi/2)) [V3 1 0 0, V3 0   1  1]
+                                                                [V3 1 0 0, V3 0 (-1) 1]
+   , QC.testProperty "Movement on the equator" . QC.expectFailure
+        $ \(S¹Polar φ₀) (S¹Polar φ₁) -> assertParTransportNeedleTargetFixpoint
+                 (S²Polar 0 0, Just "north pole")
+                 (S²Polar (pi/2) φ₀)
+                 (S²Polar (pi/2) φ₁)
+   , QC.testProperty "Just north of the equator"
+        $ \p@(S¹Polar φ₀) q@(S¹Polar φ₁) -> abs (p.-~!q) < 2
+            ==> assertParTransportNeedleTargetFixpoint
+                 (S²Polar 0 0, Just "north pole")
+                 (S²Polar (pi/2-1e-13) φ₀)
+                 (S²Polar (pi/2-1e-13) φ₁)
+   , QC.testProperty "Just slightly crossing the equator"
+        $ \(S¹Polar φ₀) (S¹Polar φ₁) -> assertParTransportNeedleTargetFixpoint
+                 (S²Polar 0 0, Just "north pole")
+                 (S²Polar (pi/2-1e-13) φ₀)
+                 (S²Polar (pi/2+1e-13) φ₁)
+   , QC.testProperty "Just south of the equator"
+        $ \p@(S¹Polar φ₀) q@(S¹Polar φ₁) -> abs (p.-~!q) < 2
+            ==> assertParTransportNeedleTargetFixpoint
+                 (S²Polar pi 0, Just "south pole")
+                 (S²Polar (pi/2+1e-13) φ₀)
+                 (S²Polar (pi/2+1e-13) φ₁)
+   , QC.testProperty "Movement on the zero meridian"
+        $ \(S¹Polar θ₀) (S¹Polar θ₁) -> assertParTransportNeedleTargetFixpoint
+                 (S²Polar (pi/2) (pi/2), Nothing)
+                 (S²Polar (abs θ₀) (if θ₀>0 then 0 else pi))
+                 (S²Polar (abs θ₁) (if θ₁>0 then 0 else pi))
+   , QC.testProperty "Rotation axis – heading-vector"
+        $ \p v -> let q = p .+~^ v :: S²
+                      w = parallelTransport p v v
+                      FibreBundle pCart vCart
+                          = embed (FibreBundle p v :: TangentBundle S²) :: TangentBundle ℝ³
+                      FibreBundle qCart wCart
+                          = embed (FibreBundle q w :: TangentBundle S²) :: TangentBundle ℝ³
+                      pxv = pCart`cross3`vCart
+                      qxw = qCart`cross3`wCart
+                    in QC.counterexample
+                           ("  𝑝 = "++SP.show p++"\t ≃ "++SP.show pCart
+                        ++"\n  𝑞 = "++SP.show q++"\t ≃ "++SP.show qCart
+                        ++"\n  𝑣 = "++SP.show v++"\t = "++SP.show vCart++" @ 𝑝"
+                        ++"\n  𝑤 = "++SP.show w++"\t = "++SP.show wCart++" @ 𝑞"
+                        ++"\n𝑝×𝑣 = "++SP.show pxv    -- rotation axis
+                        ++"\n𝑞×𝑤 = "++SP.show qxw    -- rotation axis
+                             )
+                       $ pxv ≈ qxw
+   , QC.testProperty "Rotation axis – arbitrary vectors"
+        $ \p v f -> let q = p .+~^ v :: S²
+                        g = parallelTransport p v f
+                        FibreBundle pCart fCart
+                          = embed (FibreBundle p f :: TangentBundle S²) :: TangentBundle ℝ³
+                        FibreBundle qCart gCart
+                          = embed (FibreBundle q g :: TangentBundle S²) :: TangentBundle ℝ³
+                        infix 7 ×
+                        (×) = cross3
+                        pxq = pCart×qCart
+                        fㄧg = fCart ^-^ gCart
+                        ㄍ = magnitudeSq
+                    in QC.counterexample
+                           ("              𝑝 = "++SP.show p
+                        ++"\n              𝑞 = "++SP.show q
+                        ++"\n              𝑓 = "++SP.show f
+                        ++"\n              𝑔 = "++SP.show g
+                        ++"\n            𝑝×𝑞 = "++SP.show pxq  -- rotation axis
+                        ++"\n          𝑓 − 𝑔 = "++SP.show fㄧg -- movement in the rot.-plane
+                        ++"\n    (𝑝×𝑞)×(𝑓−𝑔) = "++SP.show (pxq × fㄧg)
+                        ++"\n    (𝑝×𝑞)·(𝑓−𝑔) = "++SP.show (pxq <.> fㄧg)
+                        ++"\n ‖(𝑝×𝑞)×(𝑓−𝑔)‖² = "++SP.show (ㄍ $ pxq × fㄧg)
+                        ++"\n         ‖𝑝×𝑞‖² = "++SP.show (ㄍ pxq)
+                        ++"\n         ‖𝑓−𝑔‖² = "++SP.show (ㄍ fㄧg)
+                        ++"\n  ‖𝑝×𝑞‖²·‖𝑓−𝑔‖² = "++SP.show (ㄍ pxq*ㄍ fㄧg)
+                             )
+                       $ ㄍ (pxq × fㄧg)      -- Check that 𝑝×𝑞 and 𝑓−𝑔 are orthogonal.
+                          ≈ ㄍ pxq * ㄍ fㄧg  -- (For orthogonal 𝐚 and 𝐛, we have
+                                              -- ‖𝐚×𝐛‖ = ‖𝐚‖·‖𝐛‖.)
+   ]
+  ]
+ , testGroup "Graph structure of webs"
   [ testCase "Manually-defined empty web."
     $ toList (fst $ toGraph emptyWeb) @?= []
   , testCase "Manually-defined single-point web."
@@ -439,33 +700,80 @@
 
 infix 4 ≈
 class AEq e where
+  fuzzyEq :: ℝ -> e -> e -> Bool
+  unitEpsilon :: ℝ
+  unitEpsilon = 1e-9
   (≈) :: e -> e -> Bool
+  (≈) = fuzzyEq (unitEpsilon @e)
+
+instance AEq Double where
+  fuzzyEq η x y  = x + abs x*η >= y
+          && x - abs x*η <= y
 instance (SimpleSpace v, Needle v~v, Interior v~v, Floating (Scalar v))
              => AEq (Shade' v) where
-  Shade' c₀ σ₀ ≈ Shade' c₁ σ₁
+  fuzzyEq η (Shade' c₀ σ₀) (Shade' c₁ σ₁)
     = (σ₀|$|δ) < ε && (σ₀|$|δ) < ε
      && all (is1 . (σ₀|$|)) (normSpanningSystem' σ₁)
      && all (is1 . (σ₁|$|)) (normSpanningSystem' σ₀)
    where δ = c₁ ^-^ c₀
-         ε = 1e-2
+         ε = 1e-2 + realToFrac η
          is1 x = abs (x-1) < ε
 instance ( SimpleSpace v, DualVector (Needle' v) ~ v, Interior v ~ v
          , InnerSpace (Scalar v), Scalar (Needle' v) ~ Scalar v )
               => AEq (Shade v) where
-  Shade c₀ σ₀ ≈ Shade c₁ σ₁
+  fuzzyEq η (Shade c₀ σ₀) (Shade c₁ σ₁)
     = (dualNorm σ₀|$|δ) < ε && (dualNorm σ₀|$|δ) < ε
      && all (is1 . (dualNorm σ₀|$|)) (normSpanningSystem σ₁)
      && all (is1 . (dualNorm σ₁|$|)) (normSpanningSystem σ₀)
    where δ = c₁ ^-^ c₀
-         ε = 1e-2
+         ε = 1e-2 + realToFrac η
          is1 x = abs (x-1) < ε
 instance AEq a => AEq (Maybe a) where
-  Just x ≈ Just y = x ≈ y
-  Nothing ≈ Nothing = True
-  _ ≈ _ = False
+  fuzzyEq η (Just x) (Just y) = fuzzyEq η x y
+  fuzzyEq _ Nothing Nothing = True
+  fuzzyEq _ _ _ = False
 instance (AEq (Shade y), AEq (Shade (Needle x +> Needle y)))
               => AEq (AffineModel x y) where
-  AffineModel b₀ a₀ ≈ AffineModel b₁ a₁ = b₀ ≈ b₁ && a₀ ≈ a₁
+  fuzzyEq η (AffineModel b₀ a₀) (AffineModel b₁ a₁) = fuzzyEq η b₀ b₁ && fuzzyEq η a₀ a₁
+
+instance (AEq a, AEq b) => (AEq (a,b)) where
+  fuzzyEq η (x,y) (ξ,υ) = fuzzyEq η x ξ && fuzzyEq η y υ
+instance AEq S¹ where
+  fuzzyEq η (S¹Polar φ) (S¹Polar ϕ)
+   | φ > pi/2, ϕ < -pi/2  = fuzzyEq η (S¹Polar $ φ - 2*pi) (S¹Polar ϕ)
+   | ϕ > pi/2, φ < -pi/2  = fuzzyEq η (S¹Polar φ) (S¹Polar $ ϕ - 2*pi)
+   | otherwise            = abs (φ - ϕ) < η
+instance AEq S² where
+  fuzzyEq η (S²Polar θ φ) (S²Polar ϑ ϕ)
+   | φ > pi/2, ϕ < -pi/2  = fuzzyEq η (S²Polar θ $ φ - 2*pi) (S²Polar ϑ ϕ)
+   | ϕ > pi/2, φ < -pi/2  = fuzzyEq η (S²Polar θ φ) (S²Polar ϑ $ ϕ - 2*pi)
+   | otherwise            = abs (θ - ϑ) < η && abs (φ - ϕ) * sin θ < η
+
+instance AEq ℝ² where
+  fuzzyEq η (V2 x y) (V2 ξ υ) = abs (x - ξ) <= ε && abs (y - υ) <= ε
+   where ε = (maximum @[]) (abs<$>[x,y,ξ,υ]) * η
+instance AEq ℝ³ where
+  fuzzyEq η (V3 x y z) (V3 ξ υ ζ) = (all @[]) ((ε>=) . abs) $ [x-ξ, y-υ, z-ζ]
+   where ε = (maximum @[]) (abs<$>[x,y,z,ξ,υ,ζ]) * η
+
+instance AEq ℝP⁰ where
+  fuzzyEq _ ℝPZero ℝPZero  = True
+instance AEq ℝP¹ where
+  fuzzyEq η (HemisphereℝP¹Polar θ) (HemisphereℝP¹Polar ϑ)
+   = fuzzyEq η (S¹Polar $ θ*2) (S¹Polar $ ϑ*2)
+instance AEq ℝP² where
+  fuzzyEq η (HemisphereℝP²Polar θ φ) (HemisphereℝP²Polar ϑ ϕ)
+   | φ > pi/2, ϕ < -pi/2  = fuzzyEq η (HemisphereℝP²Polar θ $ φ - 2*pi) (HemisphereℝP²Polar ϑ ϕ)
+   | ϕ > pi/2, φ < -pi/2  = fuzzyEq η (HemisphereℝP²Polar θ φ) (HemisphereℝP²Polar ϑ $ ϕ - 2*pi)
+   | θ < pi/2             = abs (θ - ϑ) < η && abs (φ - ϕ) * θ < η
+   | φ > pi/4, ϕ < -pi/4  = fuzzyEq η (HemisphereℝP²Polar (pi/2) $ φ - pi)
+                                      (HemisphereℝP²Polar (pi/2) ϕ)
+   | ϕ > pi/4, φ < -pi/4  = fuzzyEq η (HemisphereℝP²Polar (pi/2) φ)
+                                      (HemisphereℝP²Polar (pi/2) $ ϕ - pi)
+   | otherwise            = abs (φ - ϕ) < η
+
+instance (AEq (Interior m), AEq f) => AEq (FibreBundle m f) where
+  fuzzyEq η (FibreBundle p v) (FibreBundle q w) = fuzzyEq η p q && fuzzyEq η v w
                                         
 infix 1 @?≈       
 (@?≈) :: (AEq e, Show e) => e -> e -> Assertion
@@ -473,4 +781,107 @@
  | a≈b        = return ()
  | otherwise  = assertFailure $ "Expected "++show b++", but got "++show a
 
+instance QC.Arbitrary ℝ² where
+  arbitrary = (\(x,y)->V2 x y) <$> QC.arbitrary
+  shrink (V2 x y) = V2 <$> ((/12)<$>QC.shrink (x*12))
+                       <*> ((/12)<$>QC.shrink (y*12))
 
+nearlyAssociative :: ∀ m . (AEq m, Semimanifold m, Interior m ~ m)
+                         => m -> Needle m -> Needle m -> Bool
+nearlyAssociative p v w = (p .+~^ v) .+~^ w ≈ (p .+~^ (v^+^w) :: m)
+
+originCancellation :: ∀ m . (AEq m, Manifold m, Show m, Show (Needle m))
+                         => m -> m -> QC.Property
+originCancellation p q = case ( boundarylessWitness :: BoundarylessWitness m
+                              , p.-~.q ) of
+      (BoundarylessWitness, Just v)
+          -> let p' = q.+~^v
+             in QC.counterexample ("v = "++show v++", q+v = "++show p') $ p' ≈ p
+
+embeddingBackProject :: ∀ m n . ( NaturallyEmbedded n m, AEq n, SP.Show m, SP.Show n )
+       => n -> QC.Property
+embeddingBackProject p = QC.counterexample ("Embedded: "++SP.show ep
+                                          ++", back-projected: "++SP.show p')
+                           $ p' ≈ p
+ where ep = embed p :: m
+       p' = coEmbed ep
+
+embeddingTangentiality :: ∀ m n . ( Semimanifold m, Semimanifold n
+                                  , NaturallyEmbedded n m
+                                  , NaturallyEmbedded (TangentBundle n) (TangentBundle m)
+                                  , SP.Show n, AEq n
+                                  , InnerSpace (Needle n), RealFloat (Scalar (Needle n)) )
+       => Scalar (Needle n) -> Interior n -> Needle n -> QC.Property
+embeddingTangentiality consistRadius p vub
+         = QC.counterexample ("p+v = "++SP.show q++", coEmbed (embed p+v) = "++SP.show q')
+            $ fuzzyEq (unitEpsilon @n * (1+rvub^2)) q q'
+ where rvub = realToFrac $ magnitude vub
+       v = vub ^* consistRadius
+       q, q' :: n
+       q = p .+~^ v
+       q' = coEmbed $ (pEmbd .+~^ vEmbd :: m)
+       o :: TangentBundle n
+       o = FibreBundle p v
+       FibreBundle pEmbd vEmbd = embed o :: TangentBundle m
+
+nearbyTangentSpaceEmbedding :: ∀ m n
+                     . ( Semimanifold m, Semimanifold n
+                       , m ~ Interior m, n ~ Interior n
+                       , NaturallyEmbedded n m
+                       , NaturallyEmbedded (TangentBundle n) (TangentBundle m)
+                       , ParallelTransporting (->) n (Needle n)
+                       , SP.Show n, SP.Show (Needle n), AEq (Needle n)
+                       , InnerSpace (Needle n), RealFloat (Scalar (Needle n)) )
+       => Scalar (Needle n) -> Interior n -> Needle n -> Needle n -> QC.Property
+nearbyTangentSpaceEmbedding consistRadius p vub f
+         = QC.counterexample ("𝑓 embd. at 𝑝, then proj. at 𝑝+𝑣 = "++SP.show fReProj
+                              ++", 𝑓 moved by 𝑣 = "++SP.show g)
+            $ fuzzyEq (unitEpsilon @(Needle n) * (1+rvub^2)) g fReProj
+ where rvub = realToFrac $ magnitude vub
+       v = vub ^* consistRadius
+       q :: n
+       q = p .+~^ v :: n
+       qEmbd = embed q :: m
+       FibreBundle _ fReProj :: TangentBundle n
+               = coEmbed (FibreBundle qEmbd fEmbd :: TangentBundle m)
+       g = parallelTransport p v f
+       o :: TangentBundle n
+       o = FibreBundle p f
+       FibreBundle pEmbd fEmbd = embed o :: TangentBundle m
+
+parTransportAssociativity :: ∀ m
+           . ( AEq m, Manifold m, SP.Show m
+             , ParallelTransporting (->) m (Needle m)
+             , InnerSpace (Needle m), RealFloat (Scalar (Needle m)) )
+                         => m -> Needle m -> Needle m -> QC.Property
+parTransportAssociativity p v w
+ = maximum (map magnitude [v,w]) < 1000
+       -- Very vast vectors incur inevitable floating-point uncertainty
+  ==> let q, q' :: m
+          q = (p .+~^ v) .+~^ parallelTransport p v w
+          q' = p .+~^ (v^+^w)
+      in QC.counterexample ("(p+v) + 〔pTp. v〕w = "++SP.show q++", p+(v+w) = "++SP.show q')
+          $ q ≈ q'
+
+assertParTransportNeedleTargetFixpoint :: ∀ m
+     . ( AEq m, Manifold m, SP.Show m, Show (Needle m)
+       , ParallelTransporting (->) m (Needle m) )
+    => (m, Maybe String) -> m -> m -> QC.Property
+assertParTransportNeedleTargetFixpoint (q, qName) p₀ p₁
+         = let q'= p₁ .+~^ parallelTransport p₀ (p₁ .-~! p₀) (q .-~! p₀)
+           in QC.counterexample
+                 ("Should keep pointing on "++qShw++", but got "++ SP.show q')
+               $ q' ≈ q
+ where qShw = case qName of
+        Just s  -> s
+        Nothing -> SP.show q
+
+
+sphereParallelTransportTest :: S² -> S² -> [ℝ³] -> [ℝ³] -> Assertion
+sphereParallelTransportTest p q [] [] = assert True
+sphereParallelTransportTest p q (v:vs) (w:ws)
+     = (parallelTransport p (q.-~!p) vSph @?≈ wSph)
+        >> sphereParallelTransportTest p q vs ws
+ where [FibreBundle _ vSph, FibreBundle _ wSph]
+          = [ coEmbed (FibreBundle (embed o) u :: TangentBundle ℝ³) :: TangentBundle S²
+            | (o,u) <- [(p,v), (q,w)] ]