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,write_isrf_file

;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

Ncomments=3 ;default number of comments in the .DAT file (minus the ISRF used)
c=strarr(Ncomments)
c[*]='#'
c[0]='# DUSTEM_WRAP: Interstellar radiation field '

;Testing for the two current ISRF plugin scopes
ind_replace_isrf = where(scopes EQ 'REPLACE_ISRF',count_replace_isrf)
ind_add_isrf = where(scopes EQ 'ADD_ISRF',count_add_isrf)

;stop

IF count_add_isrf NE 0 OR count_replace_isrf NE 0 THEN BEGIN
    lambir=dustem_get_wavelengths()
    ;IF ctpop NE 0 or ctusrisrf NE 0 THEN BEGIN

    ;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 BEGIN
        G0_mathis = ((*!dustem_params).g0)
    ENDIF
    ; Fortran now uses negative numbers to specify Gas.dat G0
    ; If positive G0 present in Gas.dat, it takes precedence over the G0 value in GRAIN.DAT
    ; IF ((*!dustem_params).gas.g0) NE 1. then G0_mathis = float((*!dustem_params).gas.g0)
    IF ((*!dustem_params).gas.g0) gt 0. THEN BEGIN
        G0_mathis = float((*!dustem_params).gas.g0)
    ENDIF

    ;Getting the ISRF default wavelengths
    ;using this structure

    ;st=((*!dustem_params).isrf)

    final_isrf = fltarr(n_elements(st)) ;Initializing the ISRF vector

    ;=== ADD to the ISRF

    ;IF ctpop NE 0. THEN final_isrf=final_isrf+(*(*!dustem_plugin)[tstpop].spec)
    ;JPB: This line below can have an undefined plugin ISRF if the REPLACE_ISRF plugin has not been called earlier except with /scope or paramtags keywords
    IF count_replace_isrf NE 0. THEN BEGIN
        IF ptr_valid((*!dustem_plugin)[ind_replace_isrf].spec) THEN BEGIN
                c[1]='# DUSTEM_WRAP: replaced by radiation field from plugin '+(*!dustem_plugin)[ind_replace_isrf].name
                ;final_isrf=final_isrf+(*(*!dustem_plugin)[tstusrisrf].spec)

                ;final_isrf=interpol(*(*!dustem_plugin)[ind_replace_isrf].spec,lambir,st.lambisrf)
                final_isrf=*(*!dustem_plugin)[ind_replace_isrf].spec

        ENDIF ELSE BEGIN
                message,'Replacement ISRF not found. Using Mathis',/continue
                ;stop
                mathis_isrf_datfile=!dustem_soft_dir+'data/ISRF_MATHIS.DAT'
                sst=dustem_read_isrf(mathis_isrf_datfile)
                final_isrf=sst.isrf
                 ;stop
                ;LEFT blank intensionally
       ENDELSE

    ENDIF

    IF count_add_isrf NE 0 THEN BEGIN
        IF ptr_valid((*!dustem_plugin)[ind_add_isrf].spec) THEN BEGIN
                c[2]='# DUSTEM_WRAP: added radiation field from plugin '+(*!dustem_plugin)[ind_replace_isrf].name

                final_isrf=final_isrf+(*(*!dustem_plugin)[ind_add_isrf].spec)
                ;final_isrf=final_isrf+interpol(*(*!dustem_plugin)[ind_add_isrf].spec,lambir,st.lambisrf)

        ENDIF ELSE BEGIN
                ;stop
                ;LEFT blank intensionally
        ENDELSE
    ENDIF    
    ;JPB: I don't think this loop is in fact needed, as the default value (without plugins) is the Mathis ISRF and ADD_ISRF will just add to this
    ;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 ;This is because Fortran uses ISRF*G0 which is now final_isrf+G0*mathis_isrf
    ;    c=[c,'# Mathis ISRF']
    ;    Ncomments=Ncomments+1
    ;ENDIF ELSE BEGIN
    ;    ;final_isrf = final_isrf/G0_mathis ;This is because Fortran uses ISRF*G0
    ;ENDELSE
    ;=== set the ISRF in the st structure

    st.isrf = final_isrf  

ENDIF ELSE BEGIN       ;This is Mathis field only
    ;==== default comment lines for ISRF.DAT
    Ncomments=4
    c=strarr(Ncomments)
    c[1]='# Mathis ISRF'
ENDELSE

c[Ncomments-2]='# Nbr of points'
c[Ncomments-1]='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)'

;==== Write the ISRF file
write_isrf_file:

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

;print,(*!dustem_params).g0
;print,(*!dustem_params).gas.g0
;cgplot,st.lambisrf,st.isrf,/xlog,/ylog
;stop

the_end:


END