FUNCTION PLAW, x, par 
; generates a volume normalized power law x^par(0)
; returns distribution in nr of grains : dn/da
; x : grain size 
; par(0) : power law index
; par(1) : VOLUME normalization 
; par(2) : curvature parameter beta 
; par(3) : large size threshold At

 np = n_elements(x)
 y = x^par(0)
; curvature term
 if ((par(2) NE 0) AND (par(3) NE 0)) then begin 
     psgn = par(2)/ABS(par(2))
     y = y * ( 1.d0 + ABS(par(2))*x/par(3) )^psgn
 endif 
 vy = x^4 * y 
 dx = ALOG(x(1:np-1)) - ALOG(x(0:np-2)) 
 yi = TOTAL( (vy(1:np-1) + vy(0:np-2))*0.5*dx )
 y = par(1) * y / yi
RETURN, y
END


FUNCTION LOGN, x, par
; generates a volume normalized log-normal law 
; returns distribution in nr of grains : dn/da
; x : grain size 
; par(0) : centroid of log-normal
; par(1) : sigma of log-normal
; par(2) : VOLUME normalization

 np = n_elements(x)
 x0 = par(0)
 sigma = par(1)
 y = exp(- 0.5 * ( alog(x/x0) / sigma )^2 ) / x
 vy = x^4 * y
 xm = par(0) * exp( -par(1)^2 ) ; x of max in dn/da
 print,'(W) LOGN: dn/da max @',xm*1e7,' nm'
 dx = alog(x(1:np-1)) - alog(x(0:np-2)) 
 yi = TOTAL( (vy(1:np-1) + vy(0:np-2))*0.5*dx )
 y = par(2) * y / yi
RETURN, y
END 

FUNCTION CUT_OFF, x, par
; generates the large size cut-off 
; from Weingartner & Draine 2001
; x : grain size 
; par(0) : threshold for cut-off (At)
; par(1) : shape As (roughly the size where cut_off=0.5)

 y = 0.d0*x + 1.d0
 ix = WHERE( x GT par(0), cnt )
 if CNT GT 0 then begin 
     y(ix) = exp( -( (x(ix)-par(0)) / par(1))^3. )
 endif else begin 
     print,'(W) CUT_OFF: no size larger than At found'
 endelse
RETURN, y
END


FUNCTION 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 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 ''PLAW'' or ''LOGN'''
   print,'           [Default = ''PLAW'' 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 = CREATE_VDIST( [3.e-8,2.5e-5], 3.3, slaw=''plaw'', 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 = CREATE_VDIST( [3.e-8,2.5e-5], 2.25, slaw=''logn'', 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) CREATE_VDIST: you must define a size range'
     RETURN,0
 endif 
 if n_elements( RHO ) EQ 0 then begin 
     print,'(F) 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=['PLAW']
   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 'PLAW') AND (SLAW(i) NE 'LOGN') then begin
         print,'(F) CREATE_VDIST: undefined law ',slaw(i) 
         print,'                  only ''PLAW'' and ''LOGN'' '
         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 = [ 'LOGN', 'LOGN', 'PLAW' ]
     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) 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) CREATE_VDIST : cut-off applied, size and scale are ',cutof
    vdist = vdist * 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) 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) 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) CREATE_VDIST: mass distribution written in ', strtrim(fname,2)
 endif

RETURN, vdist
END