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