;
;  statistics used for revision
;


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

 ;----What to do
  readDISCO = 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


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

  hFrac=[1.00,0.75,0.50,0.25,0.00]
  lFrac=reverse(hFrac)

;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  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)

  ecStat = fltarr(nbPop,2)
  
  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


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

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

  yBot = 0.05
  ySep = 0.03
  yTop = 0.97
  yLen = (yTop-yBot-ySep)/2.; + yBot
  yTop2= 0.977

  
  xKey = 3
  xLen = 1
  yKey = cgLogGen( 5, start=2.7e5, finish=2.7e7)
  yKey2= cgLogGen( 5, start=2e5, finish=2e7)
  kPrint=0


  yKeyLine = cgLogGen( 8, start=8.7e1, finish=4.7e4)
  yKeyLine2= cgLogGen( 8, start=8e1, finish=4e4)



  yBotDisco =  4
  yTopDisco = 20 
  yTxtDisco = 12
  xDisco = 19
  
  xRang = [2,24]
  cgPlot, 0,0, /NoData, /yLog, /Normal, $
          xTitle='Absolute magnitude (H)', $
          yTitle='Cumulative number N(<H)', $
          xStyle=1  , xRange=xRang, $
          yStyle=1, yRange=[1,1e8], $
          xMinor=2, xTickInterval=2, $
          yMinor=5, yTickInterval=1, $
          position=[0.08,yBot,0.98,yBot+yLen]


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


     if kPop eq 4 then begin

        cgText, xKey+xLen*0.5, yKey2[kPrint], 'Outer Solar System'

        cgText, xDisco+xLen*1.2, yTxtDisco, 'Margin for'
        cgText, xDisco+xLen*1.2, yBotDisco, 'discovery'


        cgPlot, 0,0, /NoData, /NoErase, /yLog, /Normal, $
               yTitle='Cumulative number N(<H)', $
               xStyle=1+8, xRange=xRang, $
               yStyle=1, yRange=[1,1e8], $
;               xTickName=replicate(' ',40), $
               xMinor=2, xTickInterval=2, $
               yMinor=5, yTickInterval=1, $
               position=[0.08,yBot+yLen+ySep,0.98,yTop]
        kPrint=0

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

      cgPlot, /OverPlot, xDisco+xLen*kPrint/(4.)+[0,0], [yBotDisco,yTopDisco], color=popC[kPop]
      cgPlot, /OverPlot, xDisco+xLen*kPrint/(4.), yTopDisco, color=popC[kPop], psym='Open circle'
      cgPlot, /OverPlot, xDisco+xLen*kPrint/(4.), yBotDisco, color=popC[kPop], psym='Filled Circle'
    endif


   ;------------------------------------------------------------------------------------------------------------------
   ;--IV.2.2-- Build Synthetic CSD:
    case popName[kPop] of
       
     ;--IV.2.2/A-- 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, range=[0,24])

        
        ;--- Read Granvik NEA model
         readcol, '/data/euclid/synthetic/nea-granvik2016.csv',synH,synN,syndN,delimiter=',',/Silent
         synH+=0.25
         synCum = total(synN,/Cumulative)-320

        ;--- 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/C-- MBA
      'MB': 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, range=[0,24])

        mbaModel = ssoModelSFD( selSFD.h, 0.51, 17.4, 0.3, reference={h:13,N:2400} ) ;-Gladman09: 0.5 - then 0.3

        selSFD.cumul.model = mbaModel.cumul.model
        low = where( selSFD.h le 14.5 )
        selSFD.cumul.model[low] = selSFD.cumul.obs[low]

;        mbaModel = ssoModelSFD( selSFD.h, 0.58, reference={h:15,N:25000} ) ;-me
;        cgPlot, mbaModel.h, mbaModel.cumul.model, /over, color='gray',linestyle=2

     end

      
      
     ;--IV.2.2/D-- Trojans (from Jewitt 2000 + Yoshida2005)
      'Trojan': begin
        sel = where( strcmp(disco.orbit.main,popName[kPop]) or strcmp(disco.orbit.sub,popName[kPop]),nbSel )
        ref=[popHmin[kPop],popHmax[kPop]]
        sfdM = ssoSynthSFD(disco[sel].h,     reference=ref, binSize=binH, range=[0,24])

