; docformat = 'rst' ; ; NAME: ; ssoDiscovery_plot ; PURPOSE: ; Plot utility for the IMCCE Database of Asteroids Discovery ; ;+ ; :Description: ; Plot utility for the IMCCE Database of Asteroids Discovery ; ; :Categories: ; Database, Asteroid ; ; :Returns: ; ; :Keywords: ; ep: in, optional, type=string/double ; Epoch for which the number of discovered object is desired ; population: in, optional, type=string ; SSO population to be plotted ('*', or any type/subtype from ; http://vo.imcce.fr/webservices/skybot/?documentation) ; absoluteMagnitude: in, optional, type=boolean, default=1 ; Set to plot SSO distribution as function of absolute magnitude. ; year: in, optional, type=boolean, default=0 ; Set to plot SSO distribution as function of time. ; cumulative: in, optional, type=boolean, default=0 ; Set to plot histograms as cumulative histograms. ; SFD: in, optional, type=boolean, default=0 ; Set to overplot the SFD power law ; histoAE: in, optional, type=boolean, default=0 ; Set to plot the (a,e) distribution of SSOs ; histoAI: in, optional, type=boolean, default=0 ; Set to plot the (a,i) distribution of SSOs ; histoBoth: in, optional, type=boolean, default=0 ; Set to plot both the (a,i) and (a,e) distribution of SSOs on the same graphic ; ; language: in, optional, type=string, default='en' ; Set the language of the graphics (en|fr) ; ; hLimit: in, optional, type=float ; A limiting absolute magnitude (H) for the number of objects to be discovered. ; vLimit: in, optional, type=float ; A limiting apparent magnitude (V) for the number of objects to be discovered. ; elongation: in, optional, type=float, default=180. ; The elongation of observations for the apparent magnitude (V) limit. ; ; dump: in, optional, type=string, default='./ssoDiscovery.eps' ; Path to the exported figure. ; config: in, optional, type=string ; Path to the configuration file for catalogs, at the ; minimum this file should contain the 2 following lines to be used ; by current routine:: ; [Tables from Literature] ; disco.imcce = PATH_TO_YOUR_IMCCE_SSO_DISCOVERY_FILE ; ; :Examples: ; IDL> ssos=ssoDiscovery_read(/dump,disco='./sso.dat') ; IDL> temps = [epRead('now'), epread('2010'), epread('2006'),epread('2000')] ; IDL> ssoDiscovery_plot, disco=ssos, ep=temps,/cum,vLimit=19.5,pop='MB', elong=90. ; IDL> ssoDiscovery_plot, disco=ssos, ep='now', /Year ; ; :Uses: ; initIDL, ssoDiscovery_read, ssoDiscovery_params, ssoHtoV, epRead, exponent, cgLibrary, lineStyleID ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Written in September 2014, B. Carry (IMCCE) ; 2017 June - B. Carry (OCA) - Coded /Year keyword ; 2017 Sep. - B. Carry (OCA) - idl2 added, Keyword -> absoluteMagnitude ; 2018 Mar. - B. Carry (OCA) - Keyword -> (a,e,i) histograms ; 2020 Mar. - J. Berthier (IMCCE) - Lettre accentuee: string() ; e.g.: é = string(233B) cf https://upload.wikimedia.org/wikipedia/commons/8/81/Table_ascii_extended.png ;- pro ssoDiscovery_plot, ep=ep, population=population, dump=dump, ancient=ancient, $ absoluteMagnitude=absoluteMagnitude, year=year, $ histoAE=histoAE, histoAI=histoAI, histoBoth=histoBoth, $ cumulative=cumulative, SFD=SFD, $ hLimit=hLimit, vLimit=vLimit, $ elongation=elongation, $ aInfo=aInfo, eInfo=eInfo, iInfo=iInfo, hInfo=hInfo, yInfo=yInfo, $ discovery=discovery, config=config, language=language COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Set IDL to Local Working Environment ------------------------------------------------ if not keyword_set(discovery) then begin ;--I.1.1-- Retrieve Default Configuration if not keyword_set(CONFIG) then config = initIDL(/Catalog) $ else config = initIDL(config, /Catalog) ;--I.1.2-- Set Path to Discovery Database discovery = config.disco.imcce endif ;--I.2-- Set Default Parameters -------------------------------------------------------------- ;--I.2.1-- Epoch and Geometry if not keyword_set(EP) then ep = epRead('now') if not keyword_set(elongation) then elongation = 180. ;--I.2.2-- Type of Graph and Output if not keyword_set(absoluteMagnitude) and not keyword_set(year) and $ not keyword_set(histoAE) and not keyword_set(histoAI) and $ not keyword_set(histoBoth) then absoluteMagnitude=1 if not keyword_set(Dump) then dump='./ssoDiscovery.eps' ;--I.2.3-- Histogram Variables if keyword_set(population) then param = ssoDiscovery_params(population) $ else param = ssoDiscovery_params() ;--I.2.4-- Semi-major axis, Eccentricity, and Inclination aDefault = {range:[0.5,6.], bin:0.01, tick:1.0, minor:5} eDefault = {range:[0. ,1.], bin:0.01, tick:0.2, minor:4} iDefault = {range:[0.,50.], bin:0.5 , tick:10 , minor:5} if not keyword_set(aInfo) then aInfo=aDefault else aInfo = updateStructure(aDefault,aInfo) if not keyword_set(eInfo) then eInfo=eDefault else eInfo = updateStructure(eDefault,eInfo) if not keyword_set(iInfo) then iInfo=iDefault else iInfo = updateStructure(iDefault,iInfo) numInfo={range:[0, 1000.]} ;--I.2.5-- Absolute magnitude hRang = {histo: [0,30], $ plot: param.hPlot, $ comp: param.hComp} ; hRang = {histo: [0,25], plot: [5,24], comp: [13,15] } ;-All ; hRang = {histo: [0,25], plot: [10,20], comp: [11.5,16.0] } ;-NEAs ; hRang = {histo: [0,25], plot: [11,18], comp: [12.0,15] } ;-MB binH = 0.5 nbBinH= (hRang.histo[1]-hRang.histo[0]) / binH + 1 ;--I.2.6-- Year of Discovery yDefault = {range:[1800,2030], bin:1 , tick:50 , minor:5} if not keyword_set(yInfo) then yInfo=yDefault else yInfo = updateStructure(yDefault,yInfo) ;--I.2.7-- Language if not keyword_set(language) then language='en' if not strCmp(language,'fr',/Fold) and not strCmp(language,'en',/Fold) then begin message, /ioError, 'Invalid language ('+strTrim(language,2)+"). Only 'en' and 'fr' are allowed." return endif ;--I.3-- Read ssoDiscovery Database ---------------------------------------------------------- if size(discovery,/Type) eq 7 then disco = ssoDiscovery_read(/Dump, discovery=discovery) $ else disco = discovery ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Subset Selection and Statistics -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Select Population ------------------------------------------------------------------ if keyword_set(population) then begin ;--II.1.1-- Test Population vs Main/Sub Dynamical Types population = strtrim(population,2) selTyp = where( strcmp(population,disco.orbit.main,/fold), nbTyp ) selSub = where( strcmp(population,disco.orbit.sub,/fold), nbSub ) ;--II.1.2-- Exception if nbTyp eq 0 and nbSub eq 0 then begin message, 'No match with population "'+population+'"' stop endif ;--II.1.3-- Concatenate Results if nbTyp ne 0 then sel=selTyp if nbSub ne 0 then sel=selSub ;--II.1.4-- Select Only specified population disco = disco(sel) endif ;--II.2-- Build Epoch Array ------------------------------------------------------------------ nbEp = n_elements(EP) jdArr= dblarr(nbEp) for kEp=0, nbEp-1 do jdArr[kEp] = epRead(ep[kEp],/JD) ord = sort(jdArr) jdArr = jdArr[ord] ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Absolute Magnitude Plot -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(absoluteMagnitude) then begin ;--II.3-- Compute Discovery Histograms ------------------------------------------------------- histoH = fltarr(nbEp,nbBinH) cumulH = fltarr(nbEp,nbBinH) for kEp=0, nbEp-1 do begin cur = where( disco.when.jd le jdArr[kEp] ) if cur[0] ne -1 then begin histoH[kEp,*] = histogram( disco[cur].H, binsize=binH, min=hRang.histo[0], max=hRang.histo[1], loc=hArr ) cumulH[kEp,*] = total(histoH[kEp,*],/Cumul) endif endfor ;-----------------------------------------------------------------------------------------------; ;--III.1-- Set Graphics Environment ---------------------------------------------- cgPS_open, Filename=dump, /Metric, /Decomposed, /Portrait, $ xsize=29.7, ysize=18, language_level=2, /quiet ;-----------------------------------------------------------------------------------------------; ;--III.2-- Plot Cumulative Histogram --------------------------------------------- if keyword_set(cumulative) then begin ;--III.2.1-- Plot Last Epoch cgPlot, hArr, cumulH[nbEp-1,*], /Normal, /NoData, $ position=[0.09,0.11,0.975,0.98], $ xTitle='Absolute magnitude (H)', $ xStyle=1, xRange=hRang.plot, xTickInt=1, xMinor=2, $ yTitle='Number of object (cumulative)', $ /yLog, YTickFormat='exponent', yRange=[1,1e6] ;--III.2.2-- Plot Ancient Epochs if keyword_set(ancient) then begin hInPlot = where( hArr ge hRang.plot[0] and hArr le hRang.plot[1]) colAncient=['Light Gray', 'Medium Gray', 'Gray', 'Slate Gray', 'Dark Gray','Charcoal','Black'] colAncient=colAncient[6-(nbEp-1):6] for kEp=0, nbEp-1 do begin cgPlot, /OverPlot, hArr, cumulH[kEp,*], color=colAncient[kEp] endfor ; cgText, hRang.plot[1]-0.2, max(cumulH[nbEp-1,hInPlot]), strMid(epRead(jdArr[nbEp-1],/ISO),0,10), align=1 cgText, hRang.plot[1]-0.2, 2.5, strMid(epRead(jdArr[nbEp-1],/ISO),0,10), align=1 ;--III.2.3-- Plot Other Epochs endif else begin hInPlot = where( hArr ge hRang.plot[0] and hArr le hRang.plot[1]) for kEp=0, nbEp-1 do begin cgPlot, /OverPlot, hArr, cumulH[kEp,*] cgText, hRang.plot[1]-0.2, 2.5, strMid(epRead(jdArr[kEp],/ISO),0,10), align=1 ; cgText, hRang.plot[1]-0.2, max(cumulH[kEp,hInPlot]), strMid(epRead(jdArr[kEp],/ISO),0,10), align=1 endfor endelse ;-----------------------------------------------------------------------------------------------; ;--III.3-- Plot Regular Histogram ------------------------------------------------ endif else begin ;--III.3.1-- Plot Last Epoch cgPlot, hArr, histoH[nbEp-1,*], /Normal, /NoData, $ position=[0.09,0.11,0.975,0.98], $ xTitle='Absolute magnitude (H)', $ xStyle=1, xRange=hRang.plot, xTickInt=1, xMinor=2, $ yTitle='Number of object', $ /yLog, YTickFormat='exponent', yRange=[1,1e5] ;--III.3.2-- Plot Ancient Epochs if keyword_set(ancient) then begin hInPlot = where( hArr ge hRang.plot[0] and hArr le hRang.plot[1]) colAncient=['Light Gray', 'Medium Gray', 'Gray', 'Slate Gray', 'Dark Gray','Charcoal','Black'] colAncient=colAncient[6-(nbEp-1):6] for kEp=0, nbEp-1 do begin cgPlot, /OverPlot, hArr, histoH[kEp,*], color=colAncient[kEp] endfor cgText, hRang.plot[1]-0.2, 2.5, strMid(epRead(jdArr[nbEp-1],/ISO),0,10), align=1 ;--III.3.3-- Plot Other Epochs endif else begin for kEp=0, nbEp-1 do begin cgPlot, /OverPlot, hArr, histoH[kEp,*] endfor cgText, hRang.plot[1]-0.2, 2.5, strMid(epRead(jdArr[nbEp-1],/ISO),0,10), align=1 endelse endelse ;-----------------------------------------------------------------------------------------------; ;--III.4-- Completeness Estimates ------------------------------------------------ ;--III.4.1-- Compute Size-Frequency Power Law if keyword_set(SFD) or keyword_set(vLimit) or keyword_set(hLimit) then begin hSet = where( hArr ge hRang.comp[0] and hArr le hRang.comp[1] ) compParam = linfit( hArr[hSet], alog10(cumulH[nbEp-1,hSet]) ) compModel = 10^(compParam[0] + hArr * compParam[1]) if keyword_set(cumulative) then $ cgPlot, /OverPlot, hArr, compModel, color='Blue', thick=2 endif ;--III.4.2-- Extrapolate the Number of Objects to be Discovered if keyword_set(vLimit) then begin if keyword_set(population) then $ hLimit = round( (vLimit-ssoHtoV(pop=population, elong=elongation))/binH) * binH $ else $ hLimit = round( (vLimit-ssoHtoV(pop='MB', elong=elongation))/binH) * binH endif if keyword_set(hLimit) then begin ;--III.4.2/A-- For each Limiting Magnitude nbLimit=n_elements( hLimit ) for kLim=0, nbLimit-1 do begin ;--III.4.2/B-- Search Limiting Bin pH = where( hArr le hLimit(kLim) and hArr ge hLimit(kLim) ) if ph[0] ne -1 then begin ;--III.4.2/C-- Compute Number to be Discovered compExtrapol = 10^(compParam[0] + hLimit(kLim) * compParam[1]) curVal = cumulH[nbEp-1,pH] nbToDisco = compExtrapol - curVal if nbToDisco lt 0 then nbToDisco=0 ;--III.4.2/D-- Plot Limit and Number cgPlot, /OverPlot, hLimit[kLim]+[0,0], [1e-2,1e7], $ color='Red', linestyle=lineStyleID('--') cgText, hLimit[kLim]-0.1, compExtrapol, strTrim(string(nbToDisco,format='(I)'),2), align=1, color='Red' endif endfor endif ;-----------------------------------------------------------------------------------------------; ;--III.5-- Figure Caption -------------------------------------------------------- if keyword_set(population) then $ cgText, 0.125, 0.90, 'Size-frequency distribution for '+strtrim(population,2), /Normal $ else $ cgText, 0.125, 0.90, 'Size-frequency distribution for all SSOs', /Normal if keyword_set(SFD) then begin ; xyouts, 0.1, 0.85, 'Power Law', /NORM ; xyouts, 0.1, 0.85, 'Power Law of Index '+strtrim(string(compParam[1],format='(F4.1)'),2), /NORM endif ;--III.6-- Close Absolute Magnitude Plot ----------------------------------------- cgPS_close, /png, Delete_PS=0 endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Timeline Plot -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(year) then begin sel = where( disco.when.jd le jdArr[0], nbSel ) if nbSel ge 1 then begin ;--IV.1-- Compute Histogram of Year of Discovery ------------------------------------------- cur=disco[sel] yearArr = round(float(strMid(cur.when.iso,0,4))) cgHistoPlot, [yearArr], minInput=yInfo.range[0], maxInput=yInfo.range[1], $ binSize=yInfo.bin, location=yArr, histData=yNum wdelete nbYear = n_elements( yArr ) cumYear = total(yNum,/Cumulative) ;--IV.2-- Open Graphics -------------------------------------------------------------------- ;--IV.2.1-- EPS cgPS_open, Filename=dump, /Metric, /Decomposed, /Portrait, $ xsize=30, ysize=15, language_level=2, charsize=1.4, /quiet ;--IV.2.2-- Language if strCmp(language,'fr',/Fold_Case) then begin xTitle='Ann'+string(233B)+'e de d'+string(233B)+'couverte' yTitle='Nombre de d'+string(233B)+'couvertes' cTitle='D'+string(233B)+'couvertes cumul'+string(233B)+'es' endif else begin xTitle='Year of discovery' yTitle='Number of discoveries' nTitle='Cumulative discoveries' endelse ;--IV.2.3-- Plot Area cgPlot, yArr, yNum, /NoData, /Normal, /yLog, $ position=[0.1,0.15,0.9,0.95], $ xStyle=1, yStyle=9, yRange=[1, 1e5], $ xTitle=xTitle, yTitle=yTitle ;--IV.3-- Plot Discovery Histogram --------------------------------------------------------- barFill = 0.9 for kY=0, nbYear-1 do begin xx = yArr[kY] + yInfo.bin/2. + [-1,1,1,-1,-1]*BarFill*yInfo.bin/2. yy = yNum[kY] * [0,0,1,1,0] cgColorFill, xx, yy, color='Cornflower Blue' cgPlot, /OverPlot, xx, yy endfor ;--IV.4-- Plot Cumulative Discoveries ------------------------------------------------------ cgPlot, yArr, cumYear, /NoErase, /Normal, /yLog, $ position=[0.1,0.15,0.9,0.95], $ xStyle=1, yStyle=5, $ yRange=[1, 1e6], color='Dark Gray', thick=2 cgAxis, yAxis=1, yRange=[1,1e6], yTitle=cTitle, color='Dark Gray' cgPS_close, /png, Delete_PS=0 endif endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- (a,e) histogram -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(histoAE) then begin ;--V.1-- Select sample --------------------------------------------------------------------- sel = where( disco.when.jd le jdArr[0], nbSel ) if nbSel ge 1 then begin ;--V.2-- Build the 2-D Histogram --------------------------------------------------------- cur = disco[sel] histo = hist_2d( [cur.a], [cur.e], bin1=aInfo.bin, bin2=eInfo.bin, $ min1=aInfo.range[0], max1=aInfo.range[1], $ min2=eInfo.range[0], max2=eInfo.range[1] ) totSSO = total(histo) ;--V.3-- Display the 2-D Histogram ------------------------------------------------------- ;--V.3.1-- Open Graphics cgPS_open, Filename=dump, /Metric, /Decomposed, /Portrait, $ xsize=30, ysize=15, language_level=2, /quiet ;--V.3.2-- Plot Position posArr=[0.085, 0.13, 0.85, 0.98] posBar=[0.89 , 0.13, 0.92, 0.98] ;--V.3.3-- Scaling of the histogram ctIndex=39 scaled = 255*alog10(histo+1)/max(alog10(numInfo.range[1])) zero=where(histo eq 0) if zero[0] ne -1 then scaled[zero]=0 scaled=byte(255-scaled) ;--V.3.4-- Display Histogram cgImage, scaled, /Normal, position=posArr, ctIndex=ctIndex, /NoErase, scale=0 ;--V.4-- Display Axes --------------------------------------------------------------------; ;--V.4.1-- Select Language if strCmp(language,'fr',/Fold_Case) then begin xTitle='Demi-grand axe (UA)' yTitle='Excentricite' nTitle="Nombre d'objets par boite en (a,e)" endif else begin xTitle='Semi-major axis (AU)' yTitle='Eccentricity' nTitle='Number of objects per (a,e) bin' endelse ;--V.4.2-- Draw Plot cgPlot, 0,0, /NoData, /NoErase, /Normal, $ position=posArr, $ xStyle=1, xRange=aInfo.range, xTickInterval=aInfo.tick, xMinor=aInfo.minor, $ yStyle=1, yRange=eInfo.range, yTickInterval=eInfo.tick, yMinor=eInfo.minor, $ xTitle=xTitle, yTitle=yTitle ;--V.5-- Display Color Bar ---------------------------------------------------------------; bar=numInfo.range[0] + findgen(255)*(numInfo.range[1]-numInfo.range[0])/255. bar= transpose([[bar],[bar]]) scalBar = 255*alog10(bar+1)/max(alog10(numInfo.range[1])) zero=where(bar eq 0) if zero[0] ne -1 then scalBar[zero]=0 scalBar= byte(255-scalBar) cgImage, byte(scalBar), /Normal, position=posBar, ctIndex=ctIndex, /NoErase, scale=0 cgPlot, 0,0,/NoErase, /NoData, /Normal, position=posBar, $ xStyle=8, xRange=[0,1], xTickInt=1, xTickName=replicate(' ',10), xMinor=1, $ yStyle=8, yRange=numInfo.Range, yTickName=replicate(' ',10) cgAxis, yAxis=1, yRange=[1,numInfo.Range[1]], /yLog, yTickLen=0.4, yMinor=10 cgText, 0.88, 0.55, orient=90, align=0.5, /Normal, nTitle ;--V.6-- Display the Date ----------------------------------------------------------------; cgText, 0.83, 0.9, align=1, /Normal, strMid(date_conv(jdArr[0],'FITS'),0,10) cgText, 0.83, 0.85, align=1, /Normal, strTrim(string(totSSO,format='(I)'),2) cgPS_close, /png, Delete_PS=0 endif endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- VI -- (a,i) histogram -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(histoAI) then begin ;--VI.1-- Select sample -------------------------------------------------------------------- sel = where( disco.when.jd le jdArr[0], nbSel ) if nbSel ge 1 then begin ;--VI.2-- Build the 2-D Histogram -------------------------------------------------------- cur = disco[sel] histo = hist_2d( [cur.a], [cur.i], bin1=aInfo.bin, bin2=iInfo.bin, $ min1=aInfo.range[0], max1=aInfo.range[1], $ min2=iInfo.range[0], max2=iInfo.range[1] ) totSSO = total(histo) ;--VI.3-- Display the 2-D Histogram ------------------------------------------------------ ;--VI.3.1-- Open Graphics cgPS_open, Filename=dump, /Metric, /Decomposed, /Portrait, $ xsize=30, ysize=15, language_level=2, /quiet ;--VI.3.2-- Plot Position posArr=[0.085, 0.13, 0.85, 0.98] posBar=[0.89 , 0.13, 0.92, 0.98] ;--VI.3.3-- Scaling of the histogram ctIndex=39 scaled = 255*alog10(histo+1)/max(alog10(numInfo.range[1])) zero=where(histo eq 0) if zero[0] ne -1 then scaled[zero]=0 scaled=byte(255-scaled) ;--VI.3.4-- Display Histogram cgImage, scaled, /Normal, position=posArr, ctIndex=ctIndex, /NoErase, scale=0 ;--VI.4-- Display Axes ------------------------------------------------------------------- ;--VI.4.1-- Select Language if strCmp(language,'fr',/Fold_Case) then begin xTitle='Demi-grand axe (UA)' yTitle='Inclinaison (!Uo!N)' nTitle="Nombre d'objets par boite en (a,i)" endif else begin xTitle='Semi-major axis (AU)' yTitle='Inclination (!Uo!N)' nTitle='Number of objects per (a,i) bin' endelse ;--VI.4.2-- Draw Plot cgPlot, 0,0, /NoData, /NoErase, /Normal, $ position=posArr, $ xStyle=1, xRange=aInfo.range, xTickInterval=aInfo.tick, xMinor=aInfo.minor, $ yStyle=1, yRange=iInfo.range, yTickInterval=iInfo.tick, yMinor=iInfo.minor, $ xTitle=xTitle, yTitle=yTitle ;--VI.5-- Display Color Bar -------------------------------------------------------------- bar=numInfo.range[0] + findgen(255)*(numInfo.range[1]-numInfo.range[0])/255. bar= transpose([[bar],[bar]]) scalBar = 255*alog10(bar+1)/max(alog10(numInfo.range[1])) zero=where(bar eq 0) if zero[0] ne -1 then scalBar[zero]=0 scalBar= byte(255-scalBar) cgImage, byte(scalBar), /Normal, position=posBar, ctIndex=ctIndex, /NoErase, scale=0 cgPlot, 0,0,/NoErase, /NoData, /Normal, position=posBar, $ xStyle=8, xRange=[0,1], xTickInt=1, xTickName=replicate(' ',10), xMinor=1, $ yStyle=8, yRange=numInfo.Range, yTickName=replicate(' ',10) cgAxis, yAxis=1, yRange=[1,numInfo.Range[1]], /yLog, yTickLen=0.4, yMinor=10 cgText, 0.88, 0.55, orient=90, align=0.5, /Normal, nTitle ;--VI.6-- Display the Date --------------------------------------------------------------- cgText, 0.83, 0.9, align=1, /Normal, strMid(date_conv(jdArr[0],'FITS'),0,10) cgText, 0.83, 0.85, align=1, /Normal, strTrim(string(totSSO,format='(I)'),2) cgPS_close, /png, Delete_PS=0 endif endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- VII -- (a,i) & (a,e) histograms -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(histoBoth) then begin ;--VII.1-- Select sample ------------------------------------------------------------------- sel = where( disco.when.jd le jdArr[0], nbSel ) if nbSel ge 1 then begin ;--VII.2-- Build the 2-D Histogram ------------------------------------------------------- cur = disco[sel] histoAI = hist_2d( [cur.a], [cur.i], bin1=aInfo.bin, bin2=iInfo.bin, $ min1=aInfo.range[0], max1=aInfo.range[1], $ min2=iInfo.range[0], max2=iInfo.range[1] ) histoAE = hist_2d( [cur.a], [cur.e], bin1=aInfo.bin, bin2=eInfo.bin, $ min1=aInfo.range[0], max1=aInfo.range[1], $ min2=eInfo.range[0], max2=eInfo.range[1] ) totSSO = total(histoAE) ;--VII.3-- Display the 2-D Histogram ----------------------------------------------------- ;--VII.3.1-- Open Graphics cgPS_open, Filename=dump, /Metric, /Decomposed, /Portrait, $ xsize=30, ysize=30, language_level=2, /quiet ;--VII.3.2-- Plot Position posBot=[0.085, 0.0703, 0.85, 0.53] posTop=[0.085, 0.53 , 0.85, 0.98] posBar=[0.89 , 0.085, 0.92, 0.965] ;--VII.3.3-- Scaling of the histogram ctIndex=39 scaledAE = 255*alog10(histoAE+1)/max(alog10(numInfo.range[1])) zero=where(histoAE eq 0) if zero[0] ne -1 then scaledAE[zero]=0 scaledAE=byte(255-scaledAE) scaledAI = 255*alog10(histoAI+1)/max(alog10(numInfo.range[1])) zero=where(histoAI eq 0) if zero[0] ne -1 then scaledAI[zero]=0 scaledAI=byte(255-scaledAI) ;--VII.3.4-- Display Histogram cgImage, scaledAE, /Normal, position=posBot, ctIndex=ctIndex, /NoErase, scale=0 cgImage, scaledAI, /Normal, position=posTop, ctIndex=ctIndex, /NoErase, scale=0 ;--VII.4-- Display Axes ------------------------------------------------------------------ ;--V.4.1-- Select Language if strCmp(language,'fr',/Fold_Case) then begin aTitle='Demi-grand axe (UA)' eTitle='Excentricite' iTitle='Inclinaison (!Uo!N)' nTitle="Nombre d'objets par boite en (a,e) et (a,i)" endif else begin aTitle='Semi-major axis (AU)' eTitle='Eccentricity' iTitle='Inclination (!Uo!N)' nTitle='Number of objects per (a,e) and (a,i) bin' endelse ;--V.4.2-- Draw Plot cgPlot, 0,0, /NoData, /NoErase, /Normal, $ position=posBot, $ xStyle=1, xRange=aInfo.range, xTickInterval=aInfo.tick, xMinor=aInfo.minor, $ yStyle=1, yRange=eInfo.range, yTickInterval=eInfo.tick, yMinor=eInfo.minor, $ xTitle=aTitle, yTitle=eTitle cgPlot, 0,0, /NoData, /NoErase, /Normal, $ position=posTop, $ xStyle=1, xRange=aInfo.range, xTickInterval=aInfo.tick, xMinor=aInfo.minor, $ yStyle=1, yRange=iInfo.range, yTickInterval=iInfo.tick, yMinor=iInfo.minor, $ xTickName=replicate(' ',25), $ yTitle=iTitle ;--VII.5-- Display Color Bar ------------------------------------------------------------- bar=numInfo.range[0] + findgen(255)*(numInfo.range[1]-numInfo.range[0])/255. bar= transpose([[bar],[bar]]) scalBar = 255*alog10(bar+1)/max(alog10(numInfo.range[1])) zero=where(bar eq 0) if zero[0] ne -1 then scalBar[zero]=0 scalBar= byte(255-scalBar) cgImage, byte(scalBar), /Normal, position=posBar, ctIndex=ctIndex, /NoErase, scale=0 cgPlot, 0,0,/NoErase, /NoData, /Normal, position=posBar, $ xStyle=8, xRange=[0,1], xTickInt=1, xTickName=replicate(' ',10), xMinor=1, $ yStyle=8, yRange=numInfo.Range, yTickName=replicate(' ',10) cgAxis, yAxis=1, yRange=[1,numInfo.Range[1]], /yLog, yTickLen=0.4, yMinor=10 cgText, 0.88, 0.55, orient=90, align=0.5, /Normal, nTitle ;--VII.6-- Display the Date -------------------------------------------------------------- cgText, 0.83, 0.95, align=1, /Normal, strMid(date_conv(jdArr[0],'FITS'),0,10) cgText, 0.83, 0.925, align=1, /Normal, strTrim(string(totSSO,format='(I)'),2) cgPS_close, /png, Delete_PS=0 endif endif end