;+ ; NAME: ; kbosum ; PURPOSE: ; Create a summary listing of KBOs in all directories of DES data ; DESCRIPTION: ; ; CATEGORY: ; Astronomy ; CALLING SEQUENCE: ; kbosum ; INPUTS: ; OPTIONAL INPUT PARAMETERS: ; KEYWORD INPUT PARAMETERS: ; ASTPATH - path to astrometry data (see kbolist.pro for default) ; REDPATH - path to reduced data, default is an empty string ('') ; RAWPATH - path to reduced data, default is a vector: ; ['/gryll/data6/buie/rawfits', '/gryll/data7/buie/rawfits'] ; OUTPUTS: ; outfile - output file for kbo list, summary.dat ; KEYWORD OUTPUT PARAMETERS: ; COMMON BLOCKS: ; SIDE EFFECTS: ; RESTRICTIONS: ; PROCEDURE: ; MODIFICATION HISTORY: ; 01/01/04, Written by Susan Kern ; 01/04/06, DBT - Fixed bug in sort, Fixed bug in KBOName, cleand up code ;- pro kbosum, ASTPATH=astpath, REDPATH=redpath, RAWPATH=rawpath ; check to see that the parameters are valid if badpar(astpath,[0,7],0,caller='kbosum:(ASTPATH) ', $ default='/gryll/data1/buie/astrometry') then return if badpar(redpath,[0,7],0,caller='kbosum:(REDPATH) ', $ default='') then return if badpar(rawpath,[0,7],0,caller='kbosum:(RAWPATH) ', $ default=['/gryll/data6/buie/rawfits','/gryll/data7/buie/rawfits'])$ then return logfile = 'kbosum.log' logerror,logfile,'Generating master KBO list',caller='kbosum' ; get just the data directory from the entire path fdir=findfile(redpath,COUNT=nfile) ; generate a master list of all observed objects moving less ; than 10 arcsec/hr using.info files for each of the kbo survey directories nelements=0 for i=0, nfile-1 do begin kbolist,fdir[i],dirname0,rawcode0,firstnight0,rmag0,logfile=logfile, $ rmagerr0,rate0,dir0,dr0,dt0,elong0,retcount,ASTPATH=astpath if retcount gt 0 then begin if nelements eq 0 then begin dirname=dirname0 rawcode=rawcode0 firstnight=firstnight0 rmag=rmag0 rmagerr=rmagerr0 rate=rate0 dir=dir0 dr=dr0 dt=dt0 elong=elong0 endif else begin dirname=[dirname,dirname0] rawcode=[rawcode,rawcode0] firstnight=[firstnight,firstnight0] rmag=[rmag,rmag0] rmagerr=[rmagerr,rmagerr0] rate=[rate,rate0] dir=[dir,dir0] dr=[dr,dr0] dt=[dt,dt0] elong=[elong,elong0] endelse nelements=nelements+retcount logerror,logfile,caller='kbosum', $ 'Directory '+fdir[i]+', new='+strn(retcount)+ $ ', cumulative count '+strn(nelements) endif endfor if nelements eq 0 then begin logerror,logfile,caller='kbosum','No entries for the KBO list' return endif ; split up the rawcode list into three options ; looker id (tmpcode), local id (lcode), and mpc designation (desig) tmpcode=strarr(nelements) lcode=strarr(nelements) desig=strarr(nelements) firstfield=firstnight rawcode=strupcase(rawcode) for i=0,nelements-1 do begin c0=strmid(rawcode[i],0,1) c1=strmid(rawcode[i],1,1) c2=strmid(rawcode[i],2,1) c3=strmid(rawcode[i],3,1) c4=strmid(rawcode[i],4,1) if c0 ge 'A' and c0 le 'Z' and $ c1 ge 'A' and c1 le 'Z' then begin lcode[i]=rawcode[i] endif else if strlen(rawcode[i]) ge 7 and $ c0 ge '0' and c0 le '9' and $ c1 ge '0' and c1 le '9' and $ c2 ge '0' and c2 le '9' and $ c3 ge '0' and c3 le '9' and $ c4 ge '0' and c4 le '9' then begin tmpcode[i]=rawcode[i] endif else begin desig[i]=rawcode[i] endelse endfor ; generate a list of local codes wih cross-reference information logerror,logfile,caller='kbosum', $ 'Generating list of local codes' rdobjdes,'newobj.dat',ncodes,objcode,field,xref,xsec,COMMENTS=comments ; determine if any local codes have designation information available for i=0,nelements-1 do begin if lcode[i] ne '' then begin z=where(lcode[i] eq objcode,count) if count eq 1 then begin desig[i]=xref[z[0]] firstfield[i]=field[z[0]] endif endif else begin z=where(desig[i] eq xref,count) if count eq 1 then begin lcode[i]=objcode[z[0]] firstfield[i]=field[z[0]] endif endelse endfor ; Find h magnitudes and the official designation ; for each of the objects known to have a mpc designation. ; Warning - Enter at your own risk logerror,logfile,caller='kbosum', $ 'Finding H mag and official designation' ; Start up Larry's name and info fetch program spawn,'getinfo',unit=pipe ; Start up Larry's program to find designations for objects spawn,'getmfad',unit=pipe2 hmag=fltarr(nelements) for i=0,nelements-1 do begin if strlen(desig[i]) eq 0 then begin hmag[i]=99999999 continue endif kboname='E'+desig[i] printf,pipe,kboname object='' readf,pipe,object,h,g,format='(a32,1x,f6.2,1x,f5.2)' object=strtrim(object,2) object=strcompress(object) words=strsplit(object,/EXTRACT) object=words[0] if strmid(object,0,10) ne 'XXXXXXXXXX' then begin desig[i]=object hmag[i]=h continue endif kboname='A'+desig[i] printf,pipe,kboname readf,pipe,object,h,g,format='(a32,1x,f6.2,1x,f5.2)' object=strtrim(object,2) object=strcompress(object) words=strsplit(object,/EXTRACT) object=words[0] if strmid(object,0,10) ne 'XXXXXXXXXX' then begin desig[i]=object hmag[i]=h continue endif printf,pipe2,kboname readf,pipe2,object,format='(a)' object=strtrim(object,2) if object ne 'X' then begin kboname='A'+object printf,pipe,kboname readf,pipe,object,h,g,format='(a32,1x,f6.2,1x,f5.2)' object=strtrim(object,2) object=strcompress(object) words=strsplit(object,/EXTRACT) ;DBT - Bad Fix? should be words[0]??? object=words[1] desig[i]=object hmag[i]=h continue endif kboname='E'+desig[i] printf,pipe2,kboname readf,pipe2,object,format='(a)' object=strtrim(object,2) if object eq 'X' then begin desig[i]='***' hmag[i]=99999999 continue endif kboname='A'+object printf,pipe,kboname readf,pipe,object,h,g,format='(a32,1x,f6.2,1x,f5.2)' object=strtrim(object,2) object=strcompress(object) words=strsplit(object,/EXTRACT) object=words[0] desig[i]=object hmag[i]=h endfor free_lun,pipe free_lun,pipe2 ; get looker ids by reading lplast.xrft files alldir=strmid(firstfield,0,6) dirlist=alldir[uniq(alldir,sort(alldir))] ndir=n_elements(dirlist) ; Look for a cross-reference file. logerror,logfile,caller='kbosum', $ 'Reading lplast.xrft files' for i=0,ndir-1 do begin tmpdir=redpath+addslash(dirlist[i]) rdlplast, lookerid, tmpid, ncount, /trim, path=tmpdir z=where(dirlist[i] eq alldir,count) for j=0,count-1 do begin ;use lcode or desig to find tmpcode if tmpcode[z[j]] eq '' then begin zz=where(lcode[z[j]] eq tmpid or desig[z[j]] eq tmpid,count) if count eq 1 then begin tmpcode[z[j]]=lookerid[zz[0]] endif endif else begin ;; Else look up lcode or desig form lookerid zz=where(tmpcode[z[j]] eq lookerid,count) if count ne 1 then begin logerror,logfile,caller='kbosum', $ '*** lookerid '+tmpcode[z[j]]+' not found in file '+tmpdir continue endif if tmpcode[z[j]] eq tmpid[zz[0]] then begin logerror,logfile,caller='kbosum', $ '*** invalid lookerid '+tmpcode[z[j]]+' in file '+tmpdir continue endif c0=strmid(tmpid[zz[0]],0,1) c1=strmid(tmpid[zz[0]],1,1) if c0 ge 'A' and c0 le 'Z' and $ c1 ge 'A' and c1 le 'Z' then begin if strlen(tmpid[zz[0]]) ne 0 then lcode[z[j]]=tmpid[zz[0]] ;; look up desig from new lcode dtz=where(lcode[z[j]] eq objcode,count) if count eq 1 then begin if strlen(xref[dtz[0]]) ne 0 then $ desig[z[j]]=strmid(xref[dtz[0]],2) endif endif else begin if strlen(tmpid[zz[0]]) ne 0 then desig[z[j]]=tmpid[zz[0]] endelse endelse endfor endfor ;; TBD ;;Sort routine is BROKEN!!! if 0 then begin ;Routine to sort data in file and to get rid of duplicates. logerror,logfile,caller='kbosum','Sorting data' ;sort on looker code only newele=0 sortx=where(strlen(lcode) eq '' and strlen(desig) eq '',count) if count ne 0 then begin tmpcode0=tmpcode[sortx] lcode0=lcode[sortx] desig0=desig[sortx] firstfield0=firstfield[sortx] rawcode0=rawcode[sortx] rmag0=rmag[sortx] rmagerr0=rmagerr[sortx] rate0=rate[sortx] dir0=dir[sortx] dr0=dr[sortx] dt0=dt[sortx] elong0=elong[sortx] hmag0=hmag[sortx] x=sort(tmpcode0) if newele eq 0 then begin ntmpcode=tmpcode0[x] nlcode=lcode0[x] ndesig=desig0[x] nfirstfield=firstfield0[x] nrawcode=rawcode0[x] nrmag=rmag0[x] nrmagerr=rmagerr0[x] nrate=rate0[x] ndir=dir0[x] ndr=dr0[x] ndt=dt0[x] nelong=elong0[x] nhmag=hmag0[x] endif else begin ntmpcode=[ntmpcode,tmpcode0[x]] nlcode=[nlcode,lcode0[x]] ndesig=[ndesig,desig0[x]] nfirstfield=[nfirstfield,firstfield0[x]] nrawcode=[nrawcode,rawcode0[x]] nrmag=[nrmag,rmag0[x]] nrmagerr=[nrmagerr,rmagerr0[x]] nrate=[nrate,rate0[x]] ndir=[ndir,dir0[x]] ndr=[ndr,dr0[x]] ndt=[ndt,dt0[x]] nelong=[nelong,elong0[x]] nhmag=[nhmag,hmag0[x]] endelse newele=n_elements(nlcode) endif ;sort on local code x=where(strlen(lcode) ne '' and strlen(desig) eq '',count) if count ne 0 then begin tmpcode0=tmpcode[x] lcode0=lcode[x] desig0=desig[x] firstfield0=firstfield[x] rawcode0=rawcode[x] rmag0=rmag[x] rmagerr0=rmagerr[x] rate0=rate[x] dir0=dir[x] dr0=dr[x] dt0=dt[x] elong0=elong[x] hmag0=hmag[x] tmplcode=long(strmid(lcode0,2)) x=sort(tmplcode) z=where((lcode0[x])[0:count-2] ne (lcode0[x])[1:count-1]) z=[0,z+1] if newele eq 0 then begin ntmpcode=tmpcode0[x[z]] nlcode=lcode0[x[z]] ndesig=desig0[x[z]] nfirstfield=firstfield0[x[z]] nrawcode=rawcode0[x[z]] nrmag=rmag0[x[z]] nrmagerr=rmagerr0[x[z]] nrate=rate0[x[z]] ndir=dir0[x[z]] ndr=dr0[x[z]] ndt=dt0[x[z]] nelong=elong0[x[z]] nhmag=hmag0[x[z]] endif else begin ntmpcode=[ntmpcode,tmpcode0[x[z]]] nlcode=[nlcode,lcode0[x[z]]] ndesig=[ndesig,desig0[x[z]]] nfirstfield=[nfirstfield,firstfield0[x[z]]] nrawcode=[nrawcode,rawcode0[x[z]]] nrmag=[nrmag,rmag0[x[z]]] nrmagerr=[nrmagerr,rmagerr0[x[z]]] nrate=[nrate,rate0[x[z]]] ndir=[ndir,dir0[x[z]]] ndr=[ndr,dr0[x[z]]] ndt=[ndt,dt0[x[z]]] nelong=[nelong,elong0[x[z]]] nhmag=[nhmag,hmag0[x[z]]] endelse nelements=n_elements(nlcode) endif ; sort on designation newele=0 x=where(strlen(desig) ne '',count) if count ne 0 then begin tmplcode=long(strmid(lcode0[x],2)) y=sort(tmplcode) tmpcode0=tmpcode[x[y]] lcode0=lcode[x[y]] desig0=desig[x[y]] firstfield0=firstfield[x[y]] rawcode0=rawcode[x[y]] rmag0=rmag[x[y]] rmagerr0=rmagerr[x[y]] rate0=rate[x[y]] dir0=dir[x[y]] dr0=dr[x[y]] dt0=dt[x[y]] elong0=elong[x[y]] hmag0=hmag[x[y]] pack=mpcdcvt(desig0) x=where(pack eq desig0,newcount) if newcount ne 0 then begin z=sort(long(pack[x])) if newele eq 0 then begin tmpcode1=tmpcode0[x[z]] lcode1=lcode0[x[z]] desig1=desig0[x[z]] firstfield1=firstfield0[x[z]] rawcode1=rawcode0[x[z]] rmag1=rmag0[x[z]] rmagerr1=rmagerr0[x[z]] rate1=rate0[x[z]] dir1=dir0[x[z]] dr1=dr0[x[z]] dt1=dt0[x[z]] elong1=elong0[x[z]] hmag1=hmag0[x[z]] endif else begin tmpcode1=[tmpcode1,tmpcode0[x[z]]] lcode1=[lcode1,lcode0[x[z]]] desig1=[desig1,desig0[x[z]]] firstfield1=[firstfield1,firstfield0[x[z]]] rawcode1=[rawcode1,rawcode0[x[z]]] rmag1=[rmag1,rmag0[x[z]]] rmagerr1=[rmagerr1,rmagerr0[x[z]]] rate1=[rate1,rate0[x[z]]] dir1=[dir1,dir0[x[z]]] dr1=[dr1,dr0[x[z]]] dt1=[dt1,dt0[x[z]]] elong1=[elong1,elong0[x[z]]] hmag1=[hmag1,hmag0[x[z]]] endelse newele=n_elements(tmpcode1) endif x=where(pack ne desig0, newcount) if newcount ne 0 then begin z=sort(pack[x]) if newele eq 0 then begin tmpcode1=tmpcode0[x[z]] lcode1=lcode0[x[z]] desig1=desig0[x[z]] firstfield1=firstfield0[x[z]] rawcode1=rawcode0[x[z]] rmag1=rmag0[x[z]] rmagerr1=rmagerr0[x[z]] rate1=rate0[x[z]] dir1=dir0[x[z]] dr1=dr0[x[z]] dt1=dt0[x[z]] elong1=elong0[x[z]] hmag1=hmag0[x[z]] endif else begin tmpcode1=[tmpcode1,tmpcode0[x[z]]] lcode1=[lcode1,lcode0[x[z]]] desig1=[desig1,desig0[x[z]]] firstfield1=[firstfield1,firstfield0[x[z]]] rawcode1=[rawcode1,rawcode0[x[z]]] rmag1=[rmag1,rmag0[x[z]]] rmagerr1=[rmagerr1,rmagerr0[x[z]]] rate1=[rate1,rate0[x[z]]] dir1=[dir1,dir0[x[z]]] dr1=[dr1,dr0[x[z]]] dt1=[dt1,dt0[x[z]]] elong1=[elong1,elong0[x[z]]] hmag1=[hmag1,hmag0[x[z]]] endelse newele=n_elements(tmpcode1) endif z=where(desig1[0:newele-2] ne desig1[1:newele-1]) z=[0,z+1] if nelements eq 0 then begin ntmpcode=[ntmpcode,tmpcode1[z]] nlcode=[nlcode,lcode1[z]] ndesig=[ndesig,desig1[z]] nfirstfield=[nfirstfield,firstfield1[z]] nrawcode=[nrawcode,rawcode1[z]] nrmag=[nrmag,rmag1[z]] nrmagerr=[nrmagerr,rmagerr1[z]] nrate=[nrate,rate1[z]] ndir=[ndir,dir1[z]] ndr=[ndr,dr1[z]] ndt=[ndt,dt1[z]] nelong=[nelong,elong1[z]] nhmag=[nhmag,hmag1[z]] endif else begin ntmpcode=[ntmpcode,tmpcode1[z]] nlcode=[nlcode,lcode1[z]] ndesig=[ndesig,desig1[z]] nfirstfield=[nfirstfield,firstfield1[z]] nrawcode=[nrawcode,rawcode1[z]] nrmag=[nrmag,rmag1[z]] nrmagerr=[nrmagerr,rmagerr1[z]] nrate=[nrate,rate1[z]] ndir=[ndir,dir1[z]] ndr=[ndr,dr1[z]] ndt=[ndt,dt1[z]] nelong=[nelong,elong1[z]] nhmag=[nhmag,hmag1[z]] endelse nelements=n_elements(nlcode) endif endif else begin ntmpcode=tmpcode nlcode=lcode ndesig=desig nfirstfield=firstfield nrawcode=rawcode nrmag=rmag nrmagerr=rmagerr nrate=rate ndir=dir ndr=dr ndt=dt nelong=elong nhmag=hmag nelements=n_elements(nlcode) endelse ;read information from .ast file in reduced directory for each object ; to get discovery jdt and ra/dec logerror,logfile,caller='kbosum', $ 'Reading .ast files' nra = strarr(nelements) ndec = strarr(nelements) njd = strarr(nelements) tmpjd = dblarr(nelements) dirlist=strmid(nfirstfield,0,6) for i=0,nelements-1 do begin ; Fill in KBOname if strlen(ndesig[i]) ne 0 then begin kboname='A'+ndesig[i] endif else if strlen(nlcode[i]) ne 0 then begin kboname='A'+nlcode[i] endif else if strlen(ntmpcode[i]) ne 0 then begin kboname='A'+ntmpcode[i] endif else begin kboname='*undefined*' endelse tmpdir=redpath+addslash(dirlist[i]) if exists (tmpdir) then begin tmpfile=tmpdir+ntmpcode[i]+'.ast' if exists(tmpfile) then begin ; read in astrometry file rdrawast,tmpfile,fn,jd,ra,dec,mag,nobs rastr,ra,1,ra0 decstr,dec,1,dec0 jdstr,jd,1,jd0 tmpjd[i]=jd[0] njd[i]=jd0[0] nra[i]=ra0[0] ndec[i]=dec0[0] endif else begin ; initialize fields to default values logerror,logfile,caller='kbosum', $ '*** error astfile='+tmpfile+' not found, KBOName='+kboname nra[i]='**:**:**.*' ndec[i]='***:**:**.*' njd[i]='*** ***' tmpjd[i]=systime(/julian) endelse endif else begin logerror,logfile,caller='kbosum', $ '*** error astdir='+tmpdir+' not found, KBOname='+kboname nra[i]='**:**:**.*' ndec[i]='***:**:**.*' njd[i]='*** ***' tmpjd[i]=systime(/julian) endelse endfor ; get orital elements for the objects using Ted's routines. logerror,logfile,caller='kbosum', $ 'Getting orbital elements for the objects' inc = dblarr(nelements) e = dblarr(nelements) a = dblarr(nelements) epherr = dblarr(nelements) pipe = '' for i=0,nelements-1 do begin ; Fill in KBOname if strlen(ndesig[i]) ne 0 then begin kboname='A'+ndesig[i] endif else if strlen(nlcode[i]) ne 0 then begin kboname='A'+nlcode[i] endif else if strlen(ntmpcode[i]) ne 0 then begin kboname='A'+ntmpcode[i] endif else begin kboname='' logerror,logfile,caller='kbosum', $ '***** KBONAME UNDEFINED *****' endelse ; If no name dont call if strlen(kboname) ne 0 then begin ;;print, 'EPHEM, ',tmpjd[i],',500, 23, ',kboname ephem,tmpjd[i],500,23,kboname,ephem0,PIPE=pipe inc[i]=ephem0[11] e[i]=ephem0[12] a[i]=ephem0[13] epherr[i]=ephem0[20] ;;22 if inc[i] eq 0.0 then inc[i]=999999 if a[i] eq 0.0 then a[i]=999999 if e[i] eq 0.0 then e[i]=999999 if epherr[i] eq 0.0 then epherr[i]=999999 endif else begin inc[i] = 999999 e[i] = 999999 a[i] = 999999 epherr[i] = 999999 endelse endfor ;free our cached pipe free_lun,pipe pipe = '' logerror,logfile,caller='kbosum', $ 'Getting FieldID and Filter info' ; Get field id and filter infomation for each object alldir=strmid(nfirstfield,0,6) numdir=n_elements(alldir) allfiles=strmid(nfirstfield,0,10) ; Get the header decoding information loadkeys,'/gryll/data6/buie/reduced/mosaic.key',hdrlist fieldid=strarr(nelements) filter=strarr(nelements) npath=n_elements(rawpath) for i=0,numdir-1 do begin for j=0,npath-1 do begin fn=addslash(rawpath[j])+addslash(alldir[i])+allfiles[i] if exists(fn+'.fits') then ft='.fits' else ft='' if exists(fn+ft) then begin hdr=headfits(fn+ft) parsekey,hdr,hdrlist,info fieldid[i]=info.object filter[i]=info.filter break endif else begin fieldid[i]=['***'] filter[i]=['***'] endelse endfor endfor logerror,logfile,caller='kbosum', $ 'Cleaning up data' ; Get rid of fields that are not discovery fields. fieldid=strupcase(fieldid) z=where(strmid(fieldid,0,1) eq 'F' and strmid(fieldid,1,1) le '9' $ or strmid(fieldid,0,1) eq '*') ndesig=ndesig[z] nlcode=nlcode[z] ntmpcpde=ntmpcode[z] fieldid=fieldid[z] nfirstfield=nfirstfield[z] njd=njd[z] nra=nra[z] ndec=ndec[z] filter=filter[z] nrmag=nrmag[z] nrmagerr=nrmagerr[z] nrate=nrate[z] ndir=ndir[z] ndr=ndr[z] ndt=ndt[z] hmag=hmag[z] a=a[z] e=e[z] inc=inc[z] epherr=epherr[z] nelong=nelong[z] nrawcode=nrawcode[z] nelements=n_elements(ndesig) ; Get quality information, image and orbit iq=fltarr(nelements) oq=fltarr(nelements) ; Put asterisk where there is object name for desig, lcode or tmpcode for i=0,nelements-1 do begin if strlen(ndesig[i]) eq 0 then ndesig[i]='*' if strlen(nlcode[i]) eq 0 then nlcode[i]='*' if strlen(ntmpcode[i]) eq 0 then ntmpcode[i]='*' endfor ; print a file called summary.dat containing the compiled info ; about our kbo survey data outfile = '/gryll/data6/buie/reduced/Summary/discoveries.3.dat.tst' logerror,logfile,caller='kbosum', $ 'Writing output file '+outfile openw, lun, logfile, /get_lun printf, lun, 'desig lcode tmpcode fieldID firstfield date'+$ ' time ra dec filter mag err rate direct'+$ ' dr dt hmag a e i epherr elong iq oq rawcode' blanks=' ' value='(A10,1x,A7,1x,A8,1x,A8,1x,A12,1x,A22,1x,A11,1x,A11,1x,A7,1x,f4.1,1x,'+$ 'f3.1,1x,f4.1,1x,f6.1,1x,f5.1,1x,f4.1,1x,f5.2,1x,f7.2,1x,f5.2,1x,f4.2,'+$ '1x,f7.2,1x,f7.2,1x,f1,1x,f1,1x,A10)' for i=0, nelements-1 do begin printf, lun, format=value,ndesig[i]+blanks,nlcode[i]+blanks,$ ntmpcode[i]+blanks,fieldid[i]+blanks,nfirstfield[i]+blanks,$ njd[i]+blanks,nra[i]+blanks,ndec[i]+blanks,filter[i]+blanks,nrmag[i],$ nrmagerr[i],nrate[i],ndir[i],ndr[i],ndt[i],hmag[i],a[i],e[i],$ inc[i],epherr[i],nelong[i],iq[i],oq[i],nrawcode[i]+blanks endfor ; Close the pipe free_lun, lun logerror,logfile,caller='kbosum', $ 'All done' end