;
;  statistics used for submission 
;    ***** WRONG as geometry was poorly used *****
;
;;


;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;----What to do
  readDISCO = 0
  onlyYear  = 0
  splitPerY = 0


 ;---------- Directories
  spawn, 'hostname', host
  case strTrim(host) of
     'hyperion': dirEuclid = '/data/euclid/'
     'endymion': dirEuclid = '/home/bcarry/work/data/euclid/'
     else: stop
  endcase

  dirFilter = dirEuclid+'filters/'
  dirTaxo   = dirEuclid+'taxonomy/'
  dirStat = dirEuclid + 'stat/'
  dirFrac = dirEuclid + 'fraction/'
  dirPop  = dirEuclid + 'synthetic/'

;  dirStat = dirEuclid + 'gaia-stat/' ;--- for Gaia only


;--- survey limits
  limEcl = 15
  limGal = 20.
  sLimit = 17     ;-Saturation
  vLimit = 24.5   ;-- YJH -- 24.5 in V, but smearing...
  binH   = 0.25

 ;-geometry of obsevration
  elongation=90
  elongMin=87
  elongMax=110


 ;-- only do plot for discovery vs years
  if onlyYear eq 1 then goto, j2y



  now = 2016
  jdNow = date_conv( string(now,format='(I4)')+'-01-01', 'JULIAN' )
  nbY = 10

;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---  III -- SSO populations                                    -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--III.1-- Read SSO populations
  readcol, dirEuclid+'pop', delimiter=',', $
           format='(A,I1,F4.1,F4.1,F4.1,F4.1,F3.1,F3,A,A)', $
           popName, popMain, popHmin, popHmax, xMin, xMax, yMin, yMax, popC, poS, /Silent
  popName = strtrim(popName,2)
  nbPop=n_elements(popName)

;nbpop=1

  popPrint = where( popMain eq 1, nbPopPrint )
  print, 'Populations considered: ('+strtrim(string(nbPop),2)+'): ', popName

 ;--III.2-- Read SSO discovery archive
  if readDISCO eq 1 or not keyword_set(disco) then $
    disco=ssoDiscovery_read()

  res = where( strCmp( strMid(disco.orbit.sub,0,8), 'resonant', /fold ), nbRes )
  disco[res].orbit.sub = 'Resonant'

 ;--III.3-- Read LSST discoveries
  readcol, dirStat+'lsst', lsstPop, lsstTot, format='(A,F)', /Silent
    lsstPop=strtrim(lsstPop,2)
    lsstY1 = 2021

 ;--III.4-- Read Euclid fraction of known object in Wide survey
  if file_test( dirFrac+'euclid-fraction.csv-ttttt',/read) then begin
    readcol, dirFrac+'euclid-fraction.csv', /Silent, delimiter=',', $
             fpopName, ftNow, $                    ;-Pop + #now
             tSky00, tSky01, tSky02, $             ;-Q25-50-75: obs/sky
             tSky10, tSky11, tSky12, $             ;-Q25-50-75: disco/sky
             fs1, fs2, fs3, $                      ;-Q25-50-75: fraction
             tEuc00, tEuc01, tEuc02, $             ;-Q25-50-75: obs/sky
             tEuc10, tEuc11, tEuc12, $             ;-Q25-50-75: disco/sky
             format='(A,IF,F,F,F,F,F,'+$
                    'F,F,F,'+$
                    'F,F,F,F,F,F)'
    fs = fltarr(nbPop,3)
    fs[*,0]=fs2
    fs[*,1]=fs1
    fs[*,2]=fs3
  endif else begin
    fs = fltarr(nbPop,3)+1.
  endelse

;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   IV -- SSO CSD and discovery                              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;--bookeeping arrays
  discoStat=fltarr(nbPop, 2, 5)            ;-Number to be discovered
  fracI = fltarr( nbPop )                  ;-Frac below/above CUT in LATITUDE
  discoVStime = fltarr(nbPop, 2, 5, nbY)   ;-Number to be disco = f(time)
  nbAllSky = fltarr(nbPop,5)               ;-Total number on sky


 ;--IV.1-- Plot parameters
  cgPS_open, Filename=dirStat+'Euclid-Stat.eps', /metric, /decomposed, /encapsulated, $
             xSize=35, ySize=20.5, language_level=2, /quiet

  xKey = 3
  xLen = 1
  yKey = cgLogGen( nbPopPrint, start=1.5e5, finish=1.5e7)
  yKey2= cgLogGen( nbPopPrint, start=1e5, finish=1e7)
  kPrint=0

  xRang = [2,24]
