anchor.diff:     st_dump_triangles, &
      st_icosahedron, &
      st_iterate_level, &
  -   st_norm_p
  +   st_norm_p, &
  +   st_rotate
   
   
   contains
  
positive.code: end subroutine st_norm_p_single
  !-------------------------------------------------------------------------------
  
  end module sphere_tesselate
  !-------------------------------------------------------------------------------
  
negative.code:    if (allocated(t%t)) deallocate(t%t)
     allocate(t%t(t%nt))
  
     pole_in = .false.
     if (present(pole)) pole_in = pole
  
     ! First orientation is the division of the octahedron into segments with the golden
     ! ratio, phi
     if (.not.pole) then
        ! Vertices of triangles, unnormalised
        t%p(1)  = point( 0._rs,    phi,  1._rs)
        t%p(2)  = point( 0._rs,   -phi,  1._rs)
        t%p(3)  = point( 0._rs,    phi, -1._rs)
        t%p(4)  = point( 0._rs,   -phi, -1._rs)
        t%p(5)  = point( 1._rs,  0._rs,    phi)
        t%p(6)  = point(-1._rs,  0._rs,    phi)
        t%p(7)  = point( 1._rs,  0._rs,   -phi)
        t%p(8)  = point(-1._rs,  0._rs,   -phi)
        t%p(9)  = point(   phi,  1._rs,  0._rs)
        t%p(10) = point(  -phi,  1._rs,  0._rs)
        t%p(11) = point(   phi, -1._rs,  0._rs)
        t%p(12) = point(  -phi, -1._rs,  0._rs)
  
        ! Indices describing location of triangle vertices
        t%t( 1)%a =  2 ;  t%t( 1)%b =  4 ;  t%t( 1)%c = 11
        t%t( 2)%a =  5 ;  t%t( 2)%b =  2 ;  t%t( 2)%c = 11
        t%t( 3)%a =  9 ;  t%t( 3)%b =  5 ;  t%t( 3)%c = 11
        t%t( 4)%a =  7 ;  t%t( 4)%b =  9 ;  t%t( 4)%c = 11
        t%t( 5)%a = 11 ;  t%t( 5)%b =  7 ;  t%t( 5)%c =  4
        t%t( 6)%a =  4 ;  t%t( 6)%b =  2 ;  t%t( 6)%c = 12
        t%t( 7)%a =  6 ;  t%t( 7)%b = 12 ;  t%t( 7)%c =  2
        t%t( 8)%a =  2 ;  t%t( 8)%b =  5 ;  t%t( 8)%c =  6
        t%t( 9)%a =  1 ;  t%t( 9)%b =  6 ;  t%t( 9)%c =  5
        t%t(10)%a =  5 ;  t%t(10)%b =  9 ;  t%t(10)%c =  1
        t%t(11)%a =  3 ;  t%t(11)%b =  1 ;  t%t(11)%c =  9
        t%t(12)%a =  9 ;  t%t(12)%b =  7 ;  t%t(12)%c =  3
        t%t(13)%a =  8 ;  t%t(13)%b =  3 ;  t%t(13)%c =  7
        t%t(14)%a =  7 ;  t%t(14)%b =  4 ;  t%t(14)%c =  8
        t%t(15)%a = 12 ;  t%t(15)%b =  8 ;  t%t(15)%c =  4
        t%t(16)%a = 12 ;  t%t(16)%b =  6 ;  t%t(16)%c = 10
        t%t(17)%a =  6 ;  t%t(17)%b =  1 ;  t%t(17)%c = 10
        t%t(18)%a =  1 ;  t%t(18)%b =  3 ;  t%t(18)%c = 10
        t%t(19)%a =  3 ;  t%t(19)%b =  8 ;  t%t(19)%c = 10
        t%t(20)%a = 10 ;  t%t(20)%b = 12 ;  t%t(20)%c =  8
  
        ! Normalise onto unit sphere
        call st_norm_p(t%p)
  
     ! Second orientation places two points at +/- z, then places points around small
     ! circles at z = +/-a = 1/sqrt(5).
     ! Coordinates are from http://mathworld.wolfram.com/Icosahedron.html
     else
        t%p(1)  = point(           0._rs,            0._rs,  1._rs)
        t%p(2)  = point(      b*cos(pi5),       b*sin(pi5),      a)
        t%p(3)  = point(b*cos(3._rs*pi5), b*sin(3._rs*pi5),      a)
        t%p(4)  = point(b*cos(5._rs*pi5), b*sin(5._rs*pi5),      a)
        t%p(5)  = point(b*cos(7._rs*pi5), b*sin(7._rs*pi5),      a)
        t%p(6)  = point(b*cos(9._rs*pi5), b*sin(9._rs*pi5),      a)
        t%p(7)  = point(               b,            0._rs,     -a)
        t%p(8)  = point(b*cos(2._rs*pi5), b*sin(2._rs*pi5),     -a)
        t%p(9)  = point(b*cos(4._rs*pi5), b*sin(4._rs*pi5),     -a)
        t%p(10) = point(b*cos(6._rs*pi5), b*sin(6._rs*pi5),     -a)
        t%p(11) = point(b*cos(8._rs*pi5), b*sin(8._rs*pi5),     -a)
        t%p(12) = point(           0._rs,            0._rs, -1._rs)
  
        ! Indices of triangle vertices
        t%t(1)  = triangle( 6,  2,  1)
        t%t(2)  = triangle( 2,  3,  1)
        t%t(3)  = triangle( 3,  4,  1)
        t%t(4)  = triangle( 4,  5,  1)
        t%t(5)  = triangle( 5,  6,  1)
        t%t(6)  = triangle(11,  7,  6)
        t%t(7)  = triangle( 7,  2,  6)
        t%t(8)  = triangle( 7,  8,  2)
        t%t(9)  = triangle( 8,  3,  2)
        t%t(10) = triangle( 8,  9,  3)
        t%t(11) = triangle( 9,  4,  3)
        t%t(12) = triangle( 9, 10,  4)
        t%t(13) = triangle(10,  5,  4)
        t%t(14) = triangle(10, 11,  5)
        t%t(15) = triangle(11,  6,  5)
        t%t(16) = triangle(12,  7, 11)
        t%t(17) = triangle(12,  8,  7)
        t%t(18) = triangle(12,  9,  8)
        t%t(19) = triangle(12, 10,  9)
        t%t(20) = triangle(12, 11, 10)
     endif
  
  end function st_icosahedron
  !-------------------------------------------------------------------------------
  