FUNCTION DUSTEM_SERKOWSKI, x, ka=ka, xmax=xmax, pmax=pmax

;+
; NAME:
;   dustem_serkwoski
;  
; PURPOSE:
;   computes the Serkowski law (as per Draine & Fraisse 2009)
;  
; CATEGORY:
;    DustEMWrap, Distributed, LowLevel, Helper
;  
; CALLING SEQUENCE:
;   serk = DUSTEM_SERKOWSKI(X,[KA=KA,XMAX=XMAX,PMAX=PMAX])
;
; INPUTS:
;  X  -- array(n_qabs) inverse wavenumber in 1/microns'
;
;  OPTIONAL INPUT PARAMETERS:
;     KA   (I): K factor, default = 0.92
;     XMAX (I): x-position of max, default = 1.82
;     PMAX (I): max polar. fraction, default = 0.03
;
; OUTPUTS:
;
; OPTIONAL OUTPUT PARAMETERS:
;
; ACCEPTED KEY-WORDS:
;    help = print this help
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURES AND SUBROUTINES USED:
;   
; EXAMPLES
;  
; MODIFICATION HISTORY:
;    Written by LV (VG ?)
;    Initially distributed with the Fortran, incorporated into DustEMWrap 2022
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-

  
  if N_PARAMS() LT 1 or keyword_set(help) then begin 
     doc_library,'dustem_serkowski'
     yy=0
     goto, the_end
  end


  ;; IF N_PARAMS() EQ 0 THEN BEGIN  
  ;;    print,'FUNCTION SERKOWSKI, x, ka=ka, xmax=xmax, pmax=pmax '
  ;;    print,' computes the Serkowski law (Draine & Fraisse 2009)'
  ;;    print,''
  ;;    print,' X    (I): array(n_qabs) inverse wavenumber in 1/microns'
  ;;    print,' KA   (I): K factor'
  ;;    print,' XMAX (I): x-position of max'
  ;;    print,' PMAX (I): max polar. fraction'
  ;; ENDIF

 IF n_elements(KA) EQ 0 THEN ka = 0.92 ; Draine & Fraisse (2009)
 IF n_elements(XMAX) EQ 0 THEN xmax = 1.82
 IF n_elements(pMAX) EQ 0 THEN pmax = 0.03

 yy = pmax * EXP( -ka * ALOG(xmax/x)^2 )
 x0 = 1/1.39
 y0 = pmax * EXP( -ka * ALOG(xmax/x0)^2 )
 ix = WHERE( x LE x0, cx )
 IF cx GT 0 THEN  yy(ix) = y0 * (x(ix)/x0)^1.7

 
the_end:
 RETURN, yy

END