;        jtaModel = ssoModelSFD( sfdM.h, 0.38, 14.0, 0.26, reference={h:9.5,N:20.} )  ;Yoshida - doesnt fit
;        jtaModel = ssoModelSFD( sfdM.h, 0.40, reference={h:9.5,N:13.5} )  ;-faint end Jewitt
        
        hRef = 12.5
        nRef = sfdM.histo.obs[ where( sfdM.h eq hRef ) ]
          nRef=nRef[0]
        print, 'JTA', hRef, nRef
        nRef = 310
        
        jtaModelL = ssoModelSFD( sfdM.h, 0.34, reference={h:hRef,N:nRef} ) ;-Jewitt-
        jtaModelM = ssoModelSFD( sfdM.h, 0.40, reference={h:hRef,N:nRef} ) ;-Jewitt
        jtaModelH = ssoModelSFD( sfdM.h, 0.46, reference={h:hRef,N:nRef} ) ;-Jewitt+

        cgPlot, jtaModelL.h, jtaModelL.cumul.model, /over, color='gray', lineStyle=2
        cgPlot, jtaModelM.h, jtaModelM.cumul.model, /over, color='gray'
        cgPlot, jtaModelH.h, jtaModelH.cumul.model, /over, color='gray', lineStyle=2

        sfdM.cumul.model = jtaModelM.cumul.model
        low = where( sfdM.h le 13 )
        sfdM.cumul.model[low] = sfdM.cumul.obs[low]

        sfdH = sfdM
        sfdH.cumul.model = jtaModelH.cumul.model
        low = where( sfdH.h le 13 )
        sfdH.cumul.model[low] = sfdH.cumul.obs[low]

        sfdL = sfdM
        sfdL.cumul.model = jtaModelL.cumul.model
        low = where( sfdL.h le 13 )
        sfdL.cumul.model[low] = sfdL.cumul.obs[low]

        selSFD=sfdM

        
;        ppJTA=where(selSFD.H eq 10.0 )
;        print, 'JTA/Jewitt ', jtaModel.cumul.model[ppJTA]

      end


     ;--IV.2.2/E-- Centaurs (from Bauer 2013 WISE)
      'Centaur': 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, range=[0,24])

        cenModel = ssoModelSFD( selSFD.h, 0.34, reference={h:10,N:30} )  ;-Bauer2013
        
        selSFD.cumul.model = cenModel.cumul.model
        low = where( selSFD.h le 9.5 )
        selSFD.cumul.model[low] = selSFD.cumul.obs[low]

;        cgPlot, cenModel.h, cenModel.cumul.model, /over, color='gray'
;        ppCEN=where(selSFD.H eq 10.0 )
;        print, 'CEN/Jewitt ', cenModel.cumul.model[ppCEN]

      end
      
      
     ;--IV.2.2/F-- KBOs (from JM Petit - Mail 2017/2/2)
      'KBO': begin

         sel = where( strcmp(disco.orbit.main,'SDO') or strcmp(disco.orbit.sub,'SDO'),nbSel )
         ref=[4.,6.]
         selSFD = ssoSynthSFD(disco[sel].h,     reference=ref, binSize=binH, range=[0,24])


        ;----------------------------------------------------------------------------------
        ;-Shankman AJ 2016  : SDO
        ;----------------------------------------------------------------------------------
         sdoModel = ssoModelSFD( selSFD.h, 0.9, 8.3, 0.5, 5.6, reference={h:8.,N:6500.} )
;         cgPlot, sdoModel.h, sdoModel.cumul.model, /over, color='gray'
         ppSDO = where(selSFD.H eq 8.0 )

        ;----------------------------------------------------------------------------------
        ;-Petit DPS 2016  : Main Belt
        ;----------------------------------------------------------------------------------
         mcbModel = ssoModelSFD( selSFD.h, 1.02, 7.0, 0.66, reference={h:7.,N:1800.} )
;         cgPlot, mcbModel.h, mcbModel.cumul.model, /over, color='gray'
         ppMCB=where(selSFD.H eq 7.0 )



        ;----------------------------------------------------------------------------------
        ;-Volk AJ 2016: Resonant Population
        ;----------------------------------------------------------------------------------
         resModel = ssoModelSFD( selSFD.h, 0.9, reference={h:8.66,N:7500.} )
         ppRES=where(selSFD.H eq 8.75 )

;         cgPlot, resModel.h, resModel.cumul.model, /over, color='gray'




         kboModel = mcbModel.cumul.model + $
                    sdoModel.cumul.model + $
                    resModel.cumul.model
