; docformat = 'rst' ; ; NAME: ; ssoHtoV ; PURPOSE: ; Compute the typical V-H for a given geometry or for a population of SSO. ;+ ; :Description: ; Compute the typical V-H for a given geometry or for a population of SSO. ; ; :Categories: ; Database, Asteroid ; ; :Returns: The V-H for the requested geometry or population of SSO. ; ; :Params: ; dSUN: in, optional, type=float ; Distance to the sun, in ua ; dOBS: in, optional, type=float ; Distance to the observer, in ua ; phase: in, optional, type=float ; Phase angle, in degree ; ; :Keywords: ; population: in, optional, type=string ; SSO population to be plotted (any type/subtype from ; http://vo.imcce.fr/webservices/skybot/?documentation) ; elongation: in, optional, type=float, default=180. ; The Solar elongation (degree) at which the V-H is desired. ; stepDist: in, optional, type=float, default=0.25 ; The step of heliocentric distance for probability distribution (AU) ; stePhase: in, optional, type=float, default=0.5 ; The step of phase angle for probability distribution (degree) ; hObs: in, optional, type=float, default=1.0 ; Heliocentric distance of the observer (AU) ; geometry: out, optional, type=structure ; Contain the distances to sun and observer and phase angle, ; for a given population and elongation. ; eproc: in, optional, type=string, default=(from initIDL) ; The EPROC directory in the astrodata tree ; configuration: 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:: ; [Asteroid Orbits] ; disco.imcce = PATH_TO_YOUR_IMCCE_DISCO_FILE ; eproc = PATH_TO_YOUR_BASIC_EPROC_INFORMATION ; ; :Examples: ; Compute the typical V-H of a Jupiter Trojan at opposition:: ; IDL> print, ssoHtoV(pop='Trojan', elong=180.) ; ; Compute the V-H for a particular geometry: heliocentric distance ; of 3 au, range to observer, 2 au, phase angle 10 degrees:: ; IDL> print, ssoHtoV( 3., 2., 10. ) ; ; :Uses: ; initIDL, readcol ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Written in October 2014, B. Carry (IMCCE) ; 2016 Feb. - B. Carry (OCA) - idl2 compile option added ; 2016 Apr. - B. Carry (OCA) - Improved peri/helio/mean distance to full range ; 2016 Sep. - B. Carry (OCA) - Corrected bug that changed elongation upon output ; 2016 Nov. - B. Carry (OCA) - Updated KBO sub-classes ; 2017 Mar. - B. Carry (OCA) - Now use a probabilitic approach. Require ssoHtoV_build being ran first. ;- function ssoHtoV, dSun, dObs, phase, population=population, elongation=elongation, $ stePhase=stePhase, stepDist=stepDist, stepMag=stepMag, hObs=hObs, $ geometry=geometry, eproc=eproc, configuration=configuration COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- General Case of Populations ---------------------------------------------------------; if not keyword_set(dSun) and not keyword_set(dObs) then begin mode = 1 ;--I.1.1-- Validation of I/O if not keyword_set(population) then begin message, /ioError, 'Syntax: HtoV = ssoHtoV( dSun, dObs [, phase]) or HtoV = ssoHtoV( population= [, elongation=]) ' return, -1 endif else population=strTrim(population,2) ;--I.1.2-- Elongation: Set/Unset? if not keyword_set(elongation) then elong = !PI else begin dimE = size(elongation) ;--I.1.3-- Elongation: Single value if dimE[0] lt 2 then begin nbElong = 1 elong = !DTOR*elongation[0] ;--I.1.4-- Elongation: PDF endif else begin nbElong = dimE[1] elong = !DTOR*elongation[*,0] pdfElong = elongation[*,1] endelse endelse ;--I.2-- Specific Case with a given Geometry -------------------------------------------------; endif else begin mode = 0 ;--I.2.1-- Validation of I/O if not keyword_set(dSun) or not keyword_set(dObs) then begin message, /ioError, 'Syntax: HtoV = ssoHtoV( dSun, dObs [, phase]) or HtoV = ssoHtoV( population= [, elongation=]) ' return, -1 endif ;--I.2.2-- Phase Angle: Set/Unset? if not keyword_set(phase) then phasArr = 0.001 * !DTOR $ else phasArr = phase * !DTOR endelse ;--I.3-- Check Eproc Directory (for mode=1) --------------------------------------------------; if mode eq 1 then begin if not keyword_set(eproc) then begin if not keyword_set(configuration) then config = initIDL(/Catalog) $ else config = initIDL(configuration, /Catalog) dirPDF = config.eproc endif else dirPDF=eproc dirPDF += 'Heliodist/' endif ;--I.4-- Heliocentric Distance Steps, Phase Array --------------------------------------------; if not keyword_set(hOBs) then hObs = 1.00 if not keyword_set(stepDist) then stepDist= 0.25 if not keyword_set(stePhase) then stePhase= 1.0 if not keyword_set(stepMag) then stepMag = 0.1 if not keyword_Set(phasArr) then begin nbPhase = 180./stePhase phasArr = stePhase*!DTOR * findgen(nbPhase) endif else nbPhase=n_elements(phasArr) ;--I.5-- Sufixes, Bowell G -------------------------------------------------------------------; sufCSV='.csv' sufEPS='.eps' G = 0.15 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Define Parameters per Population -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if mode eq 0 then begin ;--II.1-- Compute Phase Function ----------------------------------------------------------- phi1 = exp( -3.33 *tan(phasArr/2.)^0.63 ) phi2 = exp( -1.87 *tan(phasArr/2.)^1.22 ) fPhase = -2.5*alog10( (1-G)*phi1 + G*phi2 ) ;--II.2-- Compute Distance Function -------------------------------------------------------- fDist = 2.5*alog10( (dObs*dObs) * (dSun*dSun) ) ;--II.3-- Return (V-H) if Single Value ----------------------------------------------------- return, fPhase + fDist ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Define Parameters per Population -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; endif else begin ;--III.1-- PDF of Heliocentric Distances ------------------------------------------------ file = dirPDF+population+'-'+string(stepDist,format='(F05.3)')+sufCSV if not file_test( file,/read) then begin message, 'Cannot find the PDF of Heliocentric Distances for the '+population+' population, '+$ "Please run the ssoHtoV_build command: IDL> ssoHtoV_build, '"+population+"'" return, -2 endif ;--III.2-- Read the PDF ----------------------------------------------------------------- readcol, file, dSun, pdfSun, cdfSun, delimiter=',', /Silent nbSun = n_elements(dSun) ;--III.3-- Single Elongation Provided --------------------------------------------------- pdfObs = dSun *0. pdfPhase= phasArr *0. if nbElong eq 1 then begin ;--III.3.1-- Loop Over Heliocentric Distances for kSun=0, nbSun-2 do begin ;--III.3.2-- Build Min/Max Intervals in Elongation and Heliocentric Distances vElong = elong + [-1,1]*stePhase*!DTOR vSun = dSun[ [kSun,kSun+1] ] ;--III.3.3-- PDF of the range to observer dObs = hObs*cos(vElong) + sqrt( (hObs*cos(vElong))^2. - hObs^2. + vSun^2. ) minH = min( round(dObs/stepDist), max=maxH, /Nan ) if minH ge 0 and maxH ge 0 then pdfObs[ minH:maxH ] += pdfSun[kSun] ;--III.3.4-- PDF of the phase angle comPhase = abs(asin( (hObs*sin(vElong) # (1/vSun) ) ) )/!DTOR minP = min( round(comPhase/stePhase), max=maxP, /Nan ) if minP ge 0 and maxP ge 0 then pdfPhase[ minP:maxP ] += pdfSun[kSun] endfor pdfObs /= total(pdfObs) pdfPhase /= total(pdfPhase) ;--III.4-- Elongation is a PDF itself --------------------------------------------------- endif else begin ;--III.4.1-- Loop Over Heliocentric Distances for kSun=0, nbSun-2 do begin for kElong=0, nbElong-2 do begin ;--III.4.2-- Build Min/Max Intervals in Elongation and Heliocentric Distances vElong = elong[ [kElong, kElong+1] ] vSun = dSun[ [kSun,kSun+1] ] ;--III.4.3-- PDF of the range to observer dObs = hObs*cos(vElong) + sqrt( (hObs*cos(vElong))^2. - hObs^2. + vSun^2. ) minH = min( round(dObs/stepDist), max=maxH, /Nan ) if minH ge 0 and maxH ge 0 then pdfObs[ minH:maxH ] += pdfSun[kSun] ;--III.4.4-- PDF of the phase angle comPhase = abs(asin( (hObs*sin(vElong) # (1/vSun) ) ) )/!DTOR minP = min( round(comPhase/stePhase), max=maxP, /Nan ) if minP ge 0 and maxP ge 0 then pdfPhase[ minP:maxP ] += pdfSun[kSun] endfor endfor pdfObs /= total(pdfObs) pdfPhase /= total(pdfPhase) endelse dObs = dSun nbObs=nbSun ; ; ; ;--II.3-- Compute Geometry ------------------------------------------------------------- ; sel=[1,2,3] ; sel=indgen(5) ; nbSel=n_elements(sel)^2. ; ; ;--II.3.1-- Perihelion, Aphelion, Average Distance ; p = reform( a[sel] # ( 1. - e[sel] ) , nbSel) ; q = reform( a[sel] # ( 1. + e[sel] ) , nbSel) ; m = reform( a[sel] # ( 1. + 0.5* e[sel]^2. ), nbSel) ; ; ;--II.3.2-- Heliocentric Distance and Range to Observer ; dSun = reform([[p],[m],[q]],3*nbSel) ; dObs = cos(elong) + sqrt( cos(elong)^2. -1. + dSun^2. ) ; ; ;--II.3.3-- Phase Angle ; phase = abs(asin( sin(elong)/dSun )) ; ; ;--II.4-- Exporting Observing Geometry ; nbGeo = n_elements(phase) ; geometry = replicate({dOBs:0.,dSun:0.,phase:0.},nbGeo) ; geometry.dObs = dObs ; geometry.dSun = dSun ; geometry.phase = phase/!DTOR ; endelse ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Compute the PDF of the (V-H) Index -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;-IV.1-- Define (V-H) Range ------------------------------------------------------------------; minHV = -5 maxHV = 20 nbHV = (maxHV-minHV)/stepMag HV = minHV + findgen(nbHV)*stepMag pdfMag = HV*0. ;-IV.2-- Loop over Phase, dObs, dSun ---------------------------------------------------------; minProba = 0.001 okPha = where( pdfPhase ge max(pdfPhase)*minProba, nbPhaC ) okObs = where( pdfObs ge max(pdfObs) *minProba, nbObsC ) okSun = where( pdfSun ge max(pdfSun) *minProba, nbSunC ) for kPhase=0, nbPhaC-2 do begin for kObs=0, nbObsC-2 do begin for kSun=0, nbSunC-2 do begin ;--IV.2.1-- Build Min/Max Intervals vPhase = phasArr[ okPha[[kPhase, kPhase+1]] ] vSun = dSun[ okSun[[kSun,kSun+1]] ] vObs = dObs[ okObs[[kOBs,kObs+1]] ] ;--IV.2.2-- Compute Phase Function phi1 = exp( -3.33 *tan(vPhase/2.)^0.63 ) phi2 = exp( -1.87 *tan(vPhase/2.)^1.22 ) fPhase = -2.5*alog10( (1-G)*phi1 + G*phi2 ) ;--IV.2.3-- Compute Distance Function fDist = 2.5*alog10( (vObs*vObs) # (vSun*vSun) ) ;--IV.2.4-- Range of (V-H) index index = fltarr(2,2,2) index[*,*,0] = fDist + fPhase[0] index[*,*,1] = fDist + fPhase[1] ;--IV.2.4-- Compute the PDF minI = min( index, max=maxI, /Nan ) pHV = where( HV ge minI and HV le maxI ) if pHV[0] ne -1 then pdfMag[pHV] += ( pdfSun[okSun[kSun]] * pdfObs[okObs[kObs]] * pdfPhase[okPha[kPhase]] ) endfor endfor endfor pdfMag /= total( pdfMag ) ; mr = [-2,20] ; cgPlot, HV, pdfMag, xr=mr, psym=10, yst=9 ; cgPlot, /noErase, HV, total(pdfMag,/cum)/total(pdfMag), xr=mr, yr=[0,1], color='cornflowerBlue', yst=5 ; cgAxis, yAxis=1, yStyle=1, yRange=[0,1], yTitle='Cumulative probability', color='Cornflower Blue' return, {mag:HV, pdf:pdfMag} end ;vphase = findgen(170) *!DTOR ;phi1 = exp( -3.33 *tan(vPhase/2.)^0.63 ) ;phi2 = exp( -1.87 *tan(vPhase/2.)^1.22 ) ;fPhase = -2.5*alog10( (1-G)*phi1 + G*phi2 ) ;cgplot, vphase/!DTOR, fphase, xtitle='phase', ytitle='phase function' ; ;vobs=0.1*indgen(100) ;vSun=vobs+1 ;fDist = 2.5*alog10( (vObs*vObs) # (vSun*vSun) ) ;cgContour, fDist, vObs, vSun, levels=indgen(20), xr=[0,max(vObs)], yr=[0,max(vSun)] ;cgContour, fDist, vObs, vSun, levels=indgen(5)-5, color='red',/overplot