function randomize, V, N, seed, normal=normal, uniform=uniform
  if not keyword_set(N)    then N=10000
  if not keyword_set(seed) then seed=1

  if not keyword_Set(normal) and not keyword_Set(uniform) then normal=1

  if keyword_Set(normal) then begin
    return, randomn(seed,N)*V.unc + V.val 
  endif else begin
    return, V.min + randomu(seed,N)*(V.max-V.min) 
  endelse
end



;goto, j2src


;----------------------------------------------------------------------------------------------------;
;----------------------------------------------------------------------------------------------------;
;----    I - TAG ----   General Definitions                                        ------------------;
;----------------------------------------------------------------------------------------------------;
;----------------------------------------------------------------------------------------------------;

 ;--I.1-- Directories  ------------------------------------------------------------------------------;
  dir='/data/dynamics/yarkovsky/'
  dirPPT=dir+'properties/'
  
 ;--I.2-- Files  ------------------------------------------------------------------------------------;
  diamalb = '/home/bcarry/Downloads/diamalbedos.csv'

  damit   = dirPPT+'damit-2018-10-15.csv'
  mainzer = dirPPT+'mainzer2011.txt'
  fDiams  = dirPPT+'diameters.csv'
  sDiams  = dirPPT+'diameters.src'
  fDiam   = dirPPT+'sheet-diameter.csv'
  fTaxo   = dirPPT+'sheet-taxonomy.csv'
  fPer    = dirPPT+'sheet-period.csv'
  fYarko  = dir+'yarko.all.gdr2'
  tmp = '/tmp/tmp.ppt'

 ;--I.3-- Default Values  ---------------------------------------------------------------------------;
  cutSNR = 3. ; 2.2
  nbMC = 5e4
  nbMC = 1e5

  defPeriod     = {val:0.0,   unc:0.001}
  defAlbedo     = {val:0.125, unc:0.10}
  defInertia    = {val:200, unc:40.}  ;Delbo 2007
  defInertiaBin = {val:480, unc:70.}  ;Delbo 2007

  defObliPro  = {min:00., max: 90.}
  defObliRetro= {min:90., max:180.}
;  defObliquity  = {min:170., max:180.}


  seedPeriod = 1
  seedYarko  = 2
  seedAlbedo = 3
  seedDiam   = 4
  seedInertia= 5
  seedObli   = 6


;----------------------------------------------------------------------------------------------------;
;----------------------------------------------------------------------------------------------------;
;----   II - TAG ----   Read & Load Data                                           ------------------;
;----------------------------------------------------------------------------------------------------;
;----------------------------------------------------------------------------------------------------;

 ;--II.1-- Yarkovsky Detections  --------------------------------------------------------------------;
  readcol, fYarko, num, A2,sigA2, SNR, dadt, sigdadt, rms, seld, discd, orbit, S, $;diam, $
           format='(L,F,F,F,F,F,F,F,F,A,F)', /Silent

 ;--II.2-- Select Detection with SNR>3 and S Quality Factor  ----------------------------------------;
;  sel=where( SNR ge cutSNR, nbSSO)
  sel=where( SNR ge cutSNR and $
             (num ne 1864   and $
              num ne 1980   and $
              num ne 3040   and $
              num ne 3103   and $
              num ne 3199   and $
              num ne 4953   and $
              num ne 5407   and $
              num ne 5587   and $
              num ne 6455   and $
              num ne 7467   and $
              num ne 11066  and $
              num ne 16834  and $
              num ne 21088  and $
              num ne 54789  and $
              num ne 66146  and $
              num ne 68216  and $
              num ne 85804  and $
              num ne 86829  and $
              num ne 137084 and $
              num ne 137805 and $
              num ne 143649 and $
              num ne 163693 and $
              num ne 164121 and $
              num ne 162421 ), nbSSO)
  print, nbSSO, n_elements(num)


 ;--II.3-- SSO Structure  ---------------------------------------------------------------------------;
  eVal={val:0., unc:0., min:0., max:0., src:'', set:0}
  eSSO={num:0L, name:'', $
        taxo:{scheme:'', class:'', src:''}, $
        A2:eVal, dadt:eVal, $
        SNR:0., $
        diam:eVal, $
        albedo:eVal, $
        period:eVal, $
        inertia:eVal, $
        spin:  {long:eVal, lat:eVal}, $
        orbit: {long:eVal, lat:eVal}, $
        obliquity:eVal, $
        beaming:eVal, $
        emissivity:eVal, $
        mass:eVal, $
        density:eVal, $
        dSun:0., H:0., G:0., $
        meanMotion:0. }
  sso=replicate(eSSO,nbSSO)

