corrected stellar_population & write_isrf_lv

Ilyes Choubani
1 parent dcf48925
Showing 2 changed files with 101 additions and 163 deletions Show diff stats
src/idl/dustem_plugin_stellar_population.pro
src/idl/dustem_write_isrf_lv.pro
@@ -3143,76 +3143,13 @@ IF keyword_set(key) THEN BEGIN
     ;==============================================================================================    
     ENDFOR
-      
-;============This block creates a composite (9V) stellar population structure for when the function isn't used as a plugin=================== 
-ENDIF ELSE BEGIN 
-
-;IF NOT keyword_set(scope) AND NOT keyword_set(val) THEN BEGIN ;I need this condition because I don't want the accidental definition of the system variable (@the end.) that is actually tied to any parameters outside this routine.
-
-    
-    popnumber = 5.
-    comp_pop = replicate(one_pop,popnumber)
-
-    
-    ;###BY DEFAULT ALL STELLAR POPULATIONS ARE ON THE MS###
-    ;BECAUSE LESS MASSIVE (COLDER) STARS ARE MORE ABUNDANT, THEY WILL BE CHOSEN TO BUILD THE COMPOSITE STELLAR STRUCTURE
-    ;THIS MEANS THAT THE CHOSEN SPECTRAL TYPE IS 9V (MAIN SEQUENCE)
-
-    ;CHOSEN STELLAR POPULATIONS
-    ;O9V
-    ;B9V
-    ;A9V
-    ;F9V
-    ;G9V
-    
-    comp_pop(0).popid = 'O9V_stellar_population'
-    comp_pop(0).radius = 7.51*rsun2cm 
-    comp_pop(0).temperature = 33.3e3
-    comp_pop(0).distance = 10.0
-    comp_pop(0).nstars = stellar_density*(4.*!pi/3)*comp_pop(0).distance
-    
-    comp_pop(1).popid = 'B9V_stellar_population'
-    comp_pop(1).radius = 2.49*rsun2cm  
-    comp_pop(1).temperature = 10.7e3 
-    comp_pop(1).distance = 10.0
-    comp_pop(1).nstars = stellar_density*(4.*!pi/3)*comp_pop(1).distance
-    
-    comp_pop(2).popid = 'A9V_stellar_population'
-    comp_pop(2).radius = 1.747*rsun2cm
-    comp_pop(2).temperature = 7.4e3
-    comp_pop(2).distance = 10.0
-    comp_pop(2).nstars = stellar_density*(4.*!pi/3)*comp_pop(2).distance
-    
-    comp_pop(3).popid = 'F9V_stellar_population'
-    comp_pop(3).radius = 1.167*rsun2cm
-    comp_pop(3).temperature = 6.05e3
-    comp_pop(3).distance = 10.0
-    comp_pop(3).nstars = stellar_density*(4.*!pi/3)*comp_pop(3).distance
-    
-    comp_pop(4).popid = 'G9V_stellar_population'
-    comp_pop(4).radius = 0.853*rsun2cm 
-    comp_pop(4).temperature = 5.38e3
-    comp_pop(4).distance = 10.0
-    comp_pop(4).nstars = stellar_density*(4.*!pi/3)*comp_pop(4).distance
-    
-    
-;ENDIF 
-
-ENDELSE
-
-if keyword_set(key) then begin
-print, 'comp_pop is:'
-print, comp_pop
-endif
-
-;stop
 ;============================================================================================== 
 ;========================================================NOTA BENE============================================================================
 ;THIS PROCEDURE DEVIDES BY THE VALUE OF G0. THE WRAPPER USES G0 AND GAS.G0. PLEASE LOCATE IT IN THE LINES BELOW AND CHANGE G0 TO THE ONE YOU'RE USING!!!!
 ;=============================================================================================================================================
