anchor.diff:       real(rs) :: azimuth,lon1,lat1,lon2,lat2
        real(rs) :: rlon1,rlat1,rlon2,rlat2,conversion
        logical,optional :: degrees
  -     
  +
        conversion = 1._rs
        if (present(degrees)) then
           if (degrees) conversion = pi/180._rs
        endif
  -          
  +
        rlon1 = conversion*lon1  ;  rlon2 = conversion*lon2
        rlat1 = conversion*lat1  ;  rlat2 = conversion*lat2
  -     
  +
        azimuth = atan2(sin(rlon2-rlon1)*cos(rlat2) , &
                  cos(rlat1)*sin(rlat2) - sin(rlat1)*cos(rlat2)*cos(rlon2-rlon1) )
   
        if (azimuth < 0) then
           azimuth = azimuth+2._rs*pi
        endif
  -     
  +
        azimuth = azimuth / conversion
  -     
  +
   !     write(*,*)'Azimuth',azimuth
  -     
  +
        return
  -     
  +
      end function azimuth
   !==============================================================================
   
  
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:       logical,intent(in),optional :: degrees
        real(rs) :: conversion
        integer :: i,npoints
  
        ! Default to radians; convert from degrees if necessary
        conversion = 1._rs
        if (present(degrees)) then
  