Viewing contents of file '../idllib/contrib/buie/astcol.pro'
;+
; NAME:
;  astcol
; PURPOSE:
;  Collect astrometry observations for multiple objects
; DESCRIPTION:
;  
; CATEGORY:
;  Astrometry
; CALLING SEQUENCE:
;  astcol,OBSCODE=obscode,SAVEDIR=savedir
;
; INPUTS:
;
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORD INPUT PARAMETERS:
;  NONEWCODE - Flag, if set suppresses the generation of new object id codes.
;  OBSCODE - Observatory code for observations, default=688 (Lowell Obs.)
;  SAVEDIR - Directory where final astrometry is to be placed.  The default
;              value is /gryll/data1/buie/astrometry
;
; OUTPUTS:
;
; KEYWORD OUTPUT PARAMETERS:
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;    97/03/18 - Written by Marc W. Buie, Lowell Observatory
;    97/07/09 - Added the savedir keyword
;    97/07/25, MWB, slight change to the Bowell file format output.
;    98/06/22, MWB, added call to ASTLIST at end, also newcodes are taken from
;                 first available number, not last code in file.
;    98/11/04, MWB, now filters out magnitudes fainter than 80.0
;-
PRO astcol,OBSCODE=obscode,SAVEDIR=savedir,NONEWCODE=nonewcode

   if badpar(obscode,[0,1,2,3],0,CALLER='ASTCOL: (OBSCODE) ',default=688) then return
   if badpar(savedir,[0,7],0,CALLER='ASTCOL: (SAVEDIR) ', $
                default='/gryll/data1/buie/astrometry/') then return
   if badpar(nonewcode,[0,1,2,3],0,CALLER='ASTCOL: (NONEWCODE) ',default=0) then return

   savedir=addslash(savedir)

   ; Look for a cross-reference file.
   final = exists('lplast.xrft')

   IF final THEN BEGIN

      ; Read in the cross-reference file.
      openr,lun,'lplast.xrft',/get_lun
      nxrf=0
      a=''
      b=''
      WHILE not eof(lun) DO BEGIN
         readf,lun,a,b,format='(a8,4x,a)'
         a=strtrim(a,2)
         b=strtrim(b,2)
         IF nxrf eq 0 THEN BEGIN
            tmpid  = a
            realid = b
         ENDIF ELSE BEGIN
            tmpid  = [tmpid,a]
            realid = [realid,b]
         ENDELSE
         nxrf=nxrf+1
      ENDWHILE
      free_lun,lun

      if not nonewcode then begin
         ans=''
         read,'Do you want to setup final ID codes? (def=no) ',ans
         if strmid(ans,0,1) ne 'y' then nonewcode = 1
      endif

      IF not nonewcode THEN BEGIN
         ; Locate the id code file and load it
         ncodes = 0
         codetag = '?'
         IF exists(savedir+'newobj.dat') THEN BEGIN
            openr,lun,savedir+'newobj.dat',/get_lun
            line=''
            WHILE not eof(lun) DO BEGIN
               readf,lun,line,format='(a1)'
               ncodes = ncodes+1
            ENDWHILE
            point_lun,lun,0L
            tagid=strarr(ncodes)
            lines=strarr(ncodes)
            FOR i=0,ncodes-1 DO BEGIN
               readf,lun,line,format='(a)'
               lines[i]=line
               words=str_sep(strcompress(line),' ')
               tagid[i]=words[0]
            ENDFOR
            free_lun,lun
            i=0
            REPEAT BEGIN
               IF tagid[i] ne '' THEN BEGIN
                  j=0
                  REPEAT BEGIN
                     ok=valid_num(strmid(tagid[i],j,99),code,/integer)
                     IF ok and j ne 0 THEN $
                        codetag = strmid(tagid[i],0,j)
                     j=j+1
                  ENDREP UNTIL j eq strlen(tagid[i]) or ok
               ENDIF
               i=i+1
            ENDREP UNTIL codetag ne '?' or i eq ncodes
            IF codetag eq '?' THEN BEGIN
               print,'Sorry, this file does not have a clearly defined tag code.'
               print,'Unable to continue.'
            ENDIF
            tagnum=intarr(ncodes)
            FOR i=0,ncodes-1 DO BEGIN
               tagnum[i]=fix(strmid(tagid[i],strlen(codetag),99))
            ENDFOR
            print,'codetag is ',codetag
            idx=sort(tagnum)
            tagid  = tagid[idx]
            lines  = lines[idx]
            tagnum = tagnum[idx]
         ENDIF
      ENDIF
   ENDIF

   ; Get current directory, this will identify the date of observation.
   cd,'.',current=cdir
   pos=rstrpos(cdir,'/')
   date = strmid(cdir,pos+1,99)

   spawn,'ls *.ast',filelist
   nfiles=n_elements(filelist)

   IF filelist[0] eq '' THEN BEGIN
      print,'No astrometry files to collect.  Aborting.'
      return
   ENDIF

   line=''
   blanks='                   '
   tab = STRING( byte(9) )

   outast  = savedir+date+'.ast'
   outted  = savedir+date+'.ted'
   outinfo = savedir+date+'.info'

   spawn,'rm '+outted+'*'

   openw,last,outast,/get_lun
   openw,lted,outted,/get_lun

   nobs=0
   FOR i=0,nfiles-1 DO BEGIN
      obj=str_sep(filelist[i],'.')
      obj=obj[0]
      newobj=0
      IF strmid(obj,0,1) eq 'a' THEN obj=strmid(obj,1,99)
      obj=strupcase(obj)

      ; Process the data for this object
      openr,lun,filelist[i],/get_lun

      IF final THEN BEGIN
         origobj = obj
         z=where(origobj eq tmpid,count)
         IF count ne 0 THEN BEGIN
            obj = realid[z[0]]
         ENDIF ELSE IF not nonewcode THEN BEGIN
            ; Get the next valid code number
            diff = tagnum[1:ncodes-1]-tagnum[0:ncodes-2]
            zdiff = where(diff ne 1,count)
            if count eq 0 then begin
               newcode = tagnum[ncodes-1]+1
               print,' new end code ',newcode
            endif else begin
               newcode = tagnum[zdiff[0]]+1
               print,' new code ',newcode,' between ', $
                      tagnum[zdiff[0]],tagnum[zdiff[0]+1]
            endelse

            obj = strcompress(codetag+string(newcode),/remove_all)
            tag = strmid(origobj+blanks,0,8)
            line = '*   '+obj
            repwrite,'lplast.xrft',tag,tag+line

            tagnum = [tagnum,newcode]
            idx    = sort(tagnum)
            tagnum = tagnum[idx]
            ncodes = ncodes+1

            newobj=1

         ENDIF
         print,origobj,obj,format='(a,1x,a)'
      ENDIF

      ; Read through file and add to the ast output file and collect the
      ;    magnitudes
      j=0
      WHILE not eof(lun) DO BEGIN
         readf,lun,line,format='(a)'
         words=str_sep(line,' ')
         IF j eq 0 THEN mag = float(words[4]) else mag = [mag,float(words[4])]
         IF j eq 0 THEN firstfile = words[0]
         printf,last,line,' ',strtrim(string(obscode),2),' ',obj
         nobs=nobs+1
         j=j+1
      ENDWHILE

      if final and newobj then begin
         repwrite,savedir+'newobj.dat',obj,obj+tab+firstfile,/nosort
      endif

      ; Average magnitudes and then run through the file again, this time
      ;     adding to the Bowell format file.
      zg = where(mag lt 80.0,countgood)
      IF countgood gt 3 THEN $
         robomean,mag[zg],3.0,0.5,avgmag $
      ELSE if countgood eq 2 then $
         avgmag = mean(mag[zg]) $
      else if countgood eq 1 then $
         avgmag = mag[zg[0]] $
      else $
         avgmag = 99.9

      point_lun,lun,0L
      first=1
      WHILE not eof(lun) DO BEGIN
         readf,lun,line,format='(a)'
         words=str_sep(line,' ')
         jd=double(words[1])
         ra=raparse(words[2])
         dec=decparse(words[3])
         rastr,ra,3,ras
         decstr,dec,2,decs
         ras=repchar(ras,':')
         decs=repchar(decs,':')
         caldatm,jd,year,month,day,hour,minute,second
         day=float(day)+(float(hour)+(float(minute)+float(second)/60.0)/60.0)/24.0
         if first then begin
            magstr=string(avgmag,format='(f5.1)')+'R'
            npts=n_elements(mag)
            openw,lteda,outted+strcompress(npts,/remove_all),/get_lun,/append
         endif else begin
            magstr='      '
         endelse
         printf,lted,year,month,day,ras,decs,magstr,obj+blanks,obscode, $
            format='(i4,1x,i2.2,2x,f08.5,2x,a,2x,a,1x,a6,1x,a9,1x,i3)'
         printf,lteda,year,month,day,ras,decs,magstr,obj+blanks,obscode, $
            format='(i4,1x,i2.2,2x,f08.5,2x,a,2x,a,1x,a6,1x,a9,1x,i3)'
         first=0
      ENDWHILE

      free_lun,lun,lteda
   ENDFOR

   free_lun,last,lted

   print,'Astrometry saved: ',date,', ',strtrim(string(nfiles),2), $
         ' objects and ',strtrim(string(nobs),2),' measurements.'

   astlist,outast,outinfo

END