;  xRang = [10,24]
  cgPlot, 0,0, /NoData, /yLog, /Normal, $
          xTitle='Absolute magnitude (H)', $
          yTitle='Cumulative number of SSO', $
          xStyle=1+8, xRange=xRang, $
          yStyle=1, yRange=[1,1e8], $
          xMinor=2, xTickInterval=2, $
          yMinor=5, yTickInterval=1, $
          position=[0.08,0.105,0.98,0.95]

  dRang = (1329/sqrt(0.20)) * 10^(-0.2 * xRang)
  cgAxis, xAxis=1, xRange=dRang, /xLog, xStyle=1, $
          xTickName=[' ','100','10','1','0.1']
  cgText, /Normal, 0.95,0.965, 'km'
  cgText, /Normal, 0.05,0.965, 'Diameter'

 ;--IV.2-- Loop over Populations
  for kPop=0, nbPop-1 do begin

   ;--IV.2.1-- Plot main population -- current census
    if popMain[kPop] eq 1 then begin
      cgPlot, /OverPlot, xKey+[0,xLen], [yKey[kPrint],yKey[kPrint]], color=popC[kPop]
      cgText, xKey+xLen*1.2, yKey2[kPrint], popName[kPop]

      cgPlot, /OverPlot, xKey+xLen*kPrint/(nbPopPrint-1.)+[0,0], [1e4,5e4], color=popC[kPop]
      cgPlot, /OverPlot, xKey+xLen*kPrint/(nbPopPrint-1.), 5e4, color=popC[kPop], psym='Filled Circle'
      cgPlot, /OverPlot, xKey+xLen*kPrint/(nbPopPrint-1.), 1e4, color=popC[kPop], psym='Open Circle'
    endif


   ;------------------------------------------------------------------------------------------------------------------
   ;--IV.2.2-- Build Synthetic CSD:
    case popName[kPop] of
       
     ;--IV.2.2/A-- Comet (from Snodgrass2011+)
      'Comet': begin
         selSFD = saveSFD

        ;-- SFD parameters
         albedoC = 0.04
         neckH = 17.12
         BPL = -1.92
         SPL = -0.20
         coef = 311.925 * 3.5

        ;-- Build SFD
         bigC = where( selSFD.h le neckH, comp=smallC, nbBC, ncomp=nbSC )
         selSFD.cumul.obs[bigC] = coef*( (1329./sqrt(albedoC)) * 10^(-0.2*selSFD.h[bigC]) )^BPL

         selSFD.cumul.obs[smallC] = coef*( (1329./sqrt(albedoC)) * 10^(-0.2*selSFD.h[smallC]) )^SPL
         selSFD.cumul.obs[smallC] *= selSFD.cumul.obs[bigC[nbBC-1]] / selSFD.cumul.obs[smallC[0]]

         lin = linfit( selSFD.h[bigC], alog10(selSFD.cumul.obs[bigC]) )
         selSFD.cumul.model = 10^(lin[0] + selSFD.h * lin[1])

         nbSel = 100.
         nbHigh=  55.

      end

     ;--IV.2.2/B-- NEAs (from Granvik2016+)
      'NEA': begin
         sel = where( strcmp(disco.orbit.main,popName[kPop]) or strcmp(disco.orbit.sub,popName[kPop]),nbSel )
         ref=[popHmin[kPop],popHmax[kPop]]
         selSFD = ssoSynthSFD(disco[sel].h,     reference=ref, binSize=binH)

         selLow = where( (strcmp(disco.orbit.main,popName[kPop]) or strcmp(disco.orbit.sub,popName[kPop])) and disco.i le limEcl, nbLow)
         selHigh= where( (strcmp(disco.orbit.main,popName[kPop]) or strcmp(disco.orbit.sub,popName[kPop])) and disco.i gt limEcl, nbHigh)
         
        ;--- Read Granvik NEA model
         readcol, '/data/euclid/synthetic/nea-granvik2016.csv',synH,synN,syndN,delimiter=',',/Silent
         synH+=0.25
         synCum = total(synN,/Cumulative)

        ;--- Put NEA model in a selSFD structure
         modulo = selSFD.h mod 0.5
         pSel = where( selSFD.h ge 17.5 and selSFD.h le 24 and modulo eq 0 )
         pSyn = where( synH ge 17.5 and synH le 24 )

         selSFD.cumul.model = selSFD.cumul.obs
         selSFD.cumul.model[pSel] += synCum[pSyn]
         
         pInter = where( selSFD.h ge 17.5 and selSFD.h le 24 and modulo eq 0.25 )
         interp = interpol( selSFD.cumul.model[pSel], selSFD.h[pSel], selSFD.h[pInter] )
         selSFD.cumul.model[pInter]=interp

      end

      
     ;--IV.2.2/X-- Default: from ASTORB and PowerLaw
      else: begin
        sel = where( strcmp(disco.orbit.main,popName[kPop]) or strcmp(disco.orbit.sub,popName[kPop]),nbSel )
        ref=[popHmin[kPop],popHmax[kPop]]
        selSFD = ssoSynthSFD(disco[sel].h,     reference=ref, binSize=binH)
        saveSFD = selSFD

        selLow = where( (strcmp(disco.orbit.main,popName[kPop]) or strcmp(disco.orbit.sub,popName[kPop])) and disco.i le limEcl, nbLow)
        selHigh= where( (strcmp(disco.orbit.main,popName[kPop]) or strcmp(disco.orbit.sub,popName[kPop])) and disco.i gt limEcl, nbHigh)
      end

    endcase
    fracI[kPop] = 1.*nbHigh/nbSel

   ;--IV.2.4-- OverPlot synthetic distribution
    if popMain[kPop] eq 1 then begin
      cgPlot, /OverPlot, selSFD.h, selSFD.cumul.model, color=popC[kPop], linestyle=2
      cgPlot, /OverPlot, selSFD.h, selSFD.cumul.obs,   color=popC[kPop]
    endif


   ;--IV.2.5-- Convert V limiting magnitude into H 
    minH2V = ssoHtoV(pop=popName[kpop], elong=elongMin)
    avgH2V = ssoHtoV(pop=popName[kpop], elong=elongation)
    maxH2V = ssoHtoV(pop=popName[kpop], elong=elongMax)
    joint= [minH2V, avgH2V, maxH2V]
    hLimit = round( (vLimit-joint)/binH) * binH

   ;-remove redondance
    uH = uniq( hLimit, sort(hLimit) )
    hLimit=hLimit[uH]

   ;-ultra-clean
