      subroutine find_ra_limits(ra,i,rahalf,spd_start,spd_end,
     $     ra_low_lim,ra_high_lim)
C
C Description:
C
C This routine determines the RA limits that are to be searched within
C a declination band. The limits are calculated by mapping the scan width
C into RA coordinates at a particular declination. The declination is
C chosen so that the most polar part of the band is used in the calculation.
C This choice will ensure that all stars that have a chance of selection
C are examined. 
C
C Programming Interface:
C
C Name         in/out type    structure description
C
C ra           in     real*8  scalar    The RA of the center of the scan
C i            in     integer scalar    The index of the decl region. Index
C                                       consecutive integers 0..23 represent
C                                       7.5 deg regions starting at south pole
C rahalf       in     real*8  scalar    The half width of the scan
C spd_start    in     real*8  scalar    Southernmost SPD
C spd_end      in     real*8  scalar    Northernmost SPD
C ra_low_lim   out    real*8  scalar    The lower limit value of the RA at 
C                                       the given declination
C ra_high_lim  out    real*8  scalar    The upper limit value of the RA at 
C                                       the given declination
C
C Errors:
C 
C None are trapped or reported. None should occur. There is an expression that
C contains a potential division by zero. I have attempted to detect
C when this will occur and produce the answer another way.
C
      implicit none
      include 'std_io.for'
      include 'include_.for'
      real*8  ra,rahalf,spd_start,spd_end,ra_low_lim,ra_high_lim,
     $     seg_width,spdlo,spdhi,angle
      integer i
      seg_width=pi/dble(dec_regions)
C
C Find the most polar region in a SPD catalog band and calculate
C the RA limits. If the band contains the spd limit, calculate
C the width at the limit rather than the polar-most region of the band. 
C
      spdlo=seg_width*float(i)
      spdhi=seg_width*float(i+1)
C
C In the Northern hemisphere, spdhi will give the widest search limits
C while in the South spdlo will give the widest. So handle the 
C appropriate case.
C
      if (spdlo .lt. pio2) then
C
C This is the Southern case. If the Southern limit of the search
C region is in this band, the search region will provide the limit.
C Otherwise the South edge of the band will provide the limit.
C
         if ((spd_start .gt. spdlo) .and.
     $        (spd_start .lt. spdhi)) then
            spdlo=spd_start
         end if
         if (spdlo .lt. 1.0d-8) then
            angle=pio2
         else
            angle=datan(dtan(rahalf)/dsin(spdlo))
         end if
      elseif (spdhi .gt. pio2) then
C
C This is the Northern case. If the Northern limit of the search
C region is in this band, the search region will provide the limit.
C Otherwise, the North edge of the band will provide the limit
C
         if ((spd_end .gt. spdlo) .and.
     $        (spd_end .lt. spdhi)) then
            spdhi=spd_end
         end if
         if (pi-spdhi .lt. 1.0d-8) then
            angle=pio2
         else
            angle=datan(dtan(rahalf)/dsin(spdhi))
         end if
      end if
      ra_low_lim=ra-angle
      ra_high_lim=ra+angle
      return
      end