;  sso.period.unc=0.001

 ;--II.4-- Mainzer Albedo-Class Relations  ----------------------------------------------------------;
  spawn, 'grep DeMeo '+mainzer+' > '+tmp
  readcol, tmp, mSc, mClass, trash, mMed, trash, mStd, format='(A,A,I,F,F,F)', /Silent
    mClass=strTrim(mClass,2)


;----------------------------------------------------------------------------------------------------;
;----------------------------------------------------------------------------------------------------;
;----  III - TAG ----   Fill Properties                                            ------------------;
;----------------------------------------------------------------------------------------------------;
;----------------------------------------------------------------------------------------------------;
   print, '#', 'Name', 'Tax', 'SNR', $
           'Diam','Alb', 'TI', $
           'Per', 'Lon','Lat', $
           format='(A6,",",A-12,",",A-3,",",A5,",",'+$
                   'A5,",",A5,",",A5,",",'+$
                   'A8,",",A3,",",A3)'
  for k=0, nbSSO-1 do begin
    i=sel[k]

   ;--III.1-- Yarkovsky Information  ----------------------------------------------------------------;
    sso[k].num    = num[i]
    sso[k].name   = designation(sso[k].num)
    sso[k].A2.val = A2[i]    * 1e-10
    sso[k].A2.unc = sigA2[i] * 1e-10
    sso[k].dadt.val= dadt[i]
    sso[k].dadt.unc= sigdadt[i]
    sso[k].dadt.src= 'This work'
    sso[k].SNR    = SNR[i]

;    print, sso[k].num , sso[k].name, sso[k].snr, format='(I6,",",A-12,",",F5.1)'


   ;--III.2-- Diameters from General File  ------------------------------------------------------------;
    readcol, sDiams, srcId, srcBib, delimiter=',', /Silent, format='(I,A)'
    spawn, 'rm -f '+tmp
    spawn, 'grep "'+sso[k].name+'" '+fDiams+' > '+tmp
    spawn, 'wc -l '+tmp, res
;    if strCmp(sso[k].name,'Magellan') then stop
    split=strSplit(res,' ', /Extract)
    case split[0] of
       '0': begin
         if sso[k].num eq 87684 then begin
           sso[k].albedo.val     = 0.15
           sso[k].albedo.unc     = 0.10
           sso[k].diam.val       = (1329/sqrt(0.15))*10^(-0.2*16.1)
           sso[k].diam.unc       = 0.15 *  sso[k].diam.val

         endif
       end

       '1': begin
          readcol, tmp, delimiter=',', format='(I,I,A,F,F,F,F,F,F,F,F,A,A)', /Silent, $
                   idsso, sso_num,sso_name,diameter,err_diameter,albedo,err_albedo, $
                   beaming,err_beaming,emissivity,err_emissivity,iddataset,method
          sso[k].diam.val       = diameter
          sso[k].diam.unc       = err_diameter
          sso[k].albedo.val     = albedo
          sso[k].albedo.unc     = err_albedo
          sso[k].beaming.val    = beaming
          sso[k].beaming.unc    = err_beaming
          sso[k].emissivity.val = emissivity
          sso[k].emissivity.unc = err_emissivity

          pSrc=where(iddataset eq srcId)
          sso[k].diam.src       = strTrim(srcBib[pSrc[0]],2)
          sso[k].albedo.src     = strTrim(srcBib[pSrc[0]],2)
       end

       else: begin
          readcol, tmp, delimiter=',', format='(I,I,A,F,F,F,F,F,F,F,F,A,A)', /Silent, $
                   idsso, sso_num,sso_name,diameter,err_diameter,albedo,err_albedo, $
                   beaming,err_beaming,emissivity,err_emissivity,iddataset,method
          dd = meanWithUnc( diameter, err_diameter )
          sso[k].diam.val       = dd[0]
          sso[k].diam.unc       = dd[1]
          dd = meanWithUnc( albedo, err_albedo )
          sso[k].albedo.val     = dd[0]
          sso[k].albedo.unc     = dd[1]
          dd = meanWithUnc( beaming, err_beaming )
          sso[k].beaming.val    = dd[0]
          sso[k].beaming.unc    = dd[1]
          dd = meanWithUnc( emissivity, err_emissivity )
          sso[k].emissivity.val = dd[0]
          sso[k].emissivity.unc = dd[1]

          pSrc= where(iddataset[0] eq srcId)
          dd  = strTrim(srcBib[pSrc[0]],2)
