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

;+
; 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
;-

;- HISTORY: ADDED STELLAR LUMINOSITY CLASSES

IF keyword_set(test) THEN stop

IF keyword_set(help) THEN BEGIN

  doc_library,'dustem_plugin_stellar_population'

  goto,the_end
ENDIF

IF keyword_set(scope) THEN BEGIN
  message,'Scope keyword was set',/continue
  stop

ENDIF 

IF keyword_set(paramtag) THEN BEGIN
    message,'paramtag keyword was set',/continue
    stop
;    out=0
;    goto, the_paramtag
ENDIF
scope='STELLAR_POPULATION'
stellar_isrf=0.

;==== read stellar population parameters if not already in memory
IF keyword_set(initialize) THEN BEGIN
    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
    goto,the_end

ENDIF

;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)

;================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. 

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

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

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
ENDIF

;the_scope:
;scope='STELLAR_POPULATION'

the_end:

;print,minmax(stellar_isrf)
;stop

RETURN, stellar_isrf
END