anchor.diff:           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      
  -      
  +
  +       return
  +
      end subroutine sph2cart
   !==============================================================================
   
  
positive.code:       enddo
        s = s + da
  !      write(*,*) s
        
  !  Test for in or out
        if (360._rs - abs(s) <= tol) then
           sph_poly_inout = .true.
        else
           sph_poly_inout = .false.
        endif
        
        return
     end function sph_poly_inout
  !===============================================================================
  
negative.code:          deg = degrees
           if (degrees) conversion = pi/180._rs
        endif
  
        tx = conversion*x
        ty = conversion*y
  
  !  Check for point on vertex
        if (any(x == px .and. y == py)) then
           write(0,'(a)') 'spherical_geometry: sph_poly_inout: point is on vertex.'
           stop
        endif
  
  !  Check for point on poles
        if (y == 90. .or. y == -90.) then
           write(0,'(a)') 'spherical_geometry: sph_poly_inout: point is on one of the poles.'
           stop
        endif
  
  !  Loop over sides and calculate sum of angles.  If ~360, inside.  If ~0, outside
        n = size(px)
        s = 0.
  