PRO dustem_fit_sed_orion,ps=ps,png=png dustem_define_la_common ;=== initialise dustem use_mode='G17_MODELC';'DBP90';'G17_MODELA';'THEMIS';;'USER_MODEL';'THEMIS';;'THEMIS';'G17_MODELA';'THEMIS';'NY_MODELA';'G17_MODELA';'DBP90';'USER_MODEL';'THEMIS';'NY_MODELA';'G17_MODELA';'DBP90';'COMPIEGNE_ETAL2010' ;,/pol dustem_init,mode=use_mode,kwords=['logn-chrg-spin','plaw-pol','plaw-pol'],/pol;,kwords=['plaw-ed','logn','logn'];,/pol;kwords=['plaw-ed-chrg-spin'],/pol;,kwords=['plaw-ed'];kwords=['plaw-ed-chrg-spin'];,/pol; ; ;stop !dustem_verbose=1 !dustem_show_plot=1 !dustem_nocatch=1 ;specify level lev=['300.000','400.000','500.000','600.000','800.000','1000.00','1200.00','1500.00','3000.00','5000.00','10000.0','15000.0'] indxlev=6 ; meaning level 400, 1 is for 600 etc...;THIS IS FOR THE PLOTTING OF ONE LEVEL ONLY msk=['mask1','mask2','mask3'] ;this is used to be able to plot the new data that JP generated. ;1,2,3 correspond to the North-West, Center and South-East part of the OrionA integral shape filament msk=msk(1) ;mask2 center of ORION_A plckprd=['PR3','PR4','SROLL2'] plckprd=plckprd(1) ;PR4 polangmatch=['ang5','ang10','ang15'] polangmatch=polangmatch(1) ;best match. Lowest difference. co_corr = ['noCO0','noCO1','noCO2'] co_corr = co_corr(2) ;opening corresponding sed file smooth_str='_smoothedLFI3' ; dir_seds='/Users/ilyeschoubani/Downloads/SC110012_V2p2_PRO_5000001211_2/Orion/seds/' ;new SEDS have a problem ; ;dir_seds='/Users/ilyeschoubani/Downloads/SC110012_V2p2_PRO_5000001211/Orion/seds/' dir_seds = '/Users/ilyeschoubani/Downloads/SC110012_V2p2_PRO_5000001210/Orion/seds/' ;with the mask - new batch still has the error ;file_sed=dir_seds+'extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'_'+msk+'_'+co_corr+'.xcat' ;without the mask - new batch does not have the error file_sed=dir_seds+'extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'_'+co_corr+'.xcat' ;;THIS WAS PRIOR;file_sed=dir_seds+'extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'.xcat' data_struct=read_xcat(file_sed,/silent) ;This is the new general structure ;stop ; ;===reading and applying the HCD: openr,unit,file_sed,/get_lun str='' WHILE not eof(unit) DO BEGIN readf,unit,str if strtrim(strmid(str,13,3),2) EQ '857' then *!dustem_HCD=float(strtrim(strmid(str,21),2))/5.25*1.00E+26 ENDWHILE close,unit free_lun,unit ;=================== ;stop ;generating fake errors for SED and Q/USED - here 10% ;stop data_struct.sigmaII = (double((data_struct.StokesI*0.1)^2)) > 1d-20 ;this is false. Squaring it does not make sense. The square root is the one that should be taken ;NOT THE REAL REASON NOW;because of a bigger error in the variance computation in the xcat file. CHIPASS has an error value for StokesQ and StokesU but not for their associated flux. ;CHIPASS filter isn't well defined ;I looked it up and found a bandpass of 64 MHz (in 1024 channels). data_struct(-1).sigmauu = la_undef() data_struct(-1).sigmaqq = la_undef() data_struct(-1).sigmaii = la_undef() data_struct(-1).stokesu = la_undef() data_struct(-1).stokesq = la_undef() data_struct(-1).stokesi = la_undef() data_struct(-1).psi = la_undef() data_struct(-1).sigma_psi = la_undef() ; data_struct(-1).stokesi = la_undef() ; tt=where(data_struct.sigmauu EQ -39.8637) ; ;stop ; ; data_struct[tt].sigmauu = la_undef() ; data_struct[tt].sigmaqq = la_undef() ;data_struct.sigma_smallp = (double((data_struct.smallp*0.1)^2)) > 1d-20 idsigq=where(data_struct.sigmaQQ NE la_undef()) idsigu=where(data_struct.sigmaUU NE la_undef()) ;stop ;...... ;DON'T FORGET TO UNCOMMENT THIS. ;data_struct[idsigq].sigmaQQ = (double((data_struct[idsigq].StokesQ*0.1)^2)) > 1d-20 ;data_struct[idsigu].sigmaUU = (double((data_struct[idsigu].StokesU*0.1)^2)) > 1d-20 ;#######FOR THE FITTING OF Q AND U ;locating the negative values in stokesQ data so that filters are omitted from the fit indq=where(data_struct.stokesQ gt 0 and data_struct.stokesQ ne la_undef()) rmvq=(data_struct(indq).filter) ;locating the negative values in StokesU data so that filters are omitted form the fit indu=where(data_struct.stokesU lt 0 and data_struct.stokesU ne la_undef()) rmvu=(data_struct(indu).filter) rmvpsi=[rmvq,rmvu] ;removing LFI2 manually remove,1,rmvpsi ;stop ;BECAUSE SPASS band has no errors: ;NOW IT DOES ; (data_struct(-2).sigma_largeP)=(data_struct(-3).sigma_largeP) ; (data_struct(-2).sigma_smallp)=(data_struct(-3).sigma_smallp) ; (data_struct(-2).sigmaqq)=(data_struct(-3).sigmaqq) ; (data_struct(-2).sigmauu)=(data_struct(-3).sigmauu) ; (data_struct(-2).sigma_psi)=(data_struct(-3).sigma_psi) ;stop fit_struct=data_struct ;stop ;uncomment this block afterwards - because you are now just fitting the sed. ; fit_struct=dustem_mask_data(fit_struct,omit=list(rmvq,rmvu,rmvpsi,rmvpsi,rmvpsi),data_set=['QSED','USED','PSI','POLSED','POLFRAC']) ;works ; rmvq=(data_struct(-2).filter) ; ; ; fit_struct=dustem_mask_data(fit_struct,omit=list(rmvq,rmvq,rmvq,rmvq,rmvq),data_set=['QSED','USED','PSI','POLSED','POLFRAC']) ; ; ; fit_struct=dustem_mask_data(fit_struct,omit=rmvq,data_set=['QSED']) ; fit_struct=dustem_mask_data(fit_struct,omit=rmvu,data_set=['USED']) ;stop ;since we're using q and u I will have to null polsed and polfrac. Otherwise set_data will return an error ;fit_struct.largeP=la_undef() ;fit_struct.smallp=la_undef() ;doing the same for the errors as dustem_check will notice some erros have undefined associated data points. ;fit_struct.sigma_largeP=la_undef() ;fit_struct.sigma_smallp=la_undef() ;because I am now fitting only the sed ; fit_struct.sigmaQQ = la_undef() ; fit_struct.sigmaUU = la_undef() ;stop ;show structure show_st=data_struct ;stop ; show_st.largeP=la_undef() ; show_st.sigma_largeP=la_undef() ; show_st.stokesQ=la_undef() ; show_st.sigmaQQ=la_undef() ; show_st.stokesU=la_undef() ; show_st.sigmaUU=la_undef() ;################################################# ;params_TH=[5.44E-05,5.60E-04,1.22E-04,3.0,4.13E+06,2.25E+00];,0.30,53] params_G17C=[7.50E-05,1.265E-03,5.260E-04,2.739E+02,1.11E03,6.34,1.00E+07,2.73E+00,58]; dist=6.43 ;params_DBP=[1.06E-05,7.00E-04,2.32E-03,6.43,1.00E+07,2.57E+00] ;=== Set which parameters you want to fit pd = [ $ ;'(*!dustem_params).gas.G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh', $ '(*!dustem_params).grains(1).mdust_o_mh', $ '(*!dustem_params).grains(2).mdust_o_mh', $ ;'(*!dustem_params).grains(3).mdust_o_mh', $ ;'(*!dustem_params).grains(4).mdust_o_mh', $ '(*!dustem_params).gas.Tgas',$ '(*!dustem_params).gas.nh',$ 'dustem_plugin_stellar_population_O6V3',$ ;'dustem_plugin_mbbdy_isrf_1',$ ;'dustem_plugin_mbbdy_isrf_2', $ ;'dustem_plugin_mbbdy_1',$ ;'dustem_plugin_mbbdy_2', $ ;'dustem_plugin_mbbdy_3', $ ;'dustem_plugin_mbbdy_isrf_3', $ ;'dustem_plugin_mbbdy_isrf_4', $ ;'dustem_plugin_synchrotron_1', $ ;'dustem_plugin_synchrotron_2', $ ;'dustem_plugin_synchrotron_3', $ ;'dustem_plugin_synchrotron_4', $ 'dustem_plugin_freefree_1', $ ;4--->6 'dustem_plugin_freefree_2', $ ;5 ;'dustem_plugin_freefree_3', $ ;'dustem_plugin_modify_dust_pol_1',$ 'dustem_plugin_modify_dust_pol_2'$ ;'dustem_plugin_mdfy_spinning_pol_1'$ ;'dustem_plugin_continuum_2', $ ;'dustem_plugin_continuum_1' $ ;'dustem_plugin_mdfy_spinning_pol_1', $ ;'dustem_plugin_mdfy_spinning_pol_2' $ ] iv=params_G17C Npar=n_elements(pd) ulimed=replicate(0,Npar);[0] ulims=replicate(0.0,Npar) ulimed[6]=1 ;uncomment ;ulimed[6]=1 ;ulimed[9]=1 ulims[6]=1.00E+07 ;ulims[8]=5.0E-03 ;uncomment ;ulims[6]=9.00E-01;becomes 9 ;stop llimed=replicate(1,Npar);[1] ;llims=[1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20,1.00e-20];,1.00e-20,1.00e-20] llims=replicate(5.3e-6,Npar);[10] ;llims[0]=1.0e-9 ;llims[3]=1.0 llims[5]=3.0 llims[6]=1.00E+04 llims[7]=1.00E-04 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims;,tied=tied ;stop dustem_init_plugins,pd stop ; fpd = [ $ ; '(*!dustem_params).grains(1).mdust_o_mh'$ ; ;'(*!dustem_params).grains(1).mdust_o_mh'$ ; ;'(*!dustem_params).grains(4).mdust_o_mh'$ ; ] ; fiv = [1.00E-20] ; ; dustem_init_fixed_params, fpd, fiv ;=== SET THE OBSERVATION STRUCTURE ind=where(data_struct.wave GE 0.1) st_ie=dustem_set_data(m_fit=fit_struct[ind],m_show=show_st[ind],rchi2_weight=rchi2_weight,f_HI=f_HI) ; comment if second line is to be used ;stop ;=== RUN fit tol=1.e-16 xtol=1.e-16 Nitermax=60;0; just an initialization xr=[1.,5e5] yr=[5e-8,1.00e6] loadct,13 ;!y.range=[1e-8,100] ;This is to adjust plot range from outside the routine yyy="restrng=textoidl('[MJy/sr]_{"+strmid(smooth_str,9)+'|'+plckprd+"}')" yyy=execute(yyy) title=textoidl('OrionA I_{240\mum} > ')+strtrim(lev[indxlev],2)+restrng sttr=' I_{\nu} (MJy/sr) for N_{H}='+strmid(string(*!dustem_HCD,format='(1E10.2)'),2)+'H/cm^2' ytit=textoidl(sttr) ;xtit='' t1=systime(0,/sec) res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,xtit=xtit,ytit=ytit,title=title,_extra=_extra); t2=systime(0,/sec) ;help,(*!dustem_fit),/STRUCTURE ;=== SAVE FIT RESULTS ;NB: THIS LINE IS WRONG BECAUSE I AM NOT USING THE ANGLE CORRELATION CRITERIA sedid='extracted_sed_I240gt'+strtrim(lev[indxlev],2)+smooth_str+'_'+plckprd+'_masked_'+polangmatch+'_'+msk file='/tmp/DUSTEM_lastfit_orionA_'+sedid+'.sav' dustem_save_system_variables,file message,'Saved '+file,/continue ;====================================== ;====You can exit IDL here and re-enter - THIS PART APPEARS TO BE FAULTY. NEEDS INVESTIGATION ;====================================== dustem_restore_system_variables,file ;=== Plot best fit ; win=1 ; ; window,win & win=win+1 ;=== recover best fit values res=*(*!dustem_fit).current_param_values dustemwrap_plot,res,stp,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=title ; 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 ; ; ps_file=!dustem_wrap_soft_dir+'Docs/'+'Last_dustem_fit.ps' ; dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=title,yr=yr,xr=xr,/ysty,/xsty,res=res,chi2=chi2,rchi2=rchi2,/xlog,/ylog,/pol,ps=ps_file,png=png,fpol=fpol print,'dustem_mpfit_sed executed in ',t2-t1,' sec' the_end: END