diff --git a/src/idl/dustem_check_data.pro b/src/idl/dustem_check_data.pro index bfa7503..d596662 100644 --- a/src/idl/dustem_check_data.pro +++ b/src/idl/dustem_check_data.pro @@ -1,39 +1,34 @@ -FUNCTION dustem_check_data, sed=sed,ext=ext,fpol=fpol,ist=ist,pst=pst,fpst=fpst,qst=qst,ust=ust,exst=exst,pexst=pexst,qexst=qexst,uexst=uexst - - -;THIS PROCEDURE CHECKS THE DATA PROVIDED BY THE USER AND RETURNS DIFFERENT STRUCTURES TO DEFINE (SET) THE COMPOSITE DATA FITTING STRUCTURE (!dustem_data) - - -;write help section here -;I don't think this function should return anything -;This replaces a big portion of dustem_set_data, JP was right. +FUNCTION dustem_check_data,data_sed=data_sed,data_ext=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 ;CHECKING THE FORMAT OF THE SED STRUCTURE -IF KEYWORD_SET(sed) THEN BEGIN +IF KEYWORD_SET(data_sed) THEN BEGIN + ans='' - wavs=sed.wave + wavs=data_sed.wave ;=== Impose central wavelengths for photometric channels - ind=where(sed.filter NE 'SPECTRUM',count) + ind=where(data_sed.filter NE 'SPECTRUM',count) IF count NE 0 THEN BEGIN - wavs(ind)=dustem_filter2wav(sed(ind).filter) + wavs(ind)=dustem_filter2wav(data_sed(ind).filter) ENDIF - ;=== define observations structure - obs_st={instru_names:sed.instru,filt_names:sed.filter,wav:wavs, values:sed.StokesI,sigma:sqrt(abs(sed.sigmaII))} + ;=== Initializing observation structure + ;it doesn't matter that data.StokesI is used, further test will take care of what is necessary. + obs_EM={instru_names:data_sed.instru,filt_names:data_sed.filter,wav:wavs, values:data_sed.StokesI,sigma:sqrt(abs(data_sed.sigmaII))} ;INITIALIZING OUTPUT STRUCTURES - ist = obs_st - pst = ist - fpst = ist - qst = ist - ust = ist - ;This allows the user to fit polsed, qsed and used individually. + sed = obs_EM + polsed = sed + polfrac = sed + qsed = sed + used = sed + psi_em = sed + ;Test if there is StokesI data. If not move to polsed testing. - ind = where(sed.StokesI NE la_undef() and finite(sed.StokesI) EQ 1 and sed.sigmaII NE la_undef() and finite(sed.sigmaII) EQ 1,N_sed) ;Testing if the data and errors are set + ind = where(data_sed.StokesI NE la_undef() and finite(data_sed.StokesI) EQ 1 and data_sed.sigmaII NE la_undef() and finite(data_sed.sigmaII) EQ 1,N_sed) ;Testing if the data and errors are set IF N_sed EQ 0 THEN BEGIN - ist = ptr_new() + sed = ptr_new() goto, pzone ENDIF ELSE BEGIN - ind = where(finite(sed.StokesI) EQ 0 or finite(sed.sigmaII) EQ 0, countnoki) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + ind = where(finite(data_sed.StokesI) EQ 0 or finite(data_sed.sigmaII) EQ 0, countnoki) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value IF countnoki NE 0 THEN BEGIN message, "SED format isn't met. NaN value(s) detected in StokesI data.",/continue message, 'Without automatic format fitering you cannot proceed.',/continue @@ -45,27 +40,29 @@ IF KEYWORD_SET(sed) THEN BEGIN ENDIF IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ ELSE message, 'Data format is met. Removal of undefined data points...',/continue - testi=where(sed.sigmaII NE 0. and finite(sed.sigmaII) EQ 1 and sed.sigmaII NE la_undef() and sed.StokesI NE la_undef() and finite(sed.StokesI) EQ 1,ctesti) + tdaerri=where((data_sed.sigmaII EQ la_undef() and data_sed.StokesI NE la_undef()) or (data_sed.sigmaII NE la_undef() and data_sed.StokesI EQ la_undef()),ctdaerri) + IF ctdaerri NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...' + testi=where(data_sed.sigmaII NE 0. and finite(data_sed.sigmaII) EQ 1 and data_sed.sigmaII NE la_undef() and data_sed.StokesI NE la_undef() and finite(data_sed.StokesI) EQ 1,ctesti) IF ctesti NE 0 THEN BEGIN - instru_names = ist.instru_names(testi) - filt_names = ist.filt_names(testi) - wav = ist.wav(testi) - values = ist.values(testi) - sigma = ist.sigma(testi) - ist={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + instru_names = sed.instru_names(testi) + filt_names = sed.filt_names(testi) + wav = sed.wav(testi) + values = sed.values(testi) + sigma = sed.sigma(testi) + sed={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} ENDIF ENDELSE pzone: ans='' - pst.values=sed.largeP - pst.sigma=sqrt(abs(sed.sigma_largeP)) - ;Testig if there is LargeP data. If not move to polfrac data (fpol) - ind = where(sed.largeP NE la_undef() and finite(sed.largeP) EQ 1 and sed.sigma_largeP NE la_undef() and finite(sed.sigma_largeP) EQ 1,NP) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta) + polsed.values=data_sed.largeP + polsed.sigma=sqrt(abs(data_sed.sigma_largeP)) + ;Testing if there is LargeP data. If not move to polfrac data. + ind = where(data_sed.largeP NE la_undef() and finite(data_sed.largeP) EQ 1 and data_sed.sigma_largeP NE la_undef() and finite(data_sed.sigma_largeP) EQ 1,NP) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta) IF NP EQ 0 THEN BEGIN - pst = ptr_new() - goto, fpzone + polsed = ptr_new() + goto, psiemzone ENDIF ELSE BEGIN - ind = where(finite(sed.largeP) EQ 0 or finite(sed.sigma_largeP) EQ 0, countnokp) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + ind = where(finite(data_sed.largeP) EQ 0 or finite(data_sed.sigma_largeP) EQ 0, countnokp) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value IF countnokp NE 0 THEN BEGIN message, "SED format isn't met. NaN value(s) detected in LargeP data.",/continue message, 'Without automatic format fitering you cannot proceed.',/continue @@ -77,28 +74,64 @@ IF KEYWORD_SET(sed) THEN BEGIN ENDIF IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ ELSE message, 'Data format is met. Removal of undefined data points...',/continue - testp=where(sed.sigma_largeP NE 0. and finite(sed.sigma_largeP) EQ 1 and sed.sigma_largeP NE la_undef() and sed.largeP NE la_undef() and finite(sed.largeP) EQ 1,ctestp) + tdaerrp=where((data_sed.sigma_largeP EQ la_undef() and data_sed.largeP NE la_undef()) or (data_sed.sigma_largeP NE la_undef() and data_sed.largeP EQ la_undef()),ctdaerrp) + IF ctdaerrp NE 0 THEN message, 'Data and errors do not match. Undefinded value(s) detected. Aborting...' + testp=where(data_sed.sigma_largeP NE 0. and finite(data_sed.sigma_largeP) EQ 1 and data_sed.sigma_largeP NE la_undef() and data_sed.largeP NE la_undef() and finite(data_sed.largeP) EQ 1,ctestp) IF ctestp NE 0 THEN BEGIN - instru_names = pst.instru_names(testp) - filt_names = pst.filt_names(testp) - wav = pst.wav(testp) - values = pst.values(testp) - sigma = pst.sigma(testp) - pst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + instru_names = polsed.instru_names(testp) + filt_names = polsed.filt_names(testp) + wav = polsed.wav(testp) + values = polsed.values(testp) + sigma = polsed.sigma(testp) + polsed={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} ENDIF - ENDELSE + ENDELSE ; The deprecation of the fitting of POLFRAC does not happen here. But rather in dustem_mpfit_run + psiemzone: + ans='' + psi_em.values = data_sed.psi + psi_em.sigma = sqrt(abs(data_sed.sigma_psi)) + ;Testing if there is PSI (EMISSION) data. If not move to polfrac data. + ind = where(data_sed.psi NE la_undef() and finite(data_sed.psi) EQ 1 and data_sed.sigma_psi NE la_undef() and finite(data_sed.sigma_psi) EQ 1,NPSIEM) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta) + IF NPSIEM EQ 0 THEN BEGIN + psi_em = ptr_new() + goto, fpzone + ENDIF ELSE BEGIN + ind = where(finite(data_sed.psi) EQ 0 or finite(data_sed.sigma_psi) EQ 0, countnokpsiem) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + IF countnokpsiem NE 0 THEN BEGIN + message, "SED format isn't met. NaN value(s) detected in LargeP data.",/continue + message, 'Without automatic format fitering you cannot proceed.',/continue + read, ans, prompt='Would you like automatic format filtering? (Y/N):' + IF strupcase(ans) EQ 'N' THEN BEGIN + message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue + message, 'Aborting ...' + ENDIF + ENDIF + IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ + ELSE message, 'Data format is met. Removal of undefined data points...',/continue + tdaerrpsiem=where((data_sed.sigma_psi EQ la_undef() and data_sed.psi NE la_undef()) or (data_sed.sigma_psi NE la_undef() and data_sed.psi EQ la_undef()),ctdaerrpsiem) + IF ctdaerrpsiem NE 0 THEN message, 'Data and errors do not match. Undefinded value(s) detected. Aborting...' + testpsiem=where(data_sed.sigma_psi NE 0. and finite(data_sed.sigma_psi) EQ 1 and data_sed.sigma_psi NE la_undef() and data_sed.psi NE la_undef() and finite(data_sed.psi) EQ 1,ctestpsiem) + IF ctestpsiem NE 0 THEN BEGIN + instru_names = psi_em.instru_names(testpsiem) + filt_names = psi_em.filt_names(testpsiem) + wav = psi_em.wav(testpsiem) + values = psi_em.values(testpsiem) + sigma = psi_em.sigma(testpsiem) + psi_em={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + ENDIF + ENDELSE ; The deprecation of POLFRAC does not happen here. More like in dustem_mpfit_run fpzone: ans='' - fpst.values=sed.smallp - fpst.sigma=sqrt(abs(sed.sigma_smallp)) + polfrac.values=data_sed.smallp + polfrac.sigma=sqrt(abs(data_sed.sigma_smallp)) ;Testig if there is smallp data. If not move to qsed data - ind = where(sed.smallP NE la_undef() and finite(sed.smallP) EQ 1 and sed.sigma_smallP NE la_undef() and finite(sed.sigma_smallP) EQ 1,NSP) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta) + ind = where(data_sed.smallP NE la_undef() and finite(data_sed.smallP) EQ 1 and data_sed.sigma_smallP NE la_undef() and finite(data_sed.sigma_smallP) EQ 1,NSP) ;Testing on the nan values is needed to move to testing the other data fields (because if the dataset is full of Nans but not la_undef(), the code will think it's valid dta) IF NSP EQ 0 THEN BEGIN - fpst = ptr_new() + polfrac = ptr_new() ;goto, qzone goto, qzone ENDIF ELSE BEGIN - ind = where(finite(sed.smallP) EQ 0 or finite(sed.sigma_smallP) EQ 0, countnoksp) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + ind = where(finite(data_sed.smallP) EQ 0 or finite(data_sed.sigma_smallP) EQ 0, countnoksp) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value IF countnoksp NE 0 THEN BEGIN message, "SED format isn't met. NaN value(s) detected in SmallP data.",/continue message, 'Without automatic format fitering you cannot proceed.',/continue @@ -110,27 +143,29 @@ IF KEYWORD_SET(sed) THEN BEGIN ENDIF IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ ELSE message, 'Data format is met. Removal of undefined data points...',/continue - testsp=where(sed.sigma_smallP NE 0. and finite(sed.sigma_smallP) EQ 1 and sed.sigma_smallP NE la_undef() and sed.smallP NE la_undef() and finite(sed.smallP) EQ 1,ctestsp) + tdaerrsp=where((data_sed.sigma_smallp EQ la_undef() and data_sed.smallP NE la_undef()) or (data_sed.sigma_smallp NE la_undef() and data_sed.smallP EQ la_undef()),ctdaerrsp) + IF ctdaerrsp NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...' + testsp=where(data_sed.sigma_smallP NE 0. and finite(data_sed.sigma_smallP) EQ 1 and data_sed.sigma_smallP NE la_undef() and data_sed.smallP NE la_undef() and finite(data_sed.smallP) EQ 1,ctestsp) IF ctestsp NE 0 THEN BEGIN - instru_names = fpst.instru_names(testsp) - filt_names = fpst.filt_names(testsp) - wav = fpst.wav(testsp) - values = fpst.values(testsp) - sigma = fpst.sigma(testsp) - fpst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + instru_names = polfrac.instru_names(testsp) + filt_names = polfrac.filt_names(testsp) + wav = polfrac.wav(testsp) + values = polfrac.values(testsp) + sigma = polfrac.sigma(testsp) + polfrac={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} ENDIF ENDELSE qzone: ans='' - qst.values=sed.StokesQ - qst.sigma=sqrt(abs(sed.sigmaQQ)) + qsed.values=data_sed.StokesQ + qsed.sigma=sqrt(abs(data_sed.sigmaQQ)) ;Testig if there is StokesQ data. If not move to used data - ind = where(sed.StokesQ NE la_undef() and finite(sed.StokesQ) EQ 1 and sed.sigmaQQ NE la_undef() and finite(sed.sigmaQQ) EQ 1,NQ) + ind = where(data_sed.StokesQ NE la_undef() and finite(data_sed.StokesQ) EQ 1 and data_sed.sigmaQQ NE la_undef() and finite(data_sed.sigmaQQ) EQ 1,NQ) IF NQ EQ 0 THEN BEGIN - qst = ptr_new() + qsed = ptr_new() goto, uzone ENDIF ELSE BEGIN - ind = where(finite(sed.StokesQ) EQ 0 or finite(sed.sigmaQQ) EQ 0, countnokq) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + ind = where(finite(data_sed.StokesQ) EQ 0 or finite(data_sed.sigmaQQ) EQ 0, countnokq) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value IF countnokq NE 0 THEN BEGIN message, "SED format isn't met. NaN value(s) detected in StokesQ data.",/continue message, 'Without automatic format fitering you cannot proceed.',/continue @@ -142,27 +177,29 @@ IF KEYWORD_SET(sed) THEN BEGIN ENDIF IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ ELSE message, 'Data format is met. Removal of undefined data points...',/continue - testq=where(sed.sigmaQQ NE 0. and finite(sed.sigmaQQ) EQ 1 and sed.sigmaQQ NE la_undef() and sed.StokesQ NE la_undef() and finite(sed.StokesQ) EQ 1,ctestq) + tdaerrq=where((data_sed.sigmaQQ EQ la_undef() and data_sed.StokesQ NE la_undef()) or (data_sed.sigmaQQ NE la_undef() and data_sed.StokesQ EQ la_undef()),ctdaerrq) + IF ctdaerrq NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...' + testq=where(data_sed.sigmaQQ NE 0. and finite(data_sed.sigmaQQ) EQ 1 and data_sed.sigmaQQ NE la_undef() and data_sed.StokesQ NE la_undef() and finite(data_sed.StokesQ) EQ 1,ctestq) IF ctestq NE 0 THEN BEGIN - instru_names = qst.instru_names(testq) - filt_names = qst.filt_names(testq) - wav = qst.wav(testq) - values = qst.values(testq) - sigma = qst.sigma(testq) - qst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + instru_names = qsed.instru_names(testq) + filt_names = qsed.filt_names(testq) + wav = qsed.wav(testq) + values = qsed.values(testq) + sigma = qsed.sigma(testq) + qsed={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} ENDIF ENDELSE uzone: ans='' - ust.values=sed.StokesU - ust.sigma=sqrt(abs(sed.sigmaUU)) + used.values=data_sed.StokesU + used.sigma=sqrt(abs(data_sed.sigmaUU)) ;Testig if there is StokesU data. If not return empty value and end loop. - ind = where(sed.StokesU NE la_undef() and finite(sed.StokesU) EQ 1 and sed.sigmaUU NE la_undef() and finite(sed.sigmaUU) EQ 1,NU) + ind = where(data_sed.StokesU NE la_undef() and finite(data_sed.StokesU) EQ 1 and data_sed.sigmaUU NE la_undef() and finite(data_sed.sigmaUU) EQ 1,NU) IF NU EQ 0 THEN BEGIN - ust = ptr_new() + used = ptr_new() goto, end_sed ENDIF ELSE BEGIN - ind = where(finite(sed.StokesU) EQ 0 or finite(sed.sigmaUU) EQ 0, countnoku) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + ind = where(finite(data_sed.StokesU) EQ 0 or finite(data_sed.sigmaUU) EQ 0, countnoku) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value IF countnoku NE 0 THEN BEGIN message, "SED format isn't met. NaN value(s) detected in StokesU data.",/continue message, 'Without automatic format fitering you cannot proceed.',/continue @@ -174,46 +211,46 @@ IF KEYWORD_SET(sed) THEN BEGIN ENDIF IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ ELSE message, 'Data format is met. Removal of undefined data points...',/continue - testu=where(sed.sigmaUU NE 0. and finite(sed.sigmaUU) EQ 1 and sed.sigmaUU NE la_undef() and sed.StokesU NE la_undef() and finite(sed.StokesU) EQ 1,ctestu) + tdaerru=where((data_sed.sigmaUU EQ la_undef() and data_sed.StokesU NE la_undef()) or (data_sed.sigmaUU NE la_undef() and data_sed.StokesU EQ la_undef()),ctdaerru) + IF ctdaerru NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...' + testu=where(data_sed.sigmaUU NE 0. and finite(data_sed.sigmaUU) EQ 1 and data_sed.sigmaUU NE la_undef() and data_sed.StokesU NE la_undef() and finite(data_sed.StokesU) EQ 1,ctestu) IF ctestu NE 0 THEN BEGIN - instru_names = ust.instru_names(testu) - filt_names = ust.filt_names(testu) - wav = ust.wav(testu) - values = ust.values(testu) - sigma = ust.sigma(testu) - ust={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + instru_names = used.instru_names(testu) + filt_names = used.filt_names(testu) + wav = used.wav(testu) + values = used.values(testu) + sigma = used.sigma(testu) + used={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} ENDIF ENDELSE end_sed: ENDIF - -;CHECKING THE FORMAT OF THE EXT STRUCTURE -IF KEYWORD_SET(ext) THEN begin +IF KEYWORD_SET(data_ext) THEN begin + ans='' - wavs=ext.wave + wavs=data_ext.wave ;=== Impose central wavelengths for photometric channels - ind=where(ext.filter NE 'SPECTRUM',count) + ind=where(data_ext.filter NE 'SPECTRUM',count) IF count NE 0 THEN BEGIN - wavs(ind)=dustem_filter2wav(ext(ind).filter) + wavs(ind)=dustem_filter2wav(data_ext(ind).filter) ENDIF ;=== define observations structure - obs_st={instru_names:ext.instru,filt_names:ext.filter,wav:wavs, values:ext.EXT_I,sigma:sqrt(abs(ext.sigextII))} + obs_EXT={instru_names:data_ext.instru,filt_names:data_ext.filter,wav:wavs, values:data_ext.EXT_I,sigma:sqrt(abs(data_ext.sigextII))} ;INITIALIZING OUTPUT STRUCTURES - exst = obs_st - pexst = exst - qexst = exst - uexst = exst - ;This allows the user to fit polext, qext and uext individually. - ;Test if there is extinction data. If not move to polext testing. - ;this condition was updated because you cannot have data without errors and errors without data - ind = where(ext.EXT_I NE la_undef() and finite(ext.EXT_I) EQ 1 and ext.sigextII NE la_undef() and finite(ext.sigextII) EQ 1,N_ext) ;Testing if the data and errors are set + ext = obs_EXT + polext = ext + qext = ext + uext = ext + psi_ext = ext + + ind = where(data_ext.EXT_I NE la_undef() and finite(data_ext.EXT_I) EQ 1 and data_ext.sigextII NE la_undef() and finite(data_ext.sigextII) EQ 1,N_ext) ;Testing if the data and errors are set IF N_ext EQ 0 THEN BEGIN - exst = ptr_new() + ext = ptr_new() goto, pexzone ENDIF ELSE BEGIN - ind = where(finite(ext.EXT_I) EQ 0 or finite(ext.sigextII) EQ 0, countnokex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + ind = where(finite(data_ext.EXT_I) EQ 0 or finite(data_ext.sigextII) EQ 0, countnokex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value IF countnokex NE 0 THEN BEGIN message, "SED format isn't met. NaN value(s) detected in EXT_I data.",/continue message, 'Without automatic format fitering you cannot proceed.',/continue @@ -226,29 +263,31 @@ IF KEYWORD_SET(ext) THEN begin ENDIF IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ ELSE message, 'Data format is met. Removal of undefined data points...',/continue - testex=where(ext.sigextII NE 0. and finite(ext.sigextII) EQ 1 and ext.sigextII NE la_undef() and ext.EXT_I NE la_undef() and finite(ext.EXT_I) EQ 1,ctestex) + tdaerrex=where((data_ext.sigextII EQ la_undef() and data_ext.EXT_I NE la_undef()) or (data_ext.sigextII NE la_undef() and data_ext.EXT_I EQ la_undef()),ctdaerrex) + IF ctdaerrex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...' + testex=where(data_ext.sigextII NE 0. and finite(data_ext.sigextII) EQ 1 and data_ext.sigextII NE la_undef() and data_ext.EXT_I NE la_undef() and finite(data_ext.EXT_I) EQ 1,ctestex) IF ctestex NE 0 THEN BEGIN - instru_names = exst.instru_names(testex) - filt_names = exst.filt_names(testex) - wav = exst.wav(testex) - values = exst.values(testex) - sigma = exst.sigma(testex) - exst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + instru_names = ext.instru_names(testex) + filt_names = ext.filt_names(testex) + wav = ext.wav(testex) + values = ext.values(testex) + sigma = ext.sigma(testex) + ext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} ENDIF ENDELSE pexzone: ans='' - pexst.values=ext.EXT_P - pexst.sigma=sqrt(abs(ext.sigextP)) + polext.values=data_ext.EXT_P + polext.sigma=sqrt(abs(data_ext.sigextP)) ;Testig if there is EXT_P data. If not move to EXT_Q testing - ind = where(ext.EXT_P NE la_undef() and finite(ext.EXT_P) EQ 1 and ext.sigextP NE la_undef() and finite(ext.sigextII) EQ 1,NP_ext) ;Testing if the data and errors are set + ind = where(data_ext.EXT_P NE la_undef() and finite(data_ext.EXT_P) EQ 1 and data_ext.sigextP NE la_undef() and finite(data_ext.sigextP) EQ 1,NP_ext) ;Testing if the data and errors are set IF NP_ext EQ 0 THEN BEGIN - pexst = ptr_new() + polext = ptr_new() goto, qexzone ENDIF ELSE BEGIN - ind = where(finite(ext.EXT_P) EQ 0 or finite(ext.sigextP) EQ 0, countnokpex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + ind = where(finite(data_ext.EXT_P) EQ 0 or finite(data_ext.sigextP) EQ 0, countnokpex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value IF countnokpex NE 0 THEN BEGIN - message, "SED format isn't met. NaN value(s) detected in EXT_I data.",/continue + message, "SED format isn't met. NaN value(s) detected in EXT_P data.",/continue message, 'Without automatic format fitering you cannot proceed.',/continue ans='' read, ans, prompt='Would you like automatic format filtering? (Y/N):' @@ -259,30 +298,67 @@ IF KEYWORD_SET(ext) THEN begin ENDIF IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ ELSE message, 'Data format is met. Removal of undefined data points...',/continue - testpex=where(ext.sigextP NE 0. and finite(ext.sigextP) EQ 1 and ext.sigextP NE la_undef() and ext.EXT_P NE la_undef() and finite(ext.EXT_P) EQ 1,ctestpex) + tdaerrpex=where((data_ext.sigextP EQ la_undef() and data_ext.EXT_P NE la_undef()) or (data_ext.sigextP NE la_undef() and data_ext.EXT_P EQ la_undef()),ctdaerrpex) + IF ctdaerrpex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...' + testpex=where(data_ext.sigextP NE 0. and finite(data_ext.sigextP) EQ 1 and data_ext.sigextP NE la_undef() and data_ext.EXT_P NE la_undef() and finite(data_ext.EXT_P) EQ 1,ctestpex) IF ctestpex NE 0 THEN BEGIN - instru_names = pexst.instru_names(testpex) - filt_names = pexst.filt_names(testpex) - wav = pexst.wav(testpex) - values = pexst.values(testpex) - sigma = pexst.sigma(testpex) - pexst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + instru_names = polext.instru_names(testpex) + filt_names = polext.filt_names(testpex) + wav = polext.wav(testpex) + values = polext.values(testpex) + sigma = polext.sigma(testpex) + polext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + ENDIF + ENDELSE + psiexzone: + ans='' + psi_ext.values=data_ext.psi + psi_ext.sigma=sqrt(abs(data_ext.sigextpsi)) + ;Testig if there is EXT_P data. If not move to EXT_Q testing + ind = where(data_ext.psi NE la_undef() and finite(data_ext.psi) EQ 1 and data_ext.sigextpsi NE la_undef() and finite(data_ext.sigextpsi) EQ 1,NPSIEXT) ;Testing if the data and errors are set + IF NPSIEXT EQ 0 THEN BEGIN + psi_ext = ptr_new() + goto, qexzone + ENDIF ELSE BEGIN + ind = where(finite(data_ext.psi) EQ 0 or finite(data_ext.sigextpsi) EQ 0, countnokpsiex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + IF countnokpsiex NE 0 THEN BEGIN + message, "SED format isn't met. NaN value(s) detected in PSI_EXT data.",/continue + message, 'Without automatic format fitering you cannot proceed.',/continue + ans='' + read, ans, prompt='Would you like automatic format filtering? (Y/N):' + IF strupcase(ans) EQ 'N' THEN BEGIN + message,'You have chosen not to use automatic filtering. Please adjust the SED format.',/continue + message, 'Aborting ...' + ENDIF + ENDIF + IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ + ELSE message, 'Data format is met. Removal of undefined data points...',/continue + tdaerrpsiex=where((data_ext.sigextpsi EQ la_undef() and data_ext.psi NE la_undef()) or (data_ext.sigextpsi NE la_undef() and data_ext.psi EQ la_undef()),ctdaerrpsiex) + IF ctdaerrpsiex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...' + testpsiex=where(data_ext.sigextpsi NE 0. and finite(data_ext.sigextpsi) EQ 1 and data_ext.sigextpsi NE la_undef() and data_ext.psi NE la_undef() and finite(data_ext.psi) EQ 1,ctestpsiex) + IF ctestpsiex NE 0 THEN BEGIN + instru_names = psi_ext.instru_names(testpsiex) + filt_names = psi_ext.filt_names(testpsiex) + wav = psi_ext.wav(testpsiex) + values = psi_ext.values(testpsiex) + sigma = psi_ext.sigma(testpsiex) + psi_ext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} ENDIF ENDELSE qexzone: ans='' - qexst.values=ext.EXT_Q - qexst.sigma=sqrt(abs(ext.sigextQQ)) + qext.values=data_ext.EXT_Q + qext.sigma=sqrt(abs(data_ext.sigextQQ)) ;Testig if there is EXT_P data. If not move to EXT_U testing - ind = where(ext.EXT_Q NE la_undef() and finite(ext.EXT_Q) EQ 1 and ext.sigextQQ NE la_undef() and finite(ext.sigextQQ) EQ 1,NQ_ext) ;Testing if the data and errors are set + ind = where(data_ext.EXT_Q NE la_undef() and finite(data_ext.EXT_Q) EQ 1 and data_ext.sigextQQ NE la_undef() and finite(data_ext.sigextQQ) EQ 1,NQ_ext) ;Testing if the data and errors are set ;dustem_set_data.prostop IF NQ_ext EQ 0 THEN BEGIN - qexst = ptr_new() + qext = ptr_new() goto, uexzone ENDIF ELSE BEGIN - ind = where(finite(ext.EXT_Q) EQ 0 or finite(ext.sigextQQ) EQ 0, countnokqex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + ind = where(finite(data_ext.EXT_Q) EQ 0 or finite(data_ext.sigextQQ) EQ 0, countnokqex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value IF countnokqex NE 0 THEN BEGIN - message, "SED format isn't met. NaN value(s) detected in EXT_I data.",/continue + message, "SED format isn't met. NaN value(s) detected in EXT_Q data.",/continue message, 'Without automatic format fitering you cannot proceed.',/continue ans='' read, ans, prompt='Would you like automatic format filtering? (Y/N):' @@ -293,29 +369,31 @@ IF KEYWORD_SET(ext) THEN begin ENDIF IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ ELSE message, 'Data format is met. Removal of undefined data points...',/continue - testqex=where(ext.sigextQQ NE 0. and finite(ext.sigextQQ) EQ 1 and ext.sigextQQ NE la_undef() and ext.EXT_Q NE la_undef() and finite(ext.EXT_Q) EQ 1,ctestqex) + tdaerrqex=where((data_ext.sigextQQ EQ la_undef() and data_ext.EXT_Q NE la_undef()) or (data_ext.sigextQQ NE la_undef() and data_ext.EXT_Q EQ la_undef()),ctdaerrqex) + IF ctdaerrqex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...' + testqex=where(data_ext.sigextQQ NE 0. and finite(data_ext.sigextQQ) EQ 1 and data_ext.sigextQQ NE la_undef() and data_ext.EXT_Q NE la_undef() and finite(data_ext.EXT_Q) EQ 1,ctestqex) IF ctestqex NE 0 THEN BEGIN - instru_names = qexst.instru_names(testqex) - filt_names = qexst.filt_names(testqex) - wav = qexst.wav(testqex) - values = qexst.values(testqex) - sigma = qexst.sigma(testqex) - qexst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + instru_names = qext.instru_names(testqex) + filt_names = qext.filt_names(testqex) + wav = qext.wav(testqex) + values = qext.values(testqex) + sigma = qext.sigma(testqex) + qext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} ENDIF ENDELSE uexzone: ans='' - uexst.values=ext.EXT_U - uexst.sigma=sqrt(abs(ext.sigextUU)) + uext.values=data_ext.EXT_U + uext.sigma=sqrt(abs(data_ext.sigextUU)) ;Testig if there is EXT_P data. If not return empty value and end loop. - ind = where(ext.EXT_U NE la_undef() and finite(ext.EXT_U) EQ 1 and ext.sigextUU NE la_undef() and finite(ext.sigextUU) EQ 1,NU_ext) ;Testing if the data and errors are set + ind = where(data_ext.EXT_U NE la_undef() and finite(data_ext.EXT_U) EQ 1 and data_ext.sigextUU NE la_undef() and finite(data_ext.sigextUU) EQ 1,NU_ext) ;Testing if the data and errors are set IF NU_ext EQ 0 THEN BEGIN - uexst = ptr_new() + uext = ptr_new() goto, end_ext ENDIF ELSE BEGIN - ind = where(finite(ext.EXT_U) EQ 0 or finite(ext.sigextUU) EQ 0, countnokuex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value + ind = where(finite(data_ext.EXT_U) EQ 0 or finite(data_ext.sigextUU) EQ 0, countnokuex) ;Format is not okay if there are nan values but the wrapper recognizes the la_undef() value IF countnokuex NE 0 THEN BEGIN - message, "SED format isn't met. NaN value(s) detected in EXT_I data.",/continue + message, "SED format isn't met. NaN value(s) detected in EXT_U data.",/continue message, 'Without automatic format fitering you cannot proceed.',/continue ans='' read, ans, prompt='Would you like automatic format filtering? (Y/N):' @@ -326,21 +404,22 @@ IF KEYWORD_SET(ext) THEN begin ENDIF IF strupcase(ans) EQ 'Y' THEN message, 'You have chosen automatic format filtering. Removal of undefined data points...',/continue $ ELSE message, 'Data format is met. Removal of undefined data points...',/continue - testuex=where(ext.sigextUU NE 0. and finite(ext.sigextUU) EQ 1 and ext.sigextUU NE la_undef() and ext.EXT_U NE la_undef() and finite(ext.EXT_U) EQ 1,ctestuex) + tdaerruex=where((data_ext.sigextUU EQ la_undef() and data_ext.EXT_U NE la_undef()) or (data_ext.sigextUU NE la_undef() and data_ext.EXT_U EQ la_undef()),ctdaerruex) + IF ctdaerruex NE 0 THEN message, 'Data and errors do not match. Undefinded data detected. Aborting...' + testuex=where(data_ext.sigextUU NE 0. and finite(data_ext.sigextUU) EQ 1 and data_ext.sigextUU NE la_undef() and data_ext.EXT_U NE la_undef() and finite(data_ext.EXT_U) EQ 1,ctestuex) IF ctestuex NE 0 THEN BEGIN - instru_names = uexst.instru_names(testuex) - filt_names = uexst.filt_names(testuex) - wav = uexst.wav(testuex) - values = uexst.values(testuex) - sigma = uexst.sigma(testuex) - uexst={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} + instru_names = uext.instru_names(testuex) + filt_names = uext.filt_names(testuex) + wav = uext.wav(testuex) + values = uext.values(testuex) + sigma = uext.sigma(testuex) + uext={instru_names:instru_names,filt_names:filt_names,wav:wav,values:values,sigma:sigma} ENDIF ENDELSE end_ext: ENDIF - - -RETURN, 0. ;stop +RETURN, 0. + END diff --git a/src/idl/dustem_compute_polext.pro b/src/idl/dustem_compute_polext.pro index 8017806..512de67 100644 --- a/src/idl/dustem_compute_polext.pro +++ b/src/idl/dustem_compute_polext.pro @@ -10,7 +10,6 @@ IF not keyword_set(sti) THEN BEGIN sti=dustem_run(p_dim) ENDIF - If ptr_valid(*!dustem_plugin) THEN BEGIN ;ADDING PLUGIN TO SPECTRUM---------------- ;if n_tags(!dustem_data.polext) gt 1 then begin @@ -20,7 +19,6 @@ If ptr_valid(*!dustem_plugin) THEN BEGIN ;endif ;------------------------------------------ - ; ;=============================THIS BLOCK IS WRONG YOU HAVE TO CHANGE IT============================= ; ; @@ -51,8 +49,6 @@ If ptr_valid(*!dustem_plugin) THEN BEGIN ENDIF - - dustem_polext = interpol(sti.polext.ext_tot,sti.polext.wav,(*!dustem_data.polext).wav) diff --git a/src/idl/dustem_compute_polsed.pro b/src/idl/dustem_compute_polsed.pro index b67ac1d..5f4181a 100755 --- a/src/idl/dustem_compute_polsed.pro +++ b/src/idl/dustem_compute_polsed.pro @@ -59,6 +59,9 @@ frac=P/stI tes=where(finite(frac) eq 0) frac(tes)=0. +;THIS LOOP IS WRONG AND DOES NOT MAKE SENSE + + ;IF ptr_valid(!dustem_plugin) THEN BEGIN scopes=tag_names((*!dustem_scope)) IF scopes[0] NE 'NONE' THEN BEGIN @@ -77,7 +80,7 @@ FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN ENDFOR -IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(0.,Nwaves) +IF ~isa(Q_spec) && ~isa(U_spec) THEN polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(!dustem_psi,Nwaves) FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN @@ -119,7 +122,6 @@ ENDIF ;stop - ;For spectrum data points, interpolate in log-log. ;Linear interpolation leads to wrong values, in particular where few ;wavelengths points exist in the model (long wavelengths). diff --git a/src/idl/dustem_compute_sed.pro b/src/idl/dustem_compute_sed.pro index 3b8b28c..32cce1c 100755 --- a/src/idl/dustem_compute_sed.pro +++ b/src/idl/dustem_compute_sed.pro @@ -62,7 +62,6 @@ spec = st.sed.em_tot * fact ;if n_tags(!dustem_data.sed) gt 1 then begin ; executing only when !dustem_data.sed is set meaning that this allows for the creation of fake seds - scopes=tag_names((*!dustem_scope)) IF scopes[0] NE 'NONE' THEN BEGIN ;IF ptr_valid(!dustem_plugin) THEN BEGIN @@ -70,6 +69,9 @@ IF scopes[0] NE 'NONE' THEN BEGIN if total(strsplit((*(*!dustem_scope).(i)),'+',/extract) eq 'ADD_SED') then spec+=(*(*!dustem_plugin).(i))[*,0] endfor ENDIF + + + ;endif ;------------------------------------------ @@ -113,7 +115,7 @@ the_end: ;VG : error in cc IRAS1 when no PAH j=where(dustem_sed ne dustem_sed,count) if count gt 0 then dustem_sed[j]=0. - +;stop RETURN,dustem_sed END diff --git a/src/idl/dustem_compute_stokes.pro b/src/idl/dustem_compute_stokes.pro index f93be65..6b1f508 100755 --- a/src/idl/dustem_compute_stokes.pro +++ b/src/idl/dustem_compute_stokes.pro @@ -51,7 +51,7 @@ out=0. If keyword_set(result) THEN result=st -fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(st.sed).wav)*1.e20/1.e7 +fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/(st.sed).wav)*1.e20/1.e7 stI = st.sed.em_tot * fact ;This is Total intensity I P = st.polsed.em_tot * fact ;This is Polarized intensity P @@ -63,7 +63,6 @@ tes=where(finite(frac) eq 0) frac(tes)=0. - FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_QSED') THEN BEGIN @@ -78,7 +77,11 @@ FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN ENDFOR IF ~isa(Q_spec) && ~isa(U_spec) THEN BEGIN - polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(0,Nwaves) + ;maybe add a condition that takes into account the angles of the stokes parameters. + ;I believe it is best to hard code the angle information in here + ;this also includes its extraction? + + polar_ippsi2iqu,stI,Q_spec,U_spec,frac,replicate(!dustem_psi,Nwaves) ;stop ENDIF diff --git a/src/idl/dustem_init.pro b/src/idl/dustem_init.pro index c819cc3..7c351e1 100755 --- a/src/idl/dustem_init.pro +++ b/src/idl/dustem_init.pro @@ -107,16 +107,16 @@ if keyword_set(pol) then begin sed: ptr_new(), $ ext: ptr_new(), $ polext: ptr_new(), $ - polsed: ptr_new(), $ - polfrac: ptr_new(), $ + ;polsed: ptr_new(), $ ;newly deprecated + ;polfrac: ptr_new(), $ ;newly deprecated qsed: ptr_new(), $ used: ptr_new(), $ qext: ptr_new(), $ uext: ptr_new() $ } - rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.,qsed:0.,used:0.,qext:0.,uext:0.} - + ;rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.,qsed:0.,used:0.,qext:0.,uext:0.} ;polsed and polfrac newly deprecated + rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,qsed:0.,used:0.,qext:0.,uext:0.} endif else begin defsysv, '!dustem_data', { $ ;Data to fit sed: ptr_new(), $ @@ -134,6 +134,24 @@ ENDFOR instr+='}' toto=execute(instr) + +defsysv,'!dustem_show', { $ ;Data to fit + sed: ptr_new(), $ + ext: ptr_new(), $ + polext: ptr_new(), $ + polsed: ptr_new(), $ + polfrac: ptr_new(), $ + psi_em: ptr_new(), $ + psi_ext: ptr_new(), $ + qsed: ptr_new(), $ + used: ptr_new(), $ + qext: ptr_new(), $ + uext: ptr_new() $ + } + +;defsysv, '!dustem_psi',0. ;TO BE DEPRECTED... CHECK IF U STARTED IMPLEMENTING THIS ;0 is default. Meanin No PSI value(s) detected in the data xcat file. + + ;IDL> help,*!dustem_data,/str ;** Structure <2173b20>, 4 tags, length=448, data length=448, refs=1: ; INSTRU_NAMES STRING Array[14] @@ -284,6 +302,8 @@ ENDELSE st_model=dustem_read_all(dir_in,/silent,kwords=kwords) !dustem_params=ptr_new(st_model) +;ADD A last wavelength here. + ;=== create dynamical storage if needed IF not file_test(!dustem_dat,/dir) THEN BEGIN spawn,'mkdir '+!dustem_dat diff --git a/src/idl/dustem_mpfit_run.pro b/src/idl/dustem_mpfit_run.pro index f450fa0..17a9254 100755 --- a/src/idl/dustem_mpfit_run.pro +++ b/src/idl/dustem_mpfit_run.pro @@ -58,8 +58,8 @@ st_model=dustem_read_all(!dustem_soft_dir) chi2_sed = 0 chi2_ext = 0 chi2_polext = 0 -chi2_polsed = 0 -chi2_polfrac = 0 +;chi2_polsed = 0 +;chi2_polfrac = 0 chi2_qsed =0 chi2_used =0 chi2_qext =0 @@ -67,8 +67,8 @@ chi2_uext =0 rchi2_sed = 0 rchi2_ext = 0 rchi2_polext = 0 -rchi2_polsed = 0 -rchi2_polfrac = 0 +;rchi2_polsed = 0 +;rchi2_polfrac = 0 rchi2_qsed =0 rchi2_used =0 rchi2_qext =0 @@ -76,8 +76,8 @@ rchi2_uext =0 wrchi2_sed = 0 wrchi2_ext = 0 wrchi2_polext = 0 -wrchi2_polsed = 0 -wrchi2_polfrac = 0 +;wrchi2_polsed = 0 +;wrchi2_polfrac = 0 wrchi2_qsed =0 wrich2_used =0 wrchi2_qext =0 @@ -85,8 +85,8 @@ wrich2_uext =0 n_sed = 0 n_ext = 0 n_polext = 0 -n_polsed = 0 -n_polfrac = 0 +;n_polsed = 0 +;n_polfrac = 0 n_qsed =0 n_used =0 n_qext =0 @@ -110,7 +110,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN CASE tagnames(ind_data_tag(i_tag)) OF 'SED' : BEGIN - dustem_sed = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron) + dustem_sed = dustem_compute_sed(p_dim,st) chi2_sed = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2) n_sed = n_elements((*!dustem_data.sed).values) rchi2_sed = chi2_sed / n_sed @@ -130,90 +130,13 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN (*!dustem_fit).rchi2=rchi2_sed ; Plot if needed IF !dustem_show_plot NE 0 THEN BEGIN - - ;title='DUSTEM WRAP('+strtrim(string(win),2)+')' + window,win,title='DUSTEM WRAP (SED)' ;cgwindow,WMULTI=[2,2,1],/CURRENT - ;I.C - ;undeprecated the use of polfrac here - ;********** - ;I added this block so I can plot polfrac inside the sed window - ;I will probably remove it - ;********** - ;test if polfrac is in tagnames: - ;another problem is that sed and pol data don't have the same filters. - varvar=where(tagnames eq 'POLFRAC', testfpol) - if testfpol ne 0 then 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. - ;I'm not so sure about this... - - 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 - - ;I think I will have errors here!! The thing is I'm trying to make this work without (*!dustem_data.polsed) which is hard. - - 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 - ;stop - ;stop - 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 - endif - ;stop + 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 - ;stop ENDIF dustem_out = [ dustem_out, dustem_sed ] @@ -355,11 +278,67 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN END -;USE OF POLFRAC DEPRECATED--------------------------------------- - -; 'POLFRAC': BEGIN -; -; ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_POLSED AND DUSTEM_COMPUTE_SED +; 'POLFRAC': BEGIN ;Deprecate this. Fitting POLFRAC does not make any sense. It is a biased quantity. +; 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 +; +; 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 +;OLD DEPRECATED VERSION. ; ; IF !run_lin then begin ; @@ -401,45 +380,48 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN ;----------------------------------------------------------------- - 'POLSED': BEGIN - - ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED - - IF !run_lin THEN BEGIN - - 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 +; 'POLSED': 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] ; - 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 +; 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 + '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) + stop chi2_qsed =total(((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2) n_qsed = n_elements((*!dustem_data.qsed).values) rchi2_qsed = chi2_qsed / n_qsed @@ -455,16 +437,15 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN 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 +; 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,2e4] - normpol=Q_spec*fact - normpol1=dustem_qsed*fact1 + xr=[10,1e6] + ;normpol1=(*!dustem_data.qsed).values*fact1 yr=[1e-28,1e-17] ;yr=[1e-34,1e-23] @@ -473,10 +454,14 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN ylog=1 if countq ne 0 then begin - ylog=0 - yr=[-2e-31,5e-31] +; 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)' @@ -491,46 +476,46 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN 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 + END 'USED': BEGIN toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec) - chi2_used =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2) + stop + 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 + wrchi2_used = rchi2_used * !fit_rchi2_weight.used 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 - - - 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 + print,"chi2 per wl: ",((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2 + + + 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 + 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,6e4] - normpol2=U_spec*fact + 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 + ;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' @@ -547,10 +532,10 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN ;cgwindow, 'plot', - endif - END - - + endif + END +; +; 'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION @@ -670,16 +655,16 @@ if !fit_rchi2_weight.polext ne 0 then begin 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.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 @@ -708,8 +693,8 @@ total_wchi2 = chi2_ext * !fit_rchi2_weight.ext $ + chi2_sed * !fit_rchi2_weight.sed if !run_pol then begin total_wchi2 += chi2_polext * !fit_rchi2_weight.polext $ - + chi2_polsed * !fit_rchi2_weight.polsed $ - + chi2_polfrac * !fit_rchi2_weight.polfrac $ +; + chi2_polsed * !fit_rchi2_weight.polsed $ +; + chi2_polfrac * !fit_rchi2_weight.polfrac $ + chi2_qsed * !fit_rchi2_weight.qsed $ + chi2_used * !fit_rchi2_weight.used $ + chi2_qext * !fit_rchi2_weight.qext $ @@ -720,9 +705,9 @@ message,'======================================================================' print,"Total weighted chi2 driving mpfitfun = ", total_wchi2 print,"Nb of data points = ", total_n print,"DOF = ", DOF -print,"Reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",rchi2_ext,rchi2_sed,rchi2_polext,rchi2_polsed,rchi2_polfrac -print,"Weighted reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",wrchi2_ext,wrchi2_sed,wrchi2_polext,wrchi2_polsed,wrchi2_polfrac -if !run_pol then print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed,!fit_rchi2_weight.polext,!fit_rchi2_weight.polsed,!fit_rchi2_weight.polfrac $ +print,"Reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",rchi2_ext,rchi2_sed,rchi2_polext;,rchi2_polsed,rchi2_polfrac +print,"Weighted reduced chi2 EXT/SED/POLEXT/POLSED/POLFRAC",wrchi2_ext,wrchi2_sed,wrchi2_polext;,wrchi2_polsed,wrchi2_polfrac +if !run_pol then print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed,!fit_rchi2_weight.polext $;,!fit_rchi2_weight.polsed,!fit_rchi2_weight.polfrac $ else print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF) print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF) diff --git a/src/idl/dustem_plot_fit_sed.pro b/src/idl/dustem_plot_fit_sed.pro index 3b6e34a..43da5ce 100755 --- a/src/idl/dustem_plot_fit_sed.pro +++ b/src/idl/dustem_plot_fit_sed.pro @@ -46,6 +46,9 @@ ; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems. ;- + +;NOta bene: only delete the cgoplot commands + ;stop IF keyword_set(help) THEN BEGIN @@ -136,7 +139,7 @@ IF keyword_set(ytit) THEN ytit = ytit ELSE ytit = textoidl('Brightness/N_H (MJy/ ;deffo define a title variable here. The procedures are not communicating with each other. - +;pos=cgLayout() ;############################### if keyword_set(fpol) then begin if !run_lin then begin diff --git a/src/idl/dustem_plot_polsed.pro b/src/idl/dustem_plot_polsed.pro index 82e057f..68d7c95 100644 --- a/src/idl/dustem_plot_polsed.pro +++ b/src/idl/dustem_plot_polsed.pro @@ -1,4 +1,4 @@ - PRO dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win +PRO dustem_plot_polsed, st, p_dim, dustem_polsed, aligned=aligned, win=win IF keyword_set(ps) THEN BEGIN diff --git a/src/idl/dustem_plugin_modify_dust_pol.pro b/src/idl/dustem_plugin_modify_dust_pol.pro index c6c0a21..5373e32 100644 --- a/src/idl/dustem_plugin_modify_dust_pol.pro +++ b/src/idl/dustem_plugin_modify_dust_pol.pro @@ -65,7 +65,6 @@ P=((st.polsed).em_tot)*fact I=((st.sed).em_tot)*fact - Nwaves=(size(I))[1] frac=P/I*smallp diff --git a/src/idl/dustem_plugin_stellar_population.pro b/src/idl/dustem_plugin_stellar_population.pro index c915fe3..173d560 100644 --- a/src/idl/dustem_plugin_stellar_population.pro +++ b/src/idl/dustem_plugin_stellar_population.pro @@ -58,7 +58,7 @@ countf = 0. countg = 0. ;============================== - +;REPLACE THIS WHITH A PLUGIN PARAMETER THAT CAN BE FITTED ;===========Specifying whether to add the Mathis field (and its scaling via G0) to the composite ISRF=================== ismathis=0 ;this value will have to be manually changed at the present moment defsysv, '!ismathis', ismathis @@ -3177,14 +3177,12 @@ c(Ncomments-2)='# Nbr of points' c(Ncomments-1)='# wave (microns), 4*pi*Inu (erg/cm2/s/Hz)' FOR i=0L,n_elements(comp_pop.popid)-1 DO BEGIN ; Looping over all the stellar populations - - ;The initial procedure had omega multiplied by a !pi factor. Its presence in the Planck (Astron) procedure makes for a good reason to discuss this with J.P. - + omega = ((comp_pop.radius)[i]/((comp_pop.distance)[i]*pc2cm))^2 ; Dilution factor of a stellar population - Inu = planck(wave_angstrom,(comp_pop.temperature)[i])/(4.*!pi)*(wave_angstrom)^2/c2a ; ergs/cm2/s/Hz/sr + Inu = planck(wave_angstrom,(comp_pop.temperature)[i])*(wave_angstrom)^2/c2a ; ergs/cm2/s/Hz/sr - stellar_component=stellar_component+(comp_pop.nstars)[i]*omega*Inu + stellar_component=stellar_component+(comp_pop.nstars)[i]*omega*Inu ;no *4 because of Lambert's cosine law. ;Rest of the lines of the new composite ISRF.DAT file @@ -3199,10 +3197,15 @@ FOR i=0L,n_elements(comp_pop.popid)-1 DO BEGIN ; Looping over all the stellar po ENDELSE ENDFOR +;make sure which g0 is used!!! +;alwayd devide by g0 even though it equals one. +; -IF !ismathis THEN st.isrf=ma_isrf.isrf+stellar_component/(st.gas.G0) ELSE st.isrf=stellar_component;/((*!dustem_params).gas.G0) +IF !ismathis THEN st.isrf=ma_isrf.isrf+stellar_component/(st.gas.G0) ELSE st.isrf=stellar_component;/(st.gas.G0) -;don't mind my tests I'm just making sure the ISRF file changes at all +;stop + +;don't mind my tests I'm just making sure if the ISRF file changes at all print,'stellar_component is:' print, stellar_component @@ -3229,6 +3232,10 @@ close,unit free_lun,unit out=st.isrf + +;THIS BLOCK IS FALSE AND NEEDS TO BE CORRECTED. ASK JP WHAT HE WANTS FOR A DEFAULT RUN +;(a run with no parameters for instance) + ;============This block creates a composite (9V) stellar population structure for when the function isn't used as a plugin=================== ENDIF ELSE BEGIN diff --git a/src/idl/dustem_sed_plot.pro b/src/idl/dustem_sed_plot.pro index f5ec1e7..674d35a 100755 --- a/src/idl/dustem_sed_plot.pro +++ b/src/idl/dustem_sed_plot.pro @@ -1,4 +1,5 @@ PRO dustem_sed_plot,p_min,_extra=_extra,function_name=function_name,pol=pol,legend_xpos=legend_xpos,legend_ypos=legend_ypos,ps=ps,png=png,fpol=fpol +;THIS procedure is incomplete. It needs to be updated. IF not keyword_set(function_name) THEN BEGIN ;Run dustem with as an interface to mpfitfun; what does this mean? This whole procedure will have to be changed diff --git a/src/idl/dustem_set_data.pro b/src/idl/dustem_set_data.pro index e7571de..438cbee 100755 --- a/src/idl/dustem_set_data.pro +++ b/src/idl/dustem_set_data.pro @@ -1,4 +1,8 @@ -FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,fpol=fpol,rchi2_weight=rchi2_weight,f_HI=f_HI +FUNCTION dustem_set_data,st_sed_fit=st_sed_fit,st_sed_show=st_sed_show,st_ext_fit=st_ext_fit,st_ext_show=st_ext_show,rchi2_weight=rchi2_weight,f_HI=f_HI,help=help + + +;update the help + ;+ ; NAME: @@ -8,7 +12,7 @@ FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,fpol=fpol,rchi2_weight=rchi2_ ; CATEGORY: ; Dustem ; CALLING SEQUENCE: -; dustem_set_data,spec[,help=] +; dustem_set_data(data_sed=data_sed,data_ext=data_ext,[/help]) ; INPUTS: ; spec: array of structures containing a SED ; OPTIONAL INPUT PARAMETERS: @@ -42,8 +46,8 @@ FUNCTION dustem_set_data,help=help,sed=sed,ext=ext,fpol=fpol,rchi2_weight=rchi2_ ; ; V. Guillet (2012) : dustem_set_data is turned into a function which can be used ; indifferently for sed, ext, polext -; I. Choubani (2022): Adjusting the function to meet the format criteria set by the current release. -; The function uses two main input structures: SED and EXT. Polarization data is included within each. +; I. Choubani (2022): Introducing the dustem_check_data procedure where the new format outlined by the latest release is taken into account. +; This function uses two main input structures: DATA_SED and DATA_EXT. Polarization data is included within each. ;- @@ -52,132 +56,240 @@ IF keyword_set(help) THEN BEGIN goto,the_end ENDIF -;Cleaninig data structure -!dustem_data.sed = ptr_new() -!dustem_data.ext = ptr_new() +;Cleaninig data/showing structure +for i=0L,n_tags(!dustem_data)-1 do begin + !dustem_data.(i) = ptr_new() + !dustem_show.(i) = ptr_new() +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) +;outputs whatever the input structure has, into organized seperate structures formatted for the fitting process. -;SETTING THE DATA STRUCTURES FOR THE FITTING PROCESS - SED +;If SED/EXT data is present in the structure it is fitted automatically. +;If the user wished to not fit it they will have to null the corresponding fields in the input fitting structure. +IF isa(sed) THEN !dustem_data.sed = ptr_new(sed) +IF isa(ext) THEN !dustem_data.ext = ptr_new(ext) -IF keyword_set(sed) THEN BEGIN + +;#################### +;OLD VARIABLE. KEPT IT. +; 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. +endif else f_HI = 1. +if keyword_set(f_HI) and ptr_valid(!dustem_data.ext) then begin + if f_HI gt 0 then (*!dustem_data.ext).values *= f_HI else f_HI = 1. +endif else f_HI = 1. +defsysv,'!dustem_f_HI',f_HI +;#################### + +;NOTA BENE: The user is bound to fit Q and U together. +if !run_pol then begin + if keyword_set(st_sed_fit) then begin - ;maybe to activate fpol here you need to set it to one? CHECK THIS OUT. - if keyword_set(fpol) then xx=dustem_check_data(sed=sed,ist=ist,pst=pst,fpol=fpol,fpst=fpst,qst=qst,ust=ust) else xx=dustem_check_data(sed=sed,ist=ist,pst=pst,fpst=fpst,qst=qst,ust=ust) - ;stop - ;=== fill !dustem_data - IF isa(ist) THEN !dustem_data.sed = ptr_new(ist) ;ELSE !dustem_data.sed = ptr_new() - ; 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. - endif else f_HI = 1. - defsysv,'!dustem_f_HI',f_HI - if !run_pol then begin - If isa(pst) then BEGIN - ;=== fill !dustem_data - !dustem_data.polsed = ptr_new(pst) - ; If f_HI is specified, multiply the data by f_HI - if keyword_set(f_HI) then if f_HI gt 0 then $ - (*!dustem_data.polsed).values *= f_HI - ENDIF - If isa(fpst) then BEGIN - ;=== fill !dustem_data - !dustem_data.polfrac = ptr_new(fpst) - ; If f_HI is specified, multiply the data by f_HI - if keyword_set(f_HI) then if f_HI gt 0 then $ - (*!dustem_data.polfrac).values *= f_HI - ENDIF - IF isa(qst) THEN BEGIN - ;=== fill !dustem_data - !dustem_data.qsed = ptr_new(qst) - ; If f_HI is specified, multiply the data by f_HI + teststks = isa(qsed) and isa(used) + + IF (not 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 + + 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(ust) then BEGIN - ;=== fill !dustem_data - !dustem_data.used = ptr_new(ust) - ; If f_HI is specified, multiply the data by f_HI + 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 endif -ENDIF - -;SETTING THE DATA STRUCTURES FOR THE FITTING PROCESS - EXT - -IF keyword_set(ext) THEN BEGIN - xx=dustem_check_data(ext=ext,exst=exst,pexst=pexst,qexst=qexst,uexst=uexst) - ;=== fill !dustem_data - IF isa(exst) THEN !dustem_data.ext = ptr_new(exst) ;ELSE !dustem_data.ext = ptr_new() - - ; If f_HI is specified, multiply the data by f_HI - if keyword_set(f_HI) and ptr_valid(!dustem_data.ext) then begin - if f_HI gt 0 then (*!dustem_data.ext).values *= f_HI else f_HI = 1. - endif else f_HI = 1. - defsysv,'!dustem_f_HI',f_HI - if !run_pol then begin - if isa(pexst) then BEGIN - !dustem_data.polext = ptr_new(pexst) - ; If f_HI is specified, multiply the data by f_HI + 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...' + + ;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 + if isa(polext) then BEGIN + !dustem_data.polext = ptr_new(polext) if keyword_set(f_HI) then if f_HI gt 0 then $ (*!dustem_data.polext).values *= f_HI ENDIF - IF isa(qexst) THEN BEGIN - ;=== fill !dustem_data - !dustem_data.qext = ptr_new(qexst) - ; If f_HI is specified, multiply the data by f_HI + IF isa(qext) THEN BEGIN + !dustem_data.qext = ptr_new(qext) if keyword_set(f_HI) then if f_HI gt 0 then $ (*!dustem_data.qext).values *= f_HI ENDIF - if isa(uexst) then BEGIN - ;=== fill !dustem_data - !dustem_data.uext = ptr_new(uexst) - ; If f_HI is specified, multiply the data by f_HI + if isa(uext) then BEGIN + !dustem_data.uext = ptr_new(uext) if keyword_set(f_HI) then if f_HI gt 0 then $ (*!dustem_data.uext).values *= f_HI ENDIF endif -ENDIF +endif + + +if not keyword_set(st_sed_show) then st_sed_show=st_sed_fit -!fit_rchi2_weight.sed=1. -!fit_rchi2_weight.ext=1. +;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) -If !run_pol THEN BEGIN - ;SED - IF isa(pst) THEN !fit_rchi2_weight.polsed=1. - IF isa(fpst) THEN !fit_rchi2_weight.polfrac=1. - IF isa(qst) THEN !fit_rchi2_weight.qsed=1. - IF isa(ust) THEN !fit_rchi2_weight.used=1. +IF isa(sed) THEN !dustem_show.sed = ptr_new(sed) +IF isa(ext) THEN !dustem_show.ext = ptr_new(ext) + + +;#################### +;OLD VARIABLE. KEPT IT. +; 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. +endif else f_HI = 1. +if keyword_set(f_HI) and ptr_valid(!dustem_show.ext) then begin + if f_HI gt 0 then (*!dustem_show.ext).values *= f_HI else f_HI = 1. +endif else f_HI = 1. +defsysv,'!dustem_f_HI',f_HI +;#################### + +if !run_pol then begin - ;EXT - IF isa(pexst) THEN !fit_rchi2_weight.polext=1. - IF isa(qexst) THEN !fit_rchi2_weight.qext=1. - IF isa(uexst) THEN !fit_rchi2_weight.uext=1. -ENDIF + if isa(st_sed_show) then begin ; because sometimes this keyword isn't set. + + 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) 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, 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' + + 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 $ + (*!dustem_show.polsed).values *= f_HI + ENDIF + If isa(polfrac) then BEGIN + !dustem_show.polfrac = ptr_new(polfrac) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_show.polfrac).values *= f_HI + ENDIF + IF isa(qsed) THEN BEGIN + !dustem_show.qsed = ptr_new(qsed) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_show.qsed).values *= f_HI + ENDIF + if isa(used) then BEGIN + !dustem_show.used = ptr_new(used) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_show.used).values *= f_HI + ENDIF + if isa(psi_em) then BEGIN + !dustem_show.psi_em = ptr_new(psi_em) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_show.psi_em).values *= f_HI + ENDIF + endif + 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 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, 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' + + message, 'Polarization component(s) to show: '+nowshowing ,/continue + + ;EXTINCTION + if isa(polext) then BEGIN + !dustem_show.polext = ptr_new(polext) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_show.polext).values *= f_HI + ENDIF + IF isa(qext) THEN BEGIN + !dustem_show.qext = ptr_new(qext) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_show.qext).values *= f_HI + ENDIF + if isa(uext) then BEGIN + !dustem_show.uext = ptr_new(uext) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_show.uext).values *= f_HI + ENDIF + if isa(psi_ext) then BEGIN + !dustem_show.psi_ext = ptr_new(psi_ext) + if keyword_set(f_HI) then if f_HI gt 0 then $ + (*!dustem_show.psi_ext).values *= f_HI + ENDIF + endif +endif IF keyword_set(rchi2_weight) THEN BEGIN !fit_rchi2_weight.sed = rchi2_weight.sed !fit_rchi2_weight.ext = rchi2_weight.ext IF !run_pol THEN BEGIN - IF isa(pst) THEN !fit_rchi2_weight.polsed = rchi2_weight.polsed - IF isa(fpst) THEN !fit_rchi2_weight.polfrac = rchi2_weight.polfrac - IF isa(qst) THEN !fit_rchi2_weight.qsed = rchi2_weight.qsed - IF isa(ust) THEN !fit_rchi2_weight.used = rchi2_weight.used + ;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(pexst) THEN !fit_rchi2_weight.polext = rchi2_weight.polext - IF isa(qexst) THEN !fit_rchi2_weight.qext = rchi2_weight.qext - IF isa(uexst) THEN !fit_rchi2_weight.uext = rchi2_weight.uext + 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 -ENDIF +ENDIF ELSE BEGIN + + !fit_rchi2_weight.sed=1. + !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 -return,ptr_new(obs_st) +ENDELSE the_end: +;stop +return,0. + + + + + END -- libgit2 0.21.2