;         cgPlot, mcbModel.h, kboModel, /over, color='gray', lineStyle=2

;         print, 'KBO/SDO ', sdoModel.cumul.model[ppSDO]
;         print, 'KBO/RES ', resModel.cumul.model[ppRES]
;         print, 'KBO/MCB ', mcbModel.cumul.model[ppMCB]

         selSFD.cumul.model = kboModel
         low = where( selSFD.h le 4 )
         selSFD.cumul.model[low] = selSFD.cumul.obs[low]

         
      end



     ;--IV.2.2/G-- Comet (from Snodgrass2011+)
      'Comet': begin
         sel = where( strcmp(disco.orbit.main,'NEA') or strcmp(disco.orbit.sub,'NEA'),nbSel )
         ref=[popHmin[kPop],popHmax[kPop]]
         selSFD = ssoSynthSFD(disco[sel].h,     reference=ref, binSize=binH, range=[0,24])

        ;-- 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/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, range=[0,24])
        saveSFD = selSFD

      end

    endcase




    
    if popMain[kPop] eq 1 then begin

       
     ;--IV.2.4-- OverPlot synthetic distribution
      cgPlot, /OverPlot, selSFD.h, selSFD.cumul.model, color=popC[kPop], linestyle=2
      cgPlot, /OverPlot, selSFD.h, selSFD.cumul.obs,   color=popC[kPop]



      
     ;--IV.2.5-- Compute H-V
      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]

     ;-IV.2.6-- Saturation and Limiting magnitudes
      hLimit = round( (vLimit-joint)/binH) * binH
      goodH = where( hLimit le 30 and finite(hLimit) and hLimit ge popHmax[kPop] )
      hLimit = hLimit[goodH]

      lLimit = round( (sLimit-joint)/binH) * binH
      goodL = where( lLimit le 30 and finite(lLimit) and lLimit ge 0.1 )
      lLimit = lLimit[goodL]

      hBox = dataBox( hLimit)
      lBox = dataBox( lLimit)

      
      ecSFD = selSFD.cumul.model * 0.

      


     ;---Rising from Saturation
      for kQ=0,4 do begin

        pL = where( selSFD.h eq lBox[kQ] )
        pL = pL[0]

        xQ = selSFD.h[pL]

        if kQ eq 0 then begin
          yRef = selSFD.cumul.obs[pL]
          xB = xQ
          ecSFD[pL]=0

;          print, kQ, 0, selSFD.h[pL], lBox[kQ], yRef
        endif else begin

          subL = where( selSFD.h gt xB and selSFD.h le xQ, nbSub)
          iFrac = interpol( lFrac[[kQ-1,kQ]], [xB,xQ], selSFD.h[subL] )

          for kSub=0, nbSub-1 do begin
            if popName[kPop] eq 'Centaur' and selSFD.h[subL[kSub]] ge 9.5 then begin ;and xB ge 9.5 then begin
               yNew = selSFD.cumul.model[subL[kSub]]
            endif else begin
              yNew = selSFD.cumul.obs[subL[kSub]]
            endelse
            dY = (yNew - yRef)*iFrac[kSub]
            yRef += dY
            ecSFD[subL[kSub]] = dY
;            print, kQ, kSub, selSFD.h[subL[kSub]], iFrac[kSub], dY
          endfor

          xB = xQ
        endelse     

      endfor


      
     ;----follows model SFD
      pH = where( selSFD.h gt lBox[4] and selSFD.h le hBox[0], nbH )
      for kH=0, nbH-1 do begin
        yNew = selSFD.cumul.model[pH[kH]]
        dY = (yNew - yRef)
        yRef += dY
        ecSFD[pH[kH]] = dY
      endfor
      

     ;-Drop from limiting magnitude
      for kQ=0,4 do begin

        pH = where( selSFD.h eq hBox[kQ] )
        pH = pH[0]
        xQ = selSFD.h[pH]

        if kQ eq 0 then begin
          yRef = selSFD.cumul.model[pH]
          xB = xQ
        endif else begin

          subH = where( selSFD.h gt xB and selSFD.h le xQ, nbSub)
          iFrac = interpol( hFrac[[kQ-1,kQ]], [xB,xQ], selSFD.h[subH] )

          for kSub=0, nbSub-1 do begin
            yNew = selSFD.cumul.model[subH[kSub]]
            dY = (yNew - yRef)*iFrac[kSub]
            yRef += dY
            ecSFD[subH[kSub]] = dY
          endfor
          xB = xQ

       endelse     


      endfor


     ;-Euclid observations in CUMULATIVE histo
      ecCum = total(ecSFD,/cum)

     ;-clean point for figure
      case popName[kPop] of
        'MB': begin
           pH=where(selSFD.h eq 10.7500)
           ecCum[pH]=1
        end
        else:
      endcase
      cgPlot, selSFD.h, ecCum, /OverPlot, lineStyle=3, color=popC[kPop]

      

     ;---- Limiting Quartlies and margin for discovery
      for kQ=0,4 do begin
        pH = where( selSFD.h eq hBox[kQ] )
        pH = pH[0]
        xQ = selSFD.h[pH]
        yQ = ecCum[pH]
        cgPlot, /OverPlot, xQ, yQ, psym='Open Circle', color=popC[kPop]

        if kQ eq 4 then begin
          yN = selSFD.cumul.obs[pH]
          cgPlot, /OverPlot, xQ, yN, psym='Filled Circle', color=popC[kPop]
          cgPlot, /OverPlot, [xQ,xQ], [yN,yQ], color=popC[kPop]
        endif
      endfor

