First commit of IAS routines

Annie Hughes
1 parent 11afdcd2
Showing 11 changed files with 1963 additions and 0 deletions Show diff stats
src/idl/dustem_add_inst.pro
src/idl/dustem_add_mod.pro
src/idl/dustem_band_cc.pro
src/idl/dustem_blackbody.pro
src/idl/dustem_chi2.pro
src/idl/dustem_create_rfield.pro
src/idl/dustem_create_vdist.pro
src/idl/dustem_fil_chi2.pro
src/idl/dustem_serkowski.pro
src/idl/dustem_show_fortran.pro
src/idl/dustem_str_inst.pro
@@ -0,0 +1,30 @@
+FUNCTION DUSTEM_ADD_INST, sf, tag, nsz, ntrans=ntrans
+; fills in data for existing instruments in SMDAT structure 
+; or add a new instrument to SMDAT structure 
+;
+; SF     (I): SMDAT input structure
+; TAG    (I): name of instrument to add in
+; NSZ    (I): array(2) of sizes [nr of data pts, nr of grain types] 
+; NTRANS (I): nr of points for transmission (bad flux data)
+
+  stag = TAG_NAMES(sf)
+  ntag = N_TAGS(sf)
+  tag = STRUPCASE('I_'+STRTRIM(tag,2))
+  itg = WHERE( stag EQ tag, ctg )
+  if n_elements(nsz) EQ 2 then begin
+     nx = nsz(0) & ntype = nsz(1)
+  endif else begin
+     print,'(F) DUSTEM_ADD_INST: array of sizes must be 2D'
+  endelse
+  if ctg EQ 0 then begin
+     if n_elements(ntrans) EQ 0 then st = CREATE_STRUCT( tag, DUSTEM_STR_INST(nx,n2=ntype) ) $
+     else st = CREATE_STRUCT( tag, DUSTEM_STR_INST(nx,n2=ntype,n3=ntrans) )
+     sf = CREATE_STRUCT( sf, st )     
+     stag = TAG_NAMES(sf)
+     itg = WHERE( stag EQ tag )
+     sf.(itg).isel = intarr(nx) + 1
+  endif else begin
+     print,'(F) DUSTEM_ADD_INST: INST already exists in structure'
+  endelse
+  return,sf
+END
@@ -0,0 +1,54 @@
+FUNCTION DUSTEM_ADD_MOD, sf, tag, nsz, unit=unit
+; adds a model tag to a existing SMDAT structure
+; SF     (I): input SMDAT structure 
+; TAG    (I): tag of model to be added
+; NSZ    (I): int array(2) of sizes [nr of data pts, nr of grain types] 
+
+  tag = STRLOWCASE(STRTRIM(tag,2))
+  if n_elements(nsz) GT 1 then begin
+     n1 = nsz(0) & n2 = nsz(1)
+     if n_elements(nsz) EQ 3 then n3 = nsz(2) else n3 = 0
+  endif else begin
+     print,'(F) ADD_MOD: array of sizes must > 1D'
+  endelse
+  if TAG EQ 'emis' then begin
+     if n_elements(unit) EQ 0 then unit='x(microns) SED(erg/s/cm2/sr)'
+     s1 = { UNIT : unit,          $
+            X :    dblarr(n1),    $ ; wave 
+            Y :    dblarr(n1,n2), $ ; SED(wave, grain type) (index ntype is total)
+            YP :   dblarr(n1,n2)  } ; polarized SED(wave, grain type) (index ntype is total)
+  ENDIF else if TAG EQ 'ext' then begin
+     if n_elements(unit) EQ 0 then unit='x(microns) sigma(cm2/H)'
+     s1 = {UNIT :   unit,          $
+           X :      dblarr(n1),    $ ; wave
+           Y :      dblarr(n1,n2), $ ; sigma_ext(wave, grain type) (index ntype is total)
+           ABS :    dblarr(n1,n2), $ ; sigma_abs(wave, grain type) 
+           SCA :    dblarr(n1,n2), $ ; sigma_sca(wave, grain type)
+           ABS_P :  dblarr(n1,n2), $ ; absorption sigma_pol(wave, grain type)
+           SCA_P :  dblarr(n1,n2), $ ; scattering sigma_pol(wave, grain type)
+           ALB :    dblarr(n1,n2), $ ; alb(wave, grain type)
+           XR :     0d,            $ ; ref wave
+           YR_ABS : dblarr(n2),    $ ; tau_abs/NH @ XR
+           YR_SCA : dblarr(n2),    $ ; tau_sca/NH @ XR
+           RV :   0d               }
+;  ENDIF else if TAG EQ 'pol' then begin
+;     if n_elements(unit) EQ 0 then unit='x(microns) SED(erg/s/cm2/sr)'
+;     s1 = {UNIT :    unit,          $
+;           X:        dblarr(n1),    $ ; wave
+;           Y:        dblarr(n1,n2), $ ; polarized SED(wave, grain type) (index ntype is total)
+;           POL:      dblarr(n1,n2)  } ; sigma_pol(wave, grain type) (ABS + SCAT)
+  ENDIF else if tag EQ 'sdist' then begin
+     if n3 EQ 0 then begin 
+        print,'(F) ADD_MOD: 3d dimension missing '
+        return,0
+     endif
+     if n_elements(unit) EQ 0 then unit='x(cm) a^4*dn/da(cm3/H)'     
+     s1 = {UNIT :    unit,          $
+           XTOT :    dblarr(n1),    $ ; grain size
+           YTOT :    dblarr(n1,n2), $ ; size distribution
+           XI :      dblarr(n3,n2), $ ; grain size per type
+           YI :      dblarr(n3,n2)  } ; grain size dist per type
+  ENDIF
+  sm = CREATE_STRUCT( CREATE_STRUCT('M_'+tag,s1), sf )
+  return, sm
+END
@@ -0,0 +1,101 @@
+FUNCTION dustem_BAND_CC, w0, xs, ys, xt, yt, y0=y0, cc=cc, rr=rr, nint=nint
+; returns flux in given instrument band (erg/cm2/sr)
+; W0    (I): centroid of filter (any unit)
+; XS    (I): input x-grid for SED (same unit as W0)
+; YS    (I): input SED in erg/cm2/sr
+; XT  (I/O): x-grid of filter (same unit as W0)
+; YT  (I/O): filter transmission (arb. unit)
+; RR    (O): undersampling ratio: dx(SED)/dx(filt), usually 1. 
+;            If RR<1 filter oversampled in flux integration
+; Y0    (O): SED at band center
+; CC    (O): color correction
+; NINT  (O): nr of points in band integration
+
+ vlow = 1d-5
+
+; remove bad and multiple values in filter
+  xtt=xt & ytt=yt
+  iq = UNIQ(xtt,sort(xtt))
+  xtt = xtt(iq) & ytt = ytt(iq)
+  in = where( (xtt GT 0d) AND (ytt GT 0d), cn )
+  if cn GT 0 then begin
+     xtt = xtt(in)
+     ytt = ytt(in)
+  endif else begin 
+; if no good point in filter, band flux and the rest is 0 
+     sed=0d & rr=0d & y0=0d & cc=0d & nint=n_elements(xt)
+     return,sed
+  endelse
+  nfilt = n_elements(xtt)
+
+; sort arrays  
+  ix = SORT(xtt)
+  xtt = xtt(ix) & ytt = ytt(ix)
+  ytt = ytt / MAX(ytt)
+  ix = SORT(xs)
+  xss = xs(ix) & yss = ys(ix)
+  xt=xtt & yt=ytt
+
+; shrink SED array
+  ix = WHERE( xss GE MIN(xtt) AND xss LE MAX(xtt), cx)
+  if cx GT 0 then begin
+; add points at edges for interpolation
+     i1 = 0 
+     if ix(0) GT i1 then i1 = ix(0)-1
+     i2 = n_elements(xss)-1
+     if ix(cx-1) LT i2 then i2 = ix(cx-1)+1
+     xss = xss(i1:i2)
+     yss = yss(i1:i2)
+  endif else begin 
+; if no good point in SED, band flux and the rest is 0 
+     sed=0d & rr=0d & y0=0d & cc=0d & nint=n_elements(xt)
+     return,sed
+  endelse
+  nsed = n_elements(xss)
+
+; get filter sampling f_res = dx(filt)
+  ix = WHERE( ytt/MAX(ytt) GT vlow, cx )
+  f_res = ABS( MEDIAN( xtt(ix(1:cx-1))-xtt(ix(0:cx-2)) ) )
+
+; Interpolate to compute integration
+  IF NFILT GE NSED THEN BEGIN
+; finer filter grid: integrate band flux on filter x-grid
+     xi = xtt
+     rr = 1d                    ; SED sampling is filter sampling
+     yi = INTERPOL(yss,xss,xi)
+     filt = ytt
+  ENDIF ELSE BEGIN
+; finer SED grid: integrate band flux on SED x-grid
+     xi = xss & yi = yss 
+; get SED sampling: s_res = dx(SED) and sampling ratio
+     s_res = ABS( MEDIAN( xi(1:nsed-1)-xi(0:nsed-2) ) )
+     rr = s_res / f_res
+     filt = INTERPOL(ytt, xtt, xi)
+     io = WHERE( xi GT MAX(xtt) OR xi LT MIN(xtt), co)
+     if co GT 0 then filt(io) = 0d
+     in = where( filt LT 0d, cn) 
+     if cn GT 0 then filt(in) = 0d
+  ENDELSE
+  xt=xi & yt=filt
+
+; Integrate on filter bandpass (log-grid)
+; trapezium safe and accurate enough - 
+  np = n_elements(xi)
+  nint = np
+  dx = xi(1:np-1) - xi(0:np-2)
+  yint = yi*filt/xi
+;  a1 = INT_TABULATED(xi,yint, /double)
+  yint = (yint(0:np-2) + yint(1:np-1)) / 2d0
+  a1 = TOTAL( yint*dx ) 
+  yint = filt/xi
+;  a2 = INT_TABULATED(xi,yint, /double)
+  yint = (yint(0:np-2) + yint(1:np-1)) / 2d0
+  a2 = TOTAL( yint*dx )
+  
+; get color correction, color corrected band flux and sed
+  y0 = INTERPOL( yi, xi, w0)
+  cc = a1/a2/y0
+  sed = y0*cc 
+
+  RETURN, sed
+END
@@ -0,0 +1,152 @@
+FUNCTION DUSTEM_BLACKBODY, x, temp, unit=unit, wdil=wdil, check=check, g0=g0, chi=chi, $
+		     bb_to_hab=bb_to_hab, NORM=norm
+  if n_params() EQ 0 then begin 
+     print,'------------------------------------------------------------------------------------------'
+     print,'FUNCTION DUSTEM_BLACKBODY, x, temp, unit=unit, wdil=wdil, NORM=norm, check=check, g0=g0, chi=chi'
+     print,'------------------------------------------------------------------------------------------'
+     print,' computes a (or a sum of) blackbody brightness(es) B(x,T) in W/m2/(x unit)/sr '
+     print,' X	(I): blackbody abscissa in cm-1, GHz, microns or eV'
+     print,' TEMP  (I): blackbody temperature (can be an array)'
+     print,' WDIL  (I): dilution factor for each TEMP. [Default=1].'
+     print,' UNIT  (I): blackbody x unit ''K'' is B_nu (W/m2/cm-1/sr),   x in cm-1 '
+     print,'                             ''F'' is B_nu (W/m2/Hz/sr),     x in GHz [default]'
+     print,'	                         ''W'' is B_lambda (W/m2/um/sr), x in microns.' 
+     print,'	 	  	         ''E'' is B_E (W/m2/eV/sr), x in eV'
+     print,' NORM  (I): overall normalizing factor, returns NORM*BB/MAX(BB)'
+     print,' CHECK (O): check energy conservation for the BB (power is sigma*T^4)'
+     print,' CHI   (O): scaling factor @ 1000 A in ISRF unit (1.6e-3 erg/cm2/s)'
+     print,' G0    (O): scaling factor for integrated 5-13.6 eV flux in ISRF unit.'
+     print,''
+     print,' Created 2001, L Verstraete IAS'
+     print,' More units and checks, Oct 2010 LV'
+     print,'---------------------------------------------------------------------------------------'
+     return,0
+  endif
+ nx = n_elements( X )
+ nt = n_elements( TEMP )
+ nw = n_elements( WDIL )
+ if nw EQ 0 then begin 
+   wdil = dblarr(nt) + 1. 
+ endif else begin 
+   if nw LT nt then begin 
+     print,'DUSTEM_BLACKBODY (F): WDIL and TEMP have different sizes'
+     return, 0
+   endif 
+ endelse
+ 
+ if n_elements( UNIT ) EQ 0 then unit='F' $
+ else begin 
+  unit = strupcase(strcompress(unit,/rem))
+  if unit NE 'K' AND unit NE 'F' AND unit NE 'W' AND unit NE 'E' then begin 
+    print,'DUSTEM_BLACKBODY (W): not a valid BB unit --> make a B_nu'
+    unit = 'F'
+  endif
+ endelse
+
+ if n_elements( CHECK ) GT 1 then check = 0
+
+ x = double(x) 
+ temp = double(temp) 
+ wdil = double(wdil)
+
+ xpi = 3.1415926535897932385d
+ clight = 2.9979246d08          ; m/s
+ kb = 1.3806503d-23
+ hp = 6.62606876d-34
+ hck = hp*clight / kb
+ ze  = 1.60217653d-19           ; electron charge
+ hev = hp/ze                    ; h in eV unit
+
+ xly = 12.4
+ isrf0 = 1.6d-6         ; 4*pi*nu*I_nu (W/m2) in solar neighborhood
+
+;
+; loop on temperatures
+; 
+ bbt = dblarr(nx)
+ acheck = 0.
+
+ FOR j = 0, nt-1 do begin 
+
+    xi = [6., 13.6]
+
+    CASE UNIT OF
+
+       'K': begin      
+          be = exp( -(1d2*x) * hck / temp(j) )
+          bb = be * (1d2*x)^3 / (1.d0 - be)
+          bb = bb * 2d0 * hp * clight^2 
+          bb = bb * 1d2         ; m-1 to cm-1
+          xi = xi / hev/clight/1d2   
+          x1 = xly / hev/clight/1d2 
+       end
+
+       'F': begin      
+          be = exp( -(1d9*x) * hp/kb/ temp(j) )
+          bb = be * (1d9*x)^3 / (1.d0 - be)
+          bb = bb * 2d0 * hp / clight^2
+          xi = xi / hev / 1d9
+          x1 = xly / hev / 1d9
+       end
+
+       'W': begin
+          be = exp( -hp*clight/kb/ (1d-6*x) / temp(j) )
+          bb = be / (1d-6*x)^5 / (1.d0 - be)
+          bb = bb * 2d0 * hp * clight^2 
+          bb = 1d-6 * bb        ; m-1 to micron-1
+          xi = hev*clight*1d6 / xi
+          x1 = hev*clight*1d6 / xly
+       end
+
+       'E': begin
+          be = exp( -(ze*x)/kb / temp(j) )
+          bb = be * (ze*x)^3 / (1.d0 - be)
+          bb = bb * 2d0 / hp^3 / clight^2 
+          bb = bb * ze          ; J-1 to eV-1
+          xi = xi
+          x1 = xly
+       end
+
+    ENDCASE
+
+    ck = 0.
+    for jx = 1L,nx-1 do $
+       ck = ck + ( bb(jx)+bb(jx-1) )*( x(jx)-x(jx-1) ) / 2.d0
+    ck = xpi* ck / 5.6697d-8 / temp(j)^4  
+    
+    if n_elements( CHECK ) NE 0 then begin 
+       if check EQ 1 then $
+          print,'DUSTEM_BLACKBODY (W): integrated BB('+strtrim(temp(j),2)+$
+                ' K) is '+ strtrim(ck,2)+ ' times Sigma*T^4'
+    endif
+    acheck = [ acheck, ck ]
+    
+    bbt = bbt + bb*wdil(j)
+    
+ ENDFOR
+ 
+ check = acheck(1:*)
+
+; get G0
+ ig = where( x GE MIN(xi) AND x LE MAX(xi), cg)
+ if cg GT 0 then begin
+    xx=x(ig)
+    yy=xpi*x(ig)*bbt(ig)
+    s1 = TOTAL( (yy(1:cg-1)+yy(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 )
+    g0 = s1/isrf0
+ endif else g0 = 0.d0
+
+; get chi
+ ig = where( abs(2.*(x-x1)/(x+x1)) LE 1.d-1, cg) 
+ if cg GT 0 then begin 
+   tmp = min( abs(x(ig)-x1), i_min )   
+   ig = ig(i_min)
+   ig = ig(0)
+   chi = xpi * x(ig) * bbt(ig) / isrf0
+ endif else chi = 0
+
+ if n_elements( NORM ) NE 0 then bbt = norm * bbt/max(bbt)
+                                                
+the_end:
+ RETURN, BBT
+END                             ;FUNCTION BLACKBODY
@@ -0,0 +1,17 @@
+FUNCTION DUSTEM_CHI2, y, model, npar, err=err
+; returns the chi-square value of fit MODEL to data Y
+;
+; NPAR (I): nr of parameters in the fit
+; ERR  (I): error of each data point. Default is ERR=0.1*Y
+
+  ny = n_elements(y)
+  if n_elements(err) EQ 0 OR TOTAL(err) EQ 0 then begin
+     err = 0.1*y
+     print,'(W) CHI2: error missing, set to 10 %'
+  endif
+  ndof = ny - npar              ; nr of degrees of freedom
+  chi = TOTAL( ((y-model)/err)^2 )
+  if ndof GT 0 then chi = chi / ndof
+
+  return, chi
+END
@@ -0,0 +1,230 @@
+FUNCTION HABING_FIELD, x, unit=unit
+; computes the Habing interstellar radiation field SED (ISRF)
+; in W/m2 (4*!pi*nu*I_nu)
+; (from Draine & Bertoldi 1996)
+; X	(I): X-grid for the ISRF
+; UNIT	(I): unit of the ISRF. Choices are 
+;	     'W': wave in um [Default]
+;	     'F': frequency in cm-1
+;	     'E': energy in eV 
+
+ if n_elements( UNIT ) EQ 0 then unit = 'W' $
+ else unit = strupcase( strcompress( unit,/rem) )
+
+ x = double( x )
+
+ CASE unit of 
+
+   'W': x3 = 1.d1 * x 
+
+   'F': x3 = 1.d5 / x
+
+   'E': x3 = 12.4 / x
+
+ ENDCASE
+
+ field = - x3^3*4.1667 + x3^2*12.5 - x3*4.3333 
+ field = 1.d-1 * 3.d8 * 1.d-14 * field 
+
+ i_neg = where( field LT 0.d0, c_neg )
+ field( i_neg ) = 0.d0
+
+ RETURN, field
+END                             ;FUNCTION HABING_FIELD
+
+
+FUNCTION MATHIS_FIELD, x, unit=unit
+; computes the Mathis interstellar radiation field SED (ISRF)
+; in W/m2 (4*!pi*nu*I_nu)
+; from Mathis et al. 1983, A&A 128, 212
+; X	(I): X-grid for the ISRF
+; UNIT	(I): unit of the ISRF. Choices are 
+;	     'W': wave in um [Default]
+;	     'F': frequency in cm-1
+;	     'E': energy in eV
+
+ if n_elements( UNIT ) EQ 0 then unit = 'W' $
+ else unit = strupcase( strcompress( unit,/rem) )
+
+ ly_limit = 9.11267101235d-2
+
+ x = double( x )
+
+;
+; visible part 
+;
+ wdil = 4.d0 * [ 4.d-13, 1.d-13, 1.d-14]	; Mathis definition
+						; 4*!pi*Wdil*B_nu
+ field = !pi * x * BLACKBODY( x, [ 3.d3, 4.d3, 7.5d3 ], wdil=wdil, unit=unit )
+
+;
+; UV part
+;
+
+; first convert to (lambda / 1e3 AA)
+ CASE unit of 
+
+   'W': x3 = 1.d1 * x 
+
+   'F': x3 = 1.d5 / x
+
+   'E': x3 = 12.4 / x
+
+ ENDCASE
+ 
+ il = where( x3 LE 2.46, cl ) 
+ if cl GT 0 then field(il) = 0.D0
+
+ il = where( x3 GE 1d1*ly_limit AND x3 LT 1.11, cl )
+ if cl GT 0 then field(il) = 1.4817d-6 * x3(il)^(4.4172)
+ il = where( x3 GE 1.11 AND x3 LT 1.34, cl ) 
+ if cl GT 0 then field(il) = 2.0456d-6 * x3(il)
+ il = where( x3 GE 1.34 AND x3 LT 2.46, cl )
+ if cl GT 0 then field(il) = 3.3105d-6 * x3(il)^(-0.6678)
+  
+
+ RETURN, field
+END                             ;FUNCTION MATHIS_FIELD
+
+
+FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname
+
+  if N_PARAMS() LT 1 then begin 
+     print,'------------------------------------------------------------------------------------------------------'
+     print,'FUNCTION CREATE_RFIELD, temp, x=x, isrf=isrf, wdil=wdil, wcut=wcut, g0=g0, chi=chi, fname=fname'
+     print,'------------------------------------------------------------------------------------------------------'
+     print,''
+     print,' generates the radiation field for DUSTEM (ISRF.DAT)'
+     print,' in erg/cm2/s/Hz (4*!pi*I_nu)'
+     print,' RFIELD = ISRF + WDIL*PI*BB(TEMP) '
+     print,' ISRF from Mathis et al. 1983, A&A 128, 212'
+     print,' NB: if blackbody is isotropic set WDIL to 4 (to get 4*!pi factor)'
+     print,''  
+     print,' TEMP (I): blackbody temperature (can be an array). If 0 only ISRF.'
+     print,' X    (I): X-grid for the ISRF in microns. Default is 200 pts over 0.01-10^5 microns.'
+     print,'           (including WCUT point). You want to include WCUT for accurate edges.'
+     print,' ISRF (I): if set to 0 no ISRF is added, 1 is Mathis (default), 2 is Habing'
+     print,' WDIL (I): blackbody dilution factor (can be an array)'
+     print,' FNAME(I): filename for ISRF.DAT'
+     print,' WCUT (I): for wave < wcut radiation field is 0 '
+     print,' G0   (O): factor flux(6-13.6eV) wrt. Mathis field '
+     print,' CHI  (O): scaling factor at 100nm wrt. Mathis field'
+     print,''
+     print,'Example: tt = create_rfield([2d4,5d4],wdil=[1.d-14,1.d-16],x=x,fname=''ISRF.DAT'')'
+     print,''
+     print,' Created Aug. 2009, L. Verstraete, IAS'
+     print,' Force the WCUT point, May 2011, LV'
+     print,''
+     print,'------------------------------------------------------------------------------------------------------'
+     RETURN,0
+  endif
+
+; inits 
+  unit='W'
+  ly_limit =  9.11267101235d-2
+  xi = [ ly_limit, 2.07d-1 ] ; for G0 6-13.6 eV in microns
+  x1 = 1.d-1 ; for CHI
+
+  if n_elements(WCUT) EQ 0 then wcut=ly_limit  
+  if n_elements(x) EQ 0 then begin
+     nx = 199
+     xb = [ 1.d-2, 1.d5 ] ; wave boundaries in microns
+     dx = (alog(xb(1))-alog(xb(0))) / (nx-1)
+     x = alog(xb(0)) + dindgen(nx)*dx
+     x = exp(x)
+     x = [x, wcut] ; add wcut to avoid coarse edge
+     x = x(UNIQ(x,SORT(x)))
+     nx = n_elements(x)
+     if x(nx-1) NE xb(1) then x(nx-1)=xb(1) ; check rounding
+  endif else begin
+     print,'(W) CREATE-RFIELD : using input wave x'
+     x = double(x)
+     nx = n_elements(x)
+     ix = WHERE( ABS(x-wcut)/x LE 0.01, cx )
+     if cx EQ 0 then begin
+        print,'(W) CREATE_RFIELD: your x-grid does not contain wcut.'
+        print,'                   Should be included if radiation is 0 below wcut.'
+     endif 
+  endelse
+
+  rfield = dblarr(nx)
+; get Habing for normalization
+  rhabing = HABING_FIELD(x,unit=unit)
+  rhabing =  1.d3*x*rhabing/3.d14 ; erg/cm2/s/Hz
+
+; get isrf
+  rmathis = MATHIS_FIELD(x,unit=unit)
+  rmathis =  1.d3*x*rmathis/3.d14 ; erg/cm2/s/Hz
+
+  if n_elements(ISRF) EQ 0 then ISRF=1
+  if ISRF EQ 1 then begin 
+     print,'(W) CREATE_RFIELD : adding Mathis ISRF'
+     rfield = rmathis
+  endif else if ISRF EQ 2 then begin 
+     print,'(W) CREATE_RFIELD : adding Habing ISRF'
+     rfield = rhabing
+  endif
+
+; get blackbody  
+  ntemp = n_elements(TEMP)
+  if ntemp GT 0 then begin 
+     bb = BLACKBODY( x, temp, unit=unit, wdil=wdil )
+     bb = bb * 1.d3 * !pi * x^2 / 3.d14   ; erg/cm2/s/Hz
+     print,'(W) CREATE_RFIELD : adding BB with T= ',temp, format='(A38,10(1E10.4,1x))'
+     print,'                dilution factor wdil= ',wdil, format='(A38,10(1E10.4,1x))'
+  endif else bb=0.D0
+  rfield = rfield + bb
+
+; apply cut
+  ix = WHERE( x LT WCUT, cx )
+  if cx GT 0 then rfield(ix)=0.d0
+
+; get G0
+ ig = where( x GE xi(0) AND x LE xi(1), cg)
+ if cg GT 0 then begin
+    xx=x(ig)
+    yy=rfield(ig)
+    rr=rmathis(ig)
+    s1 = TOTAL( (yy(1:cg-1)+yy(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 )
+    s2 = TOTAL( (rr(1:cg-1)+rr(0:cg-2))*(xx(1:cg-1)-xx(0:cg-2))/2.D0 )
+    g0 = s1/s2
+ endif else g0 = 0.d0
+ print,'(W) CREATE_RFIELD : G0 =',g0
+
+; get chi
+ ig = where( abs(2.*(x-x1)/(x+x1)) LE 1.d-1, cg) 
+ if cg GT 0 then begin 
+   tmp = min( abs(x(ig)-x1), i_min )   
+   ig = ig(i_min)
+   ig = ig(0)
+   chi = rfield(ig) / rmathis(ig)  
+ endif else chi = 0
+ print,'(W) CREATE_RFIELD : chi =',chi
+
+; write file
+ if n_elements(FNAME) NE 0 then begin
+   OPENW, iu, fname, /get_lun
+   printf, iu, '# DUSTEM: exciting radiation field featuring '
+   if ISRF EQ 1 then printf, iu, '# Mathis ISRF'
+   if NTEMP GT 0 then begin 
+      a =  '# Blackbody with     T='
+      a1 = '# dilution factor wdil='
+      for i = 0, ntemp - 1 do begin
+         a = a   + string( format='(2x,1E11.4)', temp(i) )
+         a1 = a1 + string( format='(2x,1E11.4)', wdil(i) )
+      endfor
+      printf, iu, a
+      printf, iu, a1
+   endif
+   printf, iu, '# Nbr of points'
+   printf, iu, '# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'
+   printf, iu, nx, format='(i4)'
+   for i=0,nx-1 do begin 
+    printf, iu, x(i), rfield(i), format='(2(1E13.6,2x))'
+   endfor
+   FREE_LUN, iu
+   print,'(W) CREATE_RFIELD: radiation field written in ', strtrim(fname,2)
+ ENDIF
+
+ RETURN, RFIELD
+END
@@ -0,0 +1,224 @@
+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
@@ -0,0 +1,36 @@
+FUNCTION DUSTEM_FIL_CHI2, sf, ntype=ntype
+; fills in the Chi-square fields of an input SMDAT structure 
+; (see format in DUSTEM_GET_BAND_FLUX)
+;
+; SF    (I): input SMDAT structure 
+; NTYPE (I): index of SED to be used 
+
+  if ntype EQ 0 then begin
+     tt = SIZE(sf.sed.y)
+     ntype = tt(2)
+  endif
+  ntag = N_TAGS(sf)
+  stag = TAG_NAMES(sf)
+  itg = WHERE( STRPOS(stag,'I_') GE 0, ctg )
+  if ctg EQ 0 then begin
+     print,'(F) DUSTEM_FIL_CHI2: no instrument to fill in'
+     return,0
+  endif
+  yy=0d & mm=0d & ee=0d & ift=0
+  for k = 0, ctg-1 do begin 
+     i1 = WHERE( sf.(itg(k)).isel EQ 1 AND sf.(itg(k)).err GT 0, c1)
+     if c1 GT 0 then begin
+        sf.(itg(k)).chi2 = $
+           DUSTEM_CHI2(sf.(itg(k)).yd(i1), sf.(itg(k)).ym(i1,ntype-1), sf.(itg(k)).npar, err=sf.(itg(k)).err(i1))
+     endif
+     yy = [ yy, sf.(itg(k)).yd ]
+     ee = [ ee, sf.(itg(k)).err]
+     mm = [ mm, sf.(itg(k)).ym(*,ntype-1) ]
+     ift = [ ift, sf.(itg(k)).isel ]
+  endfor
+  yy=yy(1:*) & ee=ee(1:*) & mm=mm(1:*) & ift=ift(1:*)
+  it = WHERE( ift EQ 1 AND ee GT 0, ct)
+  if ct GT 0 then sf.chi2 = DUSTEM_CHI2( yy(it), mm(it), sf.npar, err=ee(it))
+
+  return, sf
+END
@@ -0,0 +1,22 @@
+FUNCTION DUSTEM_SERKOWSKI, x, ka=ka, xmax=xmax, pmax=pmax
+  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
+
+ RETURN, yy
+END
@@ -0,0 +1,1061 @@
+PRO dustem_show_fortran, model=model, data_dir=data_dir,dustem_dir=dustem_dir,nh=nh, fsed=fsed, sw=sw, inst=inst, smdat=smdat, com=com, wref=wref, bbpar=bbpar, wext=wext,  $
+                 xr=xr, yr=yr, tit=tit, CMB=cmb, COSMO=COSMO, DATA=data, SHOW=show, wn=wn, HARD=hard
+
+;+
+; NAME:
+;    dustem_show_fortran
+; PURPOSE:
+;    shows differents results of the DUSTEM fortran calculations,
+;    overlaid on fiducial diffuse ISM dust observables
+; CATEGORY:
+;    DustEM, Distributed, User-Example
+; CALLING SEQUENCE:
+;    dustem_show_results
+; INPUTS:
+;    DATA_DIR : path to directory containing fiducial observational
+;    SED files. Defaults to Data/EXAMPLE_OBSDATA/ subdirectory.
+;    DUSTEM_DIR : path to directory containing DustEM data files. Default
+;                 is the !dustem_dat directory (which is defined in
+;                 a user's idl_startup file)
+;    model = specifies the interstellar dust mixture used by DustEM
+;           'MC10' model from Compiegne et al 2010 (default)
+;           'DBP90' model from Desert et al 1990
+;           'DL01' model from Draine & Li 2001
+;           'WD01_RV5p5B' model from Weingartner & Draine 2002 with Rv=5.5
+;           'DL07' model from Draine & Li 2007
+;           'J13' model from Jones et al 2013, as updated in
+;                 Koehler et al 2014
+;           'G17_ModelA' model A from Guillet et al (2018). Includes
+;                 polarisation. See Tables 2 and 3 of that paper for details.
+;           'G17_ModelB' model B from Guillet et al (2018)
+;           'G17_ModelC' model C from Guillet et al (2018)
+;           'G17_ModelD' model A from Guillet et al (2018)
+;    NH     : float, fiducial column density for emission [Default is 1e20 cm-2]'
+;    FSED   : string, read SED from the file FSED 
+;    SW     : integer, width of window to SMOOTH the PAH emission 
+;    INST   : string array, instruments to be shown. Default is  ['DIRBE','FIRAS','HFI','WMAP']
+;    COM    : string to describe SMDAT structure
+;    WEXT   : normalize the UV-extinction to be 1 at wavelength WEXT (in microns). Also applied to IR-extinction
+;    WREF   : get the dust cross-section per H at wavelength WREF (in microns)
+;    BBPAR  : float tuple, T and beta for modified BB to get dust opacity from FIRAS
+;    XR,YR  : plot ranges
+;    TIT    : title for SMDAT and plots
+;    CMB    : keyword to overlay the CMB
+;    COSMO  : keyword to plot long wavelengths
+;    DATA   : if DATA=0, only the model spectrum is shown
+;             if DATA=1, only the observational SED is shown
+;             if DATA=2, both the model and data are shown
+;             both are shown (default)'
+;    SHOW   : string array, defines what to plot display. Possible
+;             options include
+;             ''emis'', ''extuv'', ''extir'', ''alb'', ''sdist'', ''polext'', ''polsed'', ''align''
+;              Default is [''emis'']. SHOW=0 no plot
+;   WN     : array of window numbers for plots
+;  HARD    : array (as show) or keyword (/HARD) for hardcopy or filename of ps files, emission'
+;             extinction albedo and size dist. 
+;             Default = [''show.ps'']
+; 
+; OPTIONAL INPUT PARAMETERS:
+;    None
+; OUTPUTS:
+; OPTIONAL OUTPUT PARAMETERS:
+; SMDAT : from DUSTEM_GET_BAND_FLUX (see format there), structure containing model projected on data points: I_INST field in SMDAT.
+;             SMDAT also contains model direct outputs in M_EMIS, M_EXT fields (see format in ADD_MOD.pro )
+; ACCEPTED KEY-WORDS:
+;    help      = If set, print this help
+;    restart   = If set, reinitialise DustEMWrap and use a new
+;                wavelength vector for integration. 
+; COMMON BLOCKS:
+;    None
+; SIDE EFFECTS:
+;    None
+; RESTRICTIONS:
+;    The DustEM fortran code must be installed
+;    The DustEMWrap IDL code must be installed
+; PROCEDURES AND SUBROUTINES USED:
+;    *** COMMENT AH --> is this really NONE? ****
+; EXAMPLES
+    ;; print,' Examples:'
+    ;; print,'  path=''/Users/lverstra/DUSTEM_LOCAL/dustem4.0/''' 
+    ;; print,'  to show SED only:                     show_dustem,path'
+    ;; print,'  SED+UV-extinction:                    show_dustem,path,show=''extuv'' '
+    ;; print,'  SED+IR extinction+size dist.:         show_dustem,path,show=[''extir'',''sdist''] '
+    ;; print,'  ps files with default names:          show_dustem,path,show=[''extir'',''alb''],/hard '
+    ;; print,'  ps files with given names:            show_dustem,path,show=[''extir'',''sdist''],hard=[''f1.ps'',''f2.ps''] '
+    ;; print,'  no plot :                             show_dustem,path,show=0'
+    ;; print,'  to get band fluxes in structure SM:   show_dustem,path,smdat=sm'
+    ;; print,'  see structure SM (overall):           help,/str,sm'
+    ;; print,'  see structure SM model field:         help,/str,sm.m_emis'
+    ;; print,'  see structure SM inst field:          help,/str,sm.i_dirbe'
+    ;; print,'  see structure SM inst flux field:     help,/str,sm.i_dirbe.flx'
+    ;; print,'                                        array(i,j+1) i:index of inst band, j:index of grain type (ntype+1 is total SED) '
+    ;; print,'  to select or add instrument field(s): show_dustem,path,smdat=sm,inst=[''spire'']  (on existing SM to add instrumentx) '
+; MODIFICATION HISTORY:
+;    Written by L Verstraete and M Compiegne, Spring 2010
+;    Modified by LV : change to cgs, add band fluxes in SMDAT, 2011
+;    Modified by LV & V Guillet : add polarization part, 2017
+;    Further evolution details on the DustEMWrap gitlab.
+;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
+
+  
+if keyword_set(help) then begin 
+  doc_library,'dustem_show_fortran'
+  goto,the_end
+END
+
+; specify the grain model
+IF keyword_set(model) THEN BEGIN
+  use_model=strupcase(model)
+ENDIF ELSE BEGIN
+  use_model='MC10'    ;Default is the Compiegne et al 2010 model
+ENDELSE
+
+; run the fortran via the wrapper once
+dustem_init,model=use_model
+st_model=dustem_read_all(!dustem_soft_dir) 
+dustem_write_all,st_model,!dustem_dat
+st=dustem_run()
+
+
+data_path = !dustem_wrap_soft_dir+'/Data/'
+fortran_path = !dustem_dat
+if keyword_set(dustem_data_dir) then fortran_path = dustem_data_dir
+if keyword_set(data_dir) then data_path = data_dir
+
+ if n_elements(NH) EQ 0 then nh = 1.d20 
+ if n_elements(INST) EQ 0 then begin 
+    inst = ['DIRBE','FIRAS','HFI','WMAP']
+ endif else begin 
+    inst = strupcase(inst)
+    inst = inst(UNIQ(inst, SORT(inst)))
+ endelse
+ n_inst = n_elements(inst)
+ if n_elements(WREF) EQ 0 then wref = 2.5d2
+ if n_elements(BBPAR) EQ 0 then BBPAR = [17.75d, 1.8d0] 
+ if n_elements(SW) EQ 0 then sw = 0
+ if n_elements(TIT) EQ 0 then tit = 'DustEM'
+ if n_elements(COM) EQ 0 then com = tit
+ if n_elements(WN) EQ 0 then wn = [0, 2, 1, 3, 6, 4, 5, 7, 8, 9, 10, 11, 12, 13] ; absolute display 
+ if n_elements( XR ) NE 2 then begin 
+    if keyword_set( COSMO ) then xr = [1., 1.e5] else xr = [ 1, 1.e3 ] 
+    endif else xr = [min(xr),max(xr)]
+ if n_elements( YR ) NE 2 then begin 
+    if keyword_set( COSMO ) then yr = [1.e-12, 1.e-3] else yr = [1e-8, 1e-4] 
+ endif else yr = [min(yr),max(yr)]
+ if n_elements( DATA ) EQ 0 then data = 0
+ if n_elements( SHOW ) EQ 0 THEN BEGIN 
+    show = ['emis'] 
+ ENDIF ELSE BEGIN
+    show = [STRLOWCASE(STRTRIM(STRING(show),2))]  
+    IF SHOW(0) NE '0' THEN BEGIN 
+       show = ['emis',show]     ; always show the SED 
+       show = show(UNIQ(show))
+    ENDIF 
+ ENDELSE
+ pgrid = [ 'emis', 'extuv', 'extir', 'alb', 'sdist','polext', 'polsed', 'align' ]
+ nplots = n_elements( SHOW ) ; nr of plots 
+ nhard = n_elements( HARD ) 
+ if nhard EQ 0 OR show(0) EQ '0' THEN hard=['0']
+ if nhard GT 0 then hard = [STRLOWCASE(STRTRIM(STRING(hard),2))]  
+ if  ((nhard GT 0) AND (nhard LT nplots)) OR (hard(0) EQ '1') then begin
+    ig = INTARR( nplots )
+    for i = 1, nplots-1 do begin
+       ig(i) = WHERE( pgrid EQ show(i))
+    endfor 
+    hard = pgrid(ig(sort(ig))) + REPLICATE('.ps',nplots)
+ endif
+ hard = strtrim(hard,2)
+ nhard = n_elements( HARD )
+
+
+; constants 
+ na = 6.02d23
+ clight = 2.9979246d10          ; cm/s
+
+;
+; plot inits
+;
+; set the rgb colors
+ red  =[0,1,1,0,0,1]
+ green=[0,1,0,1,0,1]
+ blue =[0,1,0,0,1,0]
+ tvlct,255*red, 255*green, 255*blue 
+ ls = [ [ 0,3,1,4,5 ], [ 0,3,1,4,5 ] ]
+ ihard = 0
+
+;
+; get the SED data
+;
+; New FarIR-mm SED from FBoulanger
+ RESTORE, data_path+'/EXAMPLE_OBSDATA/filters_ref_wave.xdr'
+ fact_em = 0.77                     ; to account for 20% of H to be in ionized form + 3% H2
+ to_sed_20 = 1d-17 * (1d20/1.82d18) ; MJy/sr->erg/s/cm2/sr and normalization to NH = 10^20 H/cm2
+ 
+; FIRAS
+ READCOL, data_path+'/EXAMPLE_OBSDATA/diffuse_ISM_SED.dat', wfirasf, firasf, ufirasf, skipline=5, numline=156, /sil
+ nufirasf = clight/(wfirasf*1.d-4)
+ firasf   = firasf * nufirasf * fact_em * to_sed_20 
+ ufirasf  = ufirasf * nufirasf * to_sed_20 ; error
+
+; DIRBE 60 --> 240 microns
+ READCOL, data_path+'/EXAMPLE_OBSDATA/diffuse_ISM_SED.dat', wdirbef, dirbef, udirbef, skipline=162, numline=4, /sil
+ nudirbef  = clight/ (wdirbef*1.d-4)
+ dirbef    = dirbef * nudirbef * fact_em * to_sed_20 
+ udirbef   = udirbef * nudirbef * to_sed_20 ; error
+ dwdirbef  =   dwdirbe[6:9]/2.
+
+; WMAP
+ READCOL,data_path+'/EXAMPLE_OBSDATA/diffuse_ISM_SED.dat', wwmapf, wmapf, uwmapf, skipline=167, numline=5,/sil
+ wwmapf = wwmapf[3:4]
+ wmapf  = wmapf[3:4]
+ uwmapf = uwmapf[3:4] 
+ nuwmapf = clight/(wwmapf*1.d-4)
+ wmapf   = wmapf * nuwmapf * fact_em * to_sed_20 
+ uwmapf  = uwmapf * nuwmapf * to_sed_20 ; error
+
+; DIRBE Arendt et al (1998) for |b|>25
+; numbers from Li & Draine, apj, 2001, 554, 778-802
+ Wave_ARENDT  = WDIRBE[2:*]
+ Dwave_ARENDT = DWDIRBE[2:*] / 2.
+ ARENDT       = [0.97, 1.11, 7.16, 3.57, 5.30, 18.6, 22.5, 10.1] ; 10^-26 erg/s/H/sr
+ ARENDT = ARENDT * 1.d-26 * 1.d20                                ; To erg s-1 cm-2 sr-1 for NH=10^20 H/cm2
+ err_ARENDT = 0.2 * ARENDT                                       ; from Dwek et al. 1997 DIRBE 1Sigma unc = 20%
+ nu_arendt = nudirbe[2:*] * 1.e9
+
+; Normalization of the Arendt spectrum on the Boulanger 100 microns
+ wave_arendt_midIR  = wave_arendt[0:3]
+ dwave_ARENDT_midir = Dwave_ARENDT[0:3]
+ nu_arendt_midIR    = nu_arendt[0:3]
+ correl_coeff_midIR  = [0.00183, 0.00291, 0.0462, 0.0480]               ; for Inu
+ ucorrel_coeff_midIR = [0.00001, 0.00003, 0.0001, 0.0002]               ; For Inu
+ arendt_midIR = correl_coeff_midIR * (dirbef[1]*100.*1e-4/clight*1e20)  ; in (Inu) MJy sr-1
+ err_arendt_midIR = FLTARR(N_ELEMENTS(arendt_midIR))
+ for i=0,3 do begin 
+    tmp= SQRT( ((udirbef[1]*100.*1e-4/clight*1e20)/(dirbef[1]*100.*1e-4/clight*1e20))^2. + (ucorrel_coeff_midIR[i]/correl_coeff_midIR[i])^2. )
+    err_arendt_midIR[i]  = arendt_midIR[i] * tmp
+ endfor
+ arendt_midIR = arendt_midIR / wave_arendt_midIR/1e-4*clight/1e20          ; to erg s-1 cm-2 sr-1 
+ err_arendt_midIR = err_arendt_midIR / wave_arendt_midIR/1e-4*clight/1e20  ; to erg s-1 cm-2 sr-1 
+ wdirbe_ful = [ wave_arendt_midir, wdirbef ]
+ dwdirbe_ful = [ dwave_arendt_midir, dwdirbef ]
+ dirbe_ful  = [ arendt_midir, dirbef ]
+ udirbe_ful = [ err_arendt_midir, udirbef ]
+
+; isocam galactic spectrum
+ RESTORE, data_path+'/EXAMPLE_OBSDATA/spectre_gal.xdr'                     ; wgal in microns, spec_gal in MJy/sr
+ nuisocam =   clight / wgal/1.d-4
+ norm = 0.65*0.045/filters(9)                         ; normalization for Nh=1e20 cm-2 and I12/I100 =0.045 
+ iso = 1d-17 * spec_gal*nuisocam * norm               ; galactic mid-IR SED in erg/s/cm2/sr
+ stars = 1d3 * wgal * DUSTEM_BLACKBODY(wgal, 3d3, unit='w')  ; assumed spectrum for stars in cgs
+ stars = 1.3d-6 * stars/stars(0)                      ; normalization of stellar spectrum
+ err_cvf = 0.2 
+ isocam   = iso-stars
+ smp = DUSTEM_GET_BAND_FLUX( data_path + 'FILTERS/', 'DIRBE', xs=wgal, ys=isocam )
+ isocam = isocam / smp.i_dirbe.ym(4) * arendt_midIR[2] ; factor is now 1.013027 (1.0235712 before)
+
+; Arome 3.3 microns: Giard et al 1988
+ xg = [3.3d]
+ yg = 1.1e-6 / 4./!pi * xg/0.05  ; SED erg/s/cm2/sr assumes a feature width of 0.05 mic
+ eg = 0.5e-6 / 4./!pi * xg/0.05  ; error
+
+; Modified BB overlaid to dust: default is 17.75 K and beta=1.8  (Miville-Deschênes et al 2013, Planck data)
+ np = 1000 
+ wbr = ALOG10( [20.,1.d5] )
+ dwbr = (wbr(1)-wbr(0)) / (np-1)
+ lambdaf = 10.^(wbr(0) + findgen(np)*dwbr)
+ temp1 = bbpar(0)
+ beta = bbpar(1)
+ lamb_ref = 250.
+ qf = (lambdaf/ lamb_ref)^(-beta)
+ bbm = 1d3 * lambdaf * DUSTEM_BLACKBODY(lambdaf, temp1, unit='w')
+ SP1 = bbm*qf
+ firas_nuinu = firasf 
+; normalize to FIRAS @ wave WREF < wfiras_max and get dust cross-section per H
+ wfiras_max = 8d2
+ if WREF LE wfiras_max then lamb_ref = wref 
+ tt = MIN( ABS(lambdaf-lamb_ref), ibb )
+ tt = MIN( ABS(wfirasf-lamb_ref), ifiras )
+ nw = 1
+ eps_firas = median(firasf(ifiras-nw:ifiras+nw)) / INTERPOL(sp1/qf, lambdaf, lamb_ref) / nh
+ SP1 = SP1 * median(firasf(ifiras-nw:ifiras+nw)) / INTERPOL(sp1, lambdaf, lamb_ref)
+
+; add HFI (Planck2011t, Dust in the diffuse interstellar medium and the Galactic halo)
+ nu_pep_dism = [5d3,3d3,857.,545.,353.] * 1d9 ; includes IRAS60 and 100
+ w_pep_dism = 1d4*clight / nu_pep_dism
+ dw_pep_dism = dblarr(n_elements(nu_pep_dism))
+ dw_pep_dism(0:1) = [31.,35.6] / 2.
+ dw_pep_dism(2:*) = w_pep_dism(2:*)/3. / 2.
+ y_pep_dism = [ 1.288, 6.522,5.624, 1.905, 0.465 ] * 1d-1  ; MJy/sr/1e20cm-2 
+ u_pep_dism = [ 0.436, 1.473, 1.015, 0.347, 0.100 ] * 1d-1
+ y_pep_dism = fact_em * y_pep_dism * nu_pep_dism * 1d-17 ; MJy/sr/1e20cm-2 to erg/s/cm2/sr
+ u_pep_dism = fact_em * u_pep_dism * nu_pep_dism * 1d-17
+ y_hfi_dism = y_pep_dism(0:2) & u_hfi_dism = u_pep_dism(0:2)
+
+; add cosmosomas
+ IF keyword_set( COSMO ) then begin 
+; Watson et al 2005
+;    nu=1.e9*[ 0.408, 1.42, 10.9, 12.7, 14.7, 16.3, 22.8, 33.0, 40.9, 61.3, 93.8, 1250, 2143., 3.e3 ]
+;    i_jy=[6.3,7.3,17.,18.7,26.2,32.6,42.3,40.3,33.9,34.7,77.5,9.68e4,1.31e5,6.44e4 ]
+;    e_jy=[7.8,2.,0.1,0.4,0.6,1.5,0.2,0.4,0.7,1.8,4.3,7d2,1.25d3,100.] 
+;    beam = 4.9e-4               ; steradians
+;    lref = 140d                 ; wave for normalization
+
+; nu in GHz and i_jy in Janskys
+; Planck 2011p, "New light on anomalous microwave emission from spinning dust grains"
+    nu=[0.408,0.82,1.42,10.9,12.7,14.7,16.3,22.8,28.5,33.0,40.9,44.1,61.3,70.3,93.8,100.,143.,217.,353.,545.,857.,1250,2143.,3.e3]
+    iplck = where( (nu eq 33.) or (nu eq 40.9) or (nu eq 70.3) or (nu eq 100.) or (nu eq 143.) or (nu eq 217.) )
+    nu=1d9*[0.408,0.82,1.42,10.9,12.7,14.7,16.3,22.8,28.5,33.0,40.9,44.1,61.3,70.3,93.8,100.,143.,217.,353.,545.,857.,1250,2143.,3.e3]
+    i_jy=[9.7,9.4,8.0,16.1,20.,28.4,35.8,39.8,38.1,36.7,30.8,28.5,27.4,32.2,63.,78.4,202.,1050.,3060.,1.53d4,4.87d4,9.3d4,1.17d5,5.36d4]
+    e_jy=[3.3,4.7,1.7,1.8,2.2,3.1,4.,4.2,4.,3.9,3.3,3.2,3.4,3.9,7.8,15.4,22.,128.,467.,2.09d3,6.11d3,1.29d4,1.45d4,6.67d3] 
+    fwhm = 1.81d
+    beam = !pi*(!pi/1.8d2 * fwhm/2d0)^2 ; steradians
+    lref = 140d                         ; wave for normalization
+
+    x_pep_ame = 1d4*clight/nu
+    nc = n_elements(x_pep_ame)
+    y_pep_ame = 1.e-23 * nu * i_jy / beam
+    e_pep_ame = 1.e-23 * nu * e_jy / beam
+    lref = 240d
+    iref = WHERE( wdirbe_ful EQ lref, cr)
+    if cr EQ 1 then begin
+       yic = INTERPOL(y_pep_ame,x_pep_ame,wdirbe_ful[iref])
+       y_norm = dirbe_ful[iref]/yic
+       y_norm = y_norm[0]
+    endif else y_norm = 1d
+    is = WHERE( x_pep_ame/lref GE 0.9, cs)
+    x_pep_ame=x_pep_ame(is) & y_pep_ame=y_pep_ame(is)*y_norm & e_pep_ame=e_pep_ame(is)*y_norm
+ ENDIF
+
+;
+; plot SED data
+;
+ IF SHOW(0) NE '0' THEN BEGIN
+    if HARD(0) NE '0' then begin 
+       !y.thick = 5
+       !x.thick = 5
+       !p.thick = 5
+       !p.charsize = 1.3
+       !p.charthick = 5
+       !x.ticklen = 0.04
+       set_plot,'ps'
+;       device,file=hard(0),xs=26, ys=20,/portrait,/color
+       device,file=hard(0),/portrait,/color
+       ihard = ihard + 1
+    endif else begin 
+       window, wn(0), xs=900,ys=600,tit='SHOW_DUSTEM: SED '+strtrim(wn(0),2)
+    endelse
+
+    xtit = 'Wavelength (!4l!3m)'
+    ytit = '!4m!3 I!d!4m!3!n (erg s!u-1!n cm!u-2!n sr!u-1!n)'
+    fine = 1
+;    plot_oo, INDGEN(1),/NODAT,/xs,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit,tit=tit 
+    plot_oo, INDGEN(1),/NODAT,xs=9,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit
+    axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)'
+    if DATA NE 0 then begin 
+       if HARD(0) NE '0' then begin 
+          oploterror,xg,yg,0.,eg,psym=5
+          oploterror, wfirasf, firasf, wfirasf*0d, ufirasf, errthick=0.7, /nohat
+          oploterror, wdirbe_ful(0:3), dirbe_ful(0:3), dwdirbe_ful(0:3), udirbe_ful(0:3), psym=6
+          oploterror, wdirbe_ful(6:*), dirbe_ful(6:*), dwdirbe_ful(6:*), udirbe_ful(6:*), psym=6
+;          oploterror, wdirbe_ful, dirbe_ful, dwdirbe_ful, udirbe_ful, psym=6
+          oploterror, wgal,isocam, wgal*0d, isocam*err_cvf, errthick=0.7, /nohat
+          oploterror, wwmapf, wmapf, wwmapf*0d, uwmapf, psym=4
+          oploterror, w_pep_dism, y_pep_dism, dw_pep_dism, u_pep_dism, psym=5
+       endif else begin
+          oploterror,xg,yg,0.,eg,psym=5
+          oploterror, wfirasf, firasf, wfirasf*0d, ufirasf, errthick=0.4, /nohat
+          oploterror, wdirbe_ful(0:3), dirbe_ful(0:3), dwdirbe_ful(0:3), udirbe_ful(0:3), psym=6
+          oploterror, wdirbe_ful(6:*), dirbe_ful(6:*), dwdirbe_ful(6:*), udirbe_ful(6:*), psym=6
+;          oploterror, wdirbe_ful, dirbe_ful, dwdirbe_ful, udirbe_ful, psym=6
+          oploterror, wgal,isocam, wgal*0d, isocam*err_cvf, errthick=0.4, /nohat
+          oploterror, wwmapf, wmapf, wwmapf*0d, uwmapf, psym=4
+          oploterror, w_pep_dism, y_pep_dism, dw_pep_dism, u_pep_dism, psym=5
+       endelse
+    endif
+    if DATA EQ 1 then begin
+       OPLOT, LAMBDAf, SP1, lines=1 ; FIR blackbody
+       sp1_firas = INTERPOL( sp1, lambdaf, wfirasf )
+       firas_chi2 = DUSTEM_CHI2( firasf, sp1_firas, 3, err=ufirasf )
+       print,'Grey body: Td = ',strtrim(bbpar(0),2),' and beta = ',strtrim(bbpar(1),2)
+       print,'Chi-squared of grey body to FIRAS: ', strtrim(firas_chi2,2)
+    endif
+    if keyword_set( CMB ) then begin 
+       lambda_cmb = 1.d2^( 1.d0 + dindgen(500)*0.01 )
+       sp_cmb = 1d3 * lambda_cmb * DUSTEM_BLACKBODY(LAMBDA_CMB, 2.728, unit='w') ; SED in erg/s/cm2/s/sr
+       OPLOT, LAMBDA_CMB, SP_CMB, lines=0                                 ; CMB
+    endif
+    if keyword_set( COSMO ) AND DATA NE 0 then begin 
+       oploterror, x_pep_ame, y_pep_ame, y_pep_ame*0d, e_pep_ame, ps=4
+;       oploterror, x_pep_ame(iplck), y_pep_ame(iplck), y_pep_ame(iplck)*0d, e_pep_ame(iplck), ps=5, syms=1.5, col=3
+    endif
+ ENDIF
+
+;
+; get the model SED
+;
+ fname = fortran_path + 'data/GRAIN.DAT'
+ nlines = FILE_LINES( fname )
+ OPENR, uu, fname, /get_lun
+ tmp = '#'
+ print,'(W) DUSTEM_SHOW_FORTRAN: GRAIN.DAT'
+ cnt = 0
+ WHILE STRPOS(tmp,'#') EQ 0 do begin
+    READF, uu, tmp
+    cnt = cnt + 1
+ ENDWHILE
+ r_opt = STRLOWCASE(STRTRIM(STRING(tmp),2))
+ print,r_opt
+ READF, uu, g0
+ ntype = nlines - (cnt+1)
+ nsize=intarr(ntype)
+ t_opt=strarr(ntype)
+ gtype = strarr(ntype)
+ propm = dblarr(ntype)
+ rho = dblarr(ntype)
+ for i=0,ntype-1 do begin 
+    READF,uu,tmp
+    PRINT, tmp
+    tt = strsplit(tmp, /extract)
+    gtype(i) = strtrim(tt(0))
+    nsize(i)=fix(tt(1))
+    t_opt(i)=strlowcase(strtrim(tt(2)))
+    propm(i) = double(tt(3))
+    rho(i) = double(tt(4))
+ endfor
+ close,uu
+ free_lun,uu
+ nsz_max = MAX(nsize)
+
+ IF n_elements(FSED) EQ 0 THEN fsed = fortran_path+'out/SED.RES'
+ OPENR, uu, fsed, /get_lun
+ tmp = '#'
+ WHILE (STRPOS(tmp,'#') EQ 0) do begin
+    READF, uu, tmp
+ ENDWHILE
+ tt = double( strsplit(tmp, ' ', /extract) )
+ ntype_sed = fix(tt(0))
+ if ntype_sed NE ntype then begin
+    print,'(F) DUSTEM_SHOW_FORTRAN: SED.RES & GRAIN.DAT have different NTYPE'
+    print,'                 data is not from present GRAIN.DAT'
+    return
+ endif
+ nlamb = fix(tt(1)) ; nr of wavelengths
+ x = dblarr(nlamb)
+ sedh = dblarr(nlamb,ntype+1)
+ for i=0,nlamb-1 do begin 
+    READF, uu, tmp
+    tt = double( strsplit(tmp, ' ', /extract) )
+    x(i) = tt(0) 
+    sedh(i,*) = tt(1:*)
+ endfor
+ close,uu
+ free_lun,uu
+ xlamb= x
+ sed = sedh * nh / 4. / !pi     ; in erg/s/cm2/sr
+ if keyword_set( SW ) then sed = SMOOTH( sed, sw)
+
+ if (DATA NE 1) AND (SHOW(0) NE '0') then begin 
+    dy = 10.^((ALOG10(yr(1))-ALOG10(yr(0))) / 25.)
+    dx = 10.^((ALOG10(xr(1))-ALOG10(xr(0))) / 50.)
+    yps= 10.^((ALOG10(yr(1))+ALOG10(yr(0))) / 2. + (ALOG10(yr(1))-ALOG10(yr(0))) / 2.3)
+    if keyword_set(COSMO) then begin
+       xpr=10.^((ALOG10(xr(1))+ALOG10(xr(0)))/3.) 
+       yps= 10.^((ALOG10(yr(1))+ALOG10(yr(0))) / 2. - (ALOG10(yr(1))-ALOG10(yr(0))) / 4.)
+    endif else xpr=xr(0)
+    xpr = xpr*[dx,dx^4]
+; overlay model SED in erg/s/cm2/sr
+    for i=0,ntype-1 do begin
+       oplot, x, sed(*,i), lin=ls(i+1)
+       oplot,xpr,yps*[1.,1.], lin=ls(i+1)
+       xyouts,/data,xpr(1)*1.1,yps,gtype(i),chars=1
+       yps = yps/dy
+    endfor
+    oplot, x, sed(*,ntype),lin=ls(0),col=2
+    s1 = STRMID(STRTRIM(ROUND(1e1*g0)/1e1,2), 0, 4) 
+    s2 = STRMID(STRTRIM(ROUND(1e2*alog10(nh))/1e2,2), 0, 5) 
+    xyouts,/norm,0.6,0.85,'G!d0!n='+s1+' N!dH!n=10!u'+s2+'!ncm!u-2!n'
+; overlay model polarized SED in erg/s/cm2/sr
+    ;IF (C_POLSED EQ 1) THEN oplot, x, sed_p(*,ntype),lin=ls(0),col=3
+ endif
+
+;
+; fill in SMDAT
+; 
+; first get model fluxes in instrument bands
+ unit = 'x(microns) SED(erg/s/cm2/sr)'
+ if n_elements(SMDAT) GT 0 then begin
+    smdat = DUSTEM_GET_BAND_FLUX( data_path + 'FILTERS/', inst, xs=x, ys=sed, unit=unit, smi=smdat )
+ endif else smdat = DUSTEM_GET_BAND_FLUX( data_path + 'FILTERS/', inst, xs=x, ys=sed, unit=unit )
+ smdat.com = com
+ stag = TAG_NAMES(smdat)
+ ;if c_polsed EQ 1 then smdat.m_emis.yp = sed_p
+
+; then get diffuse data available here
+ ii = WHERE( inst EQ 'DIRBE', cic)
+ if cic GT 0 then begin
+    smdat.i_dirbe.unit = smdat.i_dirbe.unit + ' YD(erg/s/cm2/sr)'
+    smdat.i_dirbe.yd = [ 0d, 0d, dirbe_ful ]
+    smdat.i_dirbe.err = [ 0d, 0d, udirbe_ful]
+    nband = n_elements( smdat.i_dirbe.x)
+    smdat.i_dirbe.isel = [ 0,0,intarr(nband-2)+1 ]
+    smdat.i_dirbe.npar = ntype
+ endif
+ 
+ ii = WHERE( inst EQ 'HFI', cic)
+ if cic GT 0 then begin 
+    smdat.i_hfi.unit = smdat.i_hfi.unit + ' YD(erg/s/cm2/sr)'
+    ibx = indgen(3)             ; only 3 first bands in PEP DISM
+    smdat.i_hfi.yd(ibx) = y_hfi_dism
+    smdat.i_hfi.err(ibx) = u_hfi_dism
+    smdat.i_hfi.isel = [1,1,1,0,0,0]
+    smdat.i_hfi.npar = 1
+ endif
+    
+ ii = WHERE( inst EQ 'WMAP', cic)
+ if cic GT 0 then begin  
+    smdat.i_wmap.unit = smdat.i_wmap.unit + ' YD(erg/s/cm2/sr)'
+    smdat.i_wmap.yd = [ 0d, 0d, 0d, wmapf ] ; only 61 and 94 GHz bands
+    smdat.i_wmap.err = [ 0d, 0d, 0d, uwmapf]
+    nband = n_elements( smdat.i_wmap.x)
+    smdat.i_wmap.isel = [ 0,0,0,intarr(nband-3)+1 ]
+    smdat.i_wmap.npar = 1
+ endif
+ 
+; put FIRAS pts in SMDAT
+ ii = WHERE( inst EQ 'FIRAS', cic)
+ if cic GT 0 then begin  
+    nband = n_elements(wfirasf)
+    smdat = DUSTEM_ADD_INST( smdat, 'FIRAS', [n_elements(firasf),ntype+1] )
+    smdat.i_firas.unit = 'x(microns) YM(erg/s/cm2/sr) FLX(MJy/sr) YD(erg/s/cm2/sr)'
+    smdat.i_firas.yd = firasf
+    smdat.i_firas.err = ufirasf
+    sd_firas = dblarr(n_elements(wfirasf),ntype+1)
+    for itp=0,ntype do sd_firas(*,itp) = INTERPOL( sed(*,itp), x, wfirasf )
+    smdat.i_firas.ym = sd_firas
+    smdat.i_firas.isel = intarr(nband)+1
+    smdat.i_firas.npar = 2
+ endif
+    
+; overlay model band fluxes
+ if (DATA NE 1) AND (SHOW(0) NE '0') then begin 
+    itg = WHERE( STRPOS(stag,'I_') GE 0 AND stag NE 'I_FIRAS', ctg )
+; then other instruments
+    if ctg GT 0 then begin
+       for k=0,ctg-1 do begin
+          oplot, smdat.(itg(k)).x, smdat.(itg(k)).ym(*,ntype), ps=6, syms=1.5, col=3
+       endfor
+    endif
+ endif
+
+;
+; get extinction
+;
+ OPENR, uu, fortran_path+'out/EXT.RES', /get_lun
+ tmp = '#'
+ WHILE STRPOS(tmp,'#') EQ 0 do begin
+    READF, uu, tmp
+ ENDWHILE
+ tt = double( strsplit(tmp, ' ', /extract) )
+ ntype_ext = fix(tt(0))
+ if ntype_ext NE ntype then begin
+    print,'(F) DUSTEM_SHOW_FORTRAN: EXT.RES & GRAIN.DAT have different NTYPE'
+    print,'                 data is not from present GRAIN.DAT'
+    return
+ endif
+ nlamb = fix(tt(1))             ; nr of wavelengths
+ x = dblarr(nlamb)
+ stmp = dblarr(nlamb,2*ntype+1)
+ ssca = dblarr(nlamb,ntype+1)
+ sabs = dblarr(nlamb,ntype+1)
+ sext = dblarr(nlamb,ntype+1)
+ alb = dblarr(nlamb,ntype+1)
+ for i=0,nlamb-1 do begin 
+    READF, uu, tmp
+    tt = double( strsplit(tmp, ' ', /extract) )
+    x(i) = tt(0)
+    stmp(i,*) = tt(1:*)
+ endfor
+ CLOSE,uu
+ FREE_LUN,uu
+ wlm = x
+ x = 1. / wlm                   ; 1/micron
+ 
+ mH = 1.67262158d-24            ; proton mass /gdust -> /gH 
+ pm = [propm,propm]
+ for i=0,ntype-1 do begin 
+    sabs(*,i) = stmp(*,i) 
+    ssca(*,i) = stmp(*,ntype+i) 
+    sext(*,i) = sabs(*,i) + ssca(*,i)
+    alb(*,i) = ssca(*,i) / sext(*,i)
+ endfor
+ for i=0,nlamb-1 do begin 
+    sabs(i,ntype) = TOTAL(sabs(i,0:ntype-1))
+    ssca(i,ntype) = TOTAL(ssca(i,0:ntype-1))  
+    sext(i,ntype) = TOTAL(sext(i,0:ntype-1))  
+    alb(i,ntype) = ssca(i,ntype)/sext(i,ntype)
+ endfor
+
+; get RV
+ av = INTERPOL( sext(*,ntype), wlm, 0.55 )
+ ab = INTERPOL( sext(*,ntype), wlm, 0.44 )
+ rv = av / (ab-av)
+
+; get 250 microns emissivity
+ eps_a = dblarr(ntype+1) & eps_e = eps_a & eps_s = eps_a
+ for i = 0,ntype-1 do begin 
+    eps_a(i)=INTERPOL(sabs(*,i), wlm, wref)
+    eps_s(i)=INTERPOL(ssca(*,i), wlm, wref)
+    eps_e(i)=INTERPOL(sext(*,i), wlm, wref)
+ endfor
+ eps_a(ntype) = INTERPOL(sabs(*,ntype), wlm, wref)
+ eps_s(ntype) = INTERPOL(ssca(*,ntype), wlm, wref)
+ eps_e(ntype) = INTERPOL(sext(*,ntype), wlm, wref)
+ print,'(W) DUSTEM_SHOW_FORTRAN: model dust cross-section @ ',STRTRIM(wref,2),' per type and total (cm2/H)', $
+       format='(A44,1E10.3,A27)'
+ print, ' Abs : ',eps_a, format='(A8,100(1E12.4))'
+ print, ' Sca : ',eps_s, format='(A8,100(1E12.4))'
+ print, ' Ext : ',eps_e, format='(A8,100(1E12.4))'
+ if WREF GT wfiras_max then begin
+    print,'(W) DUSTEM_SHOW_FORTRAN: WREF above longer FIRAS wave, using 250 microns'
+    wp = lamb_ref
+ endif else wp = wref
+ print,'(W) DUSTEM_SHOW_FORTRAN: dust cross-section @ ', STRTRIM(wp,2), ' microns from FIRAS',format='(A38,1E10.3,A20)'
+ print,'                 using T = ',bbpar(0), ' and beta = ',bbpar(1),' sigma(cm2/H) = ', eps_firas, $
+       format='(A27,1E8.2,A13,1E8.2,A17,1E12.4)'
+
+; fill in model
+ it = WHERE( stag EQ 'M_EXT', ct)
+ if ct EQ 0 then begin
+    smdat = DUSTEM_ADD_MOD( smdat, 'EXT', [n_elements(wlm),ntype+1])
+    smdat.m_ext.x = wlm
+    smdat.m_ext.y = sext
+    smdat.m_ext.abs = sabs
+    smdat.m_ext.sca = ssca
+    smdat.m_ext.alb = alb
+    smdat.m_ext.xr = wref
+    smdat.m_ext.yr_abs = eps_a
+    smdat.m_ext.yr_sca = eps_s
+    smdat.m_ext.rv = rv
+ endif
+ 
+; get standard Savage & Mathis 79 + Mathis 1990 extinction and fill in data
+ READCOL, data_path+'/EXAMPLE_OBSDATA/mean_ism_ext.dat', skip=3, x_sm, e_sm, /SIL
+ it = WHERE( stag EQ 'I_SM79', ct)
+ if ct EQ 0 then begin
+    smdat = DUSTEM_ADD_INST( smdat, 'SM79', [n_elements(x_sm),ntype+1] )
+    smdat.i_sm79.x = x_sm
+    smdat.i_sm79.ym = INTERPOL( sext(*,ntype), wlm, x_sm)
+    smdat.i_sm79.yd = e_sm
+    smdat.i_sm79.err = smdat.i_sm79.yd*0.2
+    smdat.i_sm79.unit = 'x(microns) sigma(cm2/H)'
+    smdat.i_sm79.npar = ntype
+    smdat.i_sm79.isel = intarr(n_elements(x_sm)) + 1
+ endif
+ 
+;
+; plot the vis-UV extinction
+;
+ yscl = 1d0
+ ip = WHERE( STRPOS(show,'extuv') GE 0, cp )
+ IF (SHOW(0) NE '0') AND (CP EQ 1) THEN BEGIN
+    if HARD(0) NE '0' then begin 
+       device, file=hard(ihard), /portrait,/color
+       ihard = ihard + 1
+    endif else begin 
+       window, wn(1), xs=600,ys=400, tit='DUSTEM_SHOW_FORTRAN: Extinction UV '+strtrim(wn(1),2)
+    endelse 
+    xtit = 'x (!4l!3m!u-1!n)'
+    xr = [0,11] 
+    if n_elements(WEXT) then begin 
+       xe = [x,1./wext]
+       xe = xe(SORT(XE))
+       ei = INTERPOL(sext(*,ntype),x,xe)
+       tmp = MIN( ABS(xe - 1./wext),iw  ) 
+       yr = minmax( sext(*,ntype)/ei(iw) ) * [1.,1.]
+       yscl = 1.d0 / ei(iw)
+       ytit='Normalized !4r!3!dext!n '
+       print,'(W) DUSTEM_SHOW_FORTRAN: extinction normalized to 1 at ',wext, ' microns yscl=',yscl
+    endif else begin
+       yr=[ 0, 2.5 ]
+       yscl = 1d0
+       ytit='!4r!3!dext!n (10!u-21!n cm!u2!n per H)'
+    endelse
+    plot, [1d], [1d], lin=ls(0), xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit 
+    oplot, x, sext(*,ntype)*yscl, lin=ls(0), col=2
+    oplot, 1d/x_sm, e_sm*1d21*yscl, ps=4
+;    oploterror, 1d/x_sm, e_sm*1d21, smdat.i_sm79.err*1d21, psym = 4, /nohat
+    for i=0,ntype-1 do oplot, x, (sext(*,i))*yscl, lin=ls(i+1)
+    xyouts,/norm,0.2,0.8,'R!dV!n='+STRMID(STRTRIM(ROUND(1e1*rv)/1e1,2), 0, 4)
+ ENDIF
+
+; plot IR extinction
+ ip = WHERE( STRPOS(show,'extir') GE 0, cp )
+ IF (SHOW(0) NE '0') AND (CP EQ 1) THEN BEGIN
+    if HARD(0) NE '0' then begin 
+       device, file=hard(ihard), /portrait,/color
+       ihard = ihard + 1
+    endif else begin 
+       window, wn(2), xs=600,ys=400, tit='DUSTEM_SHOW_FORTRAN: Extinction IR '+strtrim(wn(2),2)
+    endelse 
+    xtit = textoidl('\lambda (\mum)')
+    xr = [1,400] 
+    yr = [1e-6,1.]
+    if yscl NE 1d0 THEN ytit='Normalized !4r!3!dext!n' ELSE $
+       ytit='!4r!3!dext!n (10!u-21!n cm!u2!n per H)'
+    plot_oo, [1d], [1d], xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit 
+    oplot, wlm, sext(*,ntype)*yscl, lin=ls(0), col=2
+    oplot, x_sm, e_sm*1d21*yscl, ps=4
+;    oploterror, x_sm, e_sm*1d21, smdat.i_sm79.err*1d21, psym = 4, /nohat
+    for i=0,ntype-1 do oplot, wlm, sext(*,i)*yscl, lin=ls(i+1)
+    xyouts,/norm,0.8,0.85,'R!dV!n='+STRMID(STRTRIM(ROUND(1e1*rv)/1e1,2), 0, 4)
+ ENDIF
+
+;
+; and the albedo
+;
+ ip = WHERE( STRPOS(show,'alb') GE 0, cp )
+ IF (SHOW(0) NE '0') AND (CP EQ 1) THEN BEGIN
+    if HARD(0) NE '0' then begin 
+       device, file=hard(ihard), /portrait,/color
+       ihard = ihard + 1
+    endif else begin 
+       window, wn(3), xs=500,ys=350, tit='DUSTEM_SHOW_FORTRAN: Albedo '+strtrim(wn(3),2)
+    endelse 
+    xtit = 'x (!4l!3m!u-1!n)'
+    xr = [0,11] 
+    ytit='Albedo'
+    yr=[ 0, 1]
+    plot, x, alb(*,ntype), xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit 
+    oplot, x, alb(*,ntype), col=2
+    for i=0,ntype-1 do oplot, x, alb(*,i), lin=ls(i+1)
+ ENDIF
+
+
+;
+;  get the size distribution and plot
+;
+ ip = WHERE( STRPOS(show,'sdist') GE 0, cp )
+ ir = WHERE( STRPOS(r_opt,'sdist') GE 0, cr )
+ IF CR EQ 0 THEN BEGIN 
+    cp = 0
+    print,'(W) DUSTEM_SHOW_FORTRAN: SDIST keyword not set in current run'
+ ENDIF ELSE IF ((SHOW(0) NE '0') AND (CP EQ 1)) THEN BEGIN 
+    fn = fortran_path+'out/SDIST.RES'
+    OPENR, uu, fn, /get_lun
+    tt = '#'
+    WHILE STRPOS(tt,'#') EQ 0 do begin
+       READF, uu, tt
+    ENDWHILE
+    ax = dblarr(nsz_max,ntype)
+    ava = dblarr(nsz_max,ntype)
+    nsz_tot = TOTAL(nsize)
+    ava_tot = dblarr(nsz_tot,ntype+1)
+    ax_tot = 0d
+    FOR i=0,ntype-1 do begin 
+       READF,uu,tt
+       for is=0,nsize(i)-1 do begin 
+          READF, uu, tx, ty
+          ax(is,i) = tx
+          ava(is,i) = ty
+       endfor
+       ax_tot =  [ax_tot, ax(0:nsize(i)-1,i)]
+    ENDFOR
+    close,uu
+    free_lun,uu
+    ax_tot = ax_tot(1:*)
+    ax_tot = ax_tot(SORT(ax_tot))
+    FOR i=0,ntype-1 do begin 
+       tt = INTERPOL( ava(0:nsize(i)-1,i), ax(0:nsize(i)-1,i), ax_tot )
+       ix = WHERE( ax_tot LT MIN(ax(*,i)) OR ax_tot GT MAX(ax(*,i)), cx )
+       if cx GT 0 then tt(ix) = 0d ; no extrapolation
+       ava_tot(*,i) = tt
+    endfor
+    for i=0, nsz_tot-1 do ava_tot(i,ntype) = TOTAL( ava_tot(i,0:ntype-1))
+
+; fill in model
+    it = WHERE( stag EQ 'M_SDIST', ct)
+    if ct EQ 0 then begin
+       smdat = DUSTEM_ADD_MOD( smdat, 'SDIST', [nsz_tot,ntype+1,nsz_max])
+       smdat.m_sdist.xtot = ax_tot
+       smdat.m_sdist.ytot = ava_tot
+       smdat.m_sdist.xi = ax
+       smdat.m_sdist.yi = ava
+    endif
+    
+    if HARD(0) NE '0' then begin 
+       device, file=hard(nplots-1), /portrait,/color
+    endif else begin 
+       window, wn(5), xs=500,ys=350, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: Size Distribution '+strtrim(wn(5),2)
+    endelse 
+    yscl = 1.D29
+    xr=[0.1,1e4]
+    yr=[-3,3]
+    xtit='a (nm)'
+    ytit='log( 10!u29!n n!dH!u-1!n a!u4!n dn/da (cm!u3!n/H)) ' 
+    tit='DustEM'
+    plot_oi, ax(0:nsize(0)-1,0)*1.e7, alog10(ava(0:nsize(0)-1,0)*yscl), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit='DUSTEM'
+    FOR i=1,ntype-1 DO begin 
+       oplot, ax(0:nsize(i)-1,i)*1.e7, alog10(ava(0:nsize(i)-1,i)*yscl), line=ls(i+1)
+    ENDFOR
+
+ ENDIF
+
+
+; Display polarized SED
+;
+ ip = WHERE( STRPOS(show,'polsed') GE 0, c_polsed )
+ IF (SHOW(0) NE '0') AND (C_POLsed EQ 1) THEN BEGIN 
+    if HARD(0) NE '0' then begin 
+       device, file=hard(ihard), /portrait,/color
+       ihard = ihard + 1
+    endif else begin 
+       window, wn(6), xs=600,ys=400, tit='DUSTEM_SHOW_FORTRAN: Polarized SED '+strtrim(wn(6),2)
+    endelse 
+
+ ; get polarized SED
+    OPENR, uu, fortran_path+'out/SED_POL.RES', /get_lun
+    tmp = '#'
+    WHILE (STRPOS(tmp,'#') EQ 0) do begin
+       READF, uu, tmp
+    ENDWHILE
+    tt = double( strsplit(tmp, ' ', /extract) )
+    sedh_p = dblarr(nlamb,ntype+1)
+    for i=0,nlamb-1 do begin 
+       readf, uu, tmp
+       tt = double( strsplit(tmp, ' ', /extract) )
+       sedh_p(i,*) = tt(1:*)
+       x(i) = tt(0)
+    endfor
+    close,uu
+    free_lun,uu
+    sed_p = sedh_p * nh / 4. / !pi ; in erg/s/cm2/sr
+
+    yr = [1e-10,1e-5]
+    xr = [10,1e4]
+    dy = 10.^((ALOG10(yr(1))-ALOG10(yr(0))) / 25.)
+    dx = 10.^((ALOG10(xr(1))-ALOG10(xr(0))) / 50.)
+    yps= 10.^((ALOG10(yr(1))+ALOG10(yr(0))) / 2. + (ALOG10(yr(1))-ALOG10(yr(0))) / 2.3)
+    xpr=xr(0)
+    xpr = xpr*[dx,dx^4]
+
+    xtit = 'Wavelength (!4l!3m)'
+    ytit = '!4m!3 P!d!4m!3!n (erg s!u-1!n cm!u-2!n sr!u-1!n)'
+    fine = 1
+    plot_oo, INDGEN(1),/NODAT,xs=9,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit
+    axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)'
+    for i=0,ntype-1 do begin
+       oplot, x, sed_p(*,i), lin=ls(i+1)
+       oplot,xpr,yps*[1.,1.], lin=ls(i+1)
+       xyouts,/data,xpr(1)*1.1,yps,gtype(i),chars=1
+       yps = yps/dy
+    endfor
+    oplot,x,sed_p(*,ntype),lin=ls(0),col=3
+ ENDIF
+
+
+;
+; get the polarisation extinction
+;
+ ip = WHERE( STRPOS(show,'polext') GE 0, c_polext )
+ IF (SHOW(0) NE '0') AND (C_POLEXT EQ 1) THEN BEGIN 
+
+    itp = WHERE( STRPOS(t_opt,'pol') GE 0, ctp)
+    IF ctp EQ 0 THEN BEGIN 
+       c_pol = 0
+       print,'(W) DUSTEM_SHOW_FORTRAN: no POL data in current run'
+    ENDIF ELSE BEGIN
+       OPENR, uu, fortran_path+'out/EXT_POL.RES', /get_lun
+       tmp = '#'
+       WHILE STRPOS(tmp,'#') EQ 0 do begin
+          READF, uu, tmp
+       ENDWHILE
+       tt = double( strsplit(tmp, ' ', /extract) )
+       ntype_pol = fix(tt(0))
+       if ntype_pol NE ntype then begin
+          print,'(F) DUSTEM_SHOW_FORTRAN: POL.RES & GRAIN.DAT have different NTYPE'
+          print,'                 data is not from present GRAIN.DAT'
+          return
+       endif
+       nlamb = fix(tt(1))       ; nr of wavelengths
+       x = dblarr(nlamb)
+       spabs = dblarr(nlamb,ntype+1)
+       spsca = dblarr(nlamb,ntype+1)
+       for i=0,nlamb-1 do begin 
+          READF, uu, tmp
+          tt = double( strsplit(tmp, ' ', /extract) )
+          x(i) = tt(0)
+          spabs(i,0:ntype-1) = tt(1:ntype) 
+          spsca(i,0:ntype-1) = tt(ntype+1:2*ntype) 
+          spabs(i,ntype) = TOTAL(spabs(i,0:ntype-1))
+          spsca(i,ntype) = TOTAL(spsca(i,0:ntype-1))
+       endfor
+       close,uu
+       free_lun,uu
+
+; fill in model
+       smdat.m_ext.abs_p = spabs
+       smdat.m_ext.sca_p = spsca
+
+       if HARD(0) NE '0' then begin 
+          device, file=hard(ihard), /portrait,/color
+          ihard = ihard + 1
+       endif else begin 
+          window, wn(7), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: Serkowski '+strtrim(wn(7),2)
+       endelse 
+
+; plot Serkowski       
+       yscl = 1.D23
+       yscl2 = 1.D2
+       xr=[0.2,10.]
+       yr=[0.1,5]
+       xtit = textoidl('1/\lambda (\mum^{-1})')
+       ytit = textoidl('\sigma_{pol} (10^{-23}cm^2/H)')
+       tit='DUSTEM'
+       plot_oo, 1./x, yscl2*(spabs(*,0)+spsca(*,0)), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit
+       FOR i=1,ntype-1 DO oplot, 1./x, yscl2*(spabs(*,i)+spsca(*,i)), line=ls(i+1)
+       oplot, 1./x, yscl2*(spabs(*,ntype)+spsca(*,ntype)), col=3
+       ix = WHERE(x_sm LE 7d, csm) ; keep wave below 7 microns
+       if csm GT 0 then begin 
+          xx = 1./x_sm(ix)
+          xe = e_sm(ix) 
+       endif
+       ix = SORT(xx)
+       xx = xx(ix) & xe = xe(ix)
+       fpol = DUSTEM_SERKOWSKI(xx)
+       oplot,xx,yscl*fpol/5.8e21*3.1,ps=4
+;       oplot,xx,25*fpol*3.1/5.8e21/xe, ps=5
+    ENDELSE
+
+  ENDIF
+  
+;
+;  Polarization fraction in emission
+;
+ IF (SHOW(0) NE '0') AND (C_POLSED EQ 1) THEN BEGIN 
+; plot fractional polarisation       
+       if HARD(0) NE '0' then begin 
+          device, file=hard(ihard), /portrait,/color
+          ihard = ihard + 1
+       endif else begin 
+          window, wn(8), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: P/I '+strtrim(wn(8),2)
+       endelse 
+       xr=[10,1.e4]
+       yr=[0.,0.4]
+       xtit = textoidl('\lambda (\mum)')
+       ytit = textoidl('P/I')
+       tit=''
+       plot_oi, x, sed_p(*,0)/sed(*,0), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit
+       axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)'
+       FOR i=1,ntype-1 DO begin 
+          oplot, x, sed_p(*,i)/sed(*,i), line=ls(i+1)
+       ENDFOR
+       oplot, x, sed_p(*,ntype)/sed(*,ntype),col=3
+ ENDIF
+
+;
+;  Polarization fraction in extinction
+;
+ IF (SHOW(0) NE '0') AND (C_POLEXT EQ 1) THEN BEGIN 
+; plot fractional polarisation       
+       if HARD(0) NE '0' then begin 
+          device, file=hard(ihard), /portrait,/color
+          ihard = ihard + 1
+       endif else begin 
+          window, wn(9), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: p/tau '+strtrim(wn(9),2)
+       endelse 
+       yscl = 1.
+       xr=[1.e-1,1.e4]
+       yr=[0.,0.5]
+       xtit = textoidl('\lambda (\mum)')
+       ytit = textoidl('p/\tau')
+       tit=''
+       plot_oi, x, yscl*(spabs(*,0)+spsca(*,0))/sext(*,0), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit
+       FOR i=1,ntype-1 DO begin 
+          oplot, x, yscl*(spabs(*,i)+spsca(*,i))/sext(*,i), line=ls(i+1)
+       ENDFOR
+       oplot, x, yscl*(spabs(*,ntype)+spsca(*,ntype))/sext(*,ntype), col=3
+ ENDIF
+
+;
+; Alignment function
+;
+ ip = WHERE( STRPOS(show,'align') GE 0, c_align )
+ IF (SHOW(0) NE '0') AND (C_ALIGN EQ 1) THEN BEGIN 
+    fname = fortran_path + 'data/ALIGN.DAT'
+    nlines = FILE_LINES( fname )
+    OPENR, uu, fname, /get_lun
+    tmp = '#'
+    print,'(W) DUSTEM_SHOW_FORTRAN: GRAIN.DAT'
+    cnt = 0
+    WHILE STRPOS(tmp,'#') EQ 0 do begin
+       READF, uu, tmp
+       cnt = cnt + 1
+    ENDWHILE
+    flags = tmp
+    readf, uu, anisG0
+   
+    if HARD(0) NE '0' then begin 
+          device, file=hard(ihard), /portrait,/color
+          ihard = ihard + 1
+    endif else begin 
+          window, wn(10), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: f_align '+strtrim(wn(10),2)
+    endelse 
+
+    ; Parametric 
+    readf, uu, tmp
+    tt = strsplit(tmp, /extract)
+    align_type = tt(0)
+    if (align_type eq 'par') then begin
+       athresh = tt(1)
+       pstiff = tt(2)
+       plev = tt(3)
+       n = 100
+       ; Grain radius in microns
+       radmin = 1d-3
+       radmax = 10
+       radius = radmin * exp(indgen(n)*alog(radmax/radmin)/n)
+       fpol = 0.5 * plev * (1 + TANH(ALOG(radius/athresh) / pstiff ) )
+       xr=[radmin,radmax]
+       yr = [0,1.05]
+       xtit=textoidl('Grain radius (\mum)')
+       ytit='Alignment efficiency'
+       plot_oi,radius,fpol,xr=xr,yr=yr,/ys,xtit=xtit,ytit=ytit
+    endif
+       
+ ENDIF
+
+
+; get all the chi2
+ if n_elements(smdat) EQ 0 then npar = ntype
+ smdat = DUSTEM_FIL_CHI2( smdat,ntype=ntype+1 ) 
+ stag = TAG_NAMES(smdat)
+ itx = WHERE( STRPOS(stag,'I_') GE 0, ctx )
+ print,'Chi-square of model to :'
+ for ii=0,ctx-1 do print,stag(itx(ii))+': ', strtrim(smdat.(itx(ii)).chi2,2)
+ print,'Chi-square of model to all data : ', strtrim(smdat.chi2,2)
+
+;
+; reset defaults
+;
+ !x.tickname = ''
+ !y.tickname = ''
+
+ if HARD(0) NE '0' then begin
+   device, /close
+   set_plot,'x'
+   print,'(W): hardcopies in ',hard
+   !y.thick = 0
+   !x.thick = 0
+   !p.thick = 0
+   !p.charsize = 1
+   !p.charthick = 0
+   !x.ticklen = 0.02
+ endif
+
+ RETURN
+ the_end:
+ END
+
+
+
+
@@ -0,0 +1,36 @@
+FUNCTION DUSTEM_STR_INST, n1, n2=n2, n3=n3
+; returns structure containing flux and CC in band 
+; N1   (I): number of data points
+; N2   (I): nr of grain types or different models
+; N3   (I): nr of points in transmission (band flux data)
+
+  if n_elements(n2) EQ 0 then n2 = 1
+
+  if n_elements(n3) NE 0 then begin
+     strct = { NAME :  strarr(n1),       $
+               X :     dblarr(n1),       $  
+               YD :    dblarr(n1),       $
+               ERR :   dblarr(n1),       $
+               YM  :   dblarr(n1,n2),    $
+               FLX :   dblarr(n1,n2),    $
+               CC :    dblarr(n1,n2),    $                   
+               RR :    dblarr(n1),       $                   
+               ISEL :  intarr(n1)+1,     $
+               UNIT :  '',               $ 
+               NPAR :  0,                $
+               CHI2 :  0d                 }
+     s1 = CREATE_STRUCT('TRANS', {X : dblarr(n1,n3), Y : dblarr(n1,n3) })
+     strct = CREATE_STRUCT( strct, s1 )
+  endif else if n_elements(BAND) EQ 0 then begin
+        strct = { X :     dblarr(n1),           $  
+                  YD :    dblarr(n1),           $
+                  ERR :   dblarr(n1) ,          $
+                  YM :    dblarr(n1,n2),        $
+                  ISEL :  intarr(n1)+1,         $
+                  UNIT :  '',                   $ 
+                  NPAR :  0,                    $
+                  CHI2 :  0d                      }
+  endif
+
+  RETURN, strct
+END