anchor.diff:        read(20,'(1a)') comment
         if (comment == '#' .or. comment == '%' .or. comment == '>') cycle
         backspace(20)
  -      read(20,*) ray%x(npts), ray%y(npts), ray%z(npts)
  +      if (taup_input) then
  +         read(20,*) delta, r, lat, lon
  +         colat = to_rad*(90._rs - lat)
  +         colon = to_rad*lon
  +         ray%x(npts) = r * sin(colat) * cos(colon)
  +         ray%y(npts) = r * sin(colat) * sin(colon)
  +         ray%z(npts) = r * cos(colat)
  +      else
  +         read(20,*) ray%x(npts), ray%y(npts), ray%z(npts)
  +      endif
         npts = npts + 1
      enddo
      
  
positive.code:    contains
     
  !===============================================================================
     subroutine raypaths_read(file,ray)
  !===============================================================================
  !  Read in the coordinates of a ray from file.
  !  These ray files have x,y,z coordinates of ray on each line, with comment 
  !  lines containing #, % or > as the first character on the line.
     
     implicit none
     
  
negative.code:    type(raypath),intent(inout)  :: ray
     character(len=1)             :: comment
     integer                      :: npts,iostatus
     logical, intent(in), optional :: taup
     logical :: taup_input
     real(rs) :: delta,r,lon,lat,colon,colat
     
  !  Check if we have taup_input
     if (present(taup)) taup_input = taup
     
  !  Check that this is a 'new' ray
     if (allocated(ray%x) .or. allocated(ray%y) .or. allocated(ray%z)) then
  