;+
; NAME:
;  astlinks
; PURPOSE:
;  Scan for linkages among a collection of asteroid astrometric measurements
; DESCRIPTION:
;
; CATEGORY:
;  Astrometry
; CALLING SEQUENCE:
;  astlinks,otherdir
; INPUTS:
;
; OPTIONAL INPUT PARAMETERS:
;  otherdir - If specified, the local .ast files are checked against this
;               directory too.
; KEYWORD INPUT PARAMETERS:
;  NOLOCAL - Flag, if set will suppress checking for linkages within the
;               astrometry from the current directory.
;
; OUTPUTS:
;
; KEYWORD OUTPUT PARAMETERS:
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;  98/01/13, Written by Marc W. Buie, Lowell Observatory
;
;-
PRO astlinks,in_otherdir,NOLOCAL=nolocal,ONEFILE=onefile

   if badpar(in_otherdir,[0,7],0,caller='REDUCTOR: (DUMPOBJ) ', $
                                 default='.') then return
   if badpar(nolocal,[0,1,2,3],0,caller='REDUCTOR: (NOLOCAL) ', $
                                 default=0) then return
   if badpar(onefile,[0,7],[0,1],caller='REDUCTOR: (ONEFILE) ', $
                                 default='') then return

   if onefile[0] ne '' then nolocal=1

   otherdir=in_otherdir

   if onefile[0] eq '' then begin
      localfn = findfile('*.ast',count=nlocal)
      IF nlocal eq 0 THEN BEGIN
         print,'No astrometry files found in current directory.  Aborting.'
         return
      ENDIF
   endif else begin
      localfn = onefile
      nlocal  = n_elements(localfn)
   endelse

   logfile='links.log'
   openw,lunlog,logfile,/get_lun

   print,nlocal,' astrometry files in current directory.'
   printf,lunlog,nlocal,' astrometry files in current directory.'

   blanks='                                             '

   IF otherdir ne '.' THEN BEGIN
      otherfn = findfile(addslash(otherdir)+'*.ast',count=nother)
      IF nother eq 0 THEN otherdir='.'
   ENDIF

   if otherdir ne '.' then begin
      print,nother,' astrometry files from ',otherdir
      printf,lunlog,nother,' astrometry files from ',otherdir
   endif

   cr = string("15b)  ;"
   rt = string("12b)  ;"
   form='($,a,a15,1x,i5)'

   nlinks=0
   FOR i=0,nlocal-1 DO BEGIN

      print,form=form,cr,localfn[i]+blanks,nlocal-i
      rdrawast,localfn[i],fn0,jd0,ra0,dec0,mag0,nobs0

      ; Compute rate of first object
      objdir0 = atan(dec0[nobs0-1]-dec0[0],ra0[0]-ra0[nobs0-1])*!radeg
      hmotobj0=angsep(ra0[0],dec0[0],ra0[nobs0-1],dec0[nobs0-1]) $
           / (abs(jd0[nobs0-1]-jd0[0]) * 24.0 ) * !radeg * 3600.0

      IF otherdir eq '.' THEN BEGIN
         cklist = localfn[i:*]
      ENDIF ELSE if nolocal then BEGIN
         cklist = otherfn
      ENDIF ELSE BEGIN
         cklist = [localfn[i:*],otherfn]
      ENDELSE

      FOR j=0,n_elements(cklist)-1 DO BEGIN

         ckobj=cklist[j]
         loc=rstrpos(ckobj,'/')
         IF loc ne -1 THEN BEGIN
            ckobj = strmid(ckobj,loc+1,99)
         ENDIF

         IF ckobj ne localfn[i] THEN BEGIN

            rdrawast,cklist[j],fn1,jd1,ra1,dec1,mag1,nobs1
            ; Compute rate of second object
            objdir1 = atan(dec1[nobs1-1]-dec1[0],ra1[0]-ra1[nobs1-1])*!radeg
            hmotobj1=angsep(ra1[0],dec1[0],ra1[nobs1-1],dec1[nobs1-1]) $
                 / (abs(jd1[nobs1-1]-jd1[0]) * 24.0 ) * !radeg * 3600.0

            ; Compute rate assuming both 1st measures are same object
            IF jd0[0] eq jd1[0] THEN BEGIN
               hmotobj2=9999.0
            ENDIF ELSE BEGIN
               hmotobj2=angsep(ra0[0],dec0[0],ra1[0],dec1[0]) $
                    / (abs(jd1[0]-jd0[0]) * 24.0 ) * !radeg * 3600.0
            ENDELSE

            closerate = abs(hmotobj0-hmotobj1) lt 2.0 and $
                        abs(hmotobj0-hmotobj2) lt 2.0

            closedir  = abs(objdir0-objdir1) lt 5.0

            IF not closedir and objdir0 lt 0.0 THEN $
               closedir = abs((objdir0+360.0)-objdir1) lt 5.0

            IF closerate and closedir THEN BEGIN

               if nlinks eq 0 then begin
                  obj0   = localfn[i]
                  obj1   = cklist[j]
                  objcnt = 1
                  nlinks = 1
               endif else begin
                  zm=where(localfn[i] eq obj0)
                  if zm[0] eq -1 then begin
                     obj0 = [obj0,localfn[i]]
                     obj1 = [obj1,cklist[j]]
                     objcnt = [objcnt,1]
                     nlinks = nlinks + 1
                  endif else begin
                     obj1[zm[0]] = obj1[zm[0]] + ' ' + cklist[j]
                     objcnt[zm[0]] = objcnt[zm[0]]+1
                  endelse
               endelse

               print,nlinks,minmax(objcnt),n_elements(objcnt)
               print,localfn[i]+blanks,hmotobj0,objdir0,minmax(mag0), $
                  format='(2x,"[",a30,1x,f6.1,"aph ",f6.1,"d]",2(1x,f4.1))'
               print,cklist[j]+blanks,hmotobj1,objdir1,minmax(mag0),hmotobj2, $
                  format='(2x,"[",a30,1x,f6.1,"aph ",f6.1,"d]",2(1x,f4.1),2x,f6.1,"aph")'
               print,''

               printf,lunlog,''
               printf,lunlog,localfn[i]+blanks,hmotobj0,objdir0,minmax(mag0), $
                  format='(2x,"[",a30,1x,f6.1,"aph ",f6.1,"d]",2(1x,f4.1))'
               printf,lunlog,cklist[j]+blanks,hmotobj1,objdir1,minmax(mag0),hmotobj2, $
                  format='(2x,"[",a30,1x,f6.1,"aph ",f6.1,"d]",2(1x,f4.1),2x,f6.1,"aph")'
               printf,lunlog,''
            ENDIF

         ENDIF

      ENDFOR

   ENDFOR

   print,''
   free_lun,lunlog

   if nlinks gt 0 then begin
      xreffile='xref.log'
      openw,lunxref,xreffile,/get_lun
      z=where(objcnt eq 1,count1)
      if count1 gt 0 then begin
         z = z[sort(obj0[z])]
         for i=0,count1-1 do $
            printf,lunxref,obj0[z[i]],' ',obj1[z[i]]
      endif
      z=where(objcnt ne 1,count2)
      if count1 ne 0 and count2 ne 0 then printf,lunxref,'-------------------'
      if count2 gt 0 then begin
         z = z[sort(obj0[z])]
         for i=0,count2-1 do $
            printf,lunxref,obj0[z[i]],' ',obj1[z[i]]
      endif
      free_lun,lunxref
   endif

END