-st=(*!dustem_params).isrf
+st=((*!dustem_params).isrf)
 ;stop ;to check if the structure is indeed what we need. 
 c2a = 3e18 ;speed of light in ansgtroms/s (because of the Astron's PLANCK function)
 pc2cm = 3.086e18 ;cm (cgs)
@@ -3288,9 +3225,73 @@ FOR i=0L,n_waves-1 DO BEGIN
 ENDFOR
 close,unit
-free_lun,unit
+free_lun,unit 
+
+out=st.isrf     
+;============This block creates a composite (9V) stellar population structure for when the function isn't used as a plugin=================== 
+ENDIF ELSE BEGIN 
+
+;IF NOT keyword_set(scope) AND NOT keyword_set(val) THEN BEGIN ;I need this condition because I don't want the accidental definition of the system variable (@the end.) that is actually tied to any parameters outside this routine.
+
+    
+    popnumber = 5.
+    comp_pop = replicate(one_pop,popnumber)
+
+    
+    ;###BY DEFAULT ALL STELLAR POPULATIONS ARE ON THE MS###
+    ;BECAUSE LESS MASSIVE (COLDER) STARS ARE MORE ABUNDANT, THEY WILL BE CHOSEN TO BUILD THE COMPOSITE STELLAR STRUCTURE
+    ;THIS MEANS THAT THE CHOSEN SPECTRAL TYPE IS 9V (MAIN SEQUENCE)
+
+    ;CHOSEN STELLAR POPULATIONS
+    ;O9V
+    ;B9V
+    ;A9V
+    ;F9V
+    ;G9V
+    
+    comp_pop(0).popid = 'O9V_stellar_population'
+    comp_pop(0).radius = 7.51*rsun2cm 
+    comp_pop(0).temperature = 33.3e3
+    comp_pop(0).distance = 10.0
+    comp_pop(0).nstars = stellar_density*(4.*!pi/3)*comp_pop(0).distance
+    
+    comp_pop(1).popid = 'B9V_stellar_population'
+    comp_pop(1).radius = 2.49*rsun2cm  
+    comp_pop(1).temperature = 10.7e3 
+    comp_pop(1).distance = 10.0
+    comp_pop(1).nstars = stellar_density*(4.*!pi/3)*comp_pop(1).distance
+    
+    comp_pop(2).popid = 'A9V_stellar_population'
+    comp_pop(2).radius = 1.747*rsun2cm
+    comp_pop(2).temperature = 7.4e3
+    comp_pop(2).distance = 10.0
+    comp_pop(2).nstars = stellar_density*(4.*!pi/3)*comp_pop(2).distance
+    
+    comp_pop(3).popid = 'F9V_stellar_population'
+    comp_pop(3).radius = 1.167*rsun2cm
+    comp_pop(3).temperature = 6.05e3
+    comp_pop(3).distance = 10.0
+    comp_pop(3).nstars = stellar_density*(4.*!pi/3)*comp_pop(3).distance
+    
+    comp_pop(4).popid = 'G9V_stellar_population'
+    comp_pop(4).radius = 0.853*rsun2cm 
+    comp_pop(4).temperature = 5.38e3
+    comp_pop(4).distance = 10.0
+    comp_pop(4).nstars = stellar_density*(4.*!pi/3)*comp_pop(4).distance
+    
+    out=0.
+;ENDIF 
+
+ENDELSE
+
+if keyword_set(key) then begin
+print, 'comp_pop is:'
+print, comp_pop
+endif
+
+;stop
 the_end:
-return, st.isrf
+return, out
 end
 \ No newline at end of file
@@ -37,103 +37,40 @@ IF keyword_set(help) THEN BEGIN
 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
+IF ptr_valid((*!dustem_plugin).stellar_population) THEN goto, the_end
+
+;NOW YOU HAVE TO LOOP ON THE COMMENTS TO SEE IF IT IS THE COMPOSITE ISRF FILE OR THE MATHIS ONE.
+
+Ncomments=6
+c=strarr(Ncomments)
+
+;MODIFY THIS TO DISPLAY THE ACTUAL VALUES BECAUSE IT IS JUST DISPLAYING EMPTY STRINGS AT THE MOMENT
+
+;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:
 END