;          dd=iddataset[0]
          for j=1,round(float(split[0]))-1 do begin
             pSrc=where(iddataset[j] eq srcId)
             dd+=','+strTrim(srcBib[pSrc[0]],2)
          endfor
          sso[k].diam.src       = dd
          sso[k].albedo.src     = dd
       end

    endcase

   ;--III.3-- Diameters from Manual File  -------------------------------------------------------------;
    spawn, 'rm -f '+tmp
    spawn, 'grep "'+sso[k].name+'" '+fDiam+' > '+tmp
    spawn, 'wc -l '+tmp, res
    split=strSplit(res,' ', /Extract)
    if not strCmp(split[0],'0') then begin
      line=' '
      openR, 1, tmp
      readf, 1, line
      close, 1
      split=strSplit(line, ',' ,/Extract)
      if not strCmp(split[2] ,' ') then sso[k].diam.val       = float(split[2])
      if not strCmp(split[3] ,' ') then sso[k].diam.unc       = float(split[3])
      if not strCmp(split[5] ,' ') then sso[k].diam.src       = strTrim(split[5],2)
      if not strCmp(split[6] ,' ') then sso[k].albedo.val     = float(split[6])
      if not strCmp(split[7] ,' ') then sso[k].albedo.unc     = float(split[7])
      if not strCmp(split[12],' ') then sso[k].albedo.src     = strTrim(split[12],2)
      if not strCmp(split[10],' ') then begin
         sso[k].inertia.val    = float(split[10])
         sso[k].inertia.src    = strTrim(split[12],2)
      endif
      if not strCmp(split[8] ,' ') then sso[k].inertia.min    = float(split[8])
      if not strCmp(split[9] ,' ') then sso[k].inertia.max    = float(split[9])
    endif


   ;--III.4-- Absolute Magnitude and Orbit from ASTORB  -----------------------------------------------;
    elt=astElements(sso[k].num)
    sso[k].H = elt.h
    sso[k].G = 0.15 ; elt.g
    sso[k].meanMotion =  2*!PI*(elt.a^(-1.5)) / (365.2524*86400.)





    sso[k].dSun = float(elt.a) * sqrt( 1-elt.e*elt.e )
    sso[k].orbit.long.val = elt.O - 90.
    sso[k].orbit.lat.val  = 90. - elt.i



   ;--III.5-- Spin and Period from Manual File  -------------------------------------------------------;
    spawn, 'rm -f '+tmp
    spawn, 'grep "'+sso[k].name+'" '+fPer+' > '+tmp
    spawn, 'wc -l '+tmp, res
    split=strSplit(res,' ', /Extract)
    if not strCmp(split[0],'0') then begin
      line=' '
      openR, 1, tmp
      readf, 1, line
      close, 1
      split=strSplit(line, ',' ,/Extract)
      if not strCmp(split[2] ,' ') then sso[k].period.val       = float(split[2])
      if not strCmp(split[3] ,' ') then sso[k].period.unc       = float(split[3])
      if not strCmp(split[4] ,' ') then sso[k].spin.long.val    = float(split[4])
      if not strCmp(split[5] ,' ') then sso[k].spin.lat.val     = float(split[5])
      if not strCmp(split[6] ,' ') then begin
        sso[k].spin.long.unc = float(split[6])
        sso[k].spin.lat.unc  = float(split[6])
      endif
      if not strCmp(split[5] ,' ') then sso[k].spin.lat.val     = float(split[5])
      sso[k].period.src=strTrim(split[9],2)

      if sso[k].spin.long.val ne 0 then begin
        gcirc, 2, sso[k].orbit.long.val, sso[k].orbit.lat.val, $
                  sso[k].spin.long.val,  sso[k].spin.lat.val, obli
        sso[k].obliquity.val = obli/3600.

