      subroutine star_list(ndx,catalog,dec_start,dec_end,dechalf,
     $     rahalf,ra,ra_low_lim,ra_high_lim,lo_mag,hi_mag,magtype,
     $     density,list_ra,list_dec,list_mag,list_length,error)
C
C Description:
C 
C This routine examines the star database and returns a list of stars
C that meet the selection criteria.
C
C User Interface:
C
C Name         in/out type     structure  Meaning
C 
C ndx          in     integer  scalar     The logical unit number of the 
C                                         index file of the star database
C catalog      in     integer  scalar     The logical unit number of the 
C                                         catalog file of the star database
C dec_start    in     real*8   scalar     The low declination limit
C dec_end      in     real*8   scalar     The high declination limit
C dechalf      in     real*8   scalar     The half-heigt of the aperture
C rahalf       in     real*8   scalar     The half-width of the aperture
C ra           in     real*8   scalar     The RA of the center of the aperture
C ra_low_lim   in     real*8   scalar     The RA low limit
C ra_high_lim  in     real*8   scalar     The RA high limit
C lo_mag       in     real*8   scalar     The low magnitude limit
C hi_mag       in     real*8   scalar     The high magnitude limit
C magtype      in     char*1   scalar     The magnitude type
C                                         'r' is red
C                                         'b' is blue
C                                         'd' is the difference (b - r)
C density      in     real*8   scalar     The requested star density 
C                                         stars/sq-deg
C list_ra      out    integer  (100000)   The output list of star RA
C list_dec     out    integer  (100000)   The output list of star Decl
C list_mag     out    integer  (100000)   The output list of star magnitude
C list_length  out    integer  scalar     The number of stars in the list
C error        out    logical  scalar     Returned TRUE if an error occurs
C
C Errors:
C
C I/O errors are trapped in this routine. No other errors are trapped or should
C occur. If an error occurs, processing is immediately halted and the
C routine returns with 'error' set TRUE.
C
C Called routines:
C 
C err_msg
C
      implicit none
      include 'std_io.for'
      include 'include_.for'
      integer ndx,catalog,star_count,i,j,swapint4,
     $     stat,start_index,ll,ul,ira,idec,imag,recl,
     $     list_ra,list_dec,list_mag,ll2,ul2,list_length,ptr,
     $     londx,hindx,midndx,target,midra,mag,brightest,order
      real*8 ra_low_lim,ra_high_lim,dec_start,dec_end,start_ra,
     $     ra_inc,lo_mag,hi_mag,density,ra,area,dechalf,
     $     gamma,delta,arg,rahalf,rmag,bmag,cindex,rra,rdec
      logical error,new_lu,done
      character*(1) magtype
      dimension list_ra(maxlist),list_dec(maxlist),list_mag(maxlist),
     $     start_ra(ra_regions),start_index(ra_regions),
     $     star_count(ra_regions),brightest(maxlist),order(maxlist)
      ra_inc=twopi/dble(ra_regions)
      rewind(unit=ndx,err=406,iostat=stat)
      do j=1,ra_regions
         read(unit=ndx,fmt=1,err=402,iostat=stat)
     $        start_ra(j),start_index(j),star_count(j)
 1       format(f5.2,2i12)
      end do
      ll2=int(ra_low_lim/ra_inc)+1
      ul2=int(ra_high_lim/ra_inc)+1
      if (dabs(dble(ul2-1)*ra_inc - ra_high_lim) .lt. 1.0d-7) then
         ul2=ul2-1
      end if
      if (dble(ul2-1)*ra_inc .eq. ra_high_lim) then
         ul2=ul2-1
      end if
      new_lu=.true.
      recl=12
      do j=ll2,ul2
         ll=start_index(j)
         ul=start_index(j)+star_count(j)-1
C
C If the RA search region starts in the selected catalog RA segment,
C it is probably worthwhile to do a binary search to find the the
C first possible star in the catalog. It's index is returned in ll.
C Otherwise, search the whole segment.
C
         if (j .eq. ll2) then
            target=int(360000.0d0*ra_low_lim/pio180)
            londx=ll
            hindx=ul
            done=.false.
C
C An index difference of 680 is chosen in the binary search because 
C each i/o picks up about 680 stars. I think it is faster to do a linear 
C search through 680 stars than do any more I/O operations. Besides,
C we could only save 1 more io by making the search more exact.
C In a 4 MB catalog, 13 i/o ops are required to
C get within 680 stars. A linear search would require, on average,
C 250 i/o ops. If both upper and lower bounds are counted, 26 ios
C to 2500 ios or 1/10 as much i/o. In a 256 MB catalog, the difference
C is startling. 38 ios in a binary search to 16000 ios in a linear search. 
C
            do while (hindx-londx .gt. 680)
               midndx=(londx+hindx)/2
c               call read_buffer(new_lu,catalog,midndx,recl,
c     $            midra,idec,imag,error)
c               if (error) goto 400
               read(catalog,rec=midndx) midra,idec,imag
               midra=swapint4(midra)
               if (midra .gt. target) then
                  hindx=midndx
               elseif (midra .lt. target) then
                  londx=midndx
               elseif (midra .eq. target) then
                  londx=midndx
                  hindx=midndx
               end if
            end do
            ll=londx
         end if
