PRO dustem_fit_sed_polsed_readme,postcript=postcript,mode=mode,help=help ;This Readme describes how to fit total intensity and polarization SEDs jointly ;Obviously, the dustem package must have been installed succesfully ;see DustemWrap documentation for install instructions. ;+ ; NAME: ; dustem_fit_sed_polsed_readme ; PURPOSE: ; This is an example of how to fit total intensity and polarization SEDs jointly ; It is meant to be an example to follow when writing your own ; programs using the dustem IDL wrapper. ; CATEGORY: ; Dustem ; CALLING SEQUENCE: ; dustem_fit_sed_polsed_readme,postcript=postcript,mode=mode,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_sed_polsed_readme ; MODIFICATION HISTORY: ; Written by J.P. Bernard June 29th 2011 ; see evolution details on the dustem cvs maintained at CESR ; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems. ;- IF keyword_set(help) THEN BEGIN doc_library,'dustem_fit_sed_polsed_readme' goto,the_end ENDIF dustem_define_la_common ;=== initialise dustem ;dustem_init,mode='COMPIEGNE_ETAL2010',/pol use_mode='G17_MODELA' dustem_init,mode=use_mode,/pol !dustem_verbose=1 !dustem_show_plot=1 help,!run_pol ;stop ;=== make a fake one out for given filters ;filters=['DIRBE1','DIRBE2','DIRBE3','DIRBE4','IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3'] filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3'] Nfilt=n_elements(filters) sed=dustem_initialize_sed(Nfilt) sed.filter=filters sed.wave=dustem_filter2wav(filters) sed.instru=dustem_filter2instru(filters) polsed=sed ; comment if polext is to be used ;polext=sed ;comment if polsed is to be used ;UNCOMMENT THIS BLOCK AFTER THE TESTING ; ;stokes fake parameters - comment if stokes extinction parameters are to be used ; qsed=sed ; used=sed ; ;stokes extinction fake parameters - comment if stokes parameters are to be used ; qext=sed ; uext=sed st=dustem_set_data(ext=ext,polsed=polsed,polext=polext,sed=sed,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI,qsed=qsed,used=used,qext=qext,uext=uext) ;=== Set which parameters you want to fit ; pd = [ $ ; '(*!dustem_params).G0', $ ;G0 ; '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction ; '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction ; '(*!dustem_params).grains(2).mdust_o_mh',$ ;amCBEx ; 'dustem_create_synchrotron_2',$ ; 'dustem_create_stokes_1',$ ;Polarization angle psi_ref_dust ; 'dustem_create_stokes_2',$;Polarization angle psi_ref_synch ; 'dustem_create_stokes_3'$;Polarization fraction of synchrotron spectrum ; ;'dustem_create_stext_1'$;Polarization angle in extinction PSI_ext ; ] ;TESTING THE ISEF PLUGIN ; ; pd = [ $i ; ;'(*!dustem_params).G0', $ ;G0 ; ;'(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction ; ;'(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction ; ;'(*!dustem_params).grains(2).mdust_o_mh',$ ;amCBEx ; ;'dustem_create_continuum_2', $ ;intensity of NIR continuum ; 'dustem_create_stellar_population_O63',$ ; distance to the O6 stellar population ; ;'dustem_create_stellar_population_O65',$ ; number of stars of the O6 stellar population ; 'dustem_create_stellar_population_B43'$ ; distance to the B4 stellar population ; ;'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population ; ] pd = [ $ ;'(*!dustem_params).G0', $ ;G0 ;'(*!dustem_params).grains(0).mdust_o_mh', $ ;PAH0 mass fraction ;'(*!dustem_params).grains(1).mdust_o_mh', $ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil 'dustem_create_continuum_2', $ ;Intensity of NIR continuum 'dustem_create_stellar_population_O63', $ ; distance to the O6 stellar population 'dustem_create_stellar_population_O65', $ ; number of stars of the O6 stellar population 'dustem_create_stellar_population_B43', $ ; distance to the B4 stellar population 'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population ] ;p_truth = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001,13.,3.,10.,1.] ;p_truth = [5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001,13.,3.,10.,1.] p_truth = [1.8e-4,2.33e-3,8.27e-3,0.001,13.,3.,10.,1.] ;UNCOMMENT THIS WHOLE BLOCK AFTER THE TESTING ; p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,70,20,0.2] ; ;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,25] ;UNCOMMENT THIS BLOCK AFTER THE TESTING ;p_truth=[1.,7.8000E-04,7.8000E-04,2E-03,5e-3,70,20,0.2] ;TESTING THE ISRF PLUGIN ;p_truth=[1.,7.8000E-04,7.8000E-04,2E-03,2e-3,13.,3.,9.,2.] ;p_truth=[7.8000E-04,7.8000E-04,2E-03,2e-3,13.,3.];,9.,2.] ;p_truth=[7.8000E-04,7.8000E-04,2E-03,10.,8.] ;p_truth=[10.,8.] ;=== set initial value of parameters you want to fit iv=p_truth ;exact solution ;UNCOMMENT THIS BLOCK AFTER THE TESTING ; iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5,5e-2,5,5,0.2];shifted from solution ; ;iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5,5];shifted from solution ;TESTING THE ISRF PLUGIN ;iv=p_truth+[0,3.e-5,1.e-5,-2.e-5,9e-4,-3.,-2,2.,1.] ;iv=p_truth+[3.e-5,1.e-5,-2.e-5,9e-4,-3.,-2];,2.,1.] ;iv=p_truth+[3.e-5,1.e-5,-2.e-5,-2,-2] ;iv=p_truth+[-2,-2] ;iv=p_truth+[0.,5.4e-5, 5.4e-5,1.8e-5,2.33e-4,8.27e-4,1.0e-4,-2.,-2,3.,3.] ;iv=p_truth+[5.4e-5, 5.4e-5,1.8e-5,2.33e-4,8.27e-4,1.0e-4,-2.,-2,3.,3.] iv=p_truth+[1.8e-5,2.33e-4,8.27e-4,1.0e-4,-2.,-2,3.,3.] ;UNCOMMENT THIS BLOCK AFTER THE TESTING ; ; ulimed=[0 ,0 ,0 ,0 ,0 ,0 , 0 ,0] ;Polarization angles exist in the range [0-180] degress ; ; ;ulimed=[0 ,0 ,0 ,0 ,0] ;Polarization angles exist in the range [0-180] degress ; ; ; llimed=[1 , 1 ,1 ,1 ,1 ,1 ,1 ,1] ; ;llimed=[1 , 1 ,1 ,1 ,1] ; ; llims =[0. ,0. ,0. ,0. , 0. ,0. ,0. ,0.] ; ;llims =[0. ,0. ,0. ,0. ,0.] ;TESTING THE ISRF PLUGIN Npar=n_elements(pd) ulimed=replicate(0,Npar) llimed=replicate(1,Npar) ;llims =[1e-5 ,0. ,0. ,0. ,0. ,1. ,0. ,1. ,0. ] ;llims =[0. ,0. ,0. ,0. ,1. ,0. ];,1. ,0. ] ;llims =[0. ,0. ,0. ,1. ,1. ] ;llims =[1. ,1. ] ;llims = [1.0e-5 ,0. ,0. ,0. ,0. ,0. ,0. ,1. ,1. ,1. ,1. ] ;llims = [0. ,0. ,0. ,0. ,0. ,0. ,1. ,1. ,1. ,1. ] llims = [0. ,0. ,0. ,0. ,1. ,1. ,1. ,1. ] ;== 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 ;dustem_init_fixed_params,fpd,fiv ;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin' ;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin' ;(*!dustem_params).grains[4].TYPE_KEYWORDS='pol-plaw-ed' dustem_init_plugins, pd help,*!dustem_fit help,!run_pol sed.spec=dustem_compute_sed(p_truth,sst,_extra=extra) ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code sed.error=(sed.spec*0.1);>0.1 ;UNCOMMENT THIS BLOCK AFTER THE TESTING ;polsed.spec=dustem_compute_polsed(p_truth,sstp,_extra=extra,dustem_qsed,dustem_used) ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code polsed.spec=dustem_compute_polsed(p_truth,sstp,_extra=extra) polsed.error=(polsed.spec*0.1) ; ;creating fake polext ; ; polext.spec=dustem_compute_polext(p_truth,ssttp,dustem_qext,dustem_uext) ; polext.error=(polext.spec*0.1) ;UNCOMMENT THIS BLOCK AFTER THE TESTING ; ;SET qsed and used data - comment if extinction polarization data is to be used ; qsed.spec=dustem_qsed ; qsed.error=qsed.spec*0.1 ; ; ; used.spec=dustem_used ; used.error=used.spec*0.1 ; ;SET qext and uext data - comment if polarizarion data is to be used ; qext.spec=dustem_qext ; qext.error=qext.spec*0.1 ; ; ; uext.spec=dustem_uext ; uext.error=uext.spec*0.1 ; ; ;=== SET THE OBSERVATION STRUCTURE ind=where(sed.wave GE 0.1) ;UNCOMMENT THIS BLOCK AFTER THE TESTING ; st=dustem_set_data(polsed=polsed[ind],sed=sed[ind],qsed=qsed[ind],used=used[ind]) ; comment if second line is to be used ; ;st=dustem_set_data(sed=sed[ind],polext=polext[ind],qext=qext[ind],uext=uext[ind]) st=dustem_set_data(polsed=polsed[ind],sed=sed[ind]) ;=== RUN fit tol=1.e-16 xtol=1.e-16 Nitermax=30; just an initialization ;for i=0L,n_tags(!dustem_data)-1 do begin ;if ptr_valid(!dustem_data.(i)) then Nitermax+=1 ;endfor ;Nitermax=4 ;maximum number of iteration. This is the criterion which will stop the fit procedure xrange=[1.,2.e4] yrange=[1e-4,1.e3] loadct,13 ;!y.range=[1e-8,100] ;This is to ajust plot range from outside the routine t1=systime(0,/sec) res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []') t2=systime(0,/sec) print,(res-p_truth)/p_truth*100. ;help,(*!dustem_fit),/STRUCTURE ;=== SAVE FIT RESULTS file_out='/tmp/DUSTEM_polsed_fit_example.sav' dustem_save_system_variables,file_out message,'Saved '+file_out,/continue stop ;====================================== ;====You can exit IDL here and re-enter - THIS PART APPEARS TO BE FAULTY. NEEDS INVESTIGATION ;====================================== file='/tmp/DUSTEM_polsed_fit_example.sav' dustem_restore_system_variables,file ;=== Plot best fit win=1 window,win & win=win+1 xrange=[1.,2.e4] yrange=[1e-4,1.e3] ;=== recover best fit values res=*(*!dustem_fit).current_param_values chi2=(*!dustem_fit).chi2 rchi2=(*!dustem_fit).rchi2 ;errors=*(*!dustem_fit).current_param_errors errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values) ;=== Plot best fit tit='DUSTEM polarization dust 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/'+'Last_dustem_spinning_fit.ps';'/Docs/Figures/'+'Last_dustem_spinning_fit.ps' device,filename=ps_file,/color ENDIF dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog,/pol 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