;        if strCmp(sso[k].name,'Bennu') or strCmp(sso[k].name,'Phaethon') then begin
;          print, sso[k].num
;          print, sso[k].orbit.long.val, sso[k].orbit.lat.val, $
;                 sso[k].spin.long.val,  sso[k].spin.lat.val
;          print, obli
;          print, sso[k].obliquity.val
;;          stop
;        endif
        if  sso[k].spin.long.unc ne 0 then sso[k].obliquity.unc = sso[k].spin.long.unc $
                                      else sso[k].obliquity.unc = 15.


      endif

    endif


   ;--III.6-- Taxonomic classes from PDS  -------------------------------------------------------------;
    taxInfo = pds_Taxonomy_getCatalog(sso[k].num)
    if size(taxInfo,/Type) eq 8 then begin
      sso[k].taxo.scheme = taxInfo.scheme
      sso[k].taxo.class  = taxInfo.taxo
      case taxInfo.scheme of
        'Bus':   sso[k].taxo.src='2002-Icarus-158-BusII'
        'DeMeo': sso[k].taxo.src='2009-Icarus-202-DeMeo'
         else:
      endcase
    endif 
    spawn, 'rm -f '+tmp
    spawn, 'grep "'+sso[k].name+'" '+fTaxo+' > '+tmp
    spawn, 'wc -l '+tmp, res
    split=strSplit(res,' ', /Extract)
    if not strCmp(split[0],'0') then begin
      line=' '
      openR, 1, tmp
      readf, 1, line
      close, 1
      split=strSplit(line, ',' ,/Extract)
      if not strCmp(split[4] ,' ') then begin
        sso[k].taxo.class      = strTrim(split[4])
        sso[k].taxo.scheme     = strTrim(split[3])
        sso[k].taxo.src        = strTrim(split[5])
      endif
    endif

    print, sso[k].num , sso[k].name, sso[k].taxo.class, sso[k].snr, $
           sso[k].diam.val, sso[k].albedo.val, sso[k].inertia.val, $
           sso[k].period.val, sso[k].spin.long.val, sso[k].spin.lat.val, $
           format='(I6,",",A-12,",",A-3,",",F5.1,",",'+$
                   'F5.2,",",F5.3,",",I5,",",'+$
                   'F8.4,",",I3,",",I3)'

;    print, sso[k].num , sso[k].name, sso[k].snr, $
;           sso[k].diam.val, sso[k].albedo.val, sso[k].diam.src, $
;           format='(I6,",",A-12,",",F5.1,",",'+$
;                   'F5.2,",",F5.1,",",A-25)'
;
  endfor


 ;--III.7-- Statistics on Missing Parameters  -------------------------------------------------------;
  noD  = where( sso.diam.val eq 0, nbNoDiam )
  noS  = where( sso.spin.long.val eq 0, nbNoSpin )
  noA  = where( sso.albedo.val eq 0, nbNoAlbedo )
  noTI = where( sso.inertia.val eq 0, nbNoInertia )
  noTax= where( strCmp(sso.taxo.class,''), nbNoTaxo )

  print, '# Diameter: ', nbSSO-nbNoDiam   , 100-100.*nbNoDiam   /nbSSO, format='(A-12,1x,I2,2x,F5.1," %")' 
  print, '# Spin    : ', nbSSO-nbNoSpin   , 100-100.*nbNoSpin   /nbSSO, format='(A-12,1x,I2,2x,F5.1," %")' 
  print, '# Albedo  : ', nbSSO-nbNoAlbedo , 100-100.*nbNoAlbedo /nbSSO, format='(A-12,1x,I2,2x,F5.1," %")' 
  print, '# Inertia : ', nbSSO-nbNoInertia, 100-100.*nbNoInertia/nbSSO, format='(A-12,1x,I2,2x,F5.1," %")' 
  print, '# Taxonomy: ', nbSSO-nbNoTaxo   , 100-100.*nbNoTaxo   /nbSSO, format='(A-12,1x,I2,2x,F5.1," %")' 


 ;--III.8-- Fill in on Missing Parameters  -----------------------------------------------------------;
 ;--III.8.1-- Period Uncertainty
  noUncPer = where( sso.period.unc eq 0, nbNoUncPer )
  sso[noUncPer].period.unc = defPeriod.unc

 ;--III.8.2-- Albedo from Taxonomy
  for kA=0, nbNoAlbedo-1 do begin
    pA = where( strCmp(sso[noA[kA]].taxo.class, mClass) )
    sso[noA[kA]].albedo.val = mMed[pA[0]]
    sso[noA[kA]].albedo.unc = mStd[pA[0]]
  endfor

 ;--III.8.3-- Thermal Inertia
  sso[noTI].inertia.val=defInertia.val
  sso[noTI].inertia.unc=defInertia.unc

