; docformat = 'rst' ; ; NAME: ; genoidPlotMag ; PURPOSE: ; Plot the distribution of magnitude difference between an asteroid and its satellite ; ;+ ; :Description: ; Plot the distribution of magnitude difference between an asteroid and its satellite ; ; :Categories: ; Disk I/O, Genoid ; ; :Params: ; file: in, required, type=string ; The path to the VOTable to be read. Fields are:: ; .JD: Date in Julian Days ; .ISO: Date in ISO format ; .refName: Name of the ; .refSys: ; .X: .Obs: X Position - Observed ; .Err: X Position - Uncertainty ; .Cal: X Position - Computed by Genoide ; .OMC: X Position - Observed minus Computed ; .Y: similar to X structure ; .timeScale: Time scale (TT, UTC...) ; .centerFrame: Reference frame center ; .typeFrame: Type of reference frame ; .coordType: Type of coordinates ; .refFrame: Reference plane (1:equator, 2:ecliptic) ; .obsIAU: IAU code of the observatory ; .method: Photometric method used to obtain x,y measurements ; .Mag: Difference in magnitude between primary and satellite ; .dMag: Standard deviation on mag ; .QC: Quality note of the measure (A|A+|B|C|D) ; .telescope: Telescope used to acquire the measure ; .camera: Camera used to acquire the measure ; .filter: Filter used to acquire the measure ; .author: Name of the person who measured the image ; .when: Epoch at which the image was measured the last time (ISO format) ; .bib: Bibliographic reference ; ; :Keywords: ; dump: in, optional, type=string, default='./genoid-obs-omc.eps' ; The path to the figure file ; binSize: in, optional, type=float, default=0.2 ; Size of the magnitude bins ; other: in, optional, type=fltarr ; An array of magnitudes to be plotted too ; unCut: in, optional, type=float ; If set, only select magnitude which uncertainty is smaller than unCut ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Written in April 2017, B. Carry (OCA) ; 2017 Dec. - B. Carry (OCA) - Adjust Y-axis wrt number of points ;- pro genoidPlotMag, file, dump=dump, binSize=binSize, other=other, $ unCut=unCut, diameter=diameter, verbose=verbose COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Syntax Validation -------------------------------------------------------------------; if not keyword_set(file) then begin message, /ioError, 'Syntax: genoidPlotObs, file [,dump=, /omc ]' return endif ;--I.2-- File Accessibility ------------------------------------------------------------------; if not file_test(file,/read) then begin message, /Info, 'File cannot be read: '+strTrim(file,2) return endif ;--I.3-- Set Defaults ------------------------------------------------------------------------; if not keyword_set(dump) then dump = './genoid-obs-mag.eps' if not keyword_set(binSize) then binSize = 0.2 ;--I.4-- Read the Observations ---------------------------------------------------------------; inObs = genoidReadObs(file) nbObs = n_elements(inObs) ;--I.5-- Trim Observations to Best Magnitudes if Requested -----------------------------------; if not keyword_set(unCut) then obs = inObs else begin valid = where( inObs.dMag le unCut, nbValid ) ; ;-daphne ONLY ; valid = where( inObs.dMag le unCut and $ ; inObs.mag le -8, nbValid ) ; print, 'Only plotting '+strTrim(string(nbValid),2)+' data points.' if valid[0] eq -1 then begin message, 'Magnitude precision threshold too low: no data to plot' return endif obs = inObs[valid] ; flag=intarr(nbObs) ; flag[valid]=1 ; forprint, inObs.iso, inObs.camera, $ ; inObs.mag, inObs.dmag, $ ; flag, $ ; format='(A-16,2x,A-10,2x,F6.2,2x,F5.2,2x,I1)' endelse ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Define Graphical Parameters -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Size of the Figure -----------------------------------------------------------------; xSize=30. ySize=20. ;--II.2-- Graphics Corners in Normalized Coordinates -----------------------------------------; xBot = 0.08 yBot = 0.10 xEnd = 0.98 yEnd = xEnd ;--II.3-- Axes Ranges and Ticks --------------------------------------------------------------; xRange = mean(obs.mag) + [-1,1] xRange = [floor(xRange[0]), ceil(xRange[1])] xTickInt = 1 xMinor = 4 yTickInt = 5 yMinor = 5 void = replicate(' ', 40) ;print, mean(obs.mag) ;print, xRange ;--II.4-- Symbol Size and Colors -------------------------------------------------------------; barFill = 0.80 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Analyze Observations and Prepare Graphics -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Compute Magnitude Histograms ------------------------------------------------------; ;--III.1.1-- Fix Number of Bins nbBin = round( (xRange[1]-xRange[0])/binSize + 1 ) ;--III.3.2-- Compute Histograms aHist = histogram( inObs.mag, min=xRange[0], max=xRange[1], binSize=binSize ) mHist = histogram( obs.mag, loc=mArr, min=xRange[0], max=xRange[1], binSize=binSize ) yRange = [0,max(mHist)] * 1.1 ;--III.3.3-- Adjust Normal Distributions mFit = mpFitPeak( mArr, mHist, mParam, nTerms=3 ) fine = xRange[0] + (xRange[1]-xRange[0])*findgen(1000)/999. mGauss=evalGauss( fine, mParam) print, 'Average delta-mag: '+string(mParam[1], format='(F6.2)') print, 'Uncertainty : '+string(mParam[2], format='(F6.2)') print, 'Weighted average: ', meanWithUnc(obs.mag, obs.dmag) if keyword_set(diameter) then begin fluxRatio = 10.^(-0.4*mParam[1]) fluxMin = 10.^(-0.4*(mParam[1]-mParam[2])) fluxMax = 10.^(-0.4*(mParam[1]+mParam[2])) fluxUnc = fluxMax-fluxMin satDiam = diameter.val * sqrt( fluxRatio ) satUnc = satDiam*sqrt( (fluxUnc/fluxRatio)^2. + (diameter.unc/diameter.val)^2. ) print, '' print, 'Diameter (km) : '+string(satDiam,format='(F7.2)') print, 'Sigma (km) : '+string(satUnc ,format='(F7.2)') ; stop endif if max(mHist) lt 5 then begin yTickInt=1 yMinor=1 endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Create Graphics -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Open EPS ---------------------------------------------------------------------------; cgPS_open, Filename=dump, /Metric, /Decomposed, /Portrait, /Encapsulated, $ xSize=xSize, ysize=ySize, language_level=2, /Quiet ;--II.2-- Plot Distribution of Magnitudes ----------------------------------------------------; ;--II.2.1-- Plot Area cgPlot, mArr, mHist, /NoData, /Normal, $ xRange=xRange, xTickInt=xTickInt, xMinor=xMinor, xStyle=1, $ yRange=yRange, yTickInt=yTickInt, yMinor=yMinor, yStyle=1, $ position=[xBot, yBot, xEnd, yEnd], $ xTitle = 'Magnitude difference', $ yTitle = 'Distribution' ;--II.2.2-- Draw Bars All Values for kB=0, nbBin-1 do begin xx = mArr[kB] + binSize/2. + [-1,1,1,-1,-1]*BarFill*binSize/2. yy = aHist[kB] * [0,0,1,1,0] cgPlot, /OverPlot, xx, yy endfor ;--II.2.4-- Draw Bars - Selected for kB=0, nbBin-1 do begin xx = mArr[kB] + binSize/2. + [-1,1,1,-1,-1]*BarFill*binSize/2. yy = mHist[kB] * [0,0,1,1,0] cgColorFill, xx, yy, color='Cornflower Blue' cgPlot, /OverPlot, xx, yy endfor ;--II.2.4-- Overplot Gaussian Approxmation cgPlot, /OverPlot, fine+binSize/2, mGauss, lineStyle=2, color='Black', lineThick=2 ;--II.4-- Optional Other Distribution --------------------------------------------------------; if keyword_set(other) then begin ;--II.4.1-- Compute Histogram oHist = histogram( other, loc=oArr, min=xRange[0], max=xRange[1], binSize=binSize ) oSize = 0.7 ;--II.4.2-- Display Bars for kB=0, nbBin-1 do begin xx = oArr[kB] + binSize/2. + [-1,1,1,-1,-1]*BarFill*oSize*binSize/2. yy = oHist[kB] * [0,0,1,1,0] cgColorFill, xx, yy, color='Red' cgPlot, /OverPlot, xx, yy endfor ;--II.4.3-- Add Captions xKey = mean(obs.mag) + 0.7 xLen = 2.6 yKey = max(mHist)-1 yShi = 1.5 xx = xKey + binSize/2. + [-1,1,1,-1,-1]*BarFill*binSize/2. yy = yKey + [0,0,1,1,0] cgColorFill, xx, yy, color='Cornflower Blue' cgPlot, /OverPlot, xx, yy cgText, xKey + xLen*BarFill*binSize/2., yKey+0.25, 'Current' xx = xKey + binSize/2. + [-1,1,1,-1,-1]*BarFill*oSize*binSize/2. yy = yKey + [0,0,1,1,0] - yShi cgColorFill, xx, yy, color='Red' cgPlot, /OverPlot, xx, yy cgText, xKey + xLen*BarFill*binSize/2., yKey+0.25 - yShi, 'Marchis2008+' endif cgPlot, mArr, mHist, /NoData, /Normal, /NoErase, $ xRange=xRange, xTickInt=xTickInt, xMinor=xMinor, xStyle=1, $ yRange=yRange, yTickInt=yTickInt, yMinor=yMinor, yStyle=1, $ position=[xBot, yBot, xEnd, yEnd], $ xTitle = 'Magnitude difference', $ yTitle = 'Distribution' ;--II.5-- Close Graphic File -----------------------------------------------------------------; cgPS_close, /png, Delete_PS=0 end