C
C If the RA search region ends in the selected catalog RA segment,
C it is probably worthwhile to do a binary search to find the the
C last possible star in the catalog. It's index is returned in ul.
C Otherwise, search the whole segment.
C
         if (j .eq. ul2) then
            target=int(360000.0d0*ra_high_lim/pio180)
            londx=ll
            hindx=ul
            done=.false.
            do while (hindx-londx .gt. 680)
               midndx=(londx+hindx)/2
c               call read_buffer(new_lu,catalog,midndx,recl,
c     $            midra,idec,imag,error)
c               if (error) goto 400
               read(catalog,rec=midndx) midra,idec,imag
               midra=swapint4(midra)
               if (midra .gt. target) then
                  hindx=midndx
               elseif (midra .lt. target) then
                  londx=midndx
               elseif (midra .eq. target) then
                  londx=midndx
                  hindx=midndx
               end if
            end do
            ul=hindx
         end if
         do i=ll,ul
c            call read_buffer(new_lu,catalog,i,recl,ira,idec,imag,error)
c            if (error) goto 400
            read(catalog,rec=i)ira,idec,imag
            ira=swapint4(ira)
            idec=swapint4(idec)
            imag=swapint4(imag)
            if (error) goto 400
            rra=pio180*dble(ira)/360000.0d0
            rdec=pio180*dble(idec)/360000.0
            bmag=dble((mod(iabs(imag),1000000)/1000))/10.0d0
            rmag=dble((mod(iabs(imag),1000)))/10.0d0
            cindex=bmag-rmag
            if (((magtype .eq. 'r') .and.
     $           (rmag .ge. lo_mag) .and.
     $           (rmag .le. hi_mag)) .or.
     $           ((magtype .eq. 'b') .and.
     $           (bmag .ge. lo_mag) .and.
     $           (bmag .le. hi_mag)) .or.
     $           ((magtype .eq. 'd') .and.
     $           (cindex .ge. lo_mag) .and.
     $           (cindex .le. hi_mag))) then
C     
C Gamma is the difference in RA angles between the center line and the star.
C
               gamma=ra-rra
C
C This test checks whether or not the star is within the RA half-width
C
               delta=rdec
               arg=dasin(dsin(gamma)*dsin(delta))
               if (rahalf .ge. dabs(arg)) then
C
C This test checks whether the star lies within the decl limits when
C the star position is projected (at right angles) back onto the center 
C line.
C
                  if ((dabs(dabs(gamma)-pio2) .lt. 1.0d-8) .and.
     $                 (dabs(dabs(delta))-pio2 .lt. 1.0d-8)) then
                     arg=pio2
                  else
                     arg=datan2(dcos(gamma)*dsin(delta),dcos(delta))
                  end if
                  if ((dabs(arg) .le. dec_end) .and.
     $                 (dabs(arg) .ge. dec_start)) then
C          
C     The star has passed all the tests so record it.
C
                     list_length=list_length+1
                     if (list_length .gt. maxlist) then
                        error=.true.
                        call err_msg('star_list',
     $ 'Number of stars found exceeds storage allocation.')
                        goto 400
                     end if
                     list_ra(list_length)=ira
                     list_dec(list_length)=idec
                     list_mag(list_length)=imag
                  end if
               end if
            end if
         end do
      end do
C
C The following code segment verifies that the star density is
C less than 'density'. It orders the stars by brightness and
C keeps only enough of the brightest stars to match the 'density' 
C criteria. 
C
      area=rahalf*dechalf*4.0d0/pio180**2
      target=density*area
      if (target .lt. list_length) then
         do i=1,list_length
            mag=mod(iabs(list_mag(i)),1000000)
            brightest(i)=min(mag/1000,mod(mag,1000))
         end do
         order(1)=1
         call int_qsort(brightest,order,list_length)
         do i=target+1,list_length
            list_mag(order(i))=0
         end do
         ptr=1
         do i=1,list_length
            if (list_mag(i) .ne. 0) then
               list_mag(ptr)=list_mag(i)
               ptr=ptr+1
            end if
         end do
         list_length=target
      end if
C
C A way to achieve the required density and achieve a rather uniform
C distribution of stars is as follows:
C
C 1. Divide the aperture into regions that should contain exactly 1 star.
C 2. Sort the stars by declination
C 3. Loop over each band of declination that is one star wide.
C 3.1. Select the stars that are within the '1 star' declination band 
C 3.2. Sort them by RA
C 3.3. Loop over each band of RA that is one star wide.
C 3.3.1. Select the stars that meet each '1 star' ra band 
C 3.3.2. Sort by magnitude. 
C 3.5.3. Pick the brightest star.
C      end loop
C    end loop 
C If no stars are present then pick the brightest two stars from the
C next '1 star' RA band. (Generalize this step so that if a star 'debt'
C exists, the algorithm will continue pick stars from ra bands that
C contain excess stars.)
C
      goto 400
 402  call err_msg('star_list','Read error in index file.')
      error=.true.
      goto 400
 406  call err_msg('star_list','Error in rewinding index file')
      error=.true.
      goto 400
 410  call err_msg('star_list','Error in reading the catalog file.')
      error=.true.
 400  continue
      if (error) then
         call err_msg('star_list','Traceback')
      end if
      return
      end
