; docformat = 'rst' ; ; NAME: ; BinarIM ; PURPOSE: ; Call a sequence of tools for binary asteroid analysis: digital ; coronography, satellite astrometry... ; ;+ ; :Description: ; Call a sequence of tools for binary asteroid analysis: digital ; coronography, satellite astrometry... ; ; :Categories: ; AstroIM, BinarIM ; ; :Params: ; project: in, required, type=string ; Name of the project, used as root directory ; ; :Keywords: ; dates: in, optional, type=string ; An array of dates ('YYYY-MON-DD' to be processed) ; target: in, optional, type=string ; A specific target to process (if not set, all targets are analyzed) ; doDC: in, optional, type=boolean, default=0 ; Run digital coronography on images if set ; doAstrometry: in, optional, type=boolean, default=0 ; Measure satellite positions (GUI or auto) if set ; doEdge: in, optional, type=boolean, default=0 ; A boolean to trigger the edge detection procedure of the deconvolved frames ; ; sciINFO: in, optional, type=structure ; Contains the parameters for the Scientific Frame cleaning:: ; .FILTER : Filter to analyze (ALL | J | K...) ; .FLAT : Prefered flat (SKY | DOME) ; .TARGET : Target for analyze (ALL | SCI | PSF | Name ...) ; .SKY : Sky evaluation (Dark - Median - Fine) ; .SERIE : Number of the First Serie ; ; dcINFO: in, optional, type=structure ; A stucture containing the parameters for Digital Coronography:: ; .MODE : Coronography Mode ; .RADIUS: Maximum search radius ; .BIN : Radial size of the DC annulii ; .SAMPLE: Number of sample to interpolate the intersection ring/pixel ; .STEP : Mini structure with steps for Primary Halo Evaluation:: ; .FLUX : Default step in flux for halo evaluation (% of primary peak) ; .RADIUS: Default step in radius for halo evaluation (pixel) ; ; workINFO: in, optional, type=structure ; A stucture containing miscellaneous information (time, operator):: ; .ROOT : Root Directory of Projects ; .LOGS : Directory of Log Files ; .OBS : Observatory-Specific File ; .EPHEM: Directory for Ephemeris ; .AUTH : User Name ; .INST : User Institute ; .DATE : Current UTC Date ; ; :Uses: ; initIDL(), date_conv() ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change history:: ; Original Version written in February 2009, B. Carry (ESO/LESIA) ; 2013-Sep. : B.Carry (IMCCE) - Soft entirely rewritten from scratch ; 2014-Aug. : B.Carry (IMCCE) - Header in idldoc format ;- pro binarIM, project, dates=dates, target=target, nameBDD=nameBDD, $ doDC=doDC, doAstrometry=doAstrometry, doEdge=doEdge, $ sciINFO=sciINFO, dcINFO=dcINFO, edgeINFO=edgeINFO, workINFO=workINFO compile_opt hidden;, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; nbMaxMoon=3 ;--I.1-- Set IDL to Local Working Environment ----------------------------------------- confAstroIM= initIDL() confRedim = initIDL(confAstroIM.soft.binary, /REDIM ) confBinary = initIDL(confAstroIM.soft.binary, /BINARY) ;--I.2-- Combine Default & Cmd-Line Parameters ---------------------------------------- ;--I.2.1-- Parameters for Scientific Image Selection ---------------------------------- sciDefault={filter: confRedim.sci.filter, $ ;-Filter to analyze (ALL | J | K...) target: confRedim.sci.target, $ ;-Target for analyze (ALL | SCI | PSF | Name ...) serie: confRedim.sci.serie } ;-Starting Serie if not keyword_set(sciINFO) then sciINFO=sciDefault else sciINFO = updateStructure(sciDefault,sciINFO) ;--I.2.2-- Parameters for Digital Coronography ---------------------------------------- dcDefault = {mode: confBinary.mode, $ ;-Coronography Mode radius: confBinary.radius, $ ;-Maximum search radius bin: confBinary.bin, $ ;-Radial size of the DC annulii sample: confBinary.sample, $ ;-Number of sample to interpolate the intersection ring/pixel step: {flux: confBinary.step.flux,$ ;-Default step in flux for halo evaluation (% of primary peak) radius: confBinary.step.radius} } ;-Default step in radius for halo evaluation (pixel) if not keyword_set(dcINFO) then dcINFO=dcDefault else dcINFO = updateStructure(dcDefault,dcINFO) ;--I.2.3-- General Information: Who, When, Where... ----------------------------------- JD=systime(/UTC, /JULIAN) currentDate = strmid(date_conv( JD, 'FITS'),0,16) workDefault={root: confAstroIM.path.root, $ ;-Root Directory of Projects logs: confAstroIM.path.logs, $ ;-Directory of Log Files obs: confAstroIM.path.obs, $ ;-Observatory-Specific File ephem: confAstroIM.path.ephem, $ ;-Directory for Ephemeris auth: confAstroIM.user.name, $ ;-User Name inst: confAstroIM.user.inst, $ ;-User Institute date: currentDate } ;-Current UTC Date if not keyword_set(workINFO) then workINFO=workDefault else workINFO = updateStructure(workDefault,workINFO) ;--I.3-- Disk I/O - File suffixes & Working Directories--------------------------------- ;--I.3.1-- Place Suffixes in a Common Block suFITS = confAstroIM.ext.fits sufJPG = confAstroIM.ext.jpg sufREDUC = confAstroIM.suffix.reduc sufSNA = confAstroIM.suffix.sna sufDC = confAstroIM.suffix.DC ;--I.3.2-- Definition of Logs and Observatories Directories project =strtrim(project,2) ;--I.3.2/A-- Logs Directory split = strSplit( workINFO.logs, '+', count=nbField, /EXTRACT ) if nbField eq 2 then begin case split[0] of '[root]' : workINFO.logs = workINFO.root+split[1] '[project]': workINFO.logs = workINFO.root+project+'/'+split[1] '[date]' : workINFO.logs = workINFO.root+project+'/YYYY-MON-DD/'+split[1] endcase endif ;--I.3.2/B-- Observatory Database split = strSplit( workINFO.obs, '+', /Extract, count=nbField ) if nbField eq 2 then begin case split[0] of '[root]' : workINFO.obs = workINFO.root+split[1] '[project]': workINFO.obs = workINFO.root+project+'/'+split[1] '[date]' : workINFO.obs = workINFO.root+project+'/YYYY-MON-DD/'+split[1] endcase endif ;--I.3.2/C-- Ephemeris Directory split = strSplit( workINFO.ephem, '+', /Extract, count=nbField ) if nbField eq 2 then begin case split[0] of '[root]' : workINFO.ephem = workINFO.root+split[1] '[project]': workINFO.ephem = workINFO.root+project+'/'+split[1] '[date]' : workINFO.ephem = workINFO.root+project+'/YYYY-MON-DD/'+split[1] endcase endif ;--I.4-- Observing Dates -------------------------------------------------------------- ;--I.4.1-- Retrieve Night Information from RedIM Processing nightINFO=astroSummary_read(workINFO.logs+'nights.dat',/NIGHT) nbDate=n_elements(nightINFO.date) ;--I.4.2-- Sub-Select Night from Cmd-line if keyword_set(dates) then begin nbDate=n_elements(dates) sel = -1 for kDate=0, nbDate-1 do begin p = where( nightINFO.date eq dates[kDate] ) if p[0] ne -1 then sel = [sel, p] endfor nightINFO = nightINFO[sel[1:nbDate]] endif ;--I.5-- Mirror, Detector, Filters for Common Observatories -------------------------- obsINFO = astroSummary_read( workINFO.obs, /observatory ) ;--I.6-- General Constants ---------------------------------------------------------- DegToMAS = 3600000.0D MJD_to_JD = 2400000.5D ;--I.7-- Information for Edge Dectection ---------------------------------------------- edgeDefault={LoG: 1., $ ;-Size of the Laplacian of Gaussian bin: 1 ,$ ;-Pixel binning for the analysis sampling: 100 } ;-Number of points sampling the profile if not keyword_set(edgeINFO) then edgeINFO=edgeDefault else edgeINFO=updateStructure(edgeDefault,edgeINFO) ;--I.X-- Loop over Observing Nights ------------------------------------------------- for kDate=0, nbDate-1 do begin ;--I.X.1-- Directory Definition for Current Night ;--I.X.1/A-- Basename, raw data, and calibration Directories dirBASE = workINFO.root+project+'/'+nightINFO[kDate].date+'/' dirRED = dirBASE+'data/reduced/' dirDEC = dirBASE+'data/deconvolved/' dirEDGE = dirBASE+'data/edges/' dirSAT = dirBASE+'satellite/' dirFIG = dirSAT+'figs/' if keyword_set(doDC) or keyword_set(doASTROMETRY) then begin if not file_test( dirSAT, /Dir ) then file_mkdir, dirSAT if not file_test( dirFIG, /Dir ) then file_mkdir, dirFIG endif if keyword_set(doEDGE) then begin if not file_test( dirEDGE, /Dir ) then file_mkdir, dirEDGE endif ;--I.X.1/B-- Log Files Directory (if Incremental from Base Directory) posDate = strpos( workINFO.logs, 'YYYY-MON-DD' ) if posDate ge 0 then begin locPath = workINFO.logs strput, locPath, nightINFO[kDate].date, posDate workINFO.logs = locPath endif ;--I.X.1/C-- Observatories DB (if Incremental from Base Directory) posDate = strpos( workINFO.obs, 'YYYY-MON-DD' ) if posDate ge 0 then begin locPath = workINFO.obs strput, locPath, nightINFO[kDate].date, posDate workINFO.obs = locPath endif ;--I.X.1/D-- Ephemeris Directory (if Incremental from Base Directory) posDate = strpos( workINFO.ephem, 'YYYY-MON-DD' ) if posDate ge 0 then begin locPath = workINFO.ephem strput, locPath, nightINFO[kDate].date, posDate workINFO.ephem = locPath endif ;--I.X.2-- Analyze the Content of the Observations ;--I.X.2/A-- Read and Structure the Log readcol, workINFO.logs+'rawfiles-'+nightINFO[kDate].date+'.list', format='(A,A,A,I,F,I,F,A,I,I,I)', $ fileArr, typeArr, objArr, $ serieArr, expArr, nExpArr, scaleArr, filtArr, $ sizeArrX, sizeArrY, sizeArrZ, /SILENT nbFile = n_elements(fileArr) fileINFO=replicate( { name: '', type: '', target: '', serie: 0, DIT: 0., NDIT: 0, $ PS: 0., filter: '', NX: 0, NY: 0, NZ: 0 }, nbFile ) fileINFO.name = strtrim(fileArr,2) fileINFO.type = strtrim(typeArr,2) fileINFO.target= strtrim(objArr,2) fileINFO.serie = serieArr fileINFO.DIT = expArr fileINFO.NDIT = nExpArr fileINFO.PS = scaleArr fileINFO.filter= strtrim(filtArr,2) fileINFO.NX = sizeArrX fileINFO.NY = sizeArrY fileINFO.NZ = sizeArrZ ;--I.X.2/B-- On Image Cubes, Neglect the Last Slice = Cube average (ESO "cube" mode) cubeInd = where( fileINFO.NZ gt 1, nbCube ) if nbCube gt 0 then fileINFO(cubeInd).NZ-- ;--I.X.2/C-- Sort Files According to Data Types sciList = where( strcmp( typeArr, 'OBJECT' ), nbSci ) fileINFO = fileINFO(sciList) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Selection of Targets and Series -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Process Everything [Default] -------------------------------------------- uniqSCI = fileINFO( uniq(fileINFO.target,sort(fileINFO.target)) ).target ;--II.2-- User Specified some Specific Targets ------------------------------------ if keyword_set(target) then begin uniqSCI = fileINFO( where( strcmp(fileINFO.type,'OBJECT') and $ strcmp(fileINFO.target, target, /FOLD_CASE)) ).target uniqSCI = uniqSCI[0] endif nbUniqSci = n_elements(uniqSCI) ;--II.3-- Selection of the Current Working Files ---------------------------------- for kSci=0, nbUniqSCI-1 do begin ;--II.3.1-- Targets ; kTask=0 ; ssoDefault = {id:0L, type:'', name:'', parentId:0L, number:0L, others: '', $ ; ra:{sex:[0,0,0.], dec:0.}, dec:{sex:[0,0,0.], dec:0.}, $ ; label:'', epoch:'', scale:0., obs:'500'} ; ssoDefault = voReadSSoDNeT(/empty) ssoINFO = voSsoDNet_resolver( uniqSCI[kSci], dump='/tmp/ssod_binarim.dat' ) ssoINFO = ssoINFO[0] ssoINFO = updateStructure(ssoINFO, {scale:0., obs:'500',label:'',epoch:''}) ; if ssoINFO.number ne 0 then ssoINFO.label = '('+string(ssoINFO.number,format='(I0)')+') '+ssoINFO.name $ ; else ssoINFO.label = ssoINFO.name ;--II.3.1/A-- Directories dirOBJ = dirRED+uniqSCI[kSci]+'/' ;--II.3.1/B-- Analyze Series nameSum = setAstroName(TARGET=uniqSCI[kSci], /SUMMARY) serieINFO = astroSummary_read( dirOBJ+nameSum, /SCIENCE ) nbSerie = n_elements(serieINFO) ;--II.3.2-- Series for kSerie=sciINFO.serie, nbSerie-1 do begin posStar = strpos( serieINFO[kSerie].name,'*') num = float(strmid(serieINFO[kSerie].name, posStar-2,2)) ssoINFO.scale = serieINFO[kSerie].PS print, serieINFO[kSerie].name, $ serieINFO[kSerie].N, $ serieINFO[kSerie].filter, $ serieINFO[kSerie].DIT, $ serieINFO[kSerie].PS, $ serieINFO[kSerie].ISO, $ format='(A-15,1x,I3,2x,A-5,2x,F6.2,2x,F7.5,2x,A-22)' ;stop ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Digital Coronography Treatment -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(doDC) then begin ;--III-- GUI -- Console Output --------------------------------------------------------- ; kTask++ ; print, '' ; print, '------------' ; print, '-- Task '+strtrim(string(kTask,format='(I)'),2)+' ----------------------------------------------' ; print, '-- '+nightINFO[kDate].date+' -- Digital Coronography' ; print, '--------------------------------------------------------' ;--III.1-- Name of the Shift'n'Add Frame -------------------------------------------- ; print, ' Running DC' nameSNA = setAstroName( TARGET=uniqSci[kSci], FILTER=serieINFO[kSerie].filter, $ SERIE=num, /SNA ) ;--III.2-- Read Shift'n'Add Frame --------------------------------------------------- image = readfits( dirOBJ+nameSNA, hImage, /SILENT ) dims = size(image) ;--III.3-- Analyze Flux Halo -------------------------------------------------------- haloINFO = findHalo( image, mode=dcINFO.mode, $ targetINFO={X:dims[1]/2., Y:dims(2)/2., head:hImage}, $ dcINFO=dcINFO, quickDC=DC, $ figINFO={dir:dirFIG, name:nameSNA} ) ;--III.4-- Self-Substract Halo ------------------------------------------------------ if ~strCmp(dcINFO.mode,'quick',/fold) then $ DC = subtractHalo( image, mode=dcINFO.mode, quickDC=DC, $ haloINFO=haloINFO, dcINFO=dcINFO ) ;--III.5-- Write Result on Disk ----------------------------------------------------- nameDC = setAstroName( TARGET=uniqSci[kSci], FILTER=serieINFO[kSerie].filter, $ SERIE=num, /DC ) sxaddpar, hImage, 'PRIM-X', haloINFO.pos.X, ' Primary X position (pix)', before='HISTORY' sxaddpar, hImage, 'PRIM-Y', haloINFO.pos.Y, ' Primary Y position (pix)', after='PRIM-X' sxaddpar, hImage, 'PRIM-DX', haloINFO.pos.dX, ' Uncertainty on primary X position (pix)', after='PRIM-Y' sxaddpar, hImage, 'PRIM-DY', haloINFO.pos.dY, ' Uncertainty on primary Y position (pix)', after='PRIM-DX' sxaddpar, hImage, 'PRIM-F', haloINFO.flux.val,' Primary total flux', after='PRIM-DY' sxaddpar, hImage, 'PRIM-dF', haloINFO.flux.unc,' Uncertainty on primary total flux', after='PRIM-F' CHECK_FITS, DC, hImage, /UPDATE, /SILENT writefits, dirSAT+nameDC, DC, hImage ;--III.6-- Special Output, Cleaned from Keywords, For AudeLa ------------------------ nameBDD = setAstroName( TARGET=uniqSci[kSci], FILTER=serieINFO[kSerie].filter, $ SERIE=num, /NoExten ) tsplit=strSplit(project,'/', /Extract, count=nbts) targetBDD=strTrim(tsplit[nbts-1],2) if not keyword_set(nameBDD) then nameBDD=uniqSCI[kSci] hBDD = genoidHeaderBDD(DC[*,*,0], hImage, targetBDD, type='IMG' ) CHECK_FITS, DC[*,*,0], hBDD, /Update, /Silent writefits, dirSAT+nameBDD+'-Image.fits', DC[*,*,0], hBDD hBDD = genoidHeaderBDD(DC[*,*,1], hImage, targetBDD, type='CORONO1' ) CHECK_FITS, DC[*,*,1], hBDD, /Update, /Silent writefits, dirSAT+nameBDD+'-Corono1.fits', DC[*,*,1], hBDD hBDD = genoidHeaderBDD(DC[*,*,2], hImage, targetBDD, type='CORONO2' ) CHECK_FITS, DC[*,*,2], hBDD, /Update, /Silent writefits, dirSAT+nameBDD+'-Corono2.fits', DC[*,*,2], hBDD hBDD = genoidHeaderBDD(DC[*,*,3], hImage, targetBDD, type='CORONO3' ) CHECK_FITS, DC[*,*,3], hBDD, /Update, /Silent writefits, dirSAT+nameBDD+'-Corono3.fits', DC[*,*,3], hBDD hBDD = genoidHeaderBDD(DC[*,*,4], hImage, targetBDD, type='MODEL' ) CHECK_FITS, DC[*,*,4], hBDD, /Update, /Silent writefits, dirSAT+nameBDD+'-Model.fits', DC[*,*,4], hBDD print, 'Fin de traitement de '+nameBDD endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Satellite Differential Astrometry -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(doAstrometry) then begin ;--IV-- GUI -- Console Output --------------------------------------------------------- ; kTask++ ; print, '' ; print, '------------' ; print, '-- Task '+strtrim(string(kTask,format='(I)'),2)+' ----------------------------------------------' ; print, '-- '+nightINFO[kDate].date+' -- Differential Astrometry' ; print, '--------------------------------------------------------' ;--IV.1-- Read the Self-Subtracted Frame ------------------------------------------- ; print, ' Astrometry' nameDC = setAstroName(TARGET=uniqSci[kSci], FILTER=serieINFO[kSerie].filter, SERIE=num, /DC) DC = readfits( dirSAT+nameDC, hDC, /SILENT ) dims = size(DC) ;--IV.2-- Moonlet Structure for Output --------------------------------------------- moonSNA=replicate({ $ ;--IV.2.A-- Moonlet Differential Astrometry astroVal: {X: 0., $ ;-Position in pixel (X) Y: 0., $ ;-Position in pixel (Y) RA: 0., $ ;-Position in angle (Right Ascension -- ") DEC: 0. }, $ ;-Position in angle (Declination -- ") astroUnc: {X: 0., $ ;-Uncertainty in pixel (X) Y: 0., $ ;-Uncertainty in pixel (Y) RA: 0., $ ;-Uncertainty in angle (RA --") DEC: 0.}, $ ;-Uncertainty in angle (DEC --") ;--IV.2.B-- Moonlet Differential Photometry photoVal: {FLUX: 0., $ ;-Flux ratio (counts) DM : 0. }, $ ;-Magnitude difference photoUnc: {FLUX: 0., $ ;-Uncertainty in Counts DM : 0. } }, $ ;-Uncertainty in Mag nbMaxMoon ) moonINFO = moonSNA if not keyword_set(primaryDC) then begin pX = sxpar( hDC, 'PRIM-X' ) pY = sxpar( hDC, 'PRIM-Y' ) dX = sxpar( hDC, 'PRIM-DX' ) dY = sxpar( hDC, 'PRIM-DY' ) f = sxpar( hDC, 'PRIM-F' ) if pX eq 0 then pX = dims[1]/2. if pY eq 0 then pY = dims(2)/2. if dX eq 0 then dX = 1 if dY eq 0 then dY = 1 if f eq 0 then f = total( DC(*,*,0)-median(DC(*,*,0)) ) xyad, hDC, pX, pY, RA, DEC primaryDC = {pos: {X:pX, Y:pY, dX:dX, dY:dY, RA:RA, DEC:DEC}, $ flux: f} endif ;stop satINFO = findSat( DC, targetINFO={name: uniqSCI[kSci], $ pos: primaryDC.pos, $ flux: primaryDC.flux, $ head: hDC}, $ detectINFO={dMax: dcINFO.radius, bin: dcINFO.bin}, $ moonINFO=moonINFO ) moonInd = where( moonINFO.astroVal.X ne 0., nbMoon ) if moonINFO[0].astroVal.X ne 0. and moonINFO[0].astroVal.Y ne 0. then begin telID = getTelescopeID(hDC) camID = getCameraID(hDC) sel = where( strcmp(obsINFO.obsName,telID,/fold) and strcmp(obsINFO.instName,camID,/fold)) IAUcode = obsINFO(sel[0]).obsCode midTimeJD = sxpar( hDC, 'EXPMIDJD' ) midTimeISO = sxpar( hDC, 'EXPMIDIS' ) sysName = strlowcase(uniqSCI[kSci]) strput, sysName, strupcase(strmid(sysName,0,1)), 0 sat = {jd: midTimeJD, iso: date_conv(midTimeJD,'FITS'), $ refName: sysName, $ refsys: sysName, $ xobs: moonINFO[0].astroVal.RA, yobs: moonINFO[0].astroVal.DEC, $ xerr: moonINFO[0].astroUnc.RA, yerr: moonINFO[0].astroUnc.DEC, $ xcal: 0., ycal: 0., $ xomc: 0., yomc: 0., $ timescale: 'UTC', $ centerFrame: 2, typeFrame: 1, $ coordType: 1, refFrame: 1, $ obsIAU: IAUcode, bib: 'B. Carry'} nameXML= setAstroName( TARGET=uniqSci[kSci], FILTER=serieINFO[kSerie].filter, $ SERIE=num, /SAT ) ; print, dirSAT, nameXML genXML = genoide_xmlopen( dirSAT+nameXML, 1 ) genoide_xmladd, genXML, sat genoide_xmlclose, genXML endif endif ; ;------------------------------------------------------------ ; ;--III.2-- Ephemeris Generation - Geometry ; ;--III.2.1-- Generate Ephemeris if not Already Done ; nameEPH = isoObs+'.earth' ; if not FILE_EXIST( dirPOS+nameEPH ) then $ ; call_ephemcc, iauNumber, jdObs, $ ; /GEOCENTRIC, /JULIAN, $ ; /NUMBER, $ ; output=dirPOS+nameEPH ; ;--III.2.2-- Read Ephemeris ; readfmt, dirPOS+nameEPH, $ ; skipline=47, $ ; numline=1, $ ; '(1x,D20.11,6x,I2,1x,I2,1x,F8.5,1x,I3.2,1x,I2,1x,F7.4,2x,F13.9,2x,F6.2,2x,F6.2,3x,F6.2)', $ ; tpsEarth, $ ;-Julian date of observations ; earthRA_h, earthRA_m, earthRA_s, $ ;-RA ; earthDEC_d, earthDEC_m, earthDEC_s, $ ;-DEC ; earthDist, $ ;-Geocentric range ; earthVmag, $ ;-Visual magnitude ; phaseAngle, /SILENT ;-Phase angle ; earthRA = ten( earthRA_h, earthRA_m, earthRA_s ) ; earthDEC= ten( earthDEC_d, earthDEC_m, earthDEC_s ) ; ; ; ; ; ; ; ;--III.3.3-- Moonlet Related Information ; ;--III.3.3-1-- Previous Knowledge (if exist) ; if file_exist(dirSAT+nameSER+'xx-moons.info' ) then begin ; ; ;--III.3.3-1-A-- Read Summary File ; readcol, dirSAT+nameSER+'xx-moons.info', $ ; moonX, moonY, moonDX, moonDY, $ ; moonRA, moonDEC, moonDRA, moonDDEC, $ ; moonFLUX, moonDFLUX, moonDM, moonDDM, /SILENT ; ; ;--III.3.3-1-B-- Fill the Moonlet Structure ; moonSNA={ $ ; ;--III.3.3-1-B.1-- Moonlet Differential Astrometry ; astroVal: {X: moonX, $ ;-Position in pixel (X) ; Y: moonY, $ ;-Position in pixel (Y) ; RA: moonRA, $ ;-Position in angle (Right Ascension -- ") ; DEC: moonDEC }, $ ;-Position in angle (Declination -- ") ; astroUnc: {X: moonDX, $ ;-Uncertainty in pixel (X) ; Y: moonDY, $ ;-Uncertainty in pixel (Y) ; RA: moonDRA, $ ;-Uncertainty in angle (RA --") ; DEC: moonDDEC}, $ ;-Uncertainty in angle (DEC --") ; ;--III.3.3-1-B.2-- Moonlet Differential Photometry ; photoVal: {FLUX: moonFLUX, $ ;-Flux ratio (counts) ; DM : moonDM }, $ ;-Magnitude difference ; photoUnc: {FLUX: moonDFLUX, $ ;-Uncertainty in Counts ; DM : moonDDM } } ;-Uncertainty in Mag ; ; ;--III.3.3-2-- Otherwise: Create an Empty Structure ; endif else begin ; moonSNA={ $ ; ;--III.3.3-2-A-- Moonlet Differential Astrometry ; astroVal: {X: fltarr(nbMaxMoon), $ ;-Position in pixel (X) ; Y: fltarr(nbMaxMoon), $ ;-Position in pixel (Y) ; RA: fltarr(nbMaxMoon), $ ;-Position in angle (Right Ascension -- ") ; DEC: fltarr(nbMaxMoon) }, $ ;-Position in angle (Declination -- ") ; astroUnc: {X: fltarr(nbMaxMoon), $ ;-Uncertainty in pixel (X) ; Y: fltarr(nbMaxMoon), $ ;-Uncertainty in pixel (Y) ; RA: fltarr(nbMaxMoon), $ ;-Uncertainty in angle (RA --") ; DEC: fltarr(nbMaxMoon)}, $ ;-Uncertainty in angle (DEC --") ; ;--III.3.3-2-B-- Moonlet Differential Photometry ; photoVal: {FLUX: fltarr(nbMaxMoon), $ ;-Flux ratio (counts) ; DM : fltarr(nbMaxMoon) }, $ ;-Magnitude difference ; photoUnc: {FLUX: fltarr(nbMaxMoon), $ ;-Uncertainty in Counts ; DM : fltarr(nbMaxMoon) } } ;-Uncertainty in Mag ; endelse ; ; ;--III.3.3-3-- Make a Working Copy of the Moonlet Informations ; moonINFO = moonSNA ; ; ; ; ; ;------------------------------------------------------------ ; ;--III.3.4-- Digital Coronography ; dMax = 85 ; satBin =.55 ; print, currentName ; ;; if not file_exist( dirSAT+nameDC+suFITS ) then begin ; if 1 eq 1 then begin ; ;; image = readfits( '/home/bcarry/data/RAP/simulations/11-primary.fits') ; DC = coronography( image, $ ; targetINFO={posX: dims[1]/2., $ ; posY: dims(2)/2., $ ; dG: earthDist[0]*AUtoKM, $ ; scale: scaleArr(kSerie), $ ; head: hImage}, $ ;; /FASTDC , $ ; /FULLDC , $ ; detectINFO={dMax: dMax, bin: satBin }, $ ; figINFO={dir:dirFIG, name: nameSNA}, $ ; contrastINFO=contrastDC, $ ; qc=qcDC ) ; CHECK_FITS, DC, hImage, /UPDATE, /SILENT ; writefits, dirSAT+nameDC+suFITS, DC, hImage ;stop ; ; help, contrastDC, qcDC, /struc ; ; forprint, contrastDC.pos[0], contrastDC.pos[1], $ ; textout=dirSAT+nameDC+sufDAT, $ ; /NOCOMMENT, /SILENT ; ; ; endif else begin ; ; DC = readfits( dirSAT+nameDC+suFITS, hImage, /SILENT ) ; ; readcol, dirSAT+nameDC+sufDAT, $ ; px, py, /silent ; contrastDC={pos: [px,py]} ; endelse ; ; ; ;------------------------------------------------------------ ; ;--III.3.5-- Satellite Detection ; detectINFO = findsat( DC, $ ; targetINFO={posX: contrastDC.pos[0], $ ; posY: contrastDC.pos[1], $ ; dG: earthDist[0]*AUtoKM, $ ; scale: scaleArr(kSerie), $ ; head: hImage}, $ ; detectINFO={dMax: dMax, bin: satBin}, $ ; figINFO={dir:dirFIG, name: nameSNA}, $ ; moonINFO=moonINFO ) ; ; moonInd = where( moonINFO.astroVal.X ne 0., nbMoon ) ; ; midTimeJD = sxpar( hImage, 'EXPMIDJD' ) ; midTimeISO = sxpar( hImage, 'EXPMIDIS' ) ; ; sysName = strlowcase(uniqSCI(kObj)) ; ; strput, sysName, strupcase(strmid(sysName,0,1)), 0 ; sat = {jd: midTimeJD, iso: date_conv(midTimeJD,'FITS'), $ ; refName: sysName, $ ; refsys: sysName, $ ; xobs: moonINFO.astroVal.RA[0], yobs: moonINFO.astroVal.DEC[0], $ ; xerr: moonINFO.astroUnc.RA[0], yerr: moonINFO.astroUnc.DEC[0], $ ; xcal: 0., ycal: 0., $ ; xomc: 0., yomc: 0., $ ; timescale: 'UTC', $ ; centerFrame: 2, typeFrame: 1, $ ; coordType: 1, refFrame: 1, $ ; obsIAU: nightINFO.obsIAU[kDate], bib: 'B. Carry'} ; ; ; genXML = genoide_xmlopen( dirSAT+nameSAT+sufXML, 1 ) ; genoide_xmladd, genXML, sat ; genoide_xmlclose, genXML ; ; forprint, midTimeISO, $ ; moonINFO.astroVal.X, $ ; moonINFO.astroVal.Y, $ ; moonINFO.astroUnc.X, $ ; moonINFO.astroUnc.Y, $ ; moonINFO.astroVal.RA, $ ; moonINFO.astroVal.DEC, $ ; moonINFO.astroUnc.RA, $ ; moonINFO.astroUnc.DEC, $ ; subset=moonInd, $ ; format='(A23,2x,F8.3,2x,F8.3,2x,F8.3,2x,F8.3,2x,'+$ ; 'F8.4,2x,F8.4,2x,F8.4,2x,F8.4)', $ ; textout=dirSAT+nameSAT+sufDAT, /SILENT, /NOCOMMENT ; ;;stop ;; ; ;; ;; ;; java -jar ~/Documents/softs/VO/stilts.jar tcat in=alauda.H-01-sat.xml alauda.H-03-sat.xml alauda.H-05-sat.xml alauda.H-08-sat.xml ifmt=votable ofmt=votable out=test.xml ;; ;; ; ;;print, '----------BINARY----------' ;;help, detectINFO.X, detectINFO.Y, /struc ;;stop ; ;; detect = findsat( DC, $ ;; targetINFO={posX: contrastDC.pos[0], $ ;; posY: contrastDC.pos[1], $ ;; dG: earthDist[0]*AUtoKM, $ ;; scale: scaleArr(kSerie), $ ;; head: hImage}, $ ;; detectINFO={dMax: dMax, bin: satBin }, $ ;; moonINFO=moonINFO, $ ;; nameEPS = dirSAT+nameSAT+sufEPS) ;; moonInd = where( moonINFO.X ne 0., nbMoon ) ;; ;; ;; CHECK_FITS, detect, hImage, /UPDATE, /SILENT ;; writefits, dirSAT+nameSAT+suFITS, detect, hImage ;; ;; ;; ;------- OUT ARRAY ;; if moonInd[0] ne -1 then begin ;; outMoon = fltarr( nbFrameArr(kSerie)+1, 8, nbMoon) ;; moonTag = strarr( nbFrameArr(kSerie)+1 ) ;; jdMoon = dblarr( nbFrameArr(kSerie)+1 ) ;; ;; moonSNA=moonINFO ;; ;; moonTag[0] = nameSER+'xx' ;; for kMoon=0, nbMoon-1 do begin ;; outMoon(0,0,kMoon) = moonINFO.X(kMoon) ;; outMoon(0,1,kMoon) = moonINFO.Y(kMoon) ;; outMoon(0,4,kMoon) = moonINFO.DX(kMoon) ;; outMoon(0,5,kMoon) = moonINFO.DY(kMoon) ;; outMoon(0,2,kMoon) = moonINFO.RA(kMoon) ;; outMoon(0,3,kMoon) = moonINFO.DEC(kMoon) ;; outMoon(0,6,kMoon) = moonINFO.DRA(kMoon) ;; outMoon(0,7,kMoon) = moonINFO.DDEC(kMoon) ;; endfor ;; ;; forprint, moonINFO.X, $ ;; moonINFO.Y, $ ;; moonINFO.DX, $ ;; moonINFO.DY, $ ;; moonINFO.RA, $ ;; moonINFO.DEC, $ ;; moonINFO.DRA, $ ;; moonINFO.DDEC, $ ;; subset=moonInd, $ ;; format='(F8.3,2x,F8.3,2x,F8.3,2x,F8.3,2x,'+$ ;; 'F8.4,2x,F8.4,2x,F8.4,2x,F8.4)', $ ;; textout=dirSAT+nameSER+'xx-moons.info', /SILENT, /NOCOMMENT ;; ;; ;; endif ;; ;; ;; endif ;--End of III.3 ;; ;; ;; ;; ;; ;------------------------------------------------------------ ;; ;--III.4-- Analysis of all the Images of the Serie ;; if moonInd[0] ne -1 then begin ;; ;; ;; ;; ;-III.3.1 - Reading the current image ;; nameIN = currentName+sufCLEAN+suFITS ;; nameSAT = currentName+sufSAT+suFITS ;; nameHEAD= currentName+sufRED+suFITS ;; ;; image = readfits( dirRED+nameIN, /SILENT ) ;; hImage = headfits( dirRED+nameHEAD, /SILENT ) ;; ;; jdMoon(kImage+1) = max( [sxpar( hImage, 'MJD-OBS'), $ ;-Keck, ESO HImages ;; sxpar( hImage, 'MJD_OBS')]) $ ;-Gemini HImages ;; + MJD_to_JD ;; print, ' BUG!!!!' ;; print, 'AJOUTER LE TEMPS D EXPOSITION' ;; ;; ;; print, ' image '+string(kImage+1,format='(I02)')+'/'+string(nbFrameArr(kSerie),format='(I02)')+$ ;; ' '+isoObs ;; ;; moonINFO=moonSNA ;; detect = findsat( image, $ ;; targetINFO={posX:sxpar(hImage,'CENT_XC'), $ ;; posY:sxpar(hImage,'CENT_YC'), $ ;; dG: earthDist[0]*AUtoKM, $ ;; scale: scaleArr(kSerie), $ ;; head:hImage}, $ ;; moonINFO=moonINFO ) ;; ;; moonTag(kImage+1) = currentName ;; for kMoon=0, nbMoon-1 do begin ;; outMoon(kImage+1,0,kMoon) = moonINFO.X(kMoon) ;; outMoon(kImage+1,1,kMoon) = moonINFO.Y(kMoon) ;; outMoon(kImage+1,2,kMoon) = moonINFO.RA(kMoon) ;; outMoon(kImage+1,3,kMoon) = moonINFO.DEC(kMoon) ;; outMoon(kImage+1,4,kMoon) = moonINFO.DX(kMoon) ;; outMoon(kImage+1,5,kMoon) = moonINFO.DY(kMoon) ;; outMoon(kImage+1,6,kMoon) = moonINFO.DRA(kMoon) ;; outMoon(kImage+1,7,kMoon) = moonINFO.DDEC(kMoon) ;; endfor ;; ;; ;; CHECK_FITS, detect, hImage, /UPDATE, /SILENT ;; writefits, dirSAT+nameSAT, detect, hImage ;; ;; endif ;; ;; ;; ;; ;; ;-LAST FRAME = WRITE the SERIE ;; if kImage eq nbFrameArr(kSerie)-1 AND $ ;; moonInd[0] ne -1 then begin ;; ;; moonDesign =['A', 'B', 'C'] ;; ;; jdMoon[0] = mean( jdMoon(1:nbFrameArr(kSerie)) ) ;; ;; ;; for kMoon=0, nbMoon-1 do begin ;; ;; ;; if nbFrameArr(kSerie) ge 2 then begin ;; meanArr = fltarr(8) ;; medianArr= fltarr(8) ;; stdArr = fltarr(8) ;; resistArr= fltarr(8) ;; std2Arr = fltarr(8) ;; ;; for kVal=0, 7, 2 do begin ;; meanArr(kVal) = mean( outMoon(*, kVal, kMoon)) ;; medianArr(kVal) = median(outMoon(*, kVal, kMoon)) ;; stdArr(kVal) = stddev(outMoon(*, kVal, kMoon)) ;; newVec = outMoon(where( abs(outMoon(*, 0, kMoon)-medianArr[0]) le 1.5*stdArr[0] ), kVal, kMoon) ;; resistArr(kVal) = mean(newVec) ;; std2Arr(kVal) = stddev(newVec) ;; endfor ;; for kVal=1, 7, 2 do begin ;; meanArr(kVal) = mean( outMoon(*, kVal, kMoon)) ;; medianArr(kVal) = median(outMoon(*, kVal, kMoon)) ;; stdArr(kVal) = stddev(outMoon(*, kVal, kMoon)) ;; newVec = outMoon(where( abs(outMoon(*, 1, kMoon)-medianArr[1]) le 1.5*stdArr[1] ), kVal, kMoon) ;; resistArr(kVal) = mean(newVec) ;; std2Arr(kVal) = stddev(newVec) ;; endfor ;; ;; ;; forprint, [moonTag, 'mean','median','std', 'resist', 'std2'], $ ;; [jdMoon, jdMoon[0], jdMoon[0], stddev(jdMoon(1:nbFrameArr(kSerie))), jdMoon[0], jdMoon[0]], $ ;; [outMoon(*,0,kMoon), meanArr[0], medianArr[0], stdArr[0], resistArr[0], std2Arr[0]] ,$ ;; [outMoon(*,1,kMoon), meanArr[1], medianArr[1], stdArr[1], resistArr[1], std2Arr[1]] ,$ ;; [outMoon(*,2,kMoon), meanArr(2), medianArr(2), stdArr(2), resistArr(2), std2Arr(2)] ,$ ;; [outMoon(*,3,kMoon), meanArr(3), medianArr(3), stdArr(3), resistArr(3), std2Arr(3)] ,$ ;; [outMoon(*,4,kMoon), meanArr(4), medianArr(4), stdArr(4), resistArr(4), std2Arr(4)] ,$ ;; [outMoon(*,5,kMoon), meanArr(5), medianArr(5), stdArr(5), resistArr(5), std2Arr(5)] ,$ ;; [outMoon(*,6,kMoon), meanArr(6), medianArr(6), stdArr(6), resistArr(6), std2Arr(6)] ,$ ;; [outMoon(*,7,kMoon), meanArr(7), medianArr(7), stdArr(7), resistArr(7), std2Arr(7)] ,$ ;; textout=dirSAT+nameSER+'xx'+sufSAT+moonDesign(kMoon)+sufDAT, $ ;; format='(-A'+string(strlen(nameSER)+3,format='(I)')+$ ;; ',2x,D16.8,2x,F8.3,2x,F8.3,2x,F8.3,2x,F8.3,2x,'+$ ;; 'F8.4,2x,F8.4,2x,F8.4,2x,F8.4)', $ ;; /NOCOMMENT, $ ;; /SILENT ;; endif else begin ;; forprint, moonTag, $ ;; jdMoon, $ ;; outMoon(*,0,kMoon),$ ;; outMoon(*,1,kMoon),$ ;; outMoon(*,2,kMoon),$ ;; outMoon(*,3,kMoon),$ ;; outMoon(*,4,kMoon),$ ;; outMoon(*,5,kMoon),$ ;; outMoon(*,6,kMoon),$ ;; outMoon(*,7,kMoon),$ ;; textout=dirSAT+nameSER+'xx'+sufSAT+moonDesign(kMoon)+sufDAT, $ ;; format='(-A'+string(strlen(nameSER)+3,format='(I)')+$ ;; ',2x,D16.8,2x,F8.3,2x,F8.3,2x,F8.3,2x,F8.3,2x,'+$ ;; 'F8.4,2x,F8.4,2x,F8.4,2x,F8.4)', $ ;; /NOCOMMENT, $ ;; /SILENT ;; ;; ;; ;; endelse ;; ;; endfor ;; ; endif ; ; ; ;j2nextSat: ; endfor ; ; ; endfor ; ; ;j2NextObject: ; endfor ;-lopp obbjet ; ; ; ; ; endfor ; ; ; ;end ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Edge Detection of Scientific Frames -----------------------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(doEDGE) then begin ;--V-- GUI -- Console Output --------------------------------------------------------- ; kTask++ ; print, '' ; print, '------------' ; print, '-- Task '+strtrim(string(kTask,format='(I)'),2)+' ----------------------------------------------' ; print, '-- '+nightINFO[kDate].date+' -- Edge Detection' ; print, '--------------------------------------------------------' ;--V.1-- Bookkeeping Variable ------------------------------------------------------------ profNameArr = strarr(serieINFO[kSerie].N) edgeList = intarr(serieINFO[kSerie].N) ;--V.2-- Loop over Individual Frames of the Serie ---------------------------------- for kFrame=0, serieINFO[kSerie].N-1 do begin ;--V.2.2-- Name of the Deconvolved Frame -------------------------------------------- nameDEC = setAstroName( TARGET=uniqSci[kSci], FILTER=serieINFO[kSerie].filter, $ SERIE=num, frame=kFrame+1, /DEC ) if file_test( dirDEC+nameDEC, /read) then begin ;--V.2.2-- Read Deconvolved Frame --------------------------------------------------- image = readfits( dirDEC+nameDEC, hImage, /SILENT ) dims = size(image) ;--V.2.3-- Extract Edges ------------------------------------------------------------ edges = edgesDetection( image, sigmaLoG=edgeINFO.LoG, binning=edgeINFO.bin, header=hImage ) dimE=size(edges) ;--V.2.4-- Write Edges on Disk ------------------------------------------------------ if dimE[0] eq 2 then begin edgeList[kFrame] = 1 nameEdge = setAstroName( TARGET=uniqSci[kSci], FILTER=serieINFO[kSerie].filter, $ SERIE=num, frame=kFrame+1, /EDGE ) sxaddpar, hImage, 'EDGE_SIG', edgeINFO.log , 'LoG Wavelet Standard Deviation (pix)' sxaddpar, hImage, 'EDGE_BIN', edgeINFO.bin , 'Edge Determination Over-sampling', after='EDGE_SIG' sxaddpar, hImage, 'HISTORY' , workINFO.date+' : Edge Determination by '+workINFO.auth+', '+workINFO.inst CHECK_FITS, edges, hImage, /UPDATE, /SILENT writefits, dirEDGE+nameEdge, edges, hImage profNameArr[kFrame] = nameEdge endif endif endfor ;--V.2--End of Loop on Frames ;--V.3-- Combine Frames into a Single Profile -------------------------------------- ;--V.3.1-- Only Use Valid Frames sel=where( edgeList eq 1, nbEdges ) ;--V.3.2-- Define the Figure Parameters nameEdge = setAstroName( TARGET=uniqSci[kSci], FILTER=serieINFO[kSerie].filter, SERIE=num, /EDGE, /EPS ) ssoINFO.epoch=serieINFO[kSerie].iso ;--V.3.3-- Combine Edges & Generate Figure if nbEdges gt 0 then begin edge = edgesMerging( dirEDGE+profNameArr(sel), sampling=edgeINFO.sampling );, figure=figINFO ) edgesDisplay, edge, dump=dirEDGE+nameEdge, ssoINFO=ssoINFO ;--V.3.4-- Write Edges on Disk nameEdge = setAstroName( TARGET=uniqSci[kSci], FILTER=serieINFO[kSerie].filter, SERIE=num, /EDGE, /DAT ) koalaWriteAO, dirEDGE+nameEdge, edge.prof, ssoINFO endif endif ;--V--End of Edge Detection ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- XXX -- Close Loops and Say ByeBye ---------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; primaryDC=0 ;stop endfor ;--II.3.2-- End of Loop over Series endfor ;--II.3.1-- End of Loop over Targets endfor ;--I.X-- End of Loop over Nights end