From 1355825cf3651fa87141046849d32ebd8d10be92 Mon Sep 17 00:00:00 2001 From: ichoubani Date: Mon, 18 Jul 2022 14:19:14 +0200 Subject: [PATCH] General update --- src/idl/dustem_activate_plugins.pro | 9 --------- src/idl/dustem_fit_sed_readme.pro | 12 ++++++------ src/idl/dustem_init.pro | 2 +- src/idl/dustem_init_plugins.pro | 6 +++++- src/idl/dustem_mpfit_data.pro | 19 ++++++++++++++++++- src/idl/dustem_mpfit_run.pro | 326 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- src/idl/dustem_plot_fit_sed.pro | 11 ++++++----- src/idl/dustem_read_all_web3p8.pro | 2 +- src/idl/dustem_set_data.pro | 90 +++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------- src/idl/dustemcgwin_dataset.pro | 182 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------- src/idl_misc/JPBLib_for_Dustemwrap/General/modify_extra.pro | 1 + 11 files changed, 355 insertions(+), 305 deletions(-) diff --git a/src/idl/dustem_activate_plugins.pro b/src/idl/dustem_activate_plugins.pro index 374a605..4ae7aa2 100755 --- a/src/idl/dustem_activate_plugins.pro +++ b/src/idl/dustem_activate_plugins.pro @@ -82,15 +82,6 @@ FOR i=0L,n_elements(p_min)-1 DO BEGIN ;==============Assigning data to the initialized plugin data arrays============== - - print, 'ftn is:' - print, ftn - print, 'key is:' - print, index - print, 'val is:' - print, value - - if strmid(ftn,0,13) eq 'dustem_plugin' then begin diff --git a/src/idl/dustem_fit_sed_readme.pro b/src/idl/dustem_fit_sed_readme.pro index 5de56a2..6fc1fce 100755 --- a/src/idl/dustem_fit_sed_readme.pro +++ b/src/idl/dustem_fit_sed_readme.pro @@ -343,7 +343,7 @@ CASE use_model OF END 'THEMIS':BEGIN pd = [ $ - '(*!dustem_params).G0', $ ;G0 + '(*!dustem_params).G0', $ h;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ '(*!dustem_params).grains(1).mdust_o_mh',$ '(*!dustem_params).grains(2).mdust_o_mh', $ @@ -358,7 +358,7 @@ CASE use_model OF ENDCASE -dustem_init,model=use_model,pol=polarization +dustem_init,model=use_model;,pol=polarization !dustem_nocatch=1 !dustem_verbose=1 @@ -382,8 +382,8 @@ ind=where(spec.instru EQ 'FIRAS',count) IF count NE 0 THEN spec(ind).sigmaII=(0.2*spec(ind).StokesI)^2 ;=== SET THE OBSERVATION STRUCTURE -st=dustem_set_data(sed=spec) - +st=dustem_set_data(st_sed_fit=spec);sed=spec) +;stop ;== SET THE FITTED PARAMETERS dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims @@ -401,8 +401,8 @@ use_Nitermax=5 ;maximum number of iteration. This is the criterium which will st IF keyword_set(itermax) THEN use_Nitermax=itermax t1=systime(0,/sec) -yr=[1e-4,100] -xr=[1,6e4] +yr=[1.00e-4,1.00E2] +xr=[1.00E0,6.00e4] ;legend_xpos=0.6 ;legend_ypos=0.8 tit='Spectral Energy Distribution (Running)' diff --git a/src/idl/dustem_init.pro b/src/idl/dustem_init.pro index 5c9714c..fcc8abf 100755 --- a/src/idl/dustem_init.pro +++ b/src/idl/dustem_init.pro @@ -205,7 +205,7 @@ defsysv, '!dustemcgwin_ncmds', { $ ;Keeping track of number of commands run. ext: {pl:0, nrm:0}, $ polext: {pl:0, nrm:0}, $ polsed: {pl:0, nrm:0}, $ - polfrac: {pl:0, nrm:0}, $ + polfrac: {pl:0}, $ psi_em: {pl:0}, $ psi_ext: {pl:0}, $ qsed: {pl:0, nrm:0}, $ diff --git a/src/idl/dustem_init_plugins.pro b/src/idl/dustem_init_plugins.pro index cf579c0..442e525 100644 --- a/src/idl/dustem_init_plugins.pro +++ b/src/idl/dustem_init_plugins.pro @@ -1,5 +1,9 @@ PRO dustem_init_plugins, pd, fpd +;NB: +;#Fixed parameters still haven't been taken care of +;#There is use of pointers when the need does not present itself. This implies modifying the rest of the code. + ;INITIALIZING THE SCOPES DATA ARRAYS-------------------- @@ -79,7 +83,7 @@ IF Nplugins NE 0 THEN BEGIN message, 'varvar is:',/continue print, varvar - defsysv, '!dustem_scope', ptr_new(varvar) + defsysv, '!dustem_scope', ptr_new(varvar) ENDIF ELSE BEGIN ;case of one plugin diff --git a/src/idl/dustem_mpfit_data.pro b/src/idl/dustem_mpfit_data.pro index 9268932..6e3c7ae 100755 --- a/src/idl/dustem_mpfit_data.pro +++ b/src/idl/dustem_mpfit_data.pro @@ -58,7 +58,23 @@ IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax ; Build array concatenating observational datas ; Order : that of !dustem_data (sed, ext, polext, polsed) -ind_data_tag=dustem_data_tag(count=count_data_tag) + +if !run_pol then begin + + ind_data_tag=dustem_data_tag(count=count_data_tag) + tagnames=tag_names(!dustem_data) + testpol = tagnames EQ 'POLSED' or tagnames EQ 'POLEXT' or tagnames EQ 'POLFRAC' or tagnames EQ 'POLSED' or tagnames EQ 'PSI_EM' or tagnames EQ 'PSI_EXT' + data_ind=where(~testpol,ctestpol) + if ctestpol ne 0 then begin + match,ind_data_tag,data_ind,ind1,ind2 ;not using match2 because when matching indices, if the first is 0 it doesn't count which is wrong. + ;I also only need the points that match. + ;testind1 = where(ind1 ne -1,ctind1) + ;if ctind1 ne 0 then ind_data_tag = ind_data_tag[testind1] + ind_data_tag = ind_data_tag[ind1] + count_data_tag = n_elements(ind_data_tag) + endif +endif else ind_data_tag=dustem_data_tag(count=count_data_tag) + fields=['wav','values','sigma'] FOR j=0,n_elements(fields)-1 DO BEGIN @@ -104,6 +120,7 @@ p_min = dustem_mpfitfun(func_name,wav,values,sigma, $ parinfo=(*!dustem_parinfo), perror=perror, yfit=yfit, $ bestnorm=bestnorm, dof=dof,quiet=quiet, status=status, xtol=xtol, niter=niter_completed, ftol=use_tol, $ maxiter=use_Nitermax,gtol=gtol,errmsg=errmsg,functargs=_extra,nocatch=nocatch) + message,errmsg,/info message,'mpfitfun Status='+strtrim(status,2),/info p_dim=p_min*(*(*!dustem_fit).param_init_values) diff --git a/src/idl/dustem_mpfit_run.pro b/src/idl/dustem_mpfit_run.pro index ba4be12..3aeb0ac 100755 --- a/src/idl/dustem_mpfit_run.pro +++ b/src/idl/dustem_mpfit_run.pro @@ -112,26 +112,12 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta (*!dustem_fit).chi2=chi2_sed (*!dustem_fit).rchi2=rchi2_sed - ; Plot if needed - IF !dustem_show_plot NE 0 THEN BEGIN - - ;window,win,title='DUSTEM WRAP (SED)' - - ;cgwindow,WMULTI=[2,2,1],/CURRENT - - - - ;stop - ;dustemwrap_plott,p_dim,st,dustem_sed,dustem_polsed,dustem_qsed,dustem_polfrac,_extra=_extra - ;dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2,fpol=fpol - - ENDIF - - dustem_out = [ dustem_out, dustem_sed ] + + if isa(!dustem_data.sed) then dustem_out = [ dustem_out, dustem_sed ] END - 'EXT' : BEGIN + 'EXT' : BEGIN ;STILL HAVEN'T DONE THE EXTINCTION PART ;ADDING PLUGIN TO SPECTRUM---------------- IF ptr_valid(*!dustem_plugin) THEN BEGIN @@ -266,67 +252,65 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN END - ; 'POLFRAC': BEGIN -; if !run_lin then begin -; dustem_polsed = dustem_compute_polsed(p_dim,st) ; this quantity hasn't been computed yet that is why this line is present here. -; ind_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt) ;locating the filters of polfrac=smallp in SED xcat -; ind_fspec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',count_fspec) ;locating the spectrum points in polfrac data. Later testing will have to be on the wavelengths -; ; -; ; ;dustem_sed and dustem_polsed might not have the same length thus this block. -; ; -; If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin -; if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin -; dustem_polsed_x = dustem_polsed -; dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified -; ;matching the filter data points -; IF count_filt ne 0 then begin -; for i=0L,count_filt-1 do begin -; j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt) -; if testfilt ne 0 then dustem_sed_x(ind_filt(i)) = dustem_sed(j(0)) -; endfor -; endif -; ;matching the spectrum data points -; IF count_fspec ne 0 then begin -; for i=0L,count_fspec-1 do begin -; j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec) -; if testspec ne 0 then dustem_sed_x((ind_fspec(i))) = dustem_sed(j(0)) -; ; endfor -; endif -; endif ELSE begin -; -; dustem_polsed_x = dustem_sed -; dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified -; ;matching the filter data points -; IF count_filt ne 0 then begin -; for i=0L,count_filt-1 do begin -; j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt) -; if testfilt ne 0 then dustem_polsed_x(ind_filt(i)) = dustem_polsed(j(0)) -; endfor -; endif -; ;matching the spectrum data points -; IF count_fspec ne 0 then begin -; for i=0L,count_fspec-1 do begin -; j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec) -; if testspec ne 0 then dustem_polsed_x((ind_fspec(i))) = dustem_polsed(j(0)) -; endfor -; endif -; -; ENDELSE -; endif ELSE BEGIN -; dustem_polsed_x = dustem_polsed -; dustem_sed_x = dustem_sed -; ENDELSE +;!!!!!!!!!!!!!YOU NEED TO RUN 'COMPUTE_SED' AGAIN because you have no idea if the sed block will be run. + + 'POLFRAC': BEGIN + if !run_lin then begin + if ~isa(dustem_sed) then dustem_sed = dustem_compute_sed(p_dim,st) + dustem_polsed = dustem_compute_polsed(p_dim,st) ; this quantity hasn't been computed yet that is why this line is present here. + ind_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt) ;locating the filters of polfrac=smallp in SED xcat + ind_fspec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',count_fspec) ;locating the spectrum points in polfrac data. Later testing will have to be on the wavelengths ; -; dustem_polfrac = dustem_polsed_x/dustem_sed_x -; chi2_polfrac = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2) -; n_polfrac = n_elements((*!dustem_data.polfrac).values) -; rchi2_polfrac = chi2_polfrac / n_polfrac -; wrchi2_polfrac = rchi2_polfrac * !fit_rchi2_weight.polfrac -; message,"chi2 POLFRAC = "+strtrim(chi2_polfrac,2)+" rchi2 POLFRAC = "+strtrim(rchi2_polfrac,2),/continue ;," weighted rchi2 POLFRAC=",wrchi2_polfrac -; dustem_out = [ dustem_out, dustem_polfrac ] -; ;fpol=dustem_polfrac ; problem because extra structure does not get updated so you actually need the 'extra' procedure you haven't finihed coding. -; endif - ; END +; ;dustem_sed and dustem_polsed might not have the same length thus this block. +; + If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin + if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin + dustem_polsed_x = dustem_polsed + dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified + ;matching the filter data points + IF count_filt ne 0 then begin + for i=0L,count_filt-1 do begin + j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt) + if testfilt ne 0 then dustem_sed_x(ind_filt(i)) = dustem_sed(j(0)) + endfor + endif + ;matching the spectrum data points + IF count_fspec ne 0 then begin + for i=0L,count_fspec-1 do begin + j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec) + if testspec ne 0 then dustem_sed_x((ind_fspec(i))) = dustem_sed(j(0)) + endfor + endif + endif ELSE begin + + dustem_polsed_x = dustem_sed + dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified + ;matching the filter data points + IF count_filt ne 0 then begin + for i=0L,count_filt-1 do begin + j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt) + if testfilt ne 0 then dustem_polsed_x(ind_filt(i)) = dustem_polsed(j(0)) + endfor + endif + ;matching the spectrum data points + IF count_fspec ne 0 then begin + for i=0L,count_fspec-1 do begin + j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec) + if testspec ne 0 then dustem_polsed_x((ind_fspec(i))) = dustem_polsed(j(0)) + endfor + endif + + ENDELSE + endif ELSE BEGIN + dustem_polsed_x = dustem_polsed + dustem_sed_x = dustem_sed + ENDELSE + + dustem_polfrac = dustem_polsed_x/dustem_sed_x + + ENDIF + ;stop + END ;----------------------------------------------------------------- @@ -335,40 +319,16 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED IF !run_lin THEN BEGIN - ;make polsed only used for when polsed is fitted. Otherwise it's the fitting of Q_SED AND USED that fits the polarization data + dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron) - ;chi2_polsed = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2) - ;n_polsed = n_elements((*!dustem_data.polsed).values) - ;rchi2_polsed = chi2_polsed / n_polsed - ;wrchi2_polsed = rchi2_polsed * !fit_rchi2_weight.polsed - ;message,"chi2 POLSED = "+strtrim(chi2_polsed,2)+" rchi2 POLSED = "+strtrim(rchi2_polsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed - ;print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2 - ;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma - - ;dustem_out = [ dustem_out, dustem_polsed ] ; classic command - ;----------------------------------------------- - - ; Plot if needed -; IF !dustem_show_plot NE 0 THEN BEGIN -; win+=1 -; ; -; xr=[10,8e3] -; yr=[1e-34,1e-28] -; -; dustem_plot_polsed,st,p_dim,dustem_polsed,aligned=st_model.align.grains.aligned,win=win -; ; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,/donotclose,multi=1 -; ; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,multi=2 -; -; -; -; ENDIF ENDIF END 'QSED': BEGIN toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec) + ;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON ;dustem_qsed(-1)=dustem_qsed(-2) ; testspass1 = ((*!dustem_data.qsed).filt_names)(-1) eq 'SPASS1' @@ -378,118 +338,24 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN n_qsed = n_elements((*!dustem_data.qsed).values) rchi2_qsed = chi2_qsed / n_qsed wrchi2_qsed = rchi2_qsed * !fit_rchi2_weight.qsed - dustem_out = [ dustem_out, dustem_qsed ] ; classic command + if isa(!dustem_data.qsed) then dustem_out = [ dustem_out, dustem_qsed ] ; classic command message,"chi2 QSED = "+strtrim(chi2_qsed,2)+" rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2 - ;stop - - ;stop -; if !dustem_show_plot NE 0 then begin -; fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40 -; fact1 = 1.00E-20;(3.e10/(1.e-4*((*!dustem_data.qsed).wav)))/1.d40 ; for data plotting -; -; ; fact=1/(4.*!pi)/(3.e8/1.e-6/(st.polsed).wav)*1.00e17 -; ; fact1=1.E-20 -; ; -; win+=1 -; window ,win ,title='DUSTEM WRAP (QSED)' -; titstq='Stokes Q Polarization Intensity' -; ytitstq=textoidl('Q_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2') -; xtit=textoidl('\lambda (\mum)') -; xr=[10,1e6] -; -; ;normpol1=(*!dustem_data.qsed).values*fact1 -; yr=[1e-28,1e-17] -; ;yr=[1e-34,1e-23] -; indq=where(Q_spec lt 0, countq) -; ;indq=where((*!dustem_data.qsed).values lt 0, countq) -; -; ylog=1 -; if countq ne 0 then begin -; ; ylog=0 -; ; yr=[-1E-23,-1e-19];e-31 -; Q_spec=-Q_spec -; dustem_qsed =-dustem_qsed -; ((*!dustem_data.qsed).values)=-((*!dustem_data.qsed).values) -; endif -; normpol=Q_spec*fact -; normpol1=dustem_qsed*fact1 -; -; -; cgplot,st.polsed.wav,Q_spec*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)' -; cgoplot,st.polsed.wav,Q_spec*fact,color='black' -; cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,color='Tomato',psym=16,symsize=1,thick=2 -; err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1,color=cgColor('Tomato') -; cgoplot,(*!dustem_data.qsed).wav,dustem_qsed*fact1,psym=6,color='red',symsize=2 -; -; cgplot,st.polsed.wav,Q_spec*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1 -; cgoplot,st.polsed.wav,Q_spec*fact/normpol,color='black' -; cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,color='Tomato',psym=16,symsize=1,thick=2 -; err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1/normpol1,color=cgColor('Tomato') -; -; endif END 'USED': BEGIN - - ;toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec) - ;stop + ;THIS IMPLIES THAT QSED HAS BEEN RAN BEFORE chi2_used =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2) n_used = n_elements((*!dustem_data.used).values) rchi2_used = chi2_used / n_used wrchi2_used = rchi2_used * !fit_rchi2_weight.used - dustem_out = [ dustem_out, dustem_used] ; classic command + if isa(!dustem_data.used) then dustem_out = [ dustem_out, dustem_used] ; classic command message,"chi2 USED = "+strtrim(chi2_used,2)+" rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed print,"chi2 per wl: ",((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2 - dustemwrap_plot,p_dim,st,dustem_sed,dustem_qsed,dustem_used,dustem_polsed,_extra=_extra - - ;if !dustem_show_plot ne 0 then begin - ;fact = 1.00E-20;(3.e10/(1.e-4*(st.polsed.wav)))/1.d40 - ;fact1 =1.00E-20;(3.e10/(1.e-4*((*!dustem_data.used).wav)))/1.d40 ; for data plotting - - ;win+=1 - ;window ,win ,title='DUSTEM WRAP (USED)' - ;titstu='Stokes U Polarization Intensity' - ;ytitstu=textoidl('U_\nu (W/m^2/Hz/sr) for N_H=10^{20} H/cm^2') - ;xtit=textoidl('\lambda (\mum)') - ;xr=[10,1e6] - ;normpol2=U_spec*fact - ;normpol3=dustem_used*fact1 - ;yr=[min(normpol2)-min(normpol2)/10,max(normpol2)+max(normpol2)*10] - - ;conditions here are on computed spectra and they should be on the data itself. - ;This even impairs the visualization of the data. Correct this!!!!! - - ;yr=[1e-28,1e-17];[1e-34,1e-23] - ;indu=where(U_spec lt 0, countu) - ;indu=where((*!dustem_data.used).values lt 0, countu) - ;ylog=1 - ;if countu ne 0 then begin - ;ylog=0 - ;yr=[-2e-31,5e-31] - ;endif - - ;cgplot,st.polsed.wav,U_spec*fact,xtit='',ytit=ytitstu,tit=titstu,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black' - ;cgoplot,st.polsed.wav,U_spec*fact,color='black' - ;cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,color='Teal',psym=16,symsize=1,thick=2 - ;err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1,color=cgColor('Teal') - ;cgoplot,(*!dustem_data.used).wav,dustem_used*fact1,psym=6,color='red',symsize=2 - - - - ;cgplot,st.polsed.wav,U_spec*fact/normpol2,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1 - ;cgoplot,st.polsed.wav,U_spec*fact/normpol2,color='black' - ;cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,color='Teal',psym=16,symsize=1,thick=2 - ;err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1/normpol3,color=cgColor('Teal') - ; - ;cgwindow, 'plot', - - ;endif + dustemwrap_plot,p_dim,st,dustem_sed,dustem_qsed,dustem_used,dustem_polsed,dustem_polfrac,_extra=_extra + END -; -; - 'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION @@ -603,41 +469,27 @@ if !fit_rchi2_weight.sed ne 0 then begin total_rchi2 += rchi2_sed endif if !run_pol then begin -; if !fit_rchi2_weight.polext ne 0 then begin -; total_n += n_polext -; total_chi2 += chi2_polext -; total_rchi2 += rchi2_polext -; endif -; if !fit_rchi2_weight.polsed ne 0 then begin -; total_n += n_polsed -; total_chi2 += chi2_polsed -; total_rchi2 += rchi2_polsed -; endif -; if !fit_rchi2_weight.polfrac ne 0 then begin -; total_n += n_polfrac -; total_chi2 += chi2_polfrac -; total_rchi2 += rchi2_polfrac -; endif -if !fit_rchi2_weight.qsed ne 0 then begin - total_n += n_qsed - total_chi2 += chi2_qsed - total_rchi2 += rchi2_qsed -endif -if !fit_rchi2_weight.used ne 0 then begin - total_n += n_used - total_chi2 += chi2_used - total_rchi2 += rchi2_used -endif -if !fit_rchi2_weight.qext ne 0 then begin - total_n += n_qext - total_chi2 += chi2_qext - total_rchi2 += rchi2_qext -endif -if !fit_rchi2_weight.uext ne 0 then begin - total_n += n_uext - total_chi2 += chi2_uext - total_rchi2 += rchi2_uext -endif + + if !fit_rchi2_weight.qsed ne 0 then begin + total_n += n_qsed + total_chi2 += chi2_qsed + total_rchi2 += rchi2_qsed + endif + if !fit_rchi2_weight.used ne 0 then begin + total_n += n_used + total_chi2 += chi2_used + total_rchi2 += rchi2_used + endif + if !fit_rchi2_weight.qext ne 0 then begin + total_n += n_qext + total_chi2 += chi2_qext + total_rchi2 += rchi2_qext + endif + if !fit_rchi2_weight.uext ne 0 then begin + total_n += n_uext + total_chi2 += chi2_uext + total_rchi2 += rchi2_uext + endif endif diff --git a/src/idl/dustem_plot_fit_sed.pro b/src/idl/dustem_plot_fit_sed.pro index 1e8f710..f05db46 100755 --- a/src/idl/dustem_plot_fit_sed.pro +++ b/src/idl/dustem_plot_fit_sed.pro @@ -74,7 +74,7 @@ fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7 use_col_data_filt='blue' ;use_col_sed_spec='red' -use_col_sed_spec='grey' +use_col_sed_spec='Powder Blue' IF not keyword_set(col_sed) THEN BEGIN ;use_col_sed_filt=250 ;red @@ -150,7 +150,8 @@ if keyword_set(fpol) then begin IF !d.name NE 'PS' THEN cgDisplay, 600, 500 cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=ytit,tit='',/ylog,/xlog,/ys,/xs,position=[0.12,0.25,0.96,0.76],xtickformat='(A1)',_extra=_extra,charsize=1.3,/noerase endif -endif ELSE cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=ytit,tit=title,/ylog,/xlog,/ys,/xs,position=[0.12,0.35,0.96,0.90],xtickformat='(A1)',_extra=_extra,charsize=1.3 +endif ELSE cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,ytit=ytit,tit=title,/ylog,/xlog,/ys,/xs,position=[0.12,0.35,0.96,0.90],charsize=1.3,_extra=_extra,xtickformat='(A1)' +;stop ;############################### ;this is where you will add the normalized sed @@ -167,7 +168,7 @@ IF count_spec NE 0 THEN BEGIN xx=((*!dustem_data.sed).wav)[ind_spec] yy=((*!dustem_data.sed).values)[ind_spec]/norm[ind_spec] rms=3.*((*!dustem_data.sed).sigma)[ind_spec]/2./norm[ind_spec] - cgoplot,xx,yy,psym=8,syms=0.5,color=use_col_sed_spec + cgoplot,xx,yy,psym=8,color=use_col_sed_spec IF not keyword_set(no_spec_error) THEN BEGIN cgerrplot,xx,yy-rms,yy+rms,color=use_col_sed_spec ;err_bar,((*!dustem_data.sed).wav)(ind_spec),((*!dustem_data.sed).values)(ind_spec)/norm(ind_spec),yrms=3.*((*!dustem_data.sed).sigma)(ind_spec)/2./norm(ind_spec) @@ -183,7 +184,7 @@ IF count_filt NE 0 THEN BEGIN ;err_bar,((*!dustem_data.sed).wav)(ind_filt),((*!dustem_data.sed).values)(ind_filt)/norm(ind_filt),yrms=3.*((*!dustem_data.sed).sigma)(ind_filt)/2./norm(ind_filt),color=use_col_data_filt cgerrplot,xx,yy-rms,yy+rms,color='Dodger Blue';use_col_data_filt ENDIF -;=== Plot the computed SED - color corrected interpolates +;=== Plot the computed SED IF count_filt NE 0 THEN BEGIN plotsym,8 xx=((*!dustem_data.sed).wav)[ind_filt] @@ -194,7 +195,7 @@ IF count_spec NE 0 THEN BEGIN plotsym,0 xx=((*!dustem_data.sed).wav)[ind_spec] yy=sed[ind_spec]/norm[ind_spec] - cgoplot,xx,yy,color=use_col_sed_filt,psym=6,syms=2 + cgoplot,xx,yy,color='Indian Red',psym=7,syms=1 ENDIF diff --git a/src/idl/dustem_read_all_web3p8.pro b/src/idl/dustem_read_all_web3p8.pro index 933353d..82517ae 100755 --- a/src/idl/dustem_read_all_web3p8.pro +++ b/src/idl/dustem_read_all_web3p8.pro @@ -50,7 +50,7 @@ file=dir_in_dat+(*!dustem_inputs).grain ;st_grains=dustem_read_grain(file,silent=silent) st_grains=dustem_read_grain(file,silent=silent,key_str=key_str,G0=G0) Ngrains=n_elements(st_grains) - + file=dir_in_dat+'ISRF.DAT' st_isrf=dustem_read_isrf(file,silent=silent,Nisrf=Nisrf) diff --git a/src/idl/dustem_set_data.pro b/src/idl/dustem_set_data.pro index 438cbee..bec58fc 100755 --- a/src/idl/dustem_set_data.pro +++ b/src/idl/dustem_set_data.pro @@ -63,7 +63,7 @@ for i=0L,n_tags(!dustem_data)-1 do begin endfor xx=dustem_check_data(data_sed=st_sed_fit,data_ext=st_ext_fit,sed=sed,ext=ext,polext=polext,polsed=polsed,polfrac=polfrac,psi_em=psi_em,qsed=qsed,used=used,qext=qext,uext=uext,psi_ext=psi_ext) - +;stop ;outputs whatever the input structure has, into organized seperate structures formatted for the fitting process. ;If SED/EXT data is present in the structure it is fitted automatically. @@ -91,42 +91,63 @@ if !run_pol then begin teststks = isa(qsed) and isa(used) - IF (not teststks) then message, 'Stokes parameters are not set correctly. Please check the input structure. Aborting...' + IF ~teststks then message, 'Stokes parameters are not set correctly. Please check the input structure. Aborting...' ;EMISSION - if teststks then nowfitting = 'QSED and USED' - message, 'Polarization component(s) to fit: '+nowfitting ,/continue + message, 'Polarization component(s) to fit: QSED and USED' ,/continue IF isa(qsed) THEN BEGIN !dustem_data.qsed = ptr_new(qsed) if keyword_set(f_HI) then if f_HI gt 0 then $ (*!dustem_data.qsed).values *= f_HI ENDIF + if isa(used) then BEGIN !dustem_data.used = ptr_new(used) if keyword_set(f_HI) then if f_HI gt 0 then $ (*!dustem_data.used).values *= f_HI ENDIF + + if isa(polsed) then BEGIN + + !dustem_data.polsed = ptr_new(polsed) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_data.polsed).values *= f_HI + + ENDIF + + if isa(polfrac) then BEGIN + + !dustem_data.polfrac = ptr_new(polfrac) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_data.polfrac).values *= f_HI + + ENDIF + + + if isa(psi_em) then BEGIN + + !dustem_data.psi_em = ptr_new(psi_em) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_data.psi_em).values *= f_HI + + ENDIF + + endif if keyword_set(st_ext_fit) then begin ;Two tests that determine what the user will eventually fit testexstks = isa(qext) and isa(uext) - testpexstks = (isa(qext) and isa(uext) and isa(polext)) - - IF (not testexstks and not isa(polext)) then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...' - - IF testpexstks THEN $ - message, 'Presence of both extinction polarization and stokes parameter values. Fitting structure is not correctly filtered. Aborting...' + IF ~testexstks then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...' + ;EXTINCTION - allowed for the default fitting of both because I do not know if the same formalism should be applied to the extinction data because it also implies deviding extinction data into two components. - if testexstks then nowfitting = 'QEXT and UEXT' else begin ;Serkowski is default - if isa(polext) then nowfitting = 'POLEXT' - endelse - - message, 'Extinction polarization component(s) to fit: '+nowfitting ,/continue + + message, 'Extinction polarization component(s) to fit: QEXT and UEXT',/continue + if isa(polext) then BEGIN !dustem_data.polext = ptr_new(polext) if keyword_set(f_HI) then if f_HI gt 0 then $ @@ -142,6 +163,7 @@ if !run_pol then begin if keyword_set(f_HI) then if f_HI gt 0 then $ (*!dustem_data.uext).values *= f_HI ENDIF + endif endif @@ -149,8 +171,8 @@ endif if not keyword_set(st_sed_show) then st_sed_show=st_sed_fit ;Setting of the structure that will be displayed -yy=dustem_check_data(data_sed=st_sed_show,data_ext=st_ext_show,sed=sed,ext=ext,polext=polext,polsed=polsed,polfrac=polfrac,qsed=qsed,used=used,qext=qext,uext=uext) - +yy=dustem_check_data(data_sed=st_sed_show,data_ext=st_ext_show,sed=sed,ext=ext,polext=polext,polsed=polsed,polfrac=polfrac,qsed=qsed,used=used,qext=qext,uext=uext,psi_em=psi_em,psi_ext=psi_ext) +;stop IF isa(sed) THEN !dustem_show.sed = ptr_new(sed) IF isa(ext) THEN !dustem_show.ext = ptr_new(ext) @@ -173,16 +195,16 @@ if !run_pol then begin teststks = isa(qsed) and isa(used) - IF (not teststks and not isa(polsed)) or (not teststks and not isa(polfrac)) then message, 'Stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...' + IF (~teststks and ~isa(polsed)) or (~teststks and ~isa(polfrac)) then message, 'Stokes parameters/Polarization showing data is not set correctly. Please check the input structure. Aborting...' - if teststks and isa(polsed) and not (isa(polfrac)) and not(isa(psi_em)) then nowshowing = 'QSED, USED and POLSED' - if teststks and isa(polsed) and not(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED, POLSED and PSI_EM' - if teststks and not(isa(polsed)) and isa(polfrac) and not(isa(psi_em)) then nowshowing = 'QSED, USED and POLFRAC' - if teststks and not(isa(polsed)) and isa(polfrac) and isa(psi_em) then nowshowing = 'QSED, USED, POLFRAC and PSI_EM' - if teststks and not(isa(polsed)) and not(isa(polfrac)) and not(isa(psi_em)) then nowshowing = 'QSED and USED' - if teststks and isa(polsed) and isa(polfrac) and not(isa(psi_em)) then nowshowing = 'QSED, USED, POLSED and POLFRAC' + if teststks and isa(polsed) and ~(isa(polfrac)) and ~(isa(psi_em)) then nowshowing = 'QSED, USED and POLSED' + if teststks and isa(polsed) and ~(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED, POLSED and PSI_EM' + if teststks and ~(isa(polsed)) and isa(polfrac) and ~(isa(psi_em)) then nowshowing = 'QSED, USED and POLFRAC' + if teststks and ~(isa(polsed)) and isa(polfrac) and isa(psi_em) then nowshowing = 'QSED, USED, POLFRAC and PSI_EM' + if teststks and ~(isa(polsed)) and ~(isa(polfrac)) and ~(isa(psi_em)) then nowshowing = 'QSED and USED' + if teststks and isa(polsed) and isa(polfrac) and ~(isa(psi_em)) then nowshowing = 'QSED, USED, POLSED and POLFRAC' if teststks and isa(polsed) and isa(polfrac) and isa(psi_em) then nowshowing = 'QSED, USED, POLFRAC, POLSED and PSI_EM' - if teststks and not(isa(polsed)) and not(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED and PSI_EM' + if teststks and ~(isa(polsed)) and ~(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED and PSI_EM' message, 'Polarization component(s) to show: '+nowshowing ,/continue @@ -217,12 +239,12 @@ if !run_pol then begin if isa(st_ext_show) then begin teststks = isa(qext) and isa(uext) - IF (not teststks and not isa(polext)) then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...' + IF (~teststks) then message, 'Extinction stokes parameters/Polarization data is not set correctly. Please check the input structure. Aborting...' - if teststks and isa(polext) and not(isa(psi_ext)) then nowshowing = 'QEXT, UEXT and POLEXT' + if teststks and isa(polext) and ~(isa(psi_ext)) then nowshowing = 'QEXT, UEXT and POLEXT' if teststks and isa(polext) and isa(psi_ext) then nowshowing = 'QEXT, UEXT, POLEXT and PSI_EXT' - if teststks and not(isa(polext)) and isa(psi_ext) then nowshowing = 'QEXT,UEXT and PSI_EXT' - if teststks and not(isa(polext)) and not(isa(psi_ext)) then nowshowing = 'QEXT and UEXT' + if teststks and ~(isa(polext)) and isa(psi_ext) then nowshowing = 'QEXT,UEXT and PSI_EXT' + if teststks and ~(isa(polext)) and ~(isa(psi_ext)) then nowshowing = 'QEXT and UEXT' message, 'Polarization component(s) to show: '+nowshowing ,/continue @@ -255,12 +277,10 @@ IF keyword_set(rchi2_weight) THEN BEGIN !fit_rchi2_weight.ext = rchi2_weight.ext IF !run_pol THEN BEGIN - ;IF isa(polsed) THEN !fit_rchi2_weight.polsed = rchi2_weight.polsed - ;IF isa(polfrac) THEN !fit_rchi2_weight.polfrac = rchi2_weight.polfrac + IF isa(qsed) THEN !fit_rchi2_weight.qsed = rchi2_weight.qsed IF isa(used) THEN !fit_rchi2_weight.used = rchi2_weight.used - IF isa(polext) THEN !fit_rchi2_weight.polext = rchi2_weight.polext IF isa(qext) THEN !fit_rchi2_weight.qext = rchi2_weight.qext IF isa(uext) THEN !fit_rchi2_weight.uext = rchi2_weight.uext ENDIF @@ -271,13 +291,13 @@ ENDIF ELSE BEGIN !fit_rchi2_weight.ext=1. If !run_pol THEN BEGIN - ;IF isa(polsed) THEN !fit_rchi2_weight.polsed=1. - ;IF isa(polfrac) THEN !fit_rchi2_weight.polfrac=1. + IF isa(qsed) THEN !fit_rchi2_weight.qsed=1. IF isa(used) THEN !fit_rchi2_weight.used=1. - IF isa(polext) THEN !fit_rchi2_weight.polext=1. + IF isa(qext) THEN !fit_rchi2_weight.qext=1. IF isa(uext) THEN !fit_rchi2_weight.uext=1. + ENDIF ENDELSE diff --git a/src/idl/dustemcgwin_dataset.pro b/src/idl/dustemcgwin_dataset.pro index 5bdc3cd..5f5e8c3 100644 --- a/src/idl/dustemcgwin_dataset.pro +++ b/src/idl/dustemcgwin_dataset.pro @@ -44,10 +44,14 @@ if !run_pol then begin specpol = st.polsed.em_tot * fact +; stop ;Instead of matching the non-zero values in spec and specpol (dustem outputs), ;I will only test on specpol as spec does not contain any null values FOR NOW. ;If this happens to change, these lines will have to be adjusted. + ;CHECKING IF THIS IS CORRECT FOR EACH OF THE SPECIES SPECTRA = ? + + specpolfrac = specpol *0. psi_em = specpol*0. indpolsed = where(specpol ne 0, ct_nullpolsed) ;I am not using the counter for now as it serves me no direct purpose. @@ -72,8 +76,12 @@ IF scopes[0] NE 'NONE' THEN BEGIN specq = (*(*!dustem_plugin).(i))[*,1] specu = (*(*!dustem_plugin).(i))[*,2] - ;specpol = sqrt(1.0E10*specq^2+1.0E10*specu^2) ;sqrt(((*(*!dustem_plugin).(i))[*,1])^2+(*(*!dustem_plugin).(i))[*,2]^2) + + ;STILL LACKING A WORKAROUND. specpol = (specq*1.0E10)^2+(specu*1.0E10)^2 & specpol=1.0E-10*sqrt(specpol) ;for some reason cgwindow couldn't handle these mathematical operations. It seems this is caused by values that are too close to zero + + ;specpol = (specq^2 + specu^2) & specpol = sqrt(specpol) + specpolfrac[indpolsed] = specpol[indpolsed]/spec[indpolsed] psi_em[indpolsed] = 0.5*atan(specu[indpolsed],specq[indpolsed])/degtorad ;stop @@ -293,7 +301,9 @@ if keyword_set(dataset) then begin end - 'POLSED': begin + 'POLSED': begin ;#NB: There seems to be an issue with the display of the dust species. 'Carbonaceous grains' 's emission exceeds the total emission. + ;MABE AN ISSUE WITH HOW 'SPECPOL' refreshes? + ;PROCEEDING TO POLFRAC idx_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',ct_filt) idx_spec=where((*!dustem_data.polsed).filt_names EQ 'SPECTRUM',ct_spec) @@ -426,7 +436,7 @@ if keyword_set(dataset) then begin ;Locating all the hidden data points (spectrum+filter) match2,((*!dustem_data.polsed).wav),((*!dustem_show.polsed).wav),fit_polsedpts,show_polsedpts ;only show_polsedpts is needed idx_rmv_polsed=where(show_polsedpts eq -1, ct_hdnpts) ; indices of the points to hide - + ;stop if ct_hdnpts ne 0 then begin ;Hidden data points are present ;stop ;Locating the hidden spectrum and filter data points @@ -442,7 +452,7 @@ if keyword_set(dataset) then begin cgerrplot,(((*!dustem_show.polsed).wav)[idx_rmv_polsed])(idx_spec_hdn),(((*!dustem_show.polsed).values)[idx_rmv_polsed])(idx_spec_hdn)-rms,(((*!dustem_show.polsed).values)[idx_rmv_polsed])(idx_spec_hdn)+rms,color='Black' endif - if ct_filt_hdn then begin + if ct_filt_hdn ne 0 then begin ;Plotting of hidden filter data points cgplot,(((*!dustem_show.polsed).wav)[idx_rmv_polsed])(idx_filt_hdn),(((*!dustem_show.polsed).values)[idx_rmv_polsed])(idx_filt_hdn),pos=position,/ylog,/xlog,/ys,xs=9,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=yr,psym=8,syms=0.8 @@ -466,7 +476,161 @@ if keyword_set(dataset) then begin end - 'POLFRAC': begin + 'POLFRAC': begin ; YOU NEED TO CODE THIS PART - TOU STILL HAVE'NT FINISHED. + idx_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',ct_filt) + idx_spec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',ct_spec) + + ;#1) get the plotting keywords (pertaining to each data set) @here when the _extra structure is ready + + if keyword_set(nodata) then begin ;when the data is not present + + xtit = '' + cgplot,wavs,wavs*0.,/nodata,/xlog,/ys,xs=1,pos=position,noerase=1,charsize=1.15,xtickformat='(A1)',color='Powder Blue',xr=xr,xtit='',yr=[1.00E-03,50.00],/ylog,psym=8,syms=0.8 + xyouts,pospltxt[0],pospltxt[1],textoidl('p_{\nu} (%)'),color=0,/normal,charsize=1.1 + + endif else begin ;when the data is present ; normal plot + + if keyword_set(refresh) then begin ;The data points in the plot are being refreshed + + ;##################TAKEN FROM SED + + ;Plotting of the spectra of the dust species + + ;I THINK I WILL HAVE AN ERROR HERE BECAUSE I WILL HAVE TO DIVIDE TWO ARRAYS FOR EACH DUST SPECIES. + ; (st.polsed.(i+1)/st.sed.(i+1))*fact + ;There is an issue with the fact that there are grain species that are not polarized + ;test and set a counter? YES FOR NOW + ; + + FOR i=0L,Ngrains-1 DO BEGIN + + testneq = where(st.polsed.(i+1) ne 0, countneq) + + vecfin = st.polsed.(i+1)*0.0 + ;SET the condition later + if countneq ne 0 then vecfin[testneq] = (st.polsed.(i+1))[testneq]/(st.sed.(i+1))[testneq] + + cgoplot,st.polsed.wav,vecfin*100,pos=position,noerase=1,color=use_cols[i],xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=[1.00E-03,50.00] + + ENDFOR + + ;Plotting of the plugins. + + ;NB: YOU HAVE YET TO TEST THIS BLOCK. DON'T FORGET TO ADD A POLARIZATION FRACTION TO FF for instance. + ; The thing is you've already had errors of that kind. + + ;PLUGINS SEEM TO HAVE ZERO NULL VALUES. THIS IS GOOD NEWS. LESS TESTING. + ;BUT if you get errors you know it has to be the division of the two arrays. + ;THE CHANGES SEEM LIKE THEY'RE GONNA WORK. + + + ;MISSING KEYWORDS? + for i=0L,n_plgns-1 do begin + IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_POLSED') THEN begin + + cgoplot,st.polsed.wav,sqrt(((*(*!dustem_plugin).(i))[*,1])^2+((*(*!dustem_plugin).(i))[*,2])^2)/((*(*!dustem_plugin).(i))[*,0])*100,color=clrs_plgns[i],pos=position,noerase=1,linestyle=2 + + ENDIF + endfor + + + ;PLotting of the interpolates corresponding to spectrum and filter points + + IF ct_spec NE 0 THEN BEGIN + + xx=((*!dustem_data.polfrac).wav)[idx_spec] + yy=prediction[idx_spec] + cgoplot,xx,yy*100,color='red',pos=position,psym=7,syms=2,noerase=1 + ENDIF + + + IF ct_filt NE 0 THEN BEGIN + xx=((*!dustem_data.polfrac).wav)[idx_filt] + yy=prediction[idx_filt] + cgoplot,xx,yy*100,color='red',pos=position,psym=6,syms=2,noerase=1 + ENDIF + + + ;Plotting of the total dust emission spectrum + cgoplot,st.polsed.wav,specpolfrac*100,pos=position,noerase=1,/xlog,/ys,/xs + + + endif else begin ;The data points in the plot that remain unchanged. + + + ;Plotting of frequency axis + ;cgaxis, xaxis=1, xlog=1, xs=1, xminor=10, xticklen=0.05, xrange=((!const.c*1E6/(_extra.(ind_xr)))*1E-9),charsize=1.15,title=textoidl('\nu (GHz)') + + + ;Spectrum points are treated before for plotting reasons. + if ct_spec ne 0 then begin + + ;Plotting of spectrum data points (to be fitted) + cgplot,((*!dustem_data.polfrac).wav)(idx_spec),((*!dustem_data.polfrac).values)(idx_spec)*100,/xlog,/ylog,/ys,xs=1,pos=position,noerase=1,xtit=' ',charsize=1.15,xtickformat='(A1)',color='Powder Blue',xr=xr,yr=[1.00E-03,50.00],psym=16,syms=0.8;,ytickformat='dstmwrp_exp' + + ;Plotting of the spectrum error points + rms=3.*((*!dustem_data.polfrac).sigma)(idx_spec)/2. + cgerrplot,((*!dustem_data.polfrac).wav)(idx_spec),(((*!dustem_data.polfrac).values)(idx_spec)-rms)*100,(((*!dustem_data.polfrac).values)(idx_spec)+rms)*100,color='Powder Blue' + + endif + + if ct_filt ne 0 then begin + + ;Plotting of filter data points (to be fitted) + testing for prior plotting + if ct_spec ne 0 then cgoplot,((*!dustem_data.polfrac).wav)(idx_filt),((*!dustem_data.polfrac).values)(idx_filt)*100,charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,pos=position,/ys,xs=1,noerase=1,xtickformat='(A1)',xr=xr,yr=[1.00E-03,50.00],/xlog,/ylog else $;,ytickformat='dstmwrp_exp' + cgplot,((*!dustem_data.polfrac).wav)(idx_filt),((*!dustem_data.polfrac).values)(idx_filt)*100,charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,pos=position,/ys,xs=1,noerase=1,xtickformat='(A1)',xtit=' ',xr=xr,yr=[1.00E-03,50.00],/xlog,/ylog;,ytickformat='dstmwrp_exp' + + + ;Plotting of the filter error points + rms=3.*((*!dustem_data.polfrac).sigma)(idx_filt)/2.;/dustem_polfrac(idx_filt) + cgerrplot,((*!dustem_data.polfrac).wav)(idx_filt),(((*!dustem_data.polfrac).values)(idx_filt)-rms)*100,(((*!dustem_data.polfrac).values)(idx_filt)+rms)*100,color='Dodger Blue' + + endif + + + + xyouts,pospltxt[0],pospltxt[1],textoidl('p_{\nu} (%)'),color=0,/normal,charsize=1.1 + + ;Locating all the hidden data points (spectrum+filter) + + ;The filtering has been done prior (in the primary routine) + match2,((*!dustem_data.polfrac).wav),((*!dustem_show.polfrac).wav),fit_polfracpts,show_polfracpts ;only show_polfracpts is needed + idx_rmv_polfrac=where(show_polfracpts eq -1, ct_hdnpts) ; indices of the points to hide + ;stop + if ct_hdnpts ne 0 then begin ;Hidden data points are present + ;stop + + ;Locating the hidden spectrum and filter data points + idx_filt_hdn = where(((*!dustem_show.polfrac).filt_names)(idx_rmv_polfrac) NE 'SPECTRUM',ct_filt_hdn) + idx_spec_hdn = where(((*!dustem_show.polfrac).filt_names)(idx_rmv_polfrac) EQ 'SPECTRUM',ct_spec_hdn) + + if ct_spec_hdn ne 0 then begin + ;Plotting of hidden spectrum data points + cgplot,(((*!dustem_show.polfrac).wav)[idx_rmv_polfrac])(idx_spec_hdn),(((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_spec_hdn)*100,pos=position,/xlog,/ylog,/ys,xs=1,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=[1.00E-03,50.00],psym=8,syms=0.8 + + ;Plotting of hidden spectrum error points + rms=3.*(((*!dustem_show.polfrac).sigma)[idx_rmv_polfrac])(idx_spec_hdn)/2. + cgerrplot,(((*!dustem_show.polfrac).wav)[idx_rmv_polfrac])(idx_spec_hdn),((((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_spec_hdn)-rms)*100,((((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_spec_hdn)+rms)*100,color='Black' + endif + + if ct_filt_hdn ne 0 then begin + ;stop + plotsym,0, /fill + ;Plotting of hidden filter data points + cgplot,(((*!dustem_show.polfrac).wav)[idx_rmv_polfrac])(idx_filt_hdn),(((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_filt_hdn)*100,pos=position,/xlog,/ylog,/ys,xs=1,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=[1.00E-03,50.00],psym=8,syms=0.8 + ;stop + ;Plotting of hidden filter error point + rms=3.*(((*!dustem_show.polfrac).sigma)[idx_rmv_polfrac])(idx_filt_hdn)/2. + cgerrplot,(((*!dustem_show.polfrac).wav)[idx_rmv_polfrac])(idx_filt_hdn),((((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_filt_hdn)-rms)*100,((((*!dustem_show.polfrac).values)[idx_rmv_polfrac])(idx_filt_hdn)+rms)*100,color='Black' + + endif + + endif + + + endelse + + endelse end @@ -694,8 +858,8 @@ if keyword_set(dataset) then begin endif else begin ;The data points in the plot that remain unchanged ;stop xtit=textoidl('\lambda (\mum)') - if !run_pol then xtit = '' - cgplot,wavs,wavs/wavs,/xlog,/ys,xs=1,pos=position,noerase=1,xtickformat='(A1)',color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0 + ;if !run_pol then xtit = '' + cgplot,wavs,wavs/wavs,/xlog,/ys,xs=1,pos=position,noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0 xyouts,pospltxt[0],pospltxt[1],textoidl('norm'),color=0,/normal,charsize=1.1 endelse @@ -1053,8 +1217,8 @@ if keyword_set(dataset) then begin endif else begin ;The data points in the plot that remain unchanged ;stop xtit=textoidl('\lambda (\mum)') - if !run_pol then xtit = '' - cgplot,wavs,wavs/wavs,/xlog,/ys,xs=1,pos=position,noerase=1,xtickformat='(A1)',color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0 + ;if !run_pol then xtit = '' + cgplot,wavs,wavs/wavs,/xlog,/ys,xs=1,pos=position,noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0 xyouts,pospltxt[0],pospltxt[1],textoidl('norm'),color=0,/normal,charsize=1.1 endelse diff --git a/src/idl_misc/JPBLib_for_Dustemwrap/General/modify_extra.pro b/src/idl_misc/JPBLib_for_Dustemwrap/General/modify_extra.pro index 6965176..29dc45d 100644 --- a/src/idl_misc/JPBLib_for_Dustemwrap/General/modify_extra.pro +++ b/src/idl_misc/JPBLib_for_Dustemwrap/General/modify_extra.pro @@ -75,6 +75,7 @@ IF keyword_set(replace) THEN BEGIN _extra_out=_extra tags=tag_names(_extra) ind=where(tags EQ tag_name,count) + ;stop IF count NE 0 THEN BEGIN _extra_out.(ind[0])=tag_value ENDIF ELSE BEGIN -- libgit2 0.21.2