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:

;    a Stellar ISRF is added based on the content of (*(*!dustem_plugin).stellar_population)

;-

IF keyword_set(help) THEN BEGIN

  doc_library,'dustem_write_isrf_lv'

  goto,the_end
ENDIF

IF tag_exist((*!dustem_scope),'STELLAR_POPULATION') && ISA(((*!dustem_plugin).stellar_population)) THEN BEGIN
    c2a = 3e18 ;speed of light in ansgtroms/s (because of the Astron's PLANCK function)
    pc2cm = 3.086e18 ;cm (cgs)
    wave_angstrom = st.lambisrf*1.e4 ;mic to Angstrom (Astron Planck's function uses wavelengths in Angstroms)

    stellar_component=fltarr(n_elements(st)) ; array of zeros to contain the new ISRF values. 
     
    Ncomments = n_elements(((*(*!dustem_plugin).stellar_population).popid))*3+4
    c = strarr(Ncomments)
    
    
    ;First and last lines of the new composite ISRF.DAT file
    c(0)='# DUSTEM: exciting radiation field featuring'
    c(1)='# Mathis ISRF'
    c(Ncomments-2)='# Nbr of points'
    c(Ncomments-1)='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'

    FOR i=0L,n_elements(((*(*!dustem_plugin).stellar_population).popid))-1 DO BEGIN ; Looping over all the stellar populations
        
        ;The initial procedure had omega multiplied by a !pi factor. Its presence in the Planck (Astron) procedure makes for a good reason to discuss this with J.P.  
        
        omega =  (((*(*!dustem_plugin).stellar_population).radius)[i]/(((*(*!dustem_plugin).stellar_population).distance)[i]*pc2cm))^2 ; Dilution factor of a stellar population
            
        Inu = planck(wave_angstrom,((*(*!dustem_plugin).stellar_population).temperature)[i])/(4.*!pi)*(wave_angstrom)^2/c2a ; ergs/cm2/s/Hz/sr
        
        ;========================================================NOTA BENE==================================================
        ;Deving the amplitudes of the stellar populations by G0 as the composite ISRF is multiplied by the latter as a whole
        ;===================================================================================================================
        
        stellar_component=stellar_component+((*(*!dustem_plugin).stellar_population).amplitude)[i]*((*(*!dustem_plugin).stellar_population).nstars)[i]*omega*Inu;/((*!dustem_params).gas.G0)
        
        ;Rest of the lines of the new composite ISRF.DAT file
        c(3*i+2) = '#'+((*(*!dustem_plugin).stellar_population).popid)[i]   
        c(3*i+3) = '# Blackbody with     T='+string(((*(*!dustem_plugin).stellar_population).temperature)[i])  
        c(3*i+4) = '# dilution factor wdil='+string(omega)
            
    ENDFOR
    
    print,'stellar_component is:'
    
    print, stellar_component
    
    print, 'isrf is:'
    print, st.isrf

    
    st.isrf=st.isrf+stellar_component/((*!dustem_params).G0)
    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
    
    
ENDIF ELSE BEGIN

    Ncomments=6
    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(4)='# Nbr of points'
    c(5)='# 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:

ENDELSE

END