anchor.diff:           a = a*R  ; b = b*R
         endif
         f = (a-b)/a
  -      
  +
   !  Constants
         L = lon2 - lon1
         u1 = atan((1-f)*tan(lat1))
  
positive.code:       implicit none
        real(rs),intent(in) :: x,y,z
        real(rs)            :: xyz2radial(3), r
        
        r = sqrt(x**2 + y**2 + z**2)
        xyz2radial(1) = x / r
        xyz2radial(2) = y / r
        xyz2radial(3) = z / r
        
        return
     end function xyz2radial
  !===============================================================================
  
negative.code:          theta = theta_in
           phi = phi_in
        endif
  
         if (theta < 0.d0 .or. theta > pi ) then
            write(*,'(a)') 'Colatitude must be in range 0--pi (0--180deg).'
            stop
         endif
  
         x = r * sin(theta) * cos(phi)
         y = r * sin(theta) * sin(phi)
         z = r * cos(theta)
  
         return
  
     end subroutine sph2cart
  !==============================================================================
  
  