root = '/data/a-types/'
  dir = root + 'mpc/'

  readDISCO=0
  if readDISCO eq 1 or not keyword_set(disco) then $
     disco=ssoDiscovery_read()


;  uM = uniq( disco.orbit.main, sort(disco.orbit.main) )
;  uS = uniq( disco.orbit.sub, sort(disco.orbit.sub) )
;  print, disco[uM].orbit.main
;  print, disco[uS].orbit.sub





  
  popName = ['Hungaria','IMB','MMB','OMB','Cybele']
  nbPop = n_elements(popName)
  pop = replicate({name:'', aMin:0.,aMax:0., hMin:0.,hMax:0., $
                   b:0.51, knee:15.25, f:0.30, h:15., N:200. }, nbPop )
  pop.name = popName

 ;--         HUN    IMB   MMB   OMD   CYB
  pop.aMin=[ 1.80, 2.05, 2.50, 2.78, 3.25]
  pop.aMax=[ 2.05, 2.50, 2.78, 3.25, 3.70]

  pop.hMin=[ 14.0, 12.0, 12.5, 12.5, 13.0 ]
  pop.hMax=[ 17.0, 15.0, 15.0, 15.0, 15.0 ]

  pop.H = [  15.0, 15.0, 15.0, 15.0, 15.0 ]
  pop.N = [   110, 4800, 8200,12000,  256 ]
  pop.f = [  0.51, 0.36, 0.36, 0.36, 0.30 ]  ;- A bit higher than Gladman in general

  binH=0.25



 ;----------------------------------------------------------------------------------------------------
 ;----------------------------------------------------------------------------------------------------
 ;----------------------------------------------------------------------------------------------------
 ;--- Main-Belt Asteroids
  xRang = [8,20]
  cgPlot, 0,0, /NoData, /yLog, /Normal, $
          xTitle='Absolute magnitude (H)', $
          yTitle='Cumulative number N(<H)', $
          xStyle=1, xRange=xRang, $
          yStyle=1, yRange=[1e3,1e6], $
          xMinor=2, xTickInterval=2, $
          yMinor=5, yTickInterval=1;, $
;          position=[0.08,yBot,0.98,yBot+yLen]

  sel = where( strcmp(disco.orbit.main,'MB') or strcmp(disco.orbit.sub,'MB',nbSel ) )
  sfdM = ssoSynthSFD(disco[sel].h, binSize=binH, range=[0,24] )

  cgPlot, sfdM.h, sfdM.cumul.obs,   /Over
  hRef = 15
  nRef = sfdM.histo.obs[ where( sfdM.h eq hRef ) ]
    nRef=nRef[0] ; + 1e4

  modL = ssoModelSFD( sfdM.h, 0.51, 15.25, 0.28, reference={h:hRef,N:nRef} )   ;-Gladman09: 0.5 - then 0.3
  modM = ssoModelSFD( sfdM.h, 0.51, 15.25, 0.30, reference={h:hRef,N:nRef} )   ;-Gladman09: 0.5 - then 0.3
  modH = ssoModelSFD( sfdM.h, 0.51, 15.25, 0.32, reference={h:hRef,N:nRef} )   ;-Gladman09: 0.5 - then 0.3

  cgPlot, sfdM.h, sfdM.cumul.obs,   /Over

  cgPlot, modL.h, modL.cumul.model, /over, color='gray', lineStyle=2
  cgPlot, modM.h, modM.cumul.model, /over, color='gray'
  cgPlot, modH.h, modH.cumul.model, /over, color='gray', lineStyle=2


  setA = where( modM.h ge 12 and modM.h le 20 )
  forprint, modM.h, sfdM.histo.obs, modL.histo.model, modM.histo.model, modH.histo.model, $
            comment='H, nbObs, nbLow, nbModel, nbHigh', $
            format='(F5.2,",",I5,",",I6,",",I6,",",I6)', $
            subset=setA, textout=dir+'CSD-MBA.csv', /Silent
;stop


 ;----------------------------------------------------------------------------------------------------
 ;----------------------------------------------------------------------------------------------------
 ;----------------------------------------------------------------------------------------------------
 ;--- Sub Population
  xRang = [8,20]
  cgPlot, 0,0, /NoData, /yLog, /Normal, $
          xTitle='Absolute magnitude (H)', $
          yTitle='Cumulative number N(<H)', $
          xStyle=1, xRange=xRang, $
          yStyle=1, yRange=[1e2,5e5], $
          xMinor=2, xTickInterval=2, $
          yMinor=5, yTickInterval=1;, $

  cgPlot, modL.h, modL.cumul.model, /over, color='gray', lineStyle=2
  cgPlot, modM.h, modM.cumul.model, /over, color='gray'
  cgPlot, modH.h, modH.cumul.model, /over, color='gray', lineStyle=2

  
  for kPop=0, nbPop-1 do begin

    print, pop[kPop].name
    sel = where( strcmp(disco.orbit.main,pop[kPop].name) or strcmp(disco.orbit.sub,pop[kPop].name),nbSel )

    refH = [ pop[kPop].hMin, pop[kPop].hMax ]
    sfdM = ssoSynthSFD(disco[sel].h, binSize=binH, range=[0,24], reference=refH )


;    hRef = 15
;    nRef = sfdM.histo.obs[ where( sfdM.h eq hRef ) ]
;      nRef=nRef[0]
;
;    modL = ssoModelSFD( sfdM.h, 0.51, 15.25, 0.28, reference={h:hRef,N:nRef} )   ;-Gladman09: 0.5 - then 0.3
;    modM = ssoModelSFD( sfdM.h, 0.51, 15.25, 0.30, reference={h:hRef,N:nRef} )   ;-Gladman09: 0.5 - then 0.3
;    modH = ssoModelSFD( sfdM.h, 0.51, 15.25, 0.32, reference={h:hRef,N:nRef} )   ;-Gladman09: 0.5 - then 0.3
;
    cgPlot, sfdM.h, sfdM.cumul.obs,   /OverPlot, thick=2
    cgPlot, sfdM.h, sfdM.cumul.model, /OverPlot, color='gray', lineStyle=2


    valid = where( sfdM.h ge refH[0] and sfdM.h le refH[1] )
    cgPlot, sfdM.h[valid], sfdM.cumul.obs[valid], /OverPlot, color='red', thick=2


;    hRef = refH[1]
;    hRef = pop[kPop].H
;    nRef = sfdM.histo.obs[ where( sfdM.h eq hRef ) ]
;    pop[kPop].N = nRef[0]
    
    modM = ssoModelSFD( sfdM.h, pop[kPop].b, pop[kPop].knee, pop[kPop].f, $
                        reference={h:pop[kPop].H, N:pop[kPop].N} ) ;-Gladman09: 0.5 - then 0.3
    cgPlot, modM.h, modM.cumul.model, /over, color='gray'


    
;
;    cgPlot, modL.h, modL.cumul.model, /over, color='gray', lineStyle=2
;    cgPlot, modM.h, modM.cumul.model, /over, color='gray'
;    cgPlot, modH.h, modH.cumul.model, /over, color='gray', lineStyle=2


;stop
    

  endfor
  
  

  


end