function redden,lambda,flux,extin,rv,nonorm=nonorm
; This function accepts arguments lambda (nm), flux(lambda), and extinction
; parameters extin = A(V) and rv = A(V)/E(B-V).  It returns the reddened
; flux distribution on the same wavelength range, computed according to
; the parameterization of Cardelli et al. 1989 ApJ 345,245.
; By default, output fluxes are normalizes so that the integral from
; 500 to 550 nm is unity.  To suppress this normalization, set the nonorm
; keyword.

if(extin eq 0.) then begin   ; bail out cheaply in common case of no extinction
  fluxo=flux
  goto,fini
endif

; constants
lammin=500.
lammax=550.
xb1=1.1                  ; boundary between infrared and optical/NIR
xb2=3.3                  ; boundary between optical/NIR and UV
xb3=8.0                  ; max value of x for UV bump
xb4=10.                  ; max value of x for Far UV
xbf=5.9                  ; min value for which fa, fb not zero
a0=[.574,1.61]           ; coefficients for IR a
b0=[-.527,1.61]          ; coefficients for IR b
yoff=1.82                ; offset to convert x to y
a1=[1.,.17699,-.50447,-.02427,.72085,.01979,-.77530,.32999]
		   ; coeffs giving a(y) for O/NIR
b1=[0.,1.41338,2.28305,1.07233,-5.38434,-.62251,5.30260,-2.09002]
		   ; coeffs giving b(y) for O/NIR
a2=[1.752,-.316,-.104,4.62,.341]
		   ; coeffs for a(y) for UV
b2=[-3.090,1.825,1.206,4.62,.263]
		   ; coeffs for b(y) for UV
fa2=[-.04473,5.9,-.009779,5.9]
		   ; coeffs for F_a(x)
fb2=[.2130,5.9,.1207,5.9]
		   ; coeffs for F_b(x)
a3=[-1.073,-.628,8.,.137,8.,-.070,8.]
		   ; coeffs for a(x) for far UV
b3=[13.670,4.257,8.,.420,8.,.374,8.]
		   ; coeffs for b(x) for far UV

; make inverse wavelength scale x (in micron^-1)
x=1000./lambda
nlam=n_elements(lambda)

; find indices of wavelength ranges
s0=where(x gt 0. and x le xb1,ns0)
s1=where(x gt xb1 and x le xb2,ns1)
s2=where(x gt xb2 and x le xb3,ns2)
s3=where(x gt xb3 and x le xb4,ns3)

; make extinction(x) (magnitudes)
ee=fltarr(nlam)
if(ns0 gt 0) then begin
  a=a0(0)*x(s0)^a0(1)
  b=b0(0)*x(s0)^b0(1)
  ee(s0)=a+b/rv
endif
if(ns1 gt 0) then begin
  y=x(s1)-yoff
  y2=y*y
  y3=y2*y
  y4=y3*y
  y5=y4*y
  y6=y5*y
  y7=y6*y
  a=a1(0)+a1(1)*y+a1(2)*y2+a1(3)*y3+a1(4)*y4+a1(5)*y5+a1(6)*y6+a1(7)*y7
  b=b1(0)+b1(1)*y+b1(2)*y2+b1(3)*y3+b1(4)*y4+b1(5)*y5+b1(6)*y6+b1(7)*y7
  ee(s1)=a+b/rv
endif
if(ns2 gt 0) then begin
  y=x(s2)
  s2f=where(y ge xbf and y le xb3,ns2f)
  fa=fa2(0)*(y-fa2(1))^2 + fa2(2)*(y-fa2(3))^3
  fb=fb2(0)*(y-fb2(1))^2 + fb2(2)*(y-fb2(3))^3
  a=a2(0)+a2(1)*y+a2(2)/((y-a2(3))^2+a2(4))
  b=b2(0)+b2(1)*y+b2(2)/((y-b2(3))^2+b2(4))
  if(ns2f gt 0) then begin
    a(s2f)=a(s2f)+fa(s2f)
    b(s2f)=b(s2f)+fb(s2f)
  endif
  ee(s2)=a+b/rv
endif
if(ns3 gt 0) then begin
  y=x(s3)
  a=a3(0)+a3(1)*(y-a3(2))+a3(3)*(y-a3(4))^2+a3(5)*(y-a3(6))^3
  b=b3(0)+b3(1)*(y-b3(2))+b3(3)*(y-b3(4))^2+b3(5)*(y-b3(6))^3
  ee(s3)=a+b/rv
endif

; correct fluxes for extinction
fluxo=flux*10.^(-0.4*extin*ee)

; renormalize
sn=where(lambda ge lammin and lambda le lammax,nsn)
if(nsn gt 0 and not keyword_set(nonorm)) then begin
  fluxo=fluxo/total(fluxo(sn))
endif

fini: return,fluxo

end