;  noUncTI = where( sso.inertia.unc eq 0, nbNoUncTI )
;  sso[noUncTI].inertia.unc = (sso[noUncTI].inertia.max-sso[noUncTI].inertia.min)/2.


 ;---- Fill with
 ; 1) Albedo from taxonomy
 ; 2) Diameter from albedo/H
 ; 3) Inertia from single/binary
 ; 4) Spin 0-180
 ; 5) Uncertainties



;----------------------------------------------------------------------------------------------------;
;----------------------------------------------------------------------------------------------------;
;----   IV - TAG ----   Fill Properties                                            ------------------;
;----------------------------------------------------------------------------------------------------;
;----------------------------------------------------------------------------------------------------;
   print, '#', 'Name', 'Tax', 'SNR', 'dadt','unc',$
           'Diam','Alb', 'TI', $
           'Per', 'Obli','Dens','dDens', $
           format='(A6,",",A12,",",A-3,", ",'+$
                   'A5,",",A5,",",A5,",   ",'+$
                   'A5,",",A5,",",A5,",   ",'+$
                   'A8,",",A3,",   ",'+$
                   'A6,",",A6)'
;  for k=8, nbSSO-1 do begin  ;phaethon
;  for k=16, nbSSO-1 do begin  ;bennu
  for k=0, nbSSO-1 do begin  ;bennu


;     print, sso[k].name
;    print , 'TEST LICK'
;    sso[k].obliquity.val = 10.
;    sso[k].obliquity.unc = 10.
;    sso[k].period.unc = 0.000002

     
;    print, 'dadt      ', sso[k].dadt   
;    print, 'period    ', sso[k].period 
;    print, 'albedo    ', sso[k].albedo 
;    print, 'diam      ', sso[k].diam   
;    print, 'obliquity ', sso[k].obliquity
;    print, 'inertia   ', sso[k].inertia     





   ;--IV.1-- Randomize Normal Variables  ------------------------------------------------------------;
    ranYarko   = randomize(sso[k].dadt     , nbMC, seedYarko,  /Normal)
    ranPeriod  = randomize(sso[k].period   , nbMC, seedPeriod, /Normal)
    ranAlbedo  = randomize(sso[k].albedo   , nbMC, seedAlbedo, /Normal)
    ranDiam    = randomize(sso[k].diam     , nbMC, seedDiam,   /Normal)

   ;--IV.2-- Randomize Uniform/Random Variables  ----------------------------------------------------;
    if sso[k].obliquity.val eq 0 then begin

      if sso[k].dadt.val gt 0 then begin
        ranObli = randomize(defObliPro    , nbMC, seedObli, /Uniform)
        sso[k].obliquity.val = (defObliPro.max+defObliPro.min)/2.
        sso[k].obliquity.unc = (defObliPro.max-defObliPro.min)/2.
        sso[k].obliquity.src = 'Assumed'
      endif else begin
        ranObli = randomize(defObliRetro    , nbMC, seedObli, /Uniform)
        sso[k].obliquity.val = (defObliRetro.max+defObliRetro.min)/2.
        sso[k].obliquity.unc = (defObliRetro.max-defObliRetro.min)/2.
        sso[k].obliquity.src = 'Assumed'
      endelse
      
      if sso[k].num eq 12711 then begin 
        toto  = {min:170., max: 180.}
        ranObli = randomize(toto    , nbMC, seedObli, /Uniform)
        sso[k].obliquity.val = (toto.max+toto.min)/2.
        sso[k].obliquity.unc = (toto.max-toto.min)/2.
        sso[k].obliquity.src = 'Assumed'
      endif

    endif else ranObli = randomize(sso[k].obliquity, nbMC, seedObli, /Normal)


   ;--IV.3-- Randomize Skwed Variables  -------------------------------------------------------------;
    if sso[k].inertia.val eq 0 then ranInertia = randomize(sso[k].inertia, nbMC, seedInertia, /Uniform) $
                               else ranInertia = randomize(sso[k].inertia, nbMC, seedInertia, /Normal)

