      subroutine make_iraf_list(list_ra,list_dec,list_mag,
     $     list_length,dec_center,ra_center,round_aper,
     $     focallength,pixelsize,ascending,sortfield,outfilename,
     $     radius,region_count,rahalf,dechalf,error)
C
C Description:
C
C This routine writes the list of "found" stars to an IRAF file/stdout
C in 'starlist' format. This format is suitable for IRAF 'mkobject' to
C simulate a CCD image.
C
C Programming Interface:
C
C name         in/out type       structure description
C
C list_ra      in     integer     (maxlist) The list of RAs for the stars 
C                                           (arcsec*100)
C list_dec     in     integer     (maxlist) The list of decls for the stars
C                                           (arcsec*100)
C list_mag     in     integer     (maxlist) The list of magnitudes for the 
C                                           stars (format is complex - see
C                                           other documentation for details
C list_length  in     integer     scalar    The number of stars in the lists
C ra_center    in     real*8      scalar    The RA center position (rad)
C dec_center   in     real*8      scalar    The decl center position (rad)  
C round_aper   in     logical     scalar    TRUE if the aperture is round
C focallength  in     real*8      scalar    Focal length of telescope. Any
C                                           unit is acceptable as long as
C                                           is the same unit used by 
C                                           pixelsize.
C pixelsize    in     real*8      scalar    Pixel size of CCD chip. Not used
C                                           in this routine. Same units as
C                                           focallength
C ascending    in     logical     scalar    If TRUE sort order is ascending.
C sortfield    in     character*(1) scalar  Selects field on which sort is 
C                                           performed
C                                           'R' is RA
C                                           'D' is Decl.
C                                           'r' is red
C                                           'b' is blue
C                                           'd' is radius
C                                           'i' is color index
C                                           'a' is angle
C outfilename  in     character*(64) scalar The prototype name of the output 
C                                           file
C radius       in     real*8      scalar    The radius of the round aperture
C                                           (rad).
C region_count in     integer     scalar    The number of regions that
C                                           have been examined so far.
C rahalf       in     real*8      scalar    The halfwidth of the CCD (rad)
C dechalf      in     real*8      scalar    The halfheight of the CCD (rad)
C error        out    logical     scalar    Set TRUE if an error occurs
C
C Errors:
C
C The usual I/O errors can occur in this routine. If such an error
C occurs, the routine will set 'error' TRUE and return immediately.
C
C Routines accessed:
C
C lnblnk
C left_justify
C int2_qsort
C dble_qsort
C int_qsort
C err_msg
C
      implicit none
      include 'std_io.for'
      include 'include_.for'
      logical error,round_aper,ascending
      real*8 ra,dec,ra_center,dec_center,list_distance,
     $     diff,list_angle,focallength,pixelsize,radius,
     $     rahalf,dechalf,xoffset,yoffset,
     $     alpha,delta,x,y,xmin,ymin,xmax,ymax
      integer list_ra,list_dec,list_mag,list_length,i,order,
     $     luout,ll,ul,step,region_count,lnblnk,mag
      integer*2 magr,magb,maxmag,minmag,cindex
      character*(1) sortfield,tab
      character*(4) digit
      character*(64) outfilename,name,left_justify,num
      dimension list_ra(maxlist),list_dec(maxlist),list_mag(maxlist),
     $     list_distance(maxlist),order(maxlist),list_angle(maxlist),
     $     magr(maxlist),magb(maxlist),cindex(maxlist)
      round_aper=.false.
      radius=0.0d0
      mag=mod(iabs(list_mag(1)),1000000)
      maxmag=mod(mag/1000,1000)
      minmag=maxmag
      do i=1,list_length
         ra=pio180*(dble(list_ra(i))/360000.0d0)
         dec=pio180*(dble(list_dec(i))/360000.0d0)
         diff=ra_center-ra
         list_distance(i)=dacos(dcos(diff)*
     $        dsin(dec_center)*dsin(dec)+
     $        dcos(dec_center)*dcos(dec))
         if (dabs(list_distance(i)) .lt. 1.0d-8) then
            list_angle(i)=0.0d0
         elseif ((dabs(dec_center-pi) .lt. 1.0d-8) .or.
     $           (dabs(dec_center) .lt. 1.0d-8)) then
            list_angle(i)=diff
         else 
            list_angle(i)=datan2(dsin(diff)*dsin(dec),
     $           (dcos(dec)-dcos(dec_center)*dcos(list_distance(i))))
         end if
         mag=mod(iabs(list_mag(i)),1000000)
         magr(i)=mod(mag,1000)
         magb(i)=mod(mag/1000,1000)
         cindex(i)=magb(i)-magr(i)
         if (magb(i) .gt. maxmag) then
            maxmag=magb(i)
         end if
         if (magb(i) .lt. minmag) then
            minmag=magb(i)
         end if
      end do
      if (list_length .eq. 1) then
         order(1)=1
      elseif (list_length .gt. 1) then
         if (sortfield .eq. 'b') then
            call int2_qsort(magb,order,list_length)
         elseif (sortfield .eq. 'd') then
            call dble_qsort(list_distance,order,list_length)
         elseif (sortfield .eq. 'a') then
            call dble_qsort(list_angle,order,list_length)
         elseif (sortfield .eq. 'R') then
            call int_qsort(list_ra,order,list_length)
         elseif (sortfield .eq. 'D') then
            call int_qsort(list_dec,order,list_length)
         elseif (sortfield .eq. 'i') then
            call int2_qsort(cindex,order,list_length)
         elseif (sortfield .eq. 'r') then
            call int2_qsort(magr,order,list_length)
         end if
      end if
      if (list_length .gt. 0) then
         if (outfilename .ne. ' ') then
            luout=204
            write(digit,43,err=403) region_count
 43         format(i4)
            digit=left_justify(digit)
            name=outfilename(1:lnblnk(outfilename))//
     $           digit(1:lnblnk(digit))
            open(unit=luout,
     $           file=name,
     $           err=401)
         else
            luout=stdout
         end if
 100     format(i32)
         tab=char(9)
         write(luout,1,err=401) tab,name
 1       format('begin',a1,a)
         write(luout,2,err=401) tab,tab,tab
 2       format(a1,'spatial',a1,a1,'uniform')
         call plate_position(ra_center,ra_center+rahalf,
     $        dec_center-pio2,dec_center+dechalf-pio2,
     $        xmax,ymax,pixelsize,focallength)
         xmax=1.05d0*(2.0d0*xmax+2.0d0)
         ymax=1.05d0*(2.0d0*ymax+2.0d0)
         xmin=0.0d0
         ymin=0.0d0
         write(luout,3,err=401) tab,tab,tab
 3       format(a1,'xmin',a1,a1,'0.')
         write(num,100,err=403) int(xmax)
         num=left_justify(num)
         write(luout,4,err=401) tab,tab,tab,num(1:lnblnk(num))
 4       format(a1,'xmax',a1,a1,a,'.')
         write(luout,5,err=401) tab,tab,tab
 5       format(a1,'ymin',a1,a1,'0.')
         write(num,100,err=403) int(ymax)
         num=left_justify(num)
         write(luout,6,err=401) tab,tab,tab,num(1:lnblnk(num))
 6       format(a1,'ymax',a1,a1,a,'.')
         write(luout,7,err=401) tab,tab,tab,tab,tab,tab,
     $        tab,tab,real(minmag)/10.0,tab,
     $        tab,tab,real(maxmag)/10.0,tab,tab,tab,tab,tab,tab,tab,
     $        tab,tab,tab,tab,tab,tab,tab,tab,tab
 7       format(a1,'luminosity',a1,'powlaw',/,
     $        a1,'power',a1,a1,'0.6',/,
     $        a1,'minmag',a1,a1,f4.1,/,
     $        a1,'maxmag',a1,a1,f4.1,/,
     $        a1,'nssample',a1,'100',/,
     $        a1,'sorder',a1,a1,'10',/,
     $        a1,'sseed',a1,a1,'1',/,
     $        a1,'nlsample',a1,'100',/,
     $        a1,'lorder',a1,a1,'10',/,
     $        a1,'lseed',a1,a1,'1')
         write(num,100,err=403) list_length
         num=left_justify(num)
         write(luout,8,err=401) tab,tab,tab,num(1:lnblnk(num))
 8       format(a1,'nstars',a1,a1,a)
         xoffset=dble(xmax-xmin+1)/2.0d0
         yoffset=dble(ymax-ymin+1)/2.0d0
         if (ascending) then
            ll=1
            ul=list_length
            step=1
         else
            ll=list_length
            ul=1
            step=-1
         end if
         do i=ll,ul,step
            alpha=pio180*list_ra(order(i))/360000.0d0
            delta=pio180*list_dec(order(i))/360000.0d0
            delta=delta-pio2
            call plate_position(ra_center,alpha,dec_center+pio2,
     $           delta,x,y,pixelsize,focallength)
            write(luout,9,err=401) tab,x+xoffset,y+yoffset,
     $           real(magb(order(i)))/10.0
 9          format(a1,f7.2,f9.2,f9.3)
         end do
      end if
      if (outfilename .ne. ' ') then
         close(unit=luout,err=402)
      end if
      goto 400
 401  error=.true.
      call err_msg('make_iraf_list','Error opening output file.')
      goto 400
 402  error=.true.
      call err_msg('make_iraf_list','Error closing output file.')
      goto 400
 403  error=.true.
      call err_msg('make_iraf_list','Error in format conversion.')
 400  continue
      if (error) then
         call err_msg('make_iraf_list','Traceback')
      end if
      return
      end
