diff --git a/src/idl/dustem_activate_plugins.pro b/src/idl/dustem_activate_plugins.pro index dac77a6..c6d387b 100755 --- a/src/idl/dustem_activate_plugins.pro +++ b/src/idl/dustem_activate_plugins.pro @@ -76,17 +76,19 @@ endelse ;IF /avoid is used the specified scopes are avoided. They're considered when /avoid is missing. ;using '**' so not to write the entire scope dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['*STELLAR*','*REPLACE*','*ISRF*'], avoid=1 -;stop + ;=== The line below makes fawlty go Segmentation fault. st is undefined at that stage ;=== it would run fine if avoid was set to 1 +;IC: I think the segmentation fault comes from the fact that the user wants to use models that have the /pol keyword in the GRAIN.DAT file without being in polarization mode. dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['*STELLAR*','*ISRF*'] ,dustem_run=1,st=st ;dustem output is available at this step. ;added this small condition to ensure that dustem is run ;(if no ISRF plugin is used) -if ~isa(st) then begin - st=dustem_run(p_dim) - !dustem_current=ptr_new(st) -endif +; if ~isa(st) then begin +; st=dustem_run(p_dim) +; !dustem_current=ptr_new(st) +; endif +;THE ABOVE CONDITION IS UNNECESSARY SINCE DUSTEM IS CALLED ABOVE WHETHER PLUGINS ARE RUN OR NOT. dustem_run_plugins, p_dim ,param_descs, param_values, param_func, ['*REPLACE*'] diff --git a/src/idl/dustem_check_data.pro b/src/idl/dustem_check_data.pro index 496fff9..5003d6d 100644 --- a/src/idl/dustem_check_data.pro +++ b/src/idl/dustem_check_data.pro @@ -1,39 +1,52 @@ -PRO dustem_check_data,data_sed,data_ext,sed,ext,polext,polsed,polfrac,psi_em,qsed,used,qext,uext,psi_ext,fpolext +PRO dustem_check_data, data_sed, $ + data_ext, $ + 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, $ + fpolext=fpolext ;+ ; NAME: ; dustem_check_data ; ; PURPOSE: -; Verify that the input data conforms to what DustEMWrap requires +; Clean and modify the input DustEMWrap structure to meet the requirements of DustEMWrap. +; This is important for the subsequent internal handling/fitting of data. ; ; CATEGORY: ; DustEMWrap, Mid-Level, Distributed, Initialization ; ; CALLING SEQUENCE: -; dustem_check_data,data_sed,data_ext,sed,ext,polext,polsed,polfrac,psi_em,qsed,used,qext,uext,psi_ext,fpolext,[/help] +; dustem_check_data,data_sed,data_ext,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,fpolext=fpolext,[/help] ; ; INPUTS: -; data_sed : Emission dustemwrap structure -; data_ext : Extinction dustemwrap structure +; data_sed : Emission DustEMWrap structure +; data_ext : Extinction DustEMWrap structure ; ; OPTIONAL INPUT PARAMETERS: ; None ; -; OUTPUTS: -; ## adapting the data structures to the format used by the fitting algorithm -; sed : emission formatted structure -; ext : extinction formatted structure -; polext : polarized extinction formatted structure -; polsed : polarized emission formatted structure -; polfrac : fractional polarized emission formatted structure -; psi_em : emission polarization angle formatted structure -; qsed : stokesQ emission formatted structure -; used : stokesU emission formatted structure -; qext : stokesQ extinction formatted structure -; uext : stokesU extinction formatted structure -; psi_ext : extinction polarization angle formatted structure -; fpolext : fractional extinction polarization formatted structure +; OUTPUTS: +; sed : modified emission structure +; ext : modified extinction structure +; polext : modified polarized extinction structure +; polsed : modified polarized emission structure +; polfrac : modified fractional polarized emission structure +; psi_em : modified emission polarization angle structure +; qsed : modified StokesQ emission structure +; used : modified StokesU emission structure +; qext : modified StokesQ extinction structure +; uext : modified StokesU extinction structure +; psi_ext : modified extinction polarization angle structure +; fpolext : modified fractional extinction polarization structure ; ; OPTIONAL OUTPUT PARAMETERS: ; None @@ -53,12 +66,15 @@ PRO dustem_check_data,data_sed,data_ext,sed,ext,polext,polsed,polfrac,psi_em,qse ; None ; ; EXAMPLES -; dustem_check_data,data_sed,data_ext,sed,ext,polext,polsed,polfrac,psi_em,qsed,used,qext,uext,psi_ext,fpolext,[/help] -; +; dustem_check_data,data_sed,data_ext,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,fpolext=fpolext,[/help] +; dustem_check_data,/help ; MODIFICATION HISTORY: -; Written I Choubani 2022 +; Written IC 2022 ; Evolution details on the DustEMWrap gitlab. ; See http://dustemwrap.irap.omp.eu/ for FAQ and help. +; COMMENTS: +; Might still need a bit of cleaning. I am unsure of the availability of covariance data in polarized extinction. +; Still regarding extinction, some of the error messages haven't been mirrored correctly with respect to the emission block. ;- @@ -202,8 +218,7 @@ IF KEYWORD_SET(data_sed) THEN BEGIN ;Generation of POLFRAC,POLSED and PSI data by default ;Testing on the presence of covariances otherwise an error is returned - ;The testing happens here because covariances are needed for the error estimation of P,p and PSI. - + ;The testing happens here because covariances are needed for the error estimation of P,p and PSI. (Also cannot test before because we're not in !run_pol mode) tstcovu = where((data_sed(testu).sigmaIU) EQ la_undef() or (data_sed(testu).sigmaQU) EQ la_undef(), ctcovu) if ctcovu NE 0 then BEGIN ;message, 'Covariance data is missing in StokesU data. Data, errors and covariances do not match. Aborting...',/continue @@ -253,11 +268,9 @@ IF KEYWORD_SET(data_sed) THEN BEGIN ENDIF -;EXTINCTION HASN'T BEEN TREATED YET - - -;SMALLL ISSUE .. POLEXTFRAC is not present. Let's create it for completness - +;EXTINCTION +; +;Covariances?... IF KEYWORD_SET(data_ext) THEN begin ans='' diff --git a/src/idl/dustem_fit_ext_example.pro b/src/idl/dustem_fit_ext_example.pro index 1636998..afc8268 100644 --- a/src/idl/dustem_fit_ext_example.pro +++ b/src/idl/dustem_fit_ext_example.pro @@ -170,7 +170,7 @@ ext=read_xcat(file,/silent) ;=== HERE WE SET THE DATA SO THAT THE COMPUTE_ FUNCTIONS BELOW HAVE ;=== THE NECESSARY FIELDS AND INFORMATION -dustem_set_data,sed,sed,ext,ext +dustem_set_data,m_fit=m_fit,m_show=m_show,x_fit=ext,x_show=ext ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT ;== AND ACTIVATE ANY PLUGINS @@ -217,9 +217,9 @@ IF keyword_set(fits_save_and_restore) THEN BEGIN ;==== plot result taken from the saved fits table res=*(*!dustem_fit).CURRENT_PARAM_VALUES IF !dustem_noobj THEN BEGIN - dustemwrap_plot_noobj,res,dustem_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)' + dustemwrap_plot_noobj,res,st=dustem_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)' ENDIF ELSE BEGIN - dustemwrap_plot,res,dustem_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)' + dustemwrap_plot,res,st=dustem_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)' ENDELSE IF keyword_set(wait) THEN BEGIN message,'Saved the results as FITS in the file: '+fits_save_and_restore+', and made a plot using the data in this file',/info diff --git a/src/idl/dustem_fit_ext_pol_example.pro b/src/idl/dustem_fit_ext_pol_example.pro index a2fa7c8..e74b387 100755 --- a/src/idl/dustem_fit_ext_pol_example.pro +++ b/src/idl/dustem_fit_ext_pol_example.pro @@ -113,7 +113,7 @@ use_polarization=0. ;default is no polarization in models use_window=2 ; default graphics window number to use for plotting the results use_verbose=0 if keyword_set(verbose) then use_verbose=1 -use_Nitermax=5 ; maximum number of iterations for the fit +use_Nitermax=1 ; maximum number of iterations for the fit IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax dustem_define_la_common @@ -191,7 +191,7 @@ ENDIF ;=== THE NECESSARY FIELDS AND INFORMATION ;== note that m_fit and m_show do not exist -- this is deliberate ;== since we are not treating emission data -dustem_set_data, m_fit,m_show,ext,ext +dustem_set_data, m_fit=m_fit,m_show=m_show,x_fit=ext,x_show=ext ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT ;== AND ACTIVATE ANY PLUGINS (HENCE WE NEED TO DO THIS BEFORE THE COMPUTE_ CALLS) @@ -206,7 +206,7 @@ dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,u if not keyword_set(ext_file) then begin ;== GENERATE THE DATA USING THE MODEL AND TRUE VALUES ext.EXT_I = dustem_compute_ext(true_vals,st=st) - toto = dustem_compute_stokext(true_vals,st=st,qext_spec=dustem_qext,uext_spec=dustem_uext) + toto = dustem_compute_stokext(true_vals,st=st) ext.EXT_Q = toto[0] ext.EXT_U = toto[1] @@ -229,7 +229,7 @@ tol=1.e-10 ;=== INFORMATION TO MAKE THE PLOTS ;=== _x/_extinction means extinction plots, _m/_emissions means emission plots xr_x=[0.01,30] -yr_fpolext=[1.00E-10,30] +yr_fpolext=[1.00E-10,50] ;yr_fpolext=[1e-4,30] tit='Dust Optical Depth' @@ -247,7 +247,7 @@ if keyword_set(wait) then begin message,'Time taken [sec]: '+sigfig(t2-t1,2,/sci),/info wait,wait end - +stop IF keyword_set(fits_save_and_restore) THEN BEGIN message,'Writing out results structure: '+fits_save_and_restore,/info dustem_write_fits_table,filename=fits_save_and_restore,help=help @@ -257,6 +257,7 @@ IF keyword_set(fits_save_and_restore) THEN BEGIN dustem_read_fits_table,filename=fits_save_and_restore,dustem_st=dustem_st ;==== plot result taken from the saved fits table res=*(*!dustem_fit).CURRENT_PARAM_VALUES + !dustem_end=-1; subsequent calls to the fit. 0 is intialization an 1 is for the last iteration. It's initialized to 0 after the creation of the fits file. IF !dustem_noobj THEN BEGIN dustemwrap_plot_noobj,res,st=dustem_st,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit+' (From Saved FITS file)' ENDIF ELSE BEGIN diff --git a/src/idl/dustem_fit_intensity_example.pro b/src/idl/dustem_fit_intensity_example.pro index b67161e..279d36a 100644 --- a/src/idl/dustem_fit_intensity_example.pro +++ b/src/idl/dustem_fit_intensity_example.pro @@ -296,7 +296,7 @@ IF count NE 0 THEN sed[ind].sigmaII=(0.2*sed[ind].StokesI)^2 ;== SET THE OBSERVATIONAL STRUCTURE ;== sed is passed twice in the call -- the first occurrence is the SED that you ;== wish to fit, the second occurrence is the SED that you wish to visualise. -dustem_set_data,sed,sed +dustem_set_data,m_fit=sed,m_show=sed ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE ;== ADJUSTED DURING THE FIT diff --git a/src/idl/dustem_fit_intensity_mbb_example.pro b/src/idl/dustem_fit_intensity_mbb_example.pro index 443ca44..604d350 100644 --- a/src/idl/dustem_fit_intensity_mbb_example.pro +++ b/src/idl/dustem_fit_intensity_mbb_example.pro @@ -194,7 +194,7 @@ IF count NE 0 THEN sed[ind].sigmaII=(0.2*sed[ind].StokesI)^2 ;== SET THE OBSERVATIONAL STRUCTURE ;== sed is passed twice -- the first occurrence is the SED that you ;== wish to fit, the second occurrence is the SED that you wish to visualise. -dustem_set_data,sed,sed +dustem_set_data,m_fit=sed,m_show=sed ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv $ diff --git a/src/idl/dustem_fit_polarization_example.pro b/src/idl/dustem_fit_polarization_example.pro index dbdca57..6545a63 100644 --- a/src/idl/dustem_fit_polarization_example.pro +++ b/src/idl/dustem_fit_polarization_example.pro @@ -164,7 +164,7 @@ sed=read_xcat(file,/silent) ;== SET THE OBSERVATIONAL STRUCTURE ;== sed is passed twice -- the first occurrence is the SED that you ;== wish to fit, the second occurrence is the SED that you wish to visualise. -dustem_set_data,sed,sed +dustem_set_data,m_fit=sed,m_show=sed diff --git a/src/idl/dustem_fit_sed_ext_pol_example.pro b/src/idl/dustem_fit_sed_ext_pol_example.pro index 3326ddf..083d2d0 100755 --- a/src/idl/dustem_fit_sed_ext_pol_example.pro +++ b/src/idl/dustem_fit_sed_ext_pol_example.pro @@ -217,7 +217,7 @@ ENDIF ;=== HERE WE SET THE DATA SO THAT THE COMPUTE_ FUNCTIONS BELOW HAVE ;=== THE NECESSARY FIELDS AND INFORMATION -dustem_set_data,sed,sed,ext,ext +dustem_set_data,m_fit=sed,m_show=sed,x_fit=ext,x_show=ext ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS FOR THE FIT ;== AND ACTIVATE ANY PLUGINS (HENCE WE NEED TO DO THIS BEFORE THE COMPUTE_ CALLS) diff --git a/src/idl/dustem_mask_data.pro b/src/idl/dustem_mask_data.pro index b5aa4f8..3f940ad 100755 --- a/src/idl/dustem_mask_data.pro +++ b/src/idl/dustem_mask_data.pro @@ -1,157 +1,191 @@ - FUNCTION dustem_mask_data,data_struct,omit=omit,data_set=data_set,help=help +FUNCTION dustem_mask_data, data_struct, $ + filters, $ + data_set=data_set, $ + default=default, $ + help=help + ;+ ; NAME: ; dustem_mask_data ; PURPOSE: -; Removing filters from a data structure. +; Set dataset values to la_undef() for specified filter/spectrum data in an input DustEmWrap data structure. ; CATEGORY: ; DUSTEM Wrapper ; CALLING SEQUENCE: -; data_struct=dustem_mask_data(data_struct,omit=,data_set=) +; data_struct=dustem_mask_data(data_struct,filters,data_set=,/default) ; INPUTS: -; data_struct : SED Data structure formatted as outlined in the DustermWrap Users' Guide -; data_set : Array of character strings containing the different data sets to be filtered. -; ie: ['SED'],['QSED'],'USED'] or just ['QSED'] for example. The use of ['STOKESI','STOKESQ','STOKESU'] instead is equivalent. -; omit : Can be an array or list. In both cases it should contain the filter names as character strings. -; ie: omit=['SPIRE1','SPIRE2'] for an array or omit=list(['SPIRE2','SPIRE1','PACS2'],['HFI5','HFI4']) for a list. -; When using a list, its shape (dimension) should be the number of elements of the data_set array. -; ie: Each dimension is an array of filters to be removed from the corresponding data set in the data_set array. +; data_struct : SED Data structure following the format explained in the DustermWrap Users' Guide +; +; filters : Can be an array or list. In both cases it should contain the filter names as character strings. +; ie: filters=['SPIRE1','SPIRE2'] for an array or filters=list(['SPIRE2','SPIRE1','PACS2'],['HFI5','HFI4']) for a list. +; When using a list, its shape (dimension) and the number of elements of the data_set array should be equal. +; Each dimension is an array of filters to be removed from the corresponding data set in the data_set array. +; ie: filters=list(['SPIRE2','SPIRE1','PACS2'],['LFI2','LFI3']) for datas_set=['STOKESI','STOKESQ'] ; OPTIONAL INPUT PARAMETERS: -; None +; data_set : Array of character strings containing the different data sets to be filtered. +; ie: ['STOKESI','STOKESQ','STOKESU'] or just ['STOKESQ'] for example. ; OUTPUTS: ; data_struct ; OPTIONAL OUTPUT PARAMETERS: ; None ; ACCEPTED KEY-WORDS: -; help = if set, prints this help +; help = if set, print this help +; +; default = if set, remove targeted filter/spectrum point(s) from datasets that depend on the filtered dataset(s). +; EAMPLE +; masked_structure = dustem_mask_data(data_struct,['IRAS1','IRAS2'],data_set=['STOKESI']) ; COMMENTS -; Both list and array entries for omit and hide keywords are allowed. +; IC: +; If a wavelength value is supplied instead of a filter name, the procedure masks for a spectrum point with that corresponding wavelength +; If the dataset array isn't supplied, all the datasets in the input data structure will be considered. +; If both filters and dataset are arrays, all mentioned filters will be masked from all mentioned datasets. +; This procedure does not change the shape of the output structure partly because dustem_check_data.pro does. ; MODIFICATION HISTORY: -; Written by IC +; Written by IC 2022 ;- IF keyword_set(help) THEN BEGIN - doc_library,'dustem_filter_data' + doc_library,'dustem_mask_data' goto,the_end ENDIF -IF keyword_set(data_struct) and keyword_set(omit) THEN BEGIN - - ;Array of filters present in dustem - filters_dustem=[((*!dustem_filters).(0).filter_names)] +IF arg_present(data_struct) and arg_present(filters) THEN BEGIN + ;Array of filters present in dustem + filters_dustem=[((*!dustem_filters).(0).filter_names)] for i=1L,n_elements(tag_names(*!dustem_filters))-1 do begin filters_dustem=[filters_dustem,((*!dustem_filters).(i).filter_names)] endfor - ;Array of filters present in data_struct filters_data = data_struct.filter - ;Tag names present in data_struct (to be able to test on the different data sets) - tagnames = tag_names(data_struct) - - - if typename(omit) EQ 'LIST' then ind0 = omit[*] else ind0=intarr(n_elements(omit)) ; this is a better condition because it is also valid for when the user wants to omit 'spectrum' (non-filter) data points. + tagnames = tag_names(data_struct) + if typename(filters) EQ 'LIST' then ind0 = filters[*] else ind0=intarr(n_elements(filters)) ;because types can be float or double but not long. - - for i=0L,n_elements(omit)-1 do begin - dim_i = n_elements(omit(i)) ;getting the dimension of each element of the list/array + for i=0L,n_elements(filters)-1 do begin + dim_i = n_elements(filters(i)) ;getting the dimension of each element of the list/array ind0[i] = intarr(dim_i) + la_undef() ;initializing to la_undef() - endfor - + endfor ;Gathering the indices of the filters to remove. - - for i=0L,n_elements(omit)-1 do begin - - dim_i=n_elements(omit(i)) ;Getting the dimension of each element of the list - - arr_i=intarr(dim_i) + la_undef() ;initializing to la_undef() - + for i=0L,n_elements(filters)-1 do begin + dim_i=n_elements(filters(i)) ;Getting the dimension of each element of the list + arr_i=intarr(dim_i) + la_undef() ;initializing to la_undef() for j=0L,dim_i-1 do begin - ;testing whether the element corresponds to a specific wavelength. - ;accepts both numbers/strings - - if valid_num((omit[i])[j]) then begin - - ind=where(float((omit[i])[j]) EQ data_struct.wave,test0) - if test0 ne 0 then arr_i[j] = ind else message, 'WAVELENGTH '+strtrim((omit[i])[j],2)+' does not correspond to any spectrum point in the supplied data strucure.' - endif else begin - - ind = where((omit[i])[j] EQ filters_dustem,test1) ;Test if the filter is in the wrapper's database - if test1 eq 0 then message, 'Filter '+strtrim((omit[i])[j],2)+' cannot be found in the DustemWrap database.' - - ind = where((omit[i])[j] EQ filters_data,test2) ;Test if the filter in in the data itself - if test2 ne 0 then arr_i[j] = ind else message, 'Filter '+strtrim((omit[i])[j],2)+' cannot be found in the data supplied.' - + ;accepts both numbers/strings + if valid_num((filters[i])[j]) then begin + ind=where(float((filters[i])[j]) EQ data_struct.wave,test0) + if test0 ne 0 and data_struct[ind].filter EQ 'SPECTRUM' then arr_i[j] = ind else message, 'WAVELENGTH '+strtrim((filters[i])[j],2)+' does not correspond to any spectrum point in the supplied data strucure.' + endif else begin + ind = where((filters[i])[j] EQ filters_dustem,test1) + if test1 eq 0 then message, "Filter "+strtrim((filters[i])[j],2)+" isn't in the DustEMWrap Filters' list" + ind = where((filters[i])[j] EQ filters_data,test2) ;Test if the filter in in the data itself + if test2 ne 0 then arr_i[j] = ind else message, 'Filter '+strtrim((filters[i])[j],2)+' cannot be found in the supplied data structure.' endelse - - endfor - + endfor ind0[i] = arr_i ;Indices of the filters to be removed in the data structure - endfor + If typename(ind0) EQ 'LIST' then indices = ind0[i] ELSE indices = ind0 + + ;If default isn't set, errors in dustem_check_data.pro can occur because of differences in tag values. A common behavior of dustem_check_data.pro + ;example: filter removed for StokesI but not for SIGMAII + IF keyword_set(default) THEN BEGIN - eq_tags = ['STOKESI','EXT_I','STOKESQ','STOKESU','LARGEP','SMALLP','PSI','EXT_Q','EXT_U','EXT_P','EXT_smallp','PSI'] - - errvec = ['SIGMAII','SIGEXTI','SIGMAQQ','SIGMAUU','SIGMA_LARGEP','SIGMA_SMALLP','SIGMA_PSI','SIGEXTQ','SIGEXTU','SIGEXTP','SIGEXTsmallp', 'SIGEXT_PSI'] - + ;Part pertainig to the /default keyword + ;bunch of tests to locate datasets depending on STOKESI, STOKESQ, STOKESU, LARGEP, SMALLP and PSI + ;is it an emission or extinction structure? + ind_m = where(tagnames EQ 'STOKESI', ctm) + ind_x = where(tagnames EQ 'EXT_I', ctx) + IF ctm NE 0 THEN BEGIN ;no counters are used in the following because it's a consequence of the test being valid. + ;Getting the indices this way because it allows for an easy modification if/when tagnames in the xcat emission/extinction file change + ind_set_StokesI = fix([where(strmatch(tagnames,'STOKESI')), $ + where(strmatch(tagnames,'SIGMAII'))]) + ind_set_StokesQ = fix([where(strmatch(tagnames,'STOKESQ')), $ + where(strmatch(tagnames,'SIGMAQQ')), $ + where(strmatch(tagnames,'SIGMAIQ')) ]) + ind_set_StokesU = fix([where(strmatch(tagnames,'STOKESU')), $ + where(strmatch(tagnames,'SIGMAUU')), $ + where(strmatch(tagnames,'SIGMAIU')) ]) + ind_set_LargeP = fix([where(strmatch(tagnames,'LARGEP')), $ + where(strmatch(tagnames,'SIGMA_LARGEP'))]) + ind_set_smallp = fix([where(strmatch(tagnames,'SMALLP')), $ + where(strmatch(tagnames,'SIGMA_SMALLP'))]) + ind_set_psi = fix([where(strmatch(tagnames,'PSI')), $ + where(strmatch(tagnames,'SIGMA_PSI'))]) + ENDIF + IF ctx NE 0 THEN BEGIN + ind_set_EXT_I = fix([where(strmatch(tagnames,'EXT_I')), $ + where(strmatch(tagnames,'SIGEXTII')) ]) + ind_set_EXT_Q = fix([where(strmatch(tagnames,'EXT_Q')), $ + where(strmatch(tagnames,'SIGEXTQQ')), $ + where(strmatch(tagnames,'SIGEXTIQ')) ]) + ind_set_EXT_U = fix([where(strmatch(tagnames,'EXT_U')), $ + where(strmatch(tagnames,'SIGEXTUU')), $ + where(strmatch(tagnames,'SIGEXTIU')) ]) + ind_set_EXT_P = fix([where(strmatch(tagnames,'EXT_P')), $ + where(strmatch(tagnames,'SIGEXTP'))]) + ind_set_EXT_smallp = fix([where(strmatch(tagnames,'EXT_SMALLP')), $ + where(strmatch(tagnames,'SIGEXTSMALLP'))]) + ind_set_psi = fix([where(strmatch(tagnames,'PSI')), $ + where(strmatch(tagnames,'SIGEXTPSI'))]) + ENDIF + ENDIF + + ;setting filters' data to la_undef() IF keyword_set(data_set) then begin - - for i=0L,n_elements(data_set)-1 do begin - ind_set = where(tagnames EQ data_set(i),ctset) - ;here using eq_tags and errvec for another purpose - ;because I need to locate the associated errors + FOR i=0L,n_elements(data_set)-1 DO BEGIN + ind_set = where(tagnames EQ data_set(i),ctset) if ctset ne 0 then begin - - inderr = where(eq_tags EQ data_set(i)) ;no need for counters this is supposedly true at all times plus the error message is below in the next test - ind_errset = where(tagnames EQ errvec[inderr[0]]) - - ;in the case of PSI for emission and for absorption - ;I need to know the different cases: - - ;EMISSION - - if data_set(i) EQ 'PSI' and tagnames[3] EQ 'STOKESI' then ind_errset = where(tagnames EQ errvec[6]) - - - ;EXTINCION - - if data_set(i) EQ 'PSI' and tagnames[3] EQ 'EXT_I' then ind_errset = where(tagnames EQ errvec[11]) - - ;stop - endif else begin - ind_set1 = where(tag_names(*!dustem_data) EQ data_set(i),ctset) - if ctset eq 0 then message, 'Aborting. The supplied data set is not found in the DustemWrap database.' - ind_set = where(tagnames EQ eq_tags(ind_set1[0])) - ind_errset = where(tagnames EQ errvec(ind_set1[0])) - ;stop - endelse - if typename(ind0) EQ 'LIST' then begin - data_struct(ind0[i]).(ind_set) = la_undef() - data_struct(ind0[i]).(ind_errset) = la_undef() - endif else begin - data_struct(ind0).(ind_set) = la_undef() - data_struct(ind0).(ind_errset) = la_undef() - endelse - - endfor - endif else begin + ;the requires the default keyword to be set + IF keyword_set(default) THEN BEGIN + ;EMISSION CASE + IF ctm NE 0 THEN BEGIN + ;testing data_set(i) against the above strings. + test_StokesI = where(ind_set_StokesI EQ ind_set, ct_StokesI) + IF ct_StokesI NE 0 THEN data_struct(indices).(ind_set_StokesI) = la_undef() + test_StokesQ = where(ind_set_StokesQ EQ ind_set, ct_StokesQ) + IF ct_StokesQ NE 0 THEN data_struct(indices).(ind_set_StokesQ) = la_undef() + test_StokesU = where(ind_set_StokesU EQ ind_set, ct_StokesU) + IF ct_StokesU NE 0 THEN data_struct(indices).(ind_set_StokesU) = la_undef() + test_LargeP = where(ind_set_LargeP EQ ind_set, ct_LargeP) + IF ct_LargeP NE 0 THEN data_struct(indices).(ind_set_LargeP) = la_undef() + test_smallp = where(ind_set_smallp EQ ind_set, ct_smallp) + IF ct_smallp NE 0 THEN data_struct(indices).(ind_set_smallp) = la_undef() + test_psi = where(ind_set_psi EQ ind_set, ct_psi) + IF ct_psi NE 0 THEN data_struct(indices).(ind_set_psi) = la_undef() + ENDIF + ;EXTINCTION CASE + IF ctx NE 0 THEN BEGIN + ;testing data_set(i) against the above strings. + test_EXT_I = where(ind_set_EXT_I EQ ind_set, ct_EXT_I) + IF ct_EXT_I NE 0 THEN data_struct(indices).(ind_set_EXT_I) = la_undef() + test_EXT_Q = where(ind_set_EXT_Q EQ ind_set, ct_EXT_Q) + IF ct_EXT_Q NE 0 THEN data_struct(indices).(ind_set_EXT_Q) = la_undef() + test_EXT_U = where(ind_set_EXT_U EQ ind_set, ct_EXT_U) + IF ct_EXT_U NE 0 THEN data_struct(indices).(ind_set_EXT_U) = la_undef() + test_EXT_P = where(ind_set_EXT_P EQ ind_set, ct_EXT_P) + IF ct_EXT_P NE 0 THEN data_struct(indices).(ind_set_EXT_P) = la_undef() + test_EXT_smallp = where(ind_set_EXT_smallp EQ ind_set, ct_EXT_smallp) + IF ct_EXT_smallp NE 0 THEN data_struct(indices).(ind_set_EXT_smallp) = la_undef() + test_psi = where(ind_set_psi EQ ind_set, ct_psi) + IF ct_psi NE 0 THEN data_struct(indices).(ind_set_psi) = la_undef() + ENDIF + ENDIF ELSE data_struct(indices).(ind_set) = la_undef() + ENDIF + ENDFOR + ENDIF ELSE BEGIN ;If no datasets are supplied, all tag values are set to la_undef() ind_tags = where(tagnames ne 'INSTRU' and tagnames ne 'FILTER' and tagnames ne 'WAVE') - for j=0L,n_elements(ind_tags)-1 do begin - data_struct(ind0).(ind_tags(j)) = la_undef() - endfor + FOR j=0L,n_elements(ind_tags)-1 DO BEGIN + data_struct(indices).(ind_tags(j)) = la_undef() + ENDFOR endelse - - ENDIF ELSE message, 'Cannot proceed. Keyword(s) missing. Refer to help for keyword options.' -return, data_struct - +RETURN, data_struct the_end: diff --git a/src/idl/dustem_myisrf_example.pro b/src/idl/dustem_myisrf_example.pro index 521f015..1533ae5 100644 --- a/src/idl/dustem_myisrf_example.pro +++ b/src/idl/dustem_myisrf_example.pro @@ -229,7 +229,7 @@ end ;== SET THE OBSERVATIONAL STRUCTURE ;== sed is passed twice -- the first occurrence is the SED that you ;== wish to fit, the second occurrence is the SED that you wish to visualise. -dustem_set_data,sed,sed +dustem_set_data,m_fit=sed,m_show=sed ; ;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE ; ;== ADJUSTED DURING THE FIT diff --git a/src/idl/dustem_set_data.pro b/src/idl/dustem_set_data.pro index 2a38f54..fa2a322 100755 --- a/src/idl/dustem_set_data.pro +++ b/src/idl/dustem_set_data.pro @@ -1,32 +1,31 @@ -PRO dustem_set_data,m_fit,m_show,x_fit,x_show,rchi2_weight=rchi2_weight,f_HI=f_HI,help=help +PRO dustem_set_data,m_fit=m_fit,m_show=m_show,x_fit=x_fit,x_show=x_show,rchi2_weight=rchi2_weight,f_HI=f_HI,help=help ;+ ; NAME: ; dustem_set_data ; ; PURPOSE: -; Initializes the dustem data structure (*!dustem_data) used in the fitting -; Initializes the !dustem_show structure which is either a copy or a more detailed version of !dustem_data -; +; Set the *!dustem data (and *!dustem_show) datasets in emission and/or extinction from input data structures +; *!dustem_data is important for the fitting of data and its general handling by DustEMWrap. ; CATEGORY: ; DustEMWrap, Mid-Level, Distributed, Initialization ; ; CALLING SEQUENCE: -; dustem_set_data[,m_fit,x_fit][,m_show,x_show][,rchi2_weight=][,f_HI=][,/help] +; dustem_set_data,[,m_fit=,x_fit=][,m_show=,x_show=],[,rchi2_weight=][,f_HI=][,/help] ; ; INPUTS: ; None ; ; OPTIONAL INPUT PARAMETERS: -; m_fit : dustemwrap structure contaning emission data to be fitted -; x_fit : dustemwrap structure containing extinction data to be fitted -; m_show : dustemwrap structure containing emission data to be displayed -; x_show : dustemwrap structure containing extinctiond data to be displayed -; rchi2_weight : dustemwrap structure containing the wieghts for chi2 calculation +; m_fit : DustEMWrap structure contaning emission data to be fitted +; m_show : DustEMWrap structure containing emission data to be displayed +; x_fit : DusteEMWrap structure containing extinction data to be fitted +; x_show : DustEMWrap structure containing extinctiond data to be displayed +; rchi2_weight : DustEMWrap structure containing the wieghts for chi2 calculation ; f_HI : multiplicative factor for the SED values (extinction and emission) ; ; OUTPUTS: -; None +; None/ALL OF THE ABOVE? ; ; OPTIONAL OUTPUT PARAMETERS: ; None @@ -52,15 +51,14 @@ PRO dustem_set_data,m_fit,m_show,x_fit,x_show,rchi2_weight=rchi2_weight,f_HI=f_H ; dir=!dustem_wrap_soft_dir+'/Data/EXAMPLE_OBSDATA/' ; file=dir+'example_SED_1.xcat' ; spec=read_xcat(file,/silent) -; dustem_set_data,spec,x_fit,spec,x_show) -; #in this example !dustem_data and !dustem_show are identical +; dustem_set_data,m_fit=spec,m_show=spec +; nb: in this example *!dustem_data and *!dustem_show end up identical ; ; MODIFICATION HISTORY: ; Written by J.-Ph. Bernard 2007 ; V. Guillet (2012) : dustem_set_data is turned into a function which can be used ; indifferently for sed, ext, polext -; I. Choubani (2022): dustem_check_data is introduced and the data structure format of the latest release is respected. -; dustem_set_data is turned into a procedure. +; I. Choubani (2022): dustem_set_data is turned into a procedure for emission/extinction with/without polarization. ; Evolution details on the DustEMWrap gitlab. ; See http://dustemwrap.irap.omp.eu/ for FAQ and help. ;- @@ -76,25 +74,18 @@ for i=0L,n_tags(*!dustem_data)-1 do begin (*!dustem_show).(i) = ptr_new() endfor -;#FITTING STRUCTURE - -;To me this is unnecessary... as dustem_check_data is always needed as it shapes the output structures (that would otherwise not exist) -;in order to prepare them for the mpfit run. -;the *!dustem_data tags are always cleaned so data will never exist if this keyword is used. -;commenting the use of nocheck_data for now - -;IF not keyword_set(nocheck_data) THEN BEGIN - dustem_check_data,m_fit,x_fit,sed,ext,polext,polsed,polfrac,psi_em,qsed,used,qext,uext,psi_ext,fpolext -;ENDIF - -;outputs whatever the input structure has, into organized seperate structures formatted for the fitting process. +;Outputs have the format required by !dustem_data and the subsequent fitting process. +dustem_check_data, m_fit,x_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,fpolext=fpolext ;FITTING STRUCTURE INITIALIZATION IF isa(sed) THEN (*!dustem_data).sed = ptr_new(sed) IF isa(ext) THEN (*!dustem_data).ext = ptr_new(ext) ;#################### -;OLD VARIABLE. KEPT IT. +;RELIC. For Annie to judge if this is necessary ; If f_HI is specified, multiply the data by f_HI If keyword_set(f_HI) and ptr_valid((*!dustem_data).sed) then begin if f_HI gt 0 then (*(*!dustem_data).sed).values *= f_HI else f_HI = 1. @@ -105,27 +96,33 @@ endif else f_HI = 1. defsysv,'!dustem_f_HI',f_HI ;#################### +;Problem here is that the user can activate !run_pol without wanting to be in polarization mode. +;reinitialze !run_pol here? + +If !run_pol THEN BEGIN + If tag_exist(*!dustem_show,'QSED') THEN BEGIN + IF ~isa((*!dustem_data).qsed) THEN !run_pol=0 + ENDIF ELSE !run_pol=0 +ENDIF + if !run_pol then begin if keyword_set(m_fit) then begin teststks = isa(qsed) and isa(used) - ;stop IF ~teststks then message, 'Stokes parameters are not set correctly. Please check the input structure. Aborting...' - - ;☆EMISSION message, 'Polarization component(s) to fit: QSED and USED' ,/continue - - ;☆Mandatory presence in polar mode IF isa(qsed) THEN BEGIN (*!dustem_data).qsed = ptr_new(qsed) + ;RELIC if keyword_set(f_HI) then if f_HI gt 0 then $ (*(*!dustem_data).qsed).values *= f_HI ENDIF - ;☆Mandatory presence in polar mode + if isa(used) then BEGIN (*!dustem_data).used = ptr_new(used) + ;RELIC if keyword_set(f_HI) then if f_HI gt 0 then $ (*(*!dustem_data).used).values *= f_HI ENDIF @@ -133,6 +130,7 @@ if !run_pol then begin if isa(polsed) then BEGIN (*!dustem_data).polsed = ptr_new(polsed) + ;RELIC if keyword_set(f_HI) then if f_HI gt 0 then $ (*(*!dustem_data).polsed).values *= f_HI @@ -141,6 +139,7 @@ if !run_pol then begin if isa(polfrac) then BEGIN (*!dustem_data).polfrac = ptr_new(polfrac) + ;RELIC if keyword_set(f_HI) then if f_HI gt 0 then $ (*(*!dustem_data).polfrac).values *= f_HI @@ -149,23 +148,22 @@ if !run_pol then begin if isa(psi_em) then BEGIN (*!dustem_data).psi_em = ptr_new(psi_em) + ;RELIC if keyword_set(f_HI) then if f_HI gt 0 then $ (*(*!dustem_data).psi_em).values *= f_HI ENDIF - endif if keyword_set(x_fit) then begin - ;Two tests that determine what the user will eventually fit testexstks = isa(qext) and isa(uext) 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. - + ;EXTINCTION - + message, 'Extinction polarization component(s) to fit: QEXT and UEXT',/continue if isa(polext) then BEGIN @@ -200,22 +198,22 @@ if !run_pol then begin endif endif -;#SHOWING STRUCTURE +;SHOWING STRUCTURE / IT IS NEEDED TO BE ABLE TO LOCATE THE HIDDEN DATA POINTS WHEN COMPARING WITH !dustem_data. +;when creating the HUGE system variable JP wants to create. We'll be able to omit it from DusteEMWrap for good. if keyword_set(m_fit) and not keyword_set(m_show) then m_show=m_fit if keyword_set(x_fit) and not keyword_set(x_show) then x_show = x_fit - -;Setting of the structure that will be displayed -;IF not keyword_set(nocheck_data) THEN BEGIN - dustem_check_data,m_show,x_show,sed,ext,polext,polsed,polfrac,psi_em,qsed,used,qext,uext,psi_ext,fpolext -;ENDIF +dustem_check_data, m_show,x_show,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,fpolext=fpolext IF isa(sed) THEN (*!dustem_show).sed = ptr_new(sed) IF isa(ext) THEN (*!dustem_show).ext = ptr_new(ext) ;#################### -;OLD VARIABLE. KEPT IT. +;RELIC ; If f_HI is specified, multiply the data by f_HI If keyword_set(f_HI) and ptr_valid((*!dustem_show).sed) then begin if f_HI gt 0 then (*(*!dustem_show).sed).values *= f_HI else f_HI = 1. @@ -234,19 +232,6 @@ if !run_pol then begin 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 ~(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 ~(isa(polsed)) and ~(isa(polfrac)) and isa(psi_em) then nowshowing = 'QSED, USED and PSI_EM' - - - message, 'Polarization component(s) to show: '+nowshowing ,/continue - - If isa(polsed) then BEGIN ;Apply changes to !dustem_data here. (*!dustem_show).polsed = ptr_new(polsed) if keyword_set(f_HI) then if f_HI gt 0 then $ @@ -277,19 +262,7 @@ if !run_pol then begin teststks = isa(qext) and isa(uext) 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 ~(isa(polextfrac)) and ~(isa(psi_ext)) then nowshowing = 'QEXT, UEXT and POLEXT' - if teststks and isa(polext) and ~(isa(polextfrac)) and isa(psi_ext) then nowshowing = 'QEXT, UEXT, POLEXT and PSI_EXT' - if teststks and ~(isa(polext)) and isa(polextfrac) and ~(isa(psi_ext)) then nowshowing = 'QEXT, UEXT and POLEXTFRAC' - if teststks and ~(isa(polext)) and isa(polextfrac) and isa(psi_ext) then nowshowing = 'QEXT, UEXT, POLEXTFRAC and PSI_EXT' - if teststks and ~(isa(polext)) and ~(isa(polextfrac)) and ~(isa(psi_ext)) then nowshowing = 'QEXT and UEXT' - if teststks and isa(polext) and isa(polextfrac) and ~(isa(psi_ext)) then nowshowing = 'QEXT, UEXT, POLEXT and POLEXTFRAC' - if teststks and isa(polext) and isa(polextfrac) and isa(psi_ext) then nowshowing = 'QEXT, UEXT, POLEXTFRAC, POLEXT and PSI_EXT' - if teststks and ~(isa(polext)) and ~(isa(polextfrac)) and isa(psi_ext) then nowshowing = 'QEXT, UEXT and PSI_EXT' - - - message, 'Polarization component(s) to show: '+nowshowing ,/continue - + ;EXTINCTION if isa(polext) then BEGIN (*!dustem_show).polext = ptr_new(polext) @@ -364,10 +337,4 @@ ENDIF the_end: -; -; return,0. - - - - END diff --git a/src/idl/dustem_stellarpopisrf_example.pro b/src/idl/dustem_stellarpopisrf_example.pro index 4865a8b..80565eb 100644 --- a/src/idl/dustem_stellarpopisrf_example.pro +++ b/src/idl/dustem_stellarpopisrf_example.pro @@ -218,7 +218,7 @@ end for i=4l,n_tags(sed)-1 do begin sed.(i) = sed.sigmaii endfor -dustem_set_data,sed,sed +dustem_set_data,m_fit=sed,m_show=sed ;;== SET INITIAL VALUES AND LIMITS OF THE PARAMETERS THAT WILL BE ;;== ADJUSTED DURING THE FIT diff --git a/src/idl/dustemwrap_plot.pro b/src/idl/dustemwrap_plot.pro index 242c06d..76e1b9c 100755 --- a/src/idl/dustemwrap_plot.pro +++ b/src/idl/dustemwrap_plot.pro @@ -154,8 +154,7 @@ if not keyword_set(dustem_ext) and isa((*!dustem_data).ext) then begin dustem_ext = dustem_compute_ext(p_dim,st=st,EXT_spec=EXT_spec) if !run_pol and tag_exist(*!dustem_show,'qext') && !run_lin then begin ;any polarization tag (qsed,used,qext,usext,polfrac) can be used for this test - - ;stop + dustem_polext = dustem_compute_polext(p_dim,st=st,POLEXT_spec=POLEXT_spec,SPEXT_spec=SPEXT_spec,dustem_fpolext=dustem_fpolext) toto = dustem_compute_stokext(p_dim,st=st,QEXT_spec=QEXT_spec,UEXT_spec=UEXT_spec,PSIEXT_spec=PSIEXT_spec,dustem_psi_ext=dustem_psi_ext) dustem_qext = toto[0] @@ -260,6 +259,14 @@ endif else begin rchi2 = la_undef();maybe not the best initialization endelse +;changing tag_exist() by isa() because the user can mistakenly activate the !run_pol keyword + +If !run_pol THEN BEGIN + If tag_exist(*!dustem_show,'QSED') THEN BEGIN + IF ~isa((*!dustem_data).qsed) THEN !run_pol=0 + ENDIF ELSE !run_pol = 0 +ENDIF + ;===========Window positionning parameters if !run_pol and tag_exist(*!dustem_show,'qsed') then BEGIN ;qsed tag is sufficient for the test wdelta_x=60 @@ -500,7 +507,7 @@ if test_m then begin dustem_interp = list(dustem_sed) ;DustEmWrap SED predictions . dustem_spec = list(SED_spec) ;DustEM spectra predictions - position = list(p_sed,np_sed);list(list(p_sed,np_sed)) ; I think I need this ... + position = list(p_sed,np_sed) ;list(list(p_sed,np_sed)) ; I think I need this ... ;position_n =list(np_sed) datasets = ['SED'] ;replace by one keyword? @@ -516,7 +523,7 @@ if test_m then begin endif else begin - ; PLOTTING OF DATA THAT IS UNCHAED + ; PLOTTING OF DATA THAT IS UNCHANGED cgwindow,'dustem_plot_dataset', st, dataset=datasets, position=position,positive_only=positiveonly,negative_only=negativeonly, /addcmd, winid=winid_m, _extra=_extra & cmdind_m+=1 ; REFRESHING DATA - to plot the refreshed data for the first time. cgwindow,'dustem_plot_dataset', st, dustem_interp,dustem_spec,extra_spec=extra_spec,dataset=datasets,positive_only=positiveonly,negative_only=negativeonly, /refresh ,position=position, /addcmd, winid=winid_m, _extra=_extra & cmdind_m+=1 -- libgit2 0.21.2