FUNCTION DUSTEM_CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par=par, cutof=cutof, $
                       norm=norm, fname=fname, C_WD=C_WD
  
;
; Doc
;
 if N_PARAMS() LT 1 then begin 
   print,'----------------------------------------------------------------------------------------------------'
   print,'FUNCTION DUSTEM_CREATE_VDIST, ar, rho, fv=fv, mfrac=mfrac, ag=ag, ns=ns, slaw=slaw, par=par, cutof=cutof, $'
   print,'                       norm=norm, fname=fname, C_WD=C_WD'
   print,'----------------------------------------------------------------------------------------------------'
   print,''
   print,'generates a dust size distribution in volume normalized to 1 or norm'
   print,'uses power law or log-normal component'
   print,''
   print,' AR   (I): array(2) size range in cm '
   print,' RHO  (I): density (g/cm3) of grain type (array if composite)'
   print,' FV   (I): volume fractions if composite'
   print,' MFRAC(I): mass fraction if composite'
   print,' AG   (O): size grid in cm (output)'
   print,' NS   (I): nr of sizes [Default = 10]'
   print,' SLAW (I): array of size dist. law of grains ''DUSTEM_PLAW_VDIST'' or ''DUSTEM_LOGN_VDIST'''
   print,'           [Default = ''DUSTEM_PLAW_VDIST'' with index -3.5]'
   print,' PAR  (I): array(nlaw,npar) parameters for the size dist law: '
   print,'           [index,integral,beta,At] for PLAW and [center,width,integral] for LOGN'
   print,' CUTOF(I): parameters for large size cut-off in microns. Default is [0.0107,0428].'
   print,' NORM (I): normalization for the size distribution'
   print,' FNAME(I): file name to write size dist. in'
   print,' C_WD (I): keyword to generate a C dust size dist WD01 style' 
   print,'           ( PAH, VSG as log-normal+power law a^(-2.54) ) ' 
   print,''
   print,' Example 1: a power-law' 
   print,' par=dblarr(1,4) & par(0,0)=-3.21 & par(0,1)=1. & par(0,2)=0.3 & par(0,3)=0.164e-4'
   print,' sd = DUSTEM_CREATE_VDIST( [3.e-8,2.5e-5], 3.3, slaw=''dustem_plaw_vdist'', par=par )'
   print,''
   print,' Example 2: a log-normal'
   print,' par=dblarr(1,3) & par(0,0)=4.d-8 & par(0,1)=0.2 & par(0,2)=0.3 '
   print,' sd = DUSTEM_CREATE_VDIST( [3.e-8,2.5e-5], 2.25, slaw=''dustem_logn_vdist'', par=par )'
   print,''
   print,' Created March 2009, L. Verstraete, IAS'
   print,''
   print,'----------------------------------------------------------------------------------------------------'   
   RETURN, 0
 endif 
 
; inits
 if n_elements( AR ) EQ 0 then begin 
     print,'(F) DUSTEM_CREATE_VDIST: you must define a size range'
     RETURN,0
 endif 
 if n_elements( RHO ) EQ 0 then begin 
     print,'(F) DUSTEM_CREATE_VDIST: you must define a grain mass density'
     RETURN,0
 endif
 if n_elements( NS ) EQ 0 then ns=10
 if n_elements(slaw) EQ 0 then begin 
   slaw=['DUSTEM_PLAW_VDIST']
   par = dblarr(1,4)
   par(0) = -3.5
   par(1) = 1.d0
 endif else begin
   slaw=[ strupcase(strtrim(slaw,2)) ]
   for i=0,n_elements( slaw )-1 do begin 
      if (SLAW(i) NE 'DUSTEM_PLAW_VDIST') AND (SLAW(i) NE 'DUSTEM_LOGN_VDIST') then begin
         print,'(F) DUSTEM_CREATE_VDIST: undefined law ',slaw(i) 
         print,'                  only ''DUSTEM_PLAW_VDIST'' and ''DUSTEM_LOGN_VDIST'' '
         return,0
      endif
   endfor
 endelse
 nmat = n_elements( RHO ) ; nr of bulk material in composite
 if n_elements( FV ) EQ 0 then begin 
     if NMAT EQ 1 then fv = [1.] else $
          fv = (fltarr(nmat)+1.) / nmat
     mfrac = fv 
 endif
 if n_elements( CUTOF ) EQ 0 then cutof = [ 0.0107, 0.428] * 1.d-4
 if n_elements( NORM ) EQ 0 then norm = 1.d0

