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