PRO dustem_write_isrf_lv,file,st,help=help

;+
; NAME:
;    dustem_write_isrf_lv
; PURPOSE:
;    Writes file ISRF.DAT used by DustemWrapper
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_write_isrf_lv,file,st[,/help]
; INPUTS:
;    file : file name 
;    st   : input ISRF structure
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set print this help

; COMMON BLOCKS: 

;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem fortran code must be installed
;    The dustem idl wrapper must be installed
; PROCEDURE:

;    Write ISRF for the DustEmWrap

;-

IF keyword_set(help) THEN BEGIN

  doc_library,'dustem_write_isrf_lv'

  goto,the_end
ENDIF

;Replacing the old block with an iteration over the plugins' scopes 
scopes = strarr(n_tags(*!dustem_plugin))
IF n_tags(*!dustem_plugin) EQ 1 and (tag_names(*!dustem_plugin))[0] EQ 'NONE' then goto, the_after
;Getting the actual scopes 
FOR i=0L, n_tags(*!dustem_plugin)-1 DO BEGIN
    scopes[i] = (*(*!dustem_plugin).(i).scope)
ENDFOR

tstpop = where(scopes EQ 'STELLAR_POPULATION',ctpop)
tstusrisrf = where(scopes EQ 'USER_ISRF',ctusrisrf)

IF ctpop NE 0 or ctusrisrf NE 0 THEN BEGIN
    
    G0_mathis = 0. ;Just an initialization
    ;The initial test was going to be on ((*!dustem_params).g0) and ((*!dustem_params).gas.g0) but it isn't enough...
    ;this is because the case where mathis ISRF is used with G0=1 oughtn't be forgotten and should be distinguished
    ;A solution is to instead test on the parameter vectors (fixed and thawed alike)

    ;Locating G0 in the free parameters vector
    ind_G0 = where(strupcase(strmid((*(*!dustem_fit).param_descs),5,/reverse)) EQ 'MS).G0', countg0)
    IF countg0 NE 0 THEN G0_mathis = (*(*!dustem_fit).param_init_values)[ind_G0] 
    
    ;Locating G0 in the fixed parameters vector
    ;under the assumption that the user will always have a free parameter vector set but not necessarily a fixed one
    IF isa((*!dustem_fit).fixed_param_descs) THEN BEGIN
        ind_G0 = where(strupcase(strmid((*(*!dustem_fit).fixed_param_descs),5,/reverse)) EQ 'MS).G0', countg0)
        IF countg0 NE 0 THEN G0_mathis = (*(*!dustem_fit).fixed_param_init_values)[ind_G0]
    ENDIF
    
    ;Locating gas.G0 in the free parameters vector
    ind_GASG0 = where(strupcase(strmid((*(*!dustem_fit).param_descs),5,/reverse)) EQ 'GAS.G0', countgasg0)
    IF countgasg0 THEN G0_mathis = (*(*!dustem_fit).param_init_values)[ind_GASG0]
    
    ;Locating gas.G0 in the fixed parameters vector
    ;under the assumption that the user will always have a free parameter vector set but not necessarily a fixed one
    IF isa((*!dustem_fit).fixed_param_descs) THEN BEGIN
        ind_GASG0 = where(strupcase(strmid((*(*!dustem_fit).fixed_param_descs),5,/reverse)) EQ 'GAS.G0', countgasg0)
        IF countgasg0 THEN G0_mathis = (*(*!dustem_fit).fixed_param_init_values)[ind_GASG0]
    ENDIF   
    ;Getting the ISRF default wavelengths
    ;using this structure
    st=((*!dustem_params).isrf)
    final_isrf = fltarr(n_elements(st)) ;Initializing the ISRF vector
    Ncomments=4 ;default number of comments in the .DAT file
    c=strarr(Ncomments)
    c(0)='# DUSTEM: exciting radiation field featuring'
    IF ctpop NE 0 THEN final_isrf+=(*(*!dustem_plugin).(tstpop).spec)
    IF ctusrisrf NE 0 THEN final_isrf+=(*(*!dustem_plugin).(tstusrisrf).spec)

    IF ctpop NE 0 and ctusrisrf EQ 0 THEN c(1) = '#STELLAR POPULATION ISRF'
    IF ctpop EQ 0 and ctusrisrf NE 0 THEN  c(1) = '#USER-DEFINED ISRF'
    
    IF ctpop NE 0 and ctusrisrf NE 0 THEN BEGIN
        Ncomments+=1
        c(1) = '#STELLAR POPULATION ISRF'
        c(2) = '#USER-DEFINED ISRF'
    ENDIF
    
    
    IF G0_mathis NE 0 THEN BEGIN ;storing the mathis isrf in variable mathis_isrf if the user wants to use it
        ;This is reading the mathis isrf from the fortran release
        ;Let's hope they keep releasing it otherwise we'll have to include it ourselves
        ma_isrf_dir=!dustem_soft_dir+'data/ISRF_MATHIS.DAT'
        ma_isrf=dustem_read_isrf(ma_isrf_dir)
        final_isrf=final_isrf/G0_mathis+ma_isrf
        Ncomments+=1
        c(Ncomments-3) = '#MATHIS ISRF'
    ENDIF
    
    ;Will this enable us to have the current 'updated' ISRF through the !dustem_params system variable?
    st.isrf = final_isrf  
    
    ;Last .DAT header comments lines
    c(Ncomments-2)='# Nbr of points'
    c(Ncomments-1)='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'

    ;internally displaying the 'G0' of the new ISRF (integral(ISRF)/integral(MATHIS_ISRF))
    ;     un=uniq(ma_isrf.isrf)
    ;     
    ;     int_mathis=INT_TABULATED((ma_isrf.lambisrf)[un],(ma_isrf.isrf)[un])
    ;     
    ;     
    ;     ;integrating the composite isrf
    ;     int_isrf=INT_TABULATED((ma_isrf.lambisrf)[un],(st.isrf)[un])
    ;     print, 'G0_stellar_population='
    ;     print, int_isrf/(int_mathis)
    
    ;WRITING THE NEW ISRF
    
    file=!dustem_dat+'data/ISRF.DAT'  
    openw,unit,file,/get_lun
        
    FOR i=0,Ncomments-1 DO BEGIN
        printf,unit,c(i)
    ENDFOR

    n_waves=n_elements(st)
    printf,unit,n_waves

    FOR i=0L,n_waves-1 DO BEGIN
        printf,unit,st(i).lambisrf,st(i).isrf
    ENDFOR

    close,unit
    free_lun,unit 

    goto, the_end

ENDIF

the_after:
;Defaut ISRF writing (MATHIS)
;This was changed because 6 comments aren't needed
;Ncomments=6
Ncomments=4

c=strarr(Ncomments)

;First lines of the ISRF.DAT file
c(0)='# DUSTEM: exciting radiation field featuring'
c(1)='# Mathis ISRF'

;c(2)='# Blackbody with     T=   0.0000E+00'
;c(3)='# dilution factor wdil=   1.0000E+00'
c(2)='# Nbr of points'
c(3)='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'

openw,unit,file,/get_lun

FOR i=0,Ncomments-1 DO BEGIN
    printf,unit,c(i)
ENDFOR

n_waves=n_elements(st)
printf,unit,n_waves

FOR i=0L,n_waves-1 DO BEGIN
    printf,unit,st(i).lambisrf,st(i).isrf
ENDFOR

close,unit
free_lun,unit

the_end:

END