;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; root = '/data/gaia/binary/' dirStat = root+'stats/' file = 'BCarry_0265_4020_Body_Sorted.txt' cutDT = 0.010 ;- second cutMag = 1.0 ;- mag cutRho = 1.0 ;- arcsec cutV = 3.0 ;- mas/s verbose = 0 showPlot= 0 makePlot= 01 showOrbit=0 ; IP IDpred idIDT run planet TransidId fov vpu ac dpx scan helio G V dmag days date OBMT dt RA IDT Dec IDT RA FM Dec FM mura mudec Vidt val vac dra ddec rho Vr_sun Vr_gaia Gal Lat ; px deg deg mag mag mag s s deg deg deg deg mas/s mas/s mas/s mas/s mas/s " " " km/s km/s fovRad = 2.5/3600. ; 2 arcsec sTxt = 0.050/3600 ; 50 mas lScan = 0.200/3600 ; 200 mas sepSig = 0.050/3600. ; 50 mas sm={x:{w:20,s:0.12/3600.}, y:{w:3,s:0.35/3600.}} cut={sep:0.4, mag:1., $ AL:0.200, AC:0.360} iC =0 empExt={num:0L, idPred:0L, idIDT:0L, run:'', name:'', $ transitId:'', fov:0, vpu:0, $ ac:0L,dpx:0L,scan:0.,helio:0.,G:0.,V:0.,dMag:0., $ days:0.d, iso:'', obmt:0.d, dt:0., $ RA_IDT:0.d, Dec_IDT:0.d, RA_FM:0.d, Dec_FM:0.d, $ muRA:0., muDEC:0.,vIDT:0., vAL:0., vAC:0., $ dRA:0., dDEC:0., rho:0.} empDet={num:0L, name:'', $ idIDT:0L, run:'', transitId:'', $ G1:0., G2:0., $ dMag:0., dRA:0., dDec:0., dAL:0., dAC:0.} kDet=0L empVal={num:0L, name:'', N:0, $ muMag:0., sMag:0., $ muAL:0., sAL:0., $ muAC:0., sAC:0. } ;---dummy theta = 2*!PI*findgen(361)/360. ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Read File from Francois Mignard -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Read Binary Ancillary Information -------------------------------------------------- ;--II.1.1-- Blunt Read from File readcol, 'binary.csv', num, id, dM, period, genoid, $ delimiter=',', format='(I,A,F,F,I)', /Silent nbBin=n_elements(num) ;--II.1.2-- Store in Structure bins=replicate({num:0L,id:'',dm:0.,period:0.,genoid:0},nbBin) bins.num=num bins.id=strTrim(id,2) bins.dM=dM bins.period=period bins.genoid=genoid ;--II.2-- Read Mignard Extraction ------------------------------------------------------------ ;--II.2.1-- Blunt Read from File readcol, root+file, $ num, idPred, idIDT, run, name, $ transitId, fov, vpu, $ ac,dpx,scan,helio,G,V,dMag, $ days, iso, obmt, dt, $ RA_IDT, Dec_IDT, RA_FM, Dec_FM, $ muRA, muDEC,vIDT, vAL, vAC, $ dRA, dDEC, rho, $ format='(L,L,L,A,A,'+$ 'A,I,I,'+$ 'L,L,F,F,F,F,F,'+$ 'D,A,D,F,'+$ 'D,D,D,D,'+$ 'F,F,F,F,F,'+$ 'F,F,F)', $ /silent nbExt = n_elements(num) ;--II.2.2-- Store in Structure ext=replicate(empExt,nbExt) ext.num = num ext.idPred = idPred ext.idIDT = idIDT ext.run = strTrim(run,2) ext.name = strTrim(name,2) ext.transitId= strTrim(transitId,2) ext.fov = fov ext.vpu = vpu ext.ac = ac ext.dpx = dpx ext.scan = scan ext.helio = helio ext.G = G ext.V = V ext.dMag = dMag ext.days = days ext.iso = strTrim(iso,2) ext.obmt = obmt ext.dt = dt ext.RA_IDT = RA_IDT ext.Dec_IDT = Dec_IDT ext.RA_FM = RA_FM ext.Dec_FM = Dec_FM ext.muRA = muRA ext.muDEC = muDEC ext.vIDT = vIDT ext.vAL = vAL ext.vAC = vAC ext.dRA = dRA ext.dDEC = dDEC ext.rho = rho ;--II.2.3-- Prepare Output Detection Structure det=replicate(empDet,nbExt) ;--II.2.4-- Prepare Output Extraction Structure omc={o:0., c:0., omc:0.} empty = {num:0L,name:'', days:0D, $ dMag:0., dt:0., rho:0., metric:0., $ vAL:0., vAC:0., vRA:0., vDec:0., $ deltaM:omc, deltaRA:omc, deltaDEC:omc, sep:omc} res = replicate(empty,nbExt) kRes=0L ; name=strTrim(name,2) ; transitid=strTrim(transitid,2) ;--II.3-- Selection Metric ------------------------------------------------------------------- metric = sqrt( (dMag/cutMAG)^2. + $ ;-magnitude difference with pred (dt/cutDT)^2. + $ ;-time difference with pred (rho/cutRHO)^2 ) ;+ $ ;-angular separation ;--II.4-- Identify Unique SSO Identifiers ---------------------------------------------------- uSSO = uniq( num, sort(num) ) nbSSO = n_elements( uSSO ) val=replicate(empVal,nbSSO) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Set Up Loops on SSOs, Satellites, Transit -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; for kSSO=10, nbSSO-1 do begin ;--III.1-- Number of Transit ---------------------------------------------------------------; obs = where( ext.num eq ext[uSSO[kSSO]].num, nbObs ) print, ext[uSSO[kSSO]].num,ext[uSSO[kSSO]].name, nbObs, format='(I6,1x,A-20,1x,I3)' ;--III.2-- Number of Unique Observations ---------------------------------------------------; curDay = ext[obs].days uDay = uniq( curDay, sort(curDay) ) nbUDay = n_elements(uDay) print, ext[uSSO[kSSO]].num,ext[uSSO[kSSO]].name, nbObs, nbUDay, format='(I6,1x,A-20,1x,I3,1x,I3)' ;--III.3-- For Each Satellite --------------------------------------------------------------; pSat=where(ext[uSSO[kSSO]].num eq bins.num, nbSat) print, nbSat for kSat=0, nbSat-1 do begin if bins[pSat[kSat]].genoid eq 1 then begin forprint, bins.num, ' '+bins.id, bins.dM, bins.period, bins.genoid, subset=pSAt[kSat] ;--III.3.1-- Satellite ID and Parameters voID = bins[pSat[kSat]].id deltaM= bins[pSat[kSat]].dM period= bins[pSat[kSat]].period ;--III.3.2-- Satellite Ephemerides forprint, ext.iso, subset=uDay, textout='/tmp/dates.dat', /Silent, /NoComment ephem = voMiriade_callEphemCC( voID, '/tmp/dates.dat', /DateList, /Ast, /Web, obs='@gaia' ) ephem=ephem[0:-2] ;--III.4-- For Each Transit --------------------------------------------------------------; for kDay=0, nbUDay-1 do begin transit = where( ext.num eq ext[uSSO[kSSO]].num and curDay[uDay[kDay]] eq days, nbTransit ) if nbTransit gt 1 then begin ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Satellite Identification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Identify The Primary Transit -----------------------------------------------; ord = sort( metric[transit] ) prim = transit[ord[0]] other= transit[ord[1:-1]] nbO = n_elements(other) ;--IV.2-- Retrieve GDR2 Stars from the Catalog ---------------------------------------; GDR2 = voVizier_query(RA_IDT[prim], DEC_IDT[prim], '5,sec' ) ;--IV.3-- Satellite RA/Dec Separation and delta MAG ----------------------------------; sep = sqrt( (ephem[kDay].dx)^2 + (ephem[kDay].dy)^2 )/3600. satRA = ext[prim].RA_IDT + ephem[kDay].dx/3600. satDec= ext[prim].DEC_IDT + ephem[kDay].dy/3600. angSep = 3600.*sqrt( (ext[other].RA_IDT-satRA)^2. + (ext[other].DEC_IDT-satDEC)^2. ) magSep = ext[other].G-ext[prim].G - deltaM ;--IV.4-- Separation AL/AC Scan ------------------------------------------------------; sepAL = fltarr(nbO) sepAC = fltarr(nbO) sepRA = ext[other].RA_IDT-satRA sepDEC= ext[other].DEC_IDT-satDEC angle = 180-ext[other].scan sepAL = (sepRA*cos(angle*!DTOR) - sepDEC*sin(angle*!DTOR))*3600. sepAC = +(sepRA*sin(angle*!DTOR) + sepDEC*cos(angle*!DTOR))*3600. ;--IV.5-- Satellite Tentative Identification -----------------------------------------; cand = where( magSep le cut.mag and $ abs(sepAL) le cut.AL and $ abs(sepAC) le cut.AC, nbCand ) if nbCand gt 0 then begin iC=[iC,other[cand]] det[kDet].num = ext[prim].num det[kDet].name = ext[prim].name det[kDet].idIDT = ext[prim].idIDT det[kDet].run = ext[prim].run det[kDet].transitID= ext[prim].transitID det[kDet].G1 = ext[prim].G det[kDet].G2 = ext[other[cand]].G det[kDet].dMag = magSep[cand] det[kDet].dRA = sepRA [cand] *3600. det[kDet].dDec = sepDEC[cand] *3600. det[kDet].dAL = sepAL [cand] det[kDet].dAC = sepAC [cand] kDet++ endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Visualization -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if nbSat eq 1 then fileName=string(ext[prim].num,format='(I04)')+'-'+ext[prim].name+'-1' else $ fileName=string(ext[prim].num,format='(I04)')+'-'+ext[prim].name+'-2' if makePlot eq 1 then begin ;--V.1-- Set Up Graphics --------------------------------------------------------------------; ;--V.1.1-- FileNames if nbCand gt 0 then suf='-sat' else suf='' ;--V.1.2-- Graphic Device cgPS_open, Filename=root+fileName+'-'+transitID[prim]+suf+'.eps', $ /Metric, /Decomposed, /Portrait, $ xsize=20, ysize=20, language_level=2, /quiet ;--V.2-- Draw Search FOV --------------------------------------------------------------------; midRA = ext[prim].RA_IDT midDec= ext[prim].Dec_IDT cgPlot, midRA, midDEC, /NoData, /Isotropic, $ xRange = midRA +[-1,1]*fovRad, $ yRange = midDEC+[-1,1]*fovRad, $ title=string(ext[prim].num,format='(I3)')+' '+ext[prim].name+' - '+string(ext[prim].days,format='(F8.3)'), $ xTickInt=1./3600., xTickFormat='(F8.4)', $ yTickInt=1./3600., yTickFormat='(F8.4)', $ /Normal, position=[0.2,0.1,0.95,0.85] ;-FOV ; cgPlot, /OverPlot, midRA+fovRad*cos(theta), midDec+fovRad*sin(theta), color='Gray' ;--V.3-- Direction of Scan & IDT Box --------------------------------------------------------; ;--V.3.1-- Along-Scan Direction for kT=0, nbTransit-1 do $ cgPlot, /OverPlot, ext[transit[kT]].RA_IDT+ lScan*[0,cos(ext[transit[kT]].scan*!DTOR)], $ ext[transit[kT]].Dec_IDT+lScan*[0,sin(ext[transit[kT]].scan*!DTOR)], color='Gray' ;--V.3.2-- IDT Box boxX = sm.x.w*sm.x.s*[-1,1,1,-1,-1] boxY = sm.y.w*sm.y.s*[-1,-1,1,1,-1] cgPlot, /OverPlot, color='dark gray', $ midRA+ boxX*cos(ext[prim].scan*!DTOR) - boxY*sin(ext[prim].scan*!DTOR), $ midDEC+ boxX*sin(ext[prim].scan*!DTOR) + boxY*cos(ext[prim].scan*!DTOR) boxX = sm.x.s*[-1,1,1,-1,-1]/2. boxY = sm.y.s*[-1,-1,1,1,-1]/2. cgPlot, /OverPlot, color='dark gray', $ midRA+ boxX*cos(ext[prim].scan*!DTOR) - boxY*sin(ext[prim].scan*!DTOR), $ midDEC+ boxX*sin(ext[prim].scan*!DTOR) + boxY*cos(ext[prim].scan*!DTOR) ;--V.4-- Prediction Satellite Position ------------------------------------------------------; ;--V.4.1-- Orbit if showOrbit eq 1 then if nbCand gt 0 then if period le 10 then begin step = strTrim(string(period/720,format='(F8.4)'),2)+',d' t0 = epRead(ext[prim].iso+'-'+strTrim(string(period/2.,format='(F8.4)'),2)+',d',/ISO) eOrb = voMiriade_callEphemCC( voID, t0, nbd=720, step=step, /Ast, /Web, obs='@gaia' ) eOrb=eOrb[0:-2] cgPlot, midRA +eOrb.dx/3600., /OverPlot, $ midDec+eOrb.dy/3600., color='Gray', thick=2 endif ;--V.4.2-- Position and Name cgPlot, /OverPlot, psym=1, color='red', satRA, satDec cgText, color='red', satRA+sTxt, satDec+sTxt, strTrim(string(ephem[0].mag+deltaM,format='(F4.1)'),2) ;-predicted position annuli ; cgPlot, /OverPlot, midRA+sep*cos(theta), midDec+sep*sin(theta), color='Gray' ; cgPlot, /OverPlot, midRA+(sep+sepSig)*cos(theta), midDec+(sep+sepSig)*sin(theta), color='Gray' ; cgPlot, /OverPlot, midRA+(sep-sepSig)*cos(theta), midDec+(sep-sepSig)*sin(theta), color='Gray' ; ; ;--V.5-- Satellite Candidate ----------------------------------------------------------------; if nbCand gt 0 then cgPlot, /OverPlot, ext[other[cand]].RA_IDT, ext[other[cand]].Dec_IDT, $ psym='Filled Circle', color='red' ;--V.6-- GDR2 Stars -------------------------------------------------------------------------; if size(GDR2,/type) eq 8 then cgPlot, /OverPlot, psym=3, color='gray', $ GDR2.ra.dec, GDR2.dec.dec ;--V.7-- All IDT Sources --------------------------------------------------------------------; cgPlot, /OverPlot, ext[transit].RA_IDT, ext[transit].Dec_IDT, psym='Open Circle' cgPlot, /OverPlot, ext[prim].RA_IDT, ext[prim].Dec_IDT, psym='Filled Circle' cgText, ext[transit].RA_IDT+sTxt, ext[transit].Dec_IDT+sTxt, string(ext[transit].G,format='(F4.1)') cgPS_close, /png, Delete_PS=0 endif ;stop ; ; ; ; ; ; ; if showPlot eq 1 then begin ; ; ; midRA = RA_IDT[ prim] ; midDec= Dec_IDT[prim] ; cgPlot, RA_IDT[transit], Dec_IDT[transit], /NoData, /Isotropic, $ ; xRange = midRA +[-1,1]*fovRad, $ ; yRange = midDEC+[-1,1]*fovRad, $ ; title=string(num[prim],format='(I3)')+' '+name[prim]+' - '+string(days[prim],format='(F8.3)') ; ; ; ;-FOV ; theta = 2*!PI*findgen(361)/360. ; ; cgPlot, /OverPlot, midRA+fovRad*cos(theta), midDec+fovRad*sin(theta), color='Gray' ; ; ;-predicted position annuli ; cgPlot, /OverPlot, midRA+sep*cos(theta), midDec+sep*sin(theta), color='Gray' ; cgPlot, /OverPlot, midRA+(sep+sepSig)*cos(theta), midDec+(sep+sepSig)*sin(theta), color='Gray' ; cgPlot, /OverPlot, midRA+(sep-sepSig)*cos(theta), midDec+(sep-sepSig)*sin(theta), color='Gray' ; ; cgPlot, /OverPlot, psym=1, color='red', satRA, satDec ; cgText, color='red', satRA+sTxt, satDec+sTxt, strTrim(string(ephem[0].mag,format='(F4.1)'),2) ; ; ;-Direction of scan ; for kT=0, nbTransit-1 do cgPlot, /OverPlot, $ ; RA_IDT[transit[kT]]+ lScan*[0,cos(scan[transit[kT]]*!DTOR)], $ ; Dec_IDT[transit[kT]]+lScan*[0,sin(scan[transit[kT]]*!DTOR)], color='Gray' ; ; ;-IDT box ; boxX = sm.x.w*sm.x.s*[-1,1,1,-1,-1] ; boxY = sm.y.w*sm.y.s*[-1,-1,1,1,-1] ; cgPlot, /OverPlot, color='dark gray', $ ; midRA+ boxX*cos(scan[prim]*!DTOR) - boxY*sin(scan[prim]*!DTOR), $ ; midDEC+ boxX*sin(scan[prim]*!DTOR) + boxY*cos(scan[prim]*!DTOR) ; ; boxX = sm.x.s*[-1,1,1,-1,-1]/2. ; boxY = sm.y.s*[-1,-1,1,1,-1]/2. ; cgPlot, /OverPlot, color='dark gray', $ ; midRA+ boxX*cos(scan[prim]*!DTOR) - boxY*sin(scan[prim]*!DTOR), $ ; midDEC+ boxX*sin(scan[prim]*!DTOR) + boxY*cos(scan[prim]*!DTOR) ; ; ; ;-candidate satellite ; if nbCand gt 0 then begin ; cgPlot, /OverPlot, RA_IDT[other[cand]], Dec_IDT[other[cand]], psym='Filled Circle', color='red' ; ; ; endif ; ; ; ;-GDR2 sources ; if size(gdr2,/type) eq 8 then cgPlot, /OverPlot, psym=3, color='gray', $ ; gdr2.ra.dec, gdr2.dec.dec ; ; ;-IDT sources ; cgPlot, /OverPlot, RA_IDT[transit], Dec_IDT[transit], psym='Open Circle' ; cgPlot, /OverPlot, RA_IDT[prim], Dec_IDT[prim], psym='Filled Circle' ; ; ;-Magnitude in G ; cgText, RA_IDT[transit]+sTxt, Dec_IDT[transit]+sTxt, string(G[transit],format='(F4.1)') ; ; endif ; ; ;-------------- endif ;-III.4: Transits endfor ;-III.4: Transits endif ;-III.3: Satellite endfor ;stop selSSO = where(det.num eq ext[uSSO[kSSO]].num, nbCur ) val[kSSO].num = ext[prim].num val[kSSO].name = ext[prim].name val[kSSO].muMag = mean(det[selSSO].dMag) val[kSSO].sMag = stddev(det[selSSO].dMag) val[kSSO].muAL = mean(det[selSSO].dAL) val[kSSO].sAL = stddev(det[selSSO].dAL) val[kSSO].muAC = mean(det[selSSO].dAC) val[kSSO].sAC = stddev(det[selSSO].dAC) if nbCur gt 2 then begin cgPS_open, FileName=dirStat+fileName+'-dMag.eps', $ /Metric, /Decomposed, /Portrait, $ xsize=20, ysize=15, language_level=2, /quiet cgHistoPlot, det[selSSO].dMag, binSize=0.25, title=string(num[prim],format='(I3)')+' '+ext[prim].name, $ xtitle='Residual in magnitude', $ /normal, position=[0.1,0.15,0.95,0.9] cgPS_close, /png, Delete_PS=0 cgPS_open, FileName=dirStat+fileName+'-AL.eps', $ /Metric, /Decomposed, /Portrait, $ xsize=20, ysize=15, language_level=2, /quiet cgHistoPlot, det[selSSO].dAL*1000, binSize=60, title=string(num[prim],format='(I3)')+' '+ext[prim].name, $ xtitle='Residuals AL (mas)', $ xtickInt=60, xminor=1, $ /normal, position=[0.1,0.15,0.95,0.9] cgPS_close, /png, Delete_PS=0 cgPS_open, FileName=dirStat+fileName+'-AC.eps', $ /Metric, /Decomposed, /Portrait, $ xsize=20, ysize=15, language_level=2, /quiet cgHistoPlot, det[selSSO].dAC*1000, binSize=120, title=string(num[prim],format='(I3)')+' '+ext[prim].name, $ xtitle='Residuals AC (mas)', $ xtickInt=120, xminor=1, $ /normal, position=[0.1,0.15,0.95,0.9] cgPS_close, /png, Delete_PS=0 endif endfor ;-III: SSO iC=iC[1:-1] det=det[0:kDet-1] forprint, val.num, val.name, $ val.muMag, val.sMAg, $ val.muAL, val.sAL, $ val.muAC, val.sAC, $ format='(I4,1x,A-10,3(2x,F5.2,1x,F5.2))', $ comment='Num, Name, muMag, sMag, muAL, sAL, muAC, sAC', $ textout=2 end