;
;; validation with phaethon
;    vMass = yarko_to_mass( period=sso[k].period.val, $
;                           diameter=sso[k].Diam.val, $
;                           albedo=sso[k].albedo.val,$
;                           inertia=sso[k].inertia.val,$
;                           yarko=sso[k].dadt.val, $
;                           dSun=sso[k].dSun, $
;                           obliquity=sso[k].Obliquity.val, $
;                           meanMotion=sso[k].meanMotion )
;    vVol = !PI*( (sso[k].diam.val*1e3)^3)/6.
;    vDens = vMass/vVol
;
;;    print, sso[k].diam.val, vVol
;;    print, vMAss, vDens
;;stop
;
;
;
;;    window, 2
;    cgHistoPlot, ranPeriod, xtitle='Period (h)', title=sso[k].name
;    cgHistoPlot, ranDiam, xtitle='Diameter (km)', title=sso[k].name
;    cgHistoPlot, ranAlbedo, xtitle='Albedo', title=sso[k].name
;    cgHistoPlot, ranInertia, xtitle='Thermal inertia (SI)', title=sso[k].name
;    cgHistoPlot, ranYarko, xtitle='Yarkovsky da/dt', title=sso[k].name
;    cgHistoPlot, ranObli, xtitle='Obliquity (deg)', title=sso[k].name






   ;--IV.43-- Compute Mass and Density by MC  --------------------------------------------------------;
    mass=fltarr(nbMC)
    dens=fltarr(nbMC)
    for kMC=0L, nbMC-1 do begin
      mass[kMC] = yarko_to_mass( period=ranPeriod[kMC], $
                                 diameter=ranDiam[kMC], $
                                 albedo=ranAlbedo[kMC],$
                                 inertia=ranInertia[kMC],$
                                 yarko=ranYarko[kMC], $
                                 dSun=sso[k].dSun, $
                                 obliquity=ranObli[kMC], $
                                 meanMotion=sso[k].meanMotion )
       vol = !PI*( (ranDiam[kMC]*1e3)^3)/6.
       dens[kMC] = mass[kMC]/vol
    endfor
    
    sso[k].mass.val = mean(mass)
    sso[k].mass.unc = stddev(mass)
    sso[k].density.val = mean(dens)
    sso[k].density.unc = stddev(dens)

   
    print, sso[k].num , sso[k].name,  sso[k].taxo.class, $
           sso[k].snr, sso[k].dadt.val,sso[k].dadt.unc, $
           sso[k].diam.val, sso[k].albedo.val, sso[k].inertia.val, $
           sso[k].period.val, sso[k].obliquity.val, $
           sso[k].density.val, sso[k].density.unc, $
           format='(I6,",",A-12,",",A-3,", ",'+$
                   'F5.1,",",F5.1,",",F5.1,",   ",'+$
                   'F5.2,",",F5.3,",",I5,",   ",'+$
                   'F8.4,",",I3,",   ",'+$
                   'I6,",",I6)'

;    cgHistoPlot, dens, xtitle='Density', title=sso[k].name
  endfor

  forprint, sso.num , sso.name, sso.taxo.class, $
            sso.diam.val, sso.diam.unc, $
            sso.albedo.val, sso.albedo.unc, $
            sso.inertia.val, sso.inertia.unc, $
            sso.period.val, sso.period.unc, $
            sso.obliquity.val, sso.obliquity.unc, $
            sso.density.val, sso.density.unc, $
            format='(I6," & ",A-12," & ",A-3," & ",'+$
                    'F5.1," & ",F5.1," & ",'+$  ; diam
                    'F5.3," & ",F5.3," & ",'+$  ; albedo
                    'I3," & ",I3," & ",'+$      ; inertia
                    'F8.4," & ",F6.4," & ",'+$  ; period
                    'I3," & ",I3," & ",'+$      ; obliquity
                    'I4," & ",I4,"\\")', $      ; density
            textout=dir+'density_simple.tab', /Silent, /NoComment



  assI = where(strCmp(sso[*].inertia.src,''))
  sso[assI].inertia.src='Assumed'
  assA = where(strCmp(sso[*].albedo.src,''))
  sso[assA].albedo.src='Assumed'