;    goodH = where( hLimit ge 10 and hLimit le 30 )
    goodH = where( hLimit le 30 and finite(hLimit) and hLimit ge popHmax[kPop] )
    hLimit = hLimit[goodH]


    
    nbLimit=n_elements( hLimit )
    n2disco=fltarr(nbLimit)
    n2past =fltarr(nbLimit,nbY)+0.
    n2obs = fltarr(nbLimit)


   ;--IV.2.6-- Compute N to discover now and in the past
    for kLim=0, nbLimit-1 do begin

     ;--IV.2.6/A-- At current epoch
      pH = where( selSFD.h eq hLimit[kLim] )
      nbToDisco = selSFD.cumul.model[pH[0]] - selSFD.cumul.obs[pH[0]]
      n2disco[kLim] = 1.*nbToDisco
      n2obs[kLim] = selSFD.cumul.model[pH[0]]
      


      if popMain[kPop] eq 1 then begin
;        cgPlot, /OverPlot, hLimit[kLim], selSFD.cumul.obs[ph[0]], psym='Open Circle', color=popC[kPop]
;        cgPlot, /OverPlot, hLimit[kLim], selSFD.cumul.model[ph[0]], psym='Filled Circle', color=popC[kPop]
;        cgPlot, /OverPlot, hLimit[[kLim,kLim]], [selSFD.cumul.obs[ph[0]], selSFD.cumul.model[ph[0]]], $
;                color=popC[kPop]

       ;--IV.2.6/B-- In the last decade
        if ~strCmp( popName[kPop],'Comet' ) and splitPerY eq 1 then begin
          for kY=0, nbY-1 do begin
            old = now-nbY+1+kY
            jdOld = date_conv( string(old,format='(I4)')+'-01-01', 'JULIAN' )
            oldPop = where( (strcmp(disco.orbit.main,popName[kPop]) or strcmp(disco.orbit.sub,popName[kPop])) and $
                            disco.when.jd le jdOld ,nbOld )
            oldSFD = ssoSynthSFD(disco[oldPop].h,reference=ref,binSize=binH)
            
            nbToDisco = selSFD.cumul.model[pH[0]] - oldSFD.cumul.obs[pH[0]]
            n2past[kLim,kY] = 1.*nbToDisco
  
