Icosahedron → Sphere

data Vector = Vector {x, y, z :: Float}
    deriving Show

len (Vector x y z) = sqrt $ (x ^ 2) + (y ^ 2) + (z ^ 2)

(Vector ax ay az) +. (Vector bx by bz) = Vector (ax + bx) (ay + by) (az + bz)
(Vector ax ay az) -. (Vector bx by bz) = Vector (ax - bx) (ay - by) (az - bz)
k *. (Vector ax ay az) = Vector (k * ax) (k * ay) (k * az)
normalize vec = (recip $ len vec) *. vec
(Vector ax ay az) `scal` (Vector bx by bz) = ax * bx + ay * by + az * bz

type Face = (Vector, Vector, Vector)
type Figure = [Face]

icosahedron radius = top ++ middleTop ++ middleBottom ++ bottom
    where north = Vector 0 0 radius
          south = Vector 0 0 (-radius)
          topv = [Vector (radius * cos a) (radius * sin a) z | a <- init [0, 2 * pi / 5 .. 2 * pi]]
          botv = [Vector (radius * cos a) (radius * sin a) (-z) | a <- init [pi / 5, 3 * pi / 5 .. 2 * pi + pi / 5]]
          z = radius * (1 + sqrt 5) / 2
          top = [(north, a, b) | (a, b) <- zip topv (cycle topv)]
          bottom = [(south, a, b) | (a, b) <- zip botv (cycle botv)]
          middleTop = zip3 topv (cycle topv) botv
          middleBottom = zip3 botv (cycle botv) (cycle topv)
          cycle list = last list : init list

cut :: Float -> Face -> [Face]
cut r f@(fa, fb, fc) = [(fa, fab, fac), (fb, fbc, fac), (fc, fac, fbc)]
    where fab = hr *. normalize (fa +. fb)
          fbc = hr *. normalize (fb +. fc)
          fac = hr *. normalize (fa +. fc)
          hr = r / 2

triangulate :: Int -> Figure -> Figure
triangulate iter src = iterate subdivide src !! iter
    where subdivide = concatMap cutR :: Figure -> Figure
          cutR = cut (len . fst3 . head $ src) :: Face -> [Face]
          fst3 (a, _, _) = a