PRO dustem_write_isrf_release,file,st,help=help

;+
; NAME:
;    dustem_write_isrf_release
;
; PURPOSE:
;    Writes file ISRF.DAT used by DustemWrap
;
; CATEGORY:
;    DustEMWrap, Distributed, Low-level, Fortran
;
; CALLING SEQUENCE:
;    dustem_write_isrf_release,file,st[,/help]
;
; INPUTS:
;    file : file name 
;    st   : input structure containing ISRF information (will be modified)
;
; OPTIONAL INPUT PARAMETERS:
;    None
;
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUT PARAMETERS:
;    None
;
; ACCEPTED KEY-WORDS:
;    help      = If set print this help
;
; COMMON BLOCKS: 
;    None
;
; SIDE EFFECTS:
;    Writes a file
;
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
;
; PROCEDURE AND SUBROUTINES
;    
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by Ilyes Choubani 2022
;    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_write_isrf_release'
  goto,the_end
ENDIF

;LOCATING THE ISRF PLUGINS

;IF PLUGINS AREN'T SET. WRITE THE DEFAULT DUSTEM ISRF.DAT FILE

Nplugins=n_elements(*!dustem_plugin)
scopes = strarr(Nplugins)
IF (*!dustem_plugin)[0].name EQ 'NONE' THEN GOTO, the_after

;IF PLUGINS ARE SET, LOOP OVER THEIR SCOPES.

;Getting the actual scopes
FOR i=0L, Nplugins-1 DO BEGIN
    ;IF isa((*!dustem_plugin).(i).scope) THEN BEGIN
        scopes[i] = (*!dustem_plugin)[i].scope
    ;ENDIF ELSE goto, the_after

ENDFOR

;TESTING FOR THE TWO CURRENT ISRF PLUGINS
tstpop = where(scopes EQ 'STELLAR_POPULATION',ctpop)

;tstusrisrf = where(scopes EQ 'USER_ISRF',ctusrisrf)
tstusrisrf = where(scopes EQ 'REPLACE_ISRF',ctusrisrf)

IF ctpop NE 0 or ctusrisrf NE 0 THEN BEGIN

    ;stop

    ;NB: Since we're using G0_mathis=1 by default now. We will test on the !dustem_params sysvar instead of the parameter description vectors
    
    G0_mathis = 1. ;Default value
    

    IF ((*!dustem_params).g0) NE 1. then G0_mathis = ((*!dustem_params).g0)
; Fortran now uses negative numbers to specify Gas.dat G0
; If positive G0 present in Gas.dat, it takes precedence over grain.dat G0   
;   IF ((*!dustem_params).gas.g0) NE 1. then G0_mathis = float((*!dustem_params).gas.g0)
    IF ((*!dustem_params).gas.g0) gt 0. then G0_mathis = float((*!dustem_params).gas.g0)

       
    ;Getting the ISRF default wavelengths
    ;using this structure
    st=((*!dustem_params).isrf)
    final_isrf = fltarr(n_elements(st)) ;Initializing the ISRF vector
    Ncomments=3 ;default number of comments in the .DAT file (minus the ISRF used)
    c=strarr(Ncomments)

    c[0]='# DUSTEM_WRAP: exciting radiation field featuring'

    ;=== ADD to the ISRF
    IF ctpop NE 0. THEN final_isrf=final_isrf+(*(*!dustem_plugin)[tstpop].spec)
    IF ctusrisrf NE 0. THEN final_isrf=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=strarr(Ncomments)

        c[0]='# DUSTEM_WRAP: exciting radiation field featuring'  
        c[1] = '# Stellar population ISRF'
        c[2] = '# User-defined ISRF'

    ENDIF
    
    ;WHEN G0 IS USED (MATHIS FIELD ON)
    IF G0_mathis GT 1E-10 THEN BEGIN  ;discussed with JP
        ;This is reading the mathis isrf from the fortran release

        mathis_isrf_datfile=!dustem_soft_dir+'data/ISRF_MATHIS.DAT'
        mathis_isrf=dustem_read_isrf(mathis_isrf_datfile)
        final_isrf=final_isrf/G0_mathis+mathis_isrf.isrf ;[0]

        Ncomments=4
        c=strarr(Ncomments)

        c[0]='# DUSTEM_WRAP: exciting radiation field featuring'

        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+=2
            c=strarr(Ncomments)

            c[0]='# DUSTEM_WRAP: exciting radiation field featuring'  
            c[1] = '# Stellar population ISRF'
            c[2] ='# User-defined ISRF'
            c[3] = '# Mathis ISRF'
        ENDIF ELSE c[Ncomments-3] = '# Mathis ISRF'

        
    ENDIF ELSE final_isrf = final_isrf/G0_mathis ;WHEN G0 ISN'T USED (MATHIS FIELD OFF) ie: lower values than 1.0E-10.
    
    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)'

    ;Computing the 'G0' of the new ISRF (integral(ISRF)/integral(MATHIS_ISRF))
;         un=uniq(ma_isrf.isrf) ;because of the zeros
;         
;         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=4
c=strarr(Ncomments)

;First lines of the ISRF.DAT file

c[0]='# DUSTEM_WRAP: exciting radiation field featuring'
c[1]='# Mathis ISRF'
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