;            print, kPop, kLim, hLimit[kLim], ky, old, nbold, n2past[kLim,kY]
  
  
          endfor
        endif
  
        if strcmp( popName[kPop],'Comet' ) and splitPerY eq 1 then begin
          for kY=0, nbY-1 do n2past[kLim,kY] = 1.*nbToDisco
        endif






        
      endif

    endfor

    if popMain[kPop] eq 1 then begin
      sizBox = 0.2



     ;---limiting magnitude
      hLim = dataBox( hLimit)
      qLabel=' '+['0', '25', '50', '75', '100']+'%'
      qLabel=['0', '25', '50', '75', '100']
      for k=0,4 do begin
        pH = where( selSFD.h eq hLim[k] )
        xQ = selSFD.h[[pH, pH]]

        if k eq 0 then begin
          yQ = [selSFD.cumul.obs[pH],selSFD.cumul.model[pH]]
          yRef = selSFD.cumul.model[pH]
          x1 = selSFD.h[pH]
       endif else begin
          dN = (selSFD.cumul.model[pH]-yRef)*(4-k)*0.25
          yQ = [selSFD.cumul.obs[pH],yRef+dN]

          subH = where( selSFD.h ge hLim[k-1] and selSFD.h le hLim[k])
;          help, subH

          x2 = selSFD.h[pH]
          curv = interpol( [yRef,yRef+dN], [x1,x2], selSFD.h[subH] )
          cgPlot, /OverPlot, selSFD.h[subH], curv, color=popC[kPop], lineStyle=4
          x1=x2

          yRef += dN
          print, yRef
       endelse     

;        yQ = [selSFD.cumul.obs[pH],selSFD.cumul.model[pH]*((4-k)*0.25)]

        if xQ[0] ge 4 and xQ[1] le 24 then begin
;          cgPlot, /OverPlot, xQ, yQ, color=popC[kPop]
;          cgPlot, /OverPlot, xQ[0], yQ[0], color=popC[kPop], psym='Filled Circle'
          cgPlot, /OverPlot, xQ[1], yQ[1], color=popC[kPop], psym='Open Circle'
        endif

        if k eq 4 then begin
          cgPlot, /OverPlot, xQ[0], yQ[0], color=popC[kPop], psym='Filled Circle'
          cgPlot, /OverPlot, xQ, yQ, color=popC[kPop]
       endif

;        if xQ[0] ge 4 and xQ[1] le 24 and yQ[1] le 1e8 then begin
;nice           cgText, xQ[0]+0.1, yQ[1], qLabel[4-k], color=popC[kPop], orient=270, charSize=1.25,align=0
;           cgText, xQ[0]-0.1, yQ[1], qLabel[4-k], color=popC[kPop], charSize=1.25,align=1
;        endif

     endfor


     ;-saturation limit
;      hLimit2 = round( (sLimit-joint)/binH) * binH
;      ok  = where( hLimit2 gt 4, nbh)
;      hLimit2=hLimit2[ok]
;
;      if nbh ge 2 then begin
;        lLim = dataBox( hLimit2 )
;       ;---- saturation limit
;        qLabel=' '+['0', '25', '50', '75', '100']+'%  '
;  ;      qLabel=['0', '25', '50', '75', '100']
;        for k=0,4 do begin
;  
;          pH = where( selSFD.h eq lLim[k] )
;          xQ = selSFD.h[[pH, pH]]
;          if k eq 0 then minQ = selSFD.cumul.obs[pH]
;          yQ = [selSFD.cumul.obs[pH],minQ]
;          cgPlot, /OverPlot, xQ, yQ, color=popC[kPop]
;          cgPlot, /OverPlot, xQ[0], yQ[0], color=popC[kPop], psym='Filled Circle'
;;          cgPlot, /OverPlot, xQ[1], yQ[1], color=popC[kPop], psym='Filled Circle'
;          cgText, xQ[0]-0.1, yQ[0], qLabel[k], color=popC[kPop], orient=90, charSize=1.25,align=1
;        endfor
;
;     endif
;
      
    endif


   ;--IV.2.7-- Compute quartiles
    g = where( n2disco ge 0, nbG )
    if nbG ge 2 then begin
      discoStat[kPop,0,*] = dataBox(n2disco[g])
    endif else begin
      discoStat[kPop,0,*] = n2disco[g]
    endelse
    discoStat[kPop,1,*] = discoStat[kPop,0,*] * fracI[kPop]
    nbAllSky[kPop,*] = dataBox( n2obs )

    if splitPerY eq 1 and  popMain[kPop] eq 1 then begin
