Merge branch 'master' of https://gitlab.irap.omp.eu/OV-GSO-DC/dustem-wrapper_idl

Jean-Philippe Bernard
2 parents b1818a4d a38807fc
Showing 17 changed files with 143 additions and 178 deletions Show diff stats
src/idl/dustem_cut_off.pro
src/idl/dustem_fit_modified_bb_readme.pro
src/idl/dustem_habing_field.pro
src/idl/dustem_logn_vdist.pro
src/idl/dustem_mathis_field.pro
src/idl/dustem_plaw_vdist.pro
src/idl_extern/ICE_for_Dustemwrap/modifications_for_pilot.txt -> src/idl_extern/ICE_for_Dustemwrap/modifications_log.txt
src/idl/dustem_fit_polextsed_readme.pro -> src/idl_retired/dustem_fit_polextsed_readme.pro
src/idl/dustem_fit_redshift_readme.pro -> src/idl_retired/dustem_fit_redshift_readme.pro
src/idl/dustem_fit_sed_extsed_readme.pro -> src/idl_retired/dustem_fit_sed_extsed_readme.pro
src/idl/dustem_fit_sed_polsed_readme.pro -> src/idl_retired/dustem_fit_sed_polsed_readme.pro
src/idl/dustem_fit_sed_qsed_used_readme.pro -> src/idl_retired/dustem_fit_sed_qsed_used_readme.pro
src/idl/dustem_fit_sed_readme.pro -> src/idl_retired/dustem_fit_sed_readme.pro
src/idl/dustem_fit_sed_spinning_readme.pro -> src/idl_retired/dustem_fit_sed_spinning_readme.pro
src/idl/dustem_fit_spin_test_readme.pro -> src/idl_retired/dustem_fit_spin_test_readme.pro
src/idl/dustem_run_readme.pro -> src/idl_retired/dustem_run_readme.pro
src/idl/dustem_variational_study_readme.pro -> src/idl_retired/dustem_variational_study_readme.pro
@@ -0,0 +1,16 @@
+FUNCTION DUSTEM_CUT_OFF, x, par
+; generates the large size cut-off 
+; from Weingartner & Draine 2001
+; x : grain size 
+; par(0) : threshold for cut-off (At)
+; par(1) : shape As (roughly the size where cut_off=0.5)
+
+ y = 0.d0*x + 1.d0
+ ix = WHERE( x GT par(0), cnt )
+ if CNT GT 0 then begin 
+     y(ix) = exp( -( (x(ix)-par(0)) / par(1))^3. )
+ endif else begin 
+     print,'(W)  DUSTEM_CUT_OFF: no size larger than At found'
+ endelse
+RETURN, y
+END
@@ -1,178 +0,0 @@
-PRO dustem_fit_modified_bb_readme,postcript=postcript,mode=mode,help=help
-
-;This Readme describes how to fit SEDs with a modified Black-Body
-;It runs on the SED stored into the file sample_SED.xcat
-;Remember that the goal here is just to illustrate the method.
-;Note that its is not necessary to use the .xcat file format, and data SED can be provided
-;manually, but the observation structure must have the structure shown below.
-;Obviously, the dustem package must have been installed succesfully
-;see DustemWrap documentation for install instructions.
-
-;+
-; NAME:
-;    dustem_fit_modified_bb_readme
-; PURPOSE:
-;    This is an example of how to fit SEDs with a modified Black-Body
-;    using the dustem wrapper.
-;    It is meant to be an example to follow when writing your own
-;    programs using the dustem IDL wrapper.
-;    The SED used here is in sample_SED.xcat
-; CATEGORY:
-;    Dustem
-; CALLING SEQUENCE:
-;    dustem_fit_modified_bb_readme,postcript=postcript,help=help
-; INPUTS:
-;    None
-; OPTIONAL INPUT PARAMETERS:
-;    None
-; OUTPUTS:
-;    None
-; OPTIONAL OUTPUT PARAMETERS:
-;    None
-; ACCEPTED KEY-WORDS:
-;    postcript = if set plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps
-;    help      = If set print this help
-; COMMON BLOCKS:
-;    None
-; SIDE EFFECTS:
-;    None
-; RESTRICTIONS:
-;    The dustem idl wrapper must be installed
-; PROCEDURE:
-;    None
-; EXAMPLES
-;    dustem_fit_modified_bb_readme
-; MODIFICATION HISTORY:
-;    Written by J.P. Bernard June 9th 2011
-;    see evolution details on the dustem cvs maintained at CESR
-;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
-;-
-
-message,'This is obsolete.See dustem_fit_intensity_mbb_example.pro',/continue
-stop
-
-IF keyword_set(help) THEN BEGIN
-  doc_library,'dustem_fit_modified_bb_readme'
-  goto,the_end
-ENDIF
-
-;=== initialise dustem
-dustem_init
-!dustem_verbose=1
-!dustem_show_plot=1
-
-;=== Read sample SED
-;=== Composite SED from Compiegne et al 2010, gathered by C. Bot
-;dir=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/SEDs/'
-;dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
-;file=dir+'Gal_composite_spectrum.xcat'
-;sed=read_xcat(file,/silent)
-;=== Or make a fake one out for given filters
-;filters=['IRAS4','HFI1','HFI2','HFI3']
-filters=['IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3']
-Nfilt=n_elements(filters)
-sed=dustem_initialize_sed(Nfilt)
-sed.filter=filters
-sed.wave=dustem_filter2wav(filters)
-sed.instru=dustem_filter2instru(filters)
-
-;VG : dustem_set_data turned into a function to be used indifferently for sed, ext, polext
-;dustem_set_data,sed  ;This is to set filters in !dustem_data.sed. Needed by dustem_compute_gb_sed
-;!dustem_data.sed=dustem_set_data(sed) ;This is to set filters in !dustem_data.sed. Needed by dustem_compute_gb_sed
-st=dustem_set_data(sed=sed)
-
-;p_truth=[0.5D0,22.D0,1.8D0]
-p_truth=[0.5,19.,1.8]
-;sed.spec=dustem_compute_gb_sed(p_truth,_extra=extra)
-sed.stokesI=dustem_compute_gb_sed(p_truth,_extra=extra)
-;sed.error=sed.spec*0.05
-sed.sigmaII=(sed.stokesI*0.05)^2
-
-;=== SET THE OBSERVATION STRUCTURE
-
-ind=where(sed.wave GE 100.)
-;VG : dustem_set_data turned into a function to be used indifferently for sed, ext, polext
-;dustem_set_data,sed(ind)
-st=dustem_set_data(sed=sed(ind))
-;!dustem_data.sed=dustem_set_data(sed(ind))
-
-;=== Set which parameters you want to fit
-pd = ['Amplitude','T','beta']
-;=== set initial value of parameters you want to fit
-;iv =   [1.,20.,2.]
-;iv =   [0.5,22.D0,1.5D0]
-iv=[0.5,22.,1.8]
-ulimed=[0   , 0   ,0]
-llimed=[1   , 1   ,1]
-llims =[0.  ,0.   ,0.]
-
-;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
-dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
-
-
-;=== RUN fit
-;tol=1.e-20
-tol=1.e-20
-xtol=1.e-10
-Nitermax=100               ;maximum number of iteration. This is the criterium which will stop the fit procedure
-
-loadct,13
-!y.range=[1e-4,10]              ;This is to ajust plot range from outside the routine
-t1=systime(0,/sec)
-res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,function_name='dustem_greybody_mpfit',/xlog,/ylog)
-t2=systime(0,/sec)
-
-;=== SAVE FIT RESULTS
-;file_out=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
-file_out='/tmp/DUSTEM_bb_fit_example.sav'
-dustem_save_system_variables,file_out
-message,'Saved '+file_out,/continue
-
-
-;stop
-
-;======================================
-;====You can exit IDL here and re-enter
-;======================================
-
-file='/tmp/DUSTEM_bb_fit_example.sav'
-dustem_restore_system_variables,file
-
-window,1
-
-;=== RESTORE FIT RESULTS
-;=== initialise dustem
-;dustem_init
-
-;=== recover best fit values
-res=*(*!dustem_fit).current_param_values
-chi2=(*!dustem_fit).chi2
-rchi2=(*!dustem_fit).rchi2
-errors=*(*!dustem_fit).current_param_errors
-
-;=== Plot best fit
-yr=[1e-3,100]
-xr=[50,4000]
-tit='DUSTEM Modified BB Example (Saved)'
-ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
-xtit=textoidl('\lambda (\mum)')
-
-loadct,13
-IF keyword_set(postcript) THEN BEGIN
-  set_plot,'PS'
-  ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_bb_fit.ps'
-  device,filename=ps_file,/color
-ENDIF
-dustem_sed_plot,*(*!dustem_fit).current_param_values, $
-                ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit',/xlog,/ylog
-IF keyword_set(postcript) THEN BEGIN
-  device,/close
-  set_plot,'X'
-  message,'wrote '+ps_file,/info
-ENDIF
-
-print,'dustem_mpfit_sed executed in ',t2-t1,' sec'
-
-the_end:
-
-END
@@ -0,0 +1,33 @@
+FUNCTION DUSTEM_HABING_FIELD, x, unit=unit
+; computes the Habing interstellar radiation field SED (ISRF)
+; in W/m2 (4*!pi*nu*I_nu)
+; (from Draine & Bertoldi 1996)
+; X	(I): X-grid for the ISRF
+; UNIT	(I): unit of the ISRF. Choices are 
+;	     'W': wave in um [Default]
+;	     'F': frequency in cm-1
+;	     'E': energy in eV 
+
+ if n_elements( UNIT ) EQ 0 then unit = 'W' $
+ else unit = strupcase( strcompress( unit,/rem) )
+
+ x = double( x )
+
+ CASE unit of 
+
+   'W': x3 = 1.d1 * x 
+
+   'F': x3 = 1.d5 / x
+
+   'E': x3 = 12.4 / x
+
+ ENDCASE
+
+ field = - x3^3*4.1667 + x3^2*12.5 - x3*4.3333 
+ field = 1.d-1 * 3.d8 * 1.d-14 * field 
+
+ i_neg = where( field LT 0.d0, c_neg )
+ field( i_neg ) = 0.d0
+
+ RETURN, field
+END                             ;FUNCTION DUSTEM_HABING_FIELD
@@ -0,0 +1,20 @@
+FUNCTION DUSTEM_LOGN_VDIST, x, par
+; generates a volume normalized log-normal law 
+; returns distribution in nr of grains : dn/da
+; x : grain size 
+; par(0) : centroid of log-normal
+; par(1) : sigma of log-normal
+; par(2) : VOLUME normalization
+
+ np = n_elements(x)
+ x0 = par(0)
+ sigma = par(1)
+ y = exp(- 0.5 * ( alog(x/x0) / sigma )^2 ) / x
+ vy = x^4 * y
+ xm = par(0) * exp( -par(1)^2 ) ; x of max in dn/da
+ print,'(W) DUSTEM_LOGN_VDIST: dn/da max @',xm*1e7,' nm'
+ dx = alog(x(1:np-1)) - alog(x(0:np-2)) 
+ yi = TOTAL( (vy(1:np-1) + vy(0:np-2))*0.5*dx )
+ y = par(2) * y / yi
+RETURN, y
+END 
@@ -0,0 +1,52 @@
+FUNCTION DUSTEM_MATHIS_FIELD, x, unit=unit
+; computes the Mathis interstellar radiation field SED (ISRF)
+; in W/m2 (4*!pi*nu*I_nu)
+; from Mathis et al. 1983, A&A 128, 212
+; X	(I): X-grid for the ISRF
+; UNIT	(I): unit of the ISRF. Choices are 
+;	     'W': wave in um [Default]
+;	     'F': frequency in cm-1
+;	     'E': energy in eV
+
+ if n_elements( UNIT ) EQ 0 then unit = 'W' $
+ else unit = strupcase( strcompress( unit,/rem) )
+
+ ly_limit = 9.11267101235d-2
+
+ x = double( x )
+
+;
+; visible part 
+;
+ wdil = 4.d0 * [ 4.d-13, 1.d-13, 1.d-14]	; Mathis definition
+						; 4*!pi*Wdil*B_nu
+ field = !pi * x * DUSTEM_BLACKBODY( x, [ 3.d3, 4.d3, 7.5d3 ], wdil=wdil, unit=unit )
+
+;
+; UV part
+;
+
+; first convert to (lambda / 1e3 AA)
+ CASE unit of 
+
+   'W': x3 = 1.d1 * x 
+
+   'F': x3 = 1.d5 / x
+
+   'E': x3 = 12.4 / x
+
+ ENDCASE
+ 
+ il = where( x3 LE 2.46, cl ) 
+ if cl GT 0 then field(il) = 0.D0
+
+ il = where( x3 GE 1d1*ly_limit AND x3 LT 1.11, cl )
+ if cl GT 0 then field(il) = 1.4817d-6 * x3(il)^(4.4172)
+ il = where( x3 GE 1.11 AND x3 LT 1.34, cl ) 
+ if cl GT 0 then field(il) = 2.0456d-6 * x3(il)
+ il = where( x3 GE 1.34 AND x3 LT 2.46, cl )
+ if cl GT 0 then field(il) = 3.3105d-6 * x3(il)^(-0.6678)
+  
+
+ RETURN, field
+END                             ;FUNCTION DUSTEM_MATHIS_FIELD
@@ -0,0 +1,22 @@
+FUNCTION DUSTEM_PLAW_VDIST, x, par 
+; generates a volume normalized power law x^par(0)
+; returns distribution in nr of grains : dn/da
+; x : grain size 
+; par(0) : power law index
+; par(1) : VOLUME normalization 
+; par(2) : curvature parameter beta 
+; par(3) : large size threshold At
+
+ np = n_elements(x)
+ y = x^par(0)
+; curvature term
+ if ((par(2) NE 0) AND (par(3) NE 0)) then begin 
+     psgn = par(2)/ABS(par(2))
+     y = y * ( 1.d0 + ABS(par(2))*x/par(3) )^psgn
+ endif 
+ vy = x^4 * y 
+ dx = ALOG(x(1:np-1)) - ALOG(x(0:np-2)) 
+ yi = TOTAL( (vy(1:np-1) + vy(0:np-2))*0.5*dx )
+ y = par(1) * y / yi
+RETURN, y
+END