FUNCTION dustem_plugin_stellar_population, key=key, val=val, scope=scope, paramtag=paramtag,help=help

;+
; NAME:
;    dustem_plugin_stellar_population
; PURPOSE:
;    replaces the default DUSTEM ISRF with a composite stellar spectrum
; CATEGORY:
;    DUSTEM Wrapper
; CALLING SEQUENCE:
;    dustem_plugin_stellar_population(key=key,val=val)
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    key  = input parameter number
;    val  = input parameter value
; 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:
;    This is a dustem plugin
;-

;================NOTA BENE========================================================================== 
;DATA IS NOW RETRIEVED FROM THIS TEXT FILE: "A Modern Mean Dwarf Stellar Color and Effective Temperature Sequence"  
;NB: intermediate spectral classes '.5' are not taken into account - this will probably have to change
;NB: Also contact we need to contact the reasercher who wrote the text file because he said so in it.
;Other luminosity classes Should be included (only MS so far). 
;BB approximation is a first degree 'bad' approximation because of the lack of radiative transfer especially at the 
;photosphere of stars. 
;REMARKS: BECAUSE WE STILL HAVEN'T SET THE DEFAULT VALUE FOR THE MAJORIY OF THE STELLAR POULATIONS,THE CORRESPONDING LINES ARE COMMENTED INSTEAD OF SETTING ARBITRARY VALUES. 
;CONSIDERED STARS: SPEC_TYPE(N=5) = OBAFG, LUM_CLASS(N=10) = IA+,IA,IAB,IB,II,III,IV,V,VI,VII 
;KM spectral types are not included because their UV part was not that important to excite the dust.
;This will help ease and shorten the fitting procedure
;If the user wants to use them without having to read the entirety of this plugin please contact the DustEmWrap team. 

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_plugin_stellar_population'
  goto,the_end
ENDIF

;print,keyword_set(scope)
;print,keyword_set(paramtag)
;print,keyword_set(key),keyword_set(val)
;stop

IF keyword_set(scope) THEN BEGIN
  message,'Scope keyword was set',/continue
  scope='ADD_ISRF'
  stellar_isrf=0.
  goto,the_end
ENDIF

;==== read stellar population parameters if not already in memory
defsysv,'!dustem_plugin_stellar_population',exist=exist
IF exist EQ 0 THEN BEGIN
    message,'Initializing stellar population parameters',/continue
    file = !dustem_wrap_soft_dir+'Data/STELLARPOPS/EEM_dwarf_UBVIJHK_colors_Teff.xcat'
    st=read_xcat(file,/silent)
    defsysv,'!dustem_plugin_stellar_population',st
ENDIF

Nstars=n_elements(!dustem_plugin_stellar_population)
paramvalues=fltarr(Nstars,2)    ;first parameter is distance to the stars in pc, second is the number of such stars
Nparam=Nstars*2
stellar_isrf=0.
spectral_types=!dustem_plugin_stellar_population.spt

IF keyword_set(paramtag) THEN BEGIN
    message,'paramtag keyword was set',/continue
    paramtag=strarr(Nparam)
    ii=0L
    FOR i=0L,Nparam-1 DO BEGIN
            ij=index2ij([i],[Nstars,2])
            IF ij[0,1] EQ 0 THEN cc='dist to ' ELSE cc='N of '
            paramtag[ii]=cc+spectral_types[ij[0,0]]
            ii=ii+1
    ENDFOR
    stellar_isrf=0.
    goto, the_end
ENDIF
stellar_isrf=0.

;stop

;===========Constants (in cgs)========== (except for the stellar population distance (in pc))
rsun2cm = 6.957e10   ;Rsun in cm
c2a = 3e18 ;speed of light in ansgtroms/s (because of the Astron's PLANCK function)
pc2cm = 3.086e18 ;cm (cgs)

IF keyword_set(key) THEN BEGIN
    ;stop
    FOR i=0L,n_elements(key)-1 DO BEGIN
        ij=index2ij([key[i]-1],[Nstars,2])
        paramvalues[ij[0,0],ij[0,1]]=val[i]
    ENDFOR
    ;stop
    ;=== compute ISRF
    isrf_st=((*!dustem_params).isrf)
    stellar_isrf=dblarr(n_elements(isrf_st)) ; array of zeros to contain the new ISRF values. 
    FOR i=0L,Nstars-1 DO BEGIN
        IF paramvalues[i,0] NE 0. THEN BEGIN
            TEFF=!dustem_plugin_stellar_population[i].TEFF      ;in K
            Rstar=!dustem_plugin_stellar_population[i].R_Rsun   ;in solar radius
            Lstar=!dustem_plugin_stellar_population[i].logL
            dist=paramvalues[i,0]*pc2cm  ; in cm
            number_of_stars=paramvalues[i,1]
            omega=(Rstar*rsun2cm/dist)^2
            print,number_of_stars,dist,Teff,omega
            this_component=number_of_stars*omega*dustem_planck_function(Teff,isrf_st.lambisrf) ;MJy
            stellar_isrf=stellar_isrf+this_component
        ENDIF
    ENDFOR
    ;pass output into ergs/cm2/s/Hz
    fact=1./1.e20*1.e7/1.e4
    stellar_isrf=stellar_isrf*fact  ;ergs/cm2/s/Hz
    stellar_isrf=stellar_isrf/(*!dustem_params).G0
ENDIF

;stop

the_end:

RETURN,stellar_isrf
END