;    if splitPerY eq 1 and  popMain[kPop] eq 1 and ~strcmp( popName[kPop],'Comet' ) then begin
      for kY=0, nbY-1 do begin
         print, n2past[*,ky]
        discoVStime[kPop,0,*,kY] = dataBox( n2past[*,kY] )
        discoVStime[kPop,1,*,kY] = discoVStime[kPop,0,*,kY] * fracI[kPop]
        print, now-nbY+1+kY, reform(discoVStime[kPop,0,*,kY])
      endfor
    endif
;    cgPlot, now-nbY+1+indgen(nbY), discoVStime[kPop,0,2,*], yRange=[1,1e6],/yLog
;    cgPlot, /OverPlot, now-nbY+1+indgen(nbY), discoVStime[kPop,0,1,*], linestyle=2
;    cgPlot, /OverPlot, now-nbY+1+indgen(nbY), discoVStime[kPop,0,3,*], linestyle=2
;    cgPlot, /OverPlot, now-nbY+1+indgen(nbY), discoVStime[kPop,0,0,*], linestyle=1
;    cgPlot, /OverPlot, now-nbY+1+indgen(nbY), discoVStime[kPop,0,4,*], linestyle=1
;    stop



   ;--IV.2.8-- GUI - Pop + Quartiles outputs
    print, popName[kPop], fracI[kPop], $
           nbAllSky[kPop,0], nbAllSky[kPop,1], $
           nbAllSky[kPop,2], $
           nbAllSky[kPop,3], nbAllSky[kPop,4], $
           format='(A-22,2x,F6.4,2x,E8.1,2x,E8.1,2x,E8.1,2x,E8.1,2x,E8.1)'


;---    stop
;---    stop
;---    stop



    
    if popMain[kPop] eq 1 then kPrint++
 endfor

  cgText, xKey+xLen*1.2, 3e4, 'Margin for'
  cgText, xKey+xLen*1.2, 1e4, 'discovery'


  cgPS_close, /PNG, delete_PS=0
stop
  forprint, popName, fracI, $
            nbAllSky[*,0], nbAllSky[*,1], nbAllSky[*,2], nbAllSky[*,3], nbAllSky[*,4], $
            format='(A-20,",",F6.4,",",E8.1,",",E8.1,",",E8.1,",",E8.1,",",E8.1)', $
            textout=dirStat+'euclid-stat.csv', /Silent, $
            comment='Population, f_i, Min, Q25, Q50, Q75, Max'

  forprint, popName, fracI, $
            discoStat[*,0,0], $
            discoStat[*,0,1], $
            discoStat[*,0,2], $
            discoStat[*,0,3], $
            discoStat[*,0,4], $
            format='(A-20,",",F6.4,",",E8.1,",",E8.1,",",E8.1,",",E8.1,",",E8.1)', $
            textout=dirStat+'euclid-disco.csv', /Silent, $
            comment='Population, f_i, Min, Q25, Q50, Q75, Max'



  read, tCont, prompt='Continue to year graphics?'
  if tCont ne 1 then stop 





;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    V -- SSO Discovery VS YEARS                             -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