j2src:

  allSrc = [sso[*].taxo.src, sso[*].diam.src, sso[*].albedo.src, sso[*].inertia.src, sso[*].period.src]
  indiv=allSrc[0]
  for i=1,n_elements(allSrc)-1 do indiv+=','+allSrc[i]


  ;----- print indivudual biblio references
;  print, indiv

  split=strSplit(indiv,',',/Extract)
  uSplit = uniq(split,sort(split))

  uniqSrc = split[uSplit]


  nbSSO=10
  taxoSrc    = strarr(nbSSO)
  inertiaSrc = strarr(nbSSO)
  periodSrc  = strarr(nbSSO)
  diamSrc    = strarr(nbSSO)
  albSrc     = strarr(nbSSO)
  alpha      =string(bindgen(1,26)+(byte('a'))[0])
  for k=0, nbSSO-1 do begin
    pp=where(strCmp(sso[k].taxo.src,uniqSrc))
    taxoSrc[k]=strTrim(string(pp[0]+1,format='(I)'),2)
;    taxoSrc[k]=alpha[pp[0]]

    pp=where(strCmp(sso[k].inertia.src,uniqSrc))
    inertiaSrc[k]=strTrim(string(pp[0]+1,format='(I)'),2)
;    inertiaSrc[k]=alpha[pp[0]]

    pp=where(strCmp(sso[k].period.src,uniqSrc))
    periodSrc[k]=strTrim(string(pp[0]+1,format='(I)'),2)
;    periodSrc[k]=alpha[pp[0]]

    split=strSplit(sso[k].diam.src,',',/Extract,count=nbRef)
    pp=where(strCmp(split[0],uniqSrc))
    diamSrc[k]=strTrim(string(pp[0]+1,format='(I)'),2)
;    diamSrc[k]=alpha[pp[0]]
    for j=1, nbRef-1 do begin
      pp=where(strCmp(split[j],uniqSrc))
      diamSrc[k]+=','+strTrim(string(pp[0]+1,format='(I)'),2)
;      diamSrc[k]+=','+alpha[pp[0]]
    endfor

    split=strSplit(sso[k].albedo.src,',',/Extract,count=nbRef)
    pp=where(strCmp(split[0],uniqSrc))
    albSrc[k]=strTrim(string(pp[0]+1,format='(I)'),2)
;    albSrc[k]=alpha[pp[0]]
    for j=1, nbRef-1 do begin
      pp=where(strCmp(split[j],uniqSrc))
      albSrc[k]+=','+strTrim(string(pp[0]+1,format='(I)'),2)
;      albSrc[k]+=','+alpha[pp[0]]
    endfor

  endfor

  forprint, sso.num, sso.name, $
            sso.taxo.class, taxoSrc, $
            sso.diam.val, sso.diam.unc, diamSrc, $
            sso.albedo.val, sso.albedo.unc, albSrc, $
            sso.inertia.val, sso.inertia.unc, inertiaSrc, $
            sso.period.val, sso.period.unc, periodSrc, $
            sso.obliquity.val, sso.obliquity.unc, periodSrc, $
            sso.density.val, sso.density.unc, $
            format='(I6," & ",A-12," & ",'+$
                    'A-3,"$^{",A2,"}$ & ",'+$
                    'F5.1," & ",F5.1,"$^{",A12,"}$ & ",'+$  ; diam
                    'F5.3," & ",F5.3,"$^{",A12,"}$ & ",'+$  ; albedo
                    'I3," & ",I3,"$^{",A2,"}$ & ",'+$      ; inertia
                    'F8.4," & ",F6.4,"$^{",A2,"}$ & ",'+$  ; period
                    'I3," & ",I3,"$^{",A2,"}$ & ",'+$      ; obliquity
                    'I4," & ",I4,"\\")', $      ; density
            textout=dir+'density.tab', /Silent, /NoComment

  forprint, indgen(n_elements(uniqsrc))+1,uniqSrc, $
            format='("$^{",I2,"}$:~\citet{",A-22,"},")', $
            textout=dir+'density.src', /Silent, /NoComment


end