; WD01 style for carbon adapted to match cirrus emission
 if keyword_set( C_WD ) then begin 
     slaw = [ 'DUSTEM_LOGN_VDIST', 'DUSTEM_LOGN_VDIST', 'DUSTEM_PLAW_VDIST' ]
     par = dblarr(3,4)
     if TOTAL(CUTOF) NE 0 then cutof = [ 0.0107, 0.170 ] * 1.d-4

; PAH log-normal 
     par(0,0) = 4.d-8
     par(0,1) = 0.2 
     par(0,2) = 0.30

; VSG log-normal 
     par(1,0) = 7.d-8
     par(1,1) = 0.4 
     par(1,2) = 0.08

; BG power law
     par(2,0) = -2.54
     par(2,1) = 2.e-5
     par(2,2) = 0.0 ;-0.165
     par(2,3) = 0.0107e-4
 endif
 nlaw = n_elements( SLAW ) 

; define size grid (log)
 lar = alog(ar)
 da = (lar(1)-lar(0)) / (ns-1)
 a = lar(0) + dindgen(ns)*da 
 a(ns-1) = lar(1)                ; roundup error
 ag = exp(a)  

; volume distribution
 if n_elements( PAR ) EQ 0 then begin 
     print,'(F) DUSTEM_CREATE_VDIST: PAR array not defined' 
     RETURN, 0
 endif
 vdist = 0.d0
 for i = 0, nlaw-1 do begin
     pp = REFORM( par(i,*), n_elements(par(i,*)) )
     vdist = vdist + ag^(4.d0) * CALL_FUNCTION( slaw(i),ag,pp )
 endfor 
 if TOTAL(CUTOF) NE 0 then begin 
    print,'(W) DUSTEM_CREATE_VDIST : cut-off applied, size and scale are ',cutof
    vdist = vdist * DUSTEM_CUT_OFF(ag, cutof)
 endif
 
; normalize
 fac = TOTAL( vdist(1:ns-1)+vdist(0:ns-2) ) * 0.5 * da 
 vdist = norm * vdist / fac
 print,'(W) DUSTEM_CREATE_VDIST: normalization factor ',fac / norm
 yr = [1.d-4,1] * max(vdist)
 plot_oo,1.e7*ag,vdist,xtit='a (nm)',yr=yr,ytit='Normalized a!u4!ndn/da', ps=-1

; effective density
 rho_eff = TOTAL( fv*rho )
 print,'(W) DUSTEM_CREATE_VDIST : effective mass density ',rho_eff,' g/cm3'

; write in fname
 if n_elements(fname) NE 0 then begin 
   OPENW, iu, fname, /get_lun
   printf, iu, '# Size distribution of grain species'
   printf, iu, '#'
   printf, iu, '# Nbr of bulk materials'
   printf, iu, '# Bulk densities in g/cm3 '
   printf, iu, '# Mass fractions for each bulk'
   printf, iu, '# Nbr of size bins'
   printf, iu, '# [ a (cm), dloga, a^4*dn/da, rho_eff, fv ] '
   printf, iu, '# fv: volume fraction of bulks, rho_eff: volume mean density'
   printf, iu, '#'
   printf, iu, nmat, format='(i2)'
   printf, iu, rho, format='(10(e11.4,1x))'
   printf, iu, mfrac, format='(10(e11.4,1x))'
   printf, iu, ns, format='(i2,1x,e11.4)'
   for i=0,ns-1 do begin 
    printf, iu, ag(i), da, vdist(i,*), rho_eff, fv, format='(20(e11.4,1x))'
   endfor
   FREE_LUN, iu
   print,'(W) DUSTEM_CREATE_VDIST: mass distribution written in ', strtrim(fname,2)
 endif

RETURN, vdist
END