j2y:

  if binH ge 0.5 then begin
     print, 'attention!!!! '
     stop
  endif

 ;--V.1-- Plot parameters
  cgPS_open, Filename=dirStat+'Euclid-Year.eps', /metric, /decomposed, /encapsulated, $
             xSize=35, ySize=18.9, language_level=2, /quiet

  yArr = now - reverse(indgen(nbY))

  cgPlot, yArr, discoVStime[0,0,2,*], /yLog, /NoData, /Normal, $
          yStyle=1, yRange=[1,1e6], $
          xStyle=1, xRange=[2005,2025], $
          xTitle='Time', $
          yTitle='Estimation of undiscovered population', $
          xMinor=5, xTickInterval=5, $
          yMinor=5, yTickInterval=1, $
          position=[0.08,0.105,0.96,0.98]

  xKey = 2017
  xLen = 1
  xShi = 3.5
  yKey = 3.5e5
  yTxt = 3e5

  cgPlot, /OverPlot, xKey+[0,xLen]     , [yKey,yKey], linestyle=0, thick=2
  cgPlot, /OverPlot, xKey+[0,xLen]+xShi, [yKey,yKey], linestyle=2, thick=2
  cgText, xKey+xLen*1.2     , yTxt, 'Nominal'
  cgText, xKey+xLen*1.2+xShi, yTxt, 'With LSST'


 ;--V.2-- Plot each Population/Year statistics
  for kPop=0, nbPop-1 do begin
    if popMain[kPop] eq 1 then begin
      yArr = now - reverse(indgen(nbY))
      case popName[kPop] of
         'Trojan':  yArr+=0.125
         'MC':      yArr-=0.125
         'KBO':     yArr+=0.125
         'Centaur': yArr-=0.125
         'Comet':   yArr+=0.325
  ;       'NEA':     yArr+=0.075
         else:
      endcase
  
  
      mode=1
      print, kpop, popName[kPop], fs[kPop,1], format='(I2,2x,A-10,2x,F4.1)'
      discoVStime[kPop,1,*,*] = discoVStime[kPop,0,*,*] * (fs[kPop,1]/100.)

      cgErrPlot, yArr, discoVStime[kPop,mode,0,*], discoVStime[kPop,mode,4,*], color=popC[kPop]
  
  
  
      if ~strcmp(popName[kPop], 'Comet') then begin
        dA = regress( yArr, reform( discoVStime[kPop,mode,2,*]), const=dB )
      endif else begin
        dB = discoVStime[kPop,mode,2,0]
        dA=0
      endelse
  
      nbTot=190
      step=0.1
      yTot = 2006+ step*findgen(nbTot)
  
      expect = dA[0]*yTot+dB
      withLSST = expect
  
     
      selLSST = where( strcmp(lsstPop, popName[kPop]) )
      lsstPerYer = lsstTot[selLSST[0]] / 10.
      lsstStat = 0.1*indgen(100) * lsstPerYer
  
      start = where( yTot eq lsstY1 )
      withLSST[start:nbTot-1] -= lsstStat[0: (nbTot-1-start)]
  
      zero=where( withLSST le 0 )
      if zero[0] ne -1 then withLSST[zero]=1
  
      cgPlot, /OverPlot, yTot, expect,   linestyle=0, color=popC[kPop], thick=2
      cgPlot, /OverPlot, yTot, withLSST, linestyle=2, color=popC[kPop], thick=2
;      cgPlot, /OverPlot, yTot, withLSST, psym='open circle', color=popC[kPop], thick=2
  
  
;      print, popName[kPop], da, db
  
      bLen = 0.2
  
      for kY=0, nbY-1 do begin
  
        if discoVStime[kPop,mode,1,kY] lt 1 then discoVStime[kPop,mode,1,kY]=1
        boxX = bLen*0.5*[-1,1,1,-1,-1] + yArr[kY]
        boxY = [ discoVStime[kPop,mode,1,kY], discoVStime[kPop,mode,1,kY], discoVStime[kPop,mode,3,kY], discoVStime[kPop,mode,3,kY], discoVStime[kPop,mode,1,kY]]
  
        cgColorFill, boxX, boxY, color=popC[kPop]
        cgPlot, /OverPlot, boxX, boxY
        cgPlot, /OverPlot, bLen*0.5*[-1,1]+yArr[kY], [discoVStime[kPop,mode,2,kY],discoVStime[kPop,mode,2,kY]]
  
        print, yArr[ky], reform(discoVStime[kPop,mode,*,kY])
      endfor
  
      xPop = 2006.7
      case popName[kPop] of
         'MC':      popShift = +3e3
         'Comet':   begin
            ref = -1.2*( dA[0]*2006.5+dB)
            popShift = ref + 1.2e3
         end
         else:      popShift = 0
      endcase
      cgText, xPop, 1.2*( dA[0]*2006.5+dB)+popShift, popName[kPop], align=1, color=popC[kPop]
    endif
  endfor


  cgPS_close, /PNG, delete_PS=0

 ;-compute fraction SSO with i>limEcl
 ;-compute fraction of time when asteroid is above i on its orbit
 ;-for each population, what's the H limit based on Euclid observing mode
 ;-a+b=c --> number of objects per population (roughly)

 ;- VIS < 24.5
 ;- NIS < 24 in Y -> V <~25
 ;              J -> V < 25.1
 ;              H -> V < 25.4 

end