;      print, popName[kPop], ' ', yQ, yQ-yN

     ;---- Limiting Quartlies for saturation
      for kQ=0,4 do begin
        pL = where( selSFD.h eq lBox[kQ] )
        pL = pL[0]
        xQ = selSFD.h[pL]
        yQ = ecCum[pL]
        cgPlot, /OverPlot, xQ, yQ, psym='Open square', color=popC[kPop]
      endfor


      ecStat[kPop,0] = max(ecCum)
      ecStat[kPop,1] = max(ecCum)-max(selSFD.cumul.obs)
      
    endif
    
   
    if popMain[kPop] eq 1 then kPrint++
  endfor

  cgText, xKey+xLen*0.5, yKey2[kPrint], 'Inner Solar System'


 ;--- Diameter range
  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, yTop2, 'km'
  cgText, /Normal, 0.05, yTop2, 'Diameter'

 
        kPrint=7
        cgText, xKey+xLen*0.5, yKeyLine2[kPrint], 'Populations', color='Dark Gray'
        kPrint--
        cgPlot, /OverPlot, xKey+[0,xLen], yKeyLine[[kPrint,kPrint]], thick=6, color='Dark Gray'
        cgText, xKey+xLen*1.2, yKeyLine2[kPrint], 'Known', color='Dark Gray'
        kPrint--
        cgPlot, /OverPlot, xKey+[0,xLen], yKeyLine[[kPrint,kPrint]], thick=6, lineStyle=2, color='Dark Gray'
        cgText, xKey+xLen*1.2, yKeyLine2[kPrint], 'Synthetic', color='Dark Gray'
        kPrint--
        cgPlot, /OverPlot, xKey+[0,xLen], yKeyLine[[kPrint,kPrint]], thick=6, lineStyle=3, color='Dark Gray'
        cgText, xKey+xLen*1.2, yKeyLine2[kPrint], 'by Euclid', color='Dark Gray'

        kPrint--
        kPrint--
        cgText, xKey+xLen*0.5, yKeyLine2[kPrint], 'Detections', color='Dark Gray'
        kPrint--
        cgPlot, /OverPlot, xKey+xLen/2., yKeyLine[[kPrint]], color='Dark Gray', psym='Open Circle'
        cgText, xKey+xLen*1.2, yKeyLine2[kPrint], 'Photon-limited', color='Dark Gray'
        kPrint--
        cgPlot, /OverPlot, xKey+xLen/2., yKeyLine[[kPrint]], color='Dark Gray', psym='Open square'
        cgText, xKey+xLen*1.2, yKeyLine2[kPrint], 'Saturated', color='Dark Gray'

  
  cgText, xDisco+xLen*1.2, yTxtDisco, 'Margin for'
  cgText, xDisco+xLen*1.2, yBotDisco, 'discovery'




  
  cgPS_close, /PNG, delete_PS=0



  forprint, popName, ecStat[*,0], ecStat[*,1], ecStat[*,0], ecStat[*,1], $
            format='(A-20,", ",I8,", ",I8,", ",E9.3,", ",E9.3)', $
            textout=dirStat+'euclid-stat.csv', /Silent, comment='Population, Nall, Ndisco, NallE, NdiscoE'

;  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'
;

end