OV-GSO-DC / dustem-wrapper_idl

Commit 3196a378f3dc9f99b9509f34b6ed0a2a0a3c27cf

Authored by Ilyes Choubani
1 parent 7817238c
Exists in master

Complementing fits handling of extiction data... Frozen data are not yet incorporated.

Showing 5 changed files with 232 additions and 115 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_fillup_systvar_from_fits.pro
  Diff comments   View file @3196a37
1 1 PRO dustem_fillup_systvar_from_fits,syst_var,str_input_SED,used_pol
2   -
3   - Nsed=n_elements(str_input_SED)
  2 +
  3 + ;IC: some way for the procedure to be aware of emission and extinction
  4 + test_m = fix(total(tag_names (str_input_sed) EQ 'STOKESI')) ;a proper boolean
  5 +
  6 + Nsed=n_elements(str_input_SED)
4 7 sed=dustem_initialize_internal_sed(Nsed)
5   -
6   - (syst_var).sed=ptr_new(sed)
7   -
8   - (*syst_var.sed).instru_names=dustem_filter2instru(str_input_SED.filter)
9   - (*syst_var.sed).filt_names=strtrim(str_input_SED.filter,2)
10   - (*syst_var.sed).wav=str_input_SED.wavelength
11   - (*syst_var.sed).values=str_input_SED.STOKESI
12   - (*syst_var.sed).sigma=la_power(str_input_SED.varianceII,0.5)
13   -
  8 +
  9 + IF test_m THEN BEGIN
  10 + (syst_var).sed=ptr_new(sed)
  11 + (*syst_var.sed).instru_names=dustem_filter2instru(str_input_SED.filter)
  12 + (*syst_var.sed).filt_names=strtrim(str_input_SED.filter,2)
  13 + (*syst_var.sed).wav=str_input_SED.wavelength
  14 + (*syst_var.sed).values=str_input_SED.STOKESI
  15 + (*syst_var.sed).sigma=la_power(str_input_SED.varianceII,0.5)
  16 + ENDIF ELSE BEGIN
  17 + (syst_var).ext=ptr_new(sed)
  18 + (*syst_var.ext).instru_names=dustem_filter2instru(str_input_SED.filter)
  19 + (*syst_var.ext).filt_names=strtrim(str_input_SED.filter,2)
  20 + (*syst_var.ext).wav=str_input_SED.wavelength
  21 + (*syst_var.ext).values=str_input_SED.STOKESI_EXT
  22 + (*syst_var.ext).sigma=la_power(str_input_SED.varianceII_EXT,0.5)
  23 + ENDELSE
  24 +
14 25 IF used_pol THEN BEGIN
15 26 ;stop
16   - ;===== note: those should include only points with polarization data, not all sed points
17   - indq=where(str_input_SED.STOKESQ NE la_undef(),Nqsed)
18   - indu=where(str_input_SED.STOKESU NE la_undef(),Nused)
19   - qsed=dustem_initialize_internal_sed(Nqsed)
20   - used=dustem_initialize_internal_sed(Nused)
21   - indpolar=where(str_input_SED.STOKESU NE la_undef() AND str_input_SED.STOKESQ NE la_undef(),Npolar)
22   -
23   - psi_sed=dustem_initialize_internal_sed(Npolar)
24   - smallp_sed=dustem_initialize_internal_sed(Npolar)
25   - largep_sed=dustem_initialize_internal_sed(Npolar)
26   -
27   - syst_var.qsed=ptr_new(qsed)
28   - syst_var.used=ptr_new(used)
29   - (*syst_var.qsed).instru_names=dustem_filter2instru(str_input_SED[indq].filter)
30   - (*syst_var.qsed).filt_names=strtrim(str_input_SED[indq].filter,2)
31   - (*syst_var.qsed).wav=str_input_SED[indq].wavelength
32   - (*syst_var.qsed).values=str_input_SED[indq].STOKESQ
33   - (*syst_var.qsed).sigma=la_power(str_input_SED[indq].varianceQQ,0.5)
34   - ;now U
35   - (*syst_var.used).instru_names=dustem_filter2instru(str_input_SED[indu].filter)
36   - (*syst_var.used).filt_names=strtrim(str_input_SED[indu].filter,2)
37   - (*syst_var.used).wav=str_input_SED[indu].wavelength
38   - (*syst_var.used).values=str_input_SED[indu].STOKESU
39   - (*syst_var.used).sigma=la_power(str_input_SED[indu].varianceUU,0.5)
40   - ;=== add up dependent fields
41   - syst_var.psi_em=ptr_new(psi_sed)
42   - syst_var.polfrac=ptr_new(smallp_sed)
43   - syst_var.polsed=ptr_new(largep_sed)
44   - polar_iqu2ippsi,str_input_SED[indpolar].STOKESI,str_input_SED[indpolar].STOKESQ,str_input_SED[indpolar].STOKESU,p,psi,largep=largep
45   - (*syst_var.psi_em).values=psi
46   - (*syst_var.polfrac).values=p
47   - (*syst_var.polsed).values=largep
48   - CIQ=str_input_SED[indpolar].VARIANCEII & CIQ[*]=0. ;because input SEDs don't include covariances, yet
49   - CIU=str_input_SED[indpolar].VARIANCEII & CIU[*]=0. ;because input SEDs don't include covariances, yet
50   - CQU=str_input_SED[indpolar].VARIANCEII & CQU[*]=0. ;because input SEDs don't include covariances, yet
51   - polar_variance_iqu2ippsi,str_input_SED[indpolar].STOKESI $
52   - ,str_input_SED[indpolar].STOKESQ $
53   - ,str_input_SED[indpolar].STOKESU $
54   - ,str_input_SED[indpolar].VARIANCEII $
55   - ,str_input_SED[indpolar].VARIANCEQQ $
56   - ,str_input_SED[indpolar].VARIANCEUU $
57   - ,CIQ,CIU,CQU,variance_smallp,variance_psi,variance_largep,/lac,snr_largeP=snr_largeP
58   -
59   - (*syst_var.psi_em).instru_names=dustem_filter2instru(str_input_SED[indpolar].filter)
60   - (*syst_var.psi_em).filt_names=strtrim(str_input_SED[indpolar].filter,2)
61   - (*syst_var.psi_em).wav=str_input_SED[indpolar].wavelength
62   -
63   - (*syst_var.polfrac).instru_names = (*syst_var.psi_em).instru_names
64   - (*syst_var.polfrac).filt_names = (*syst_var.psi_em).filt_names
65   - (*syst_var.polfrac).wav = (*syst_var.psi_em).wav
66   -
67   - (*syst_var.polsed).instru_names = (*syst_var.psi_em).instru_names
68   - (*syst_var.polsed).filt_names = (*syst_var.psi_em).filt_names
69   - (*syst_var.polsed).wav = (*syst_var.psi_em).wav
  27 + IF test_m THEN BEGIN;IC: quickest way I found to make this work
  28 + ;===== note: those should include only points with polarization data, not all sed points
  29 + indq=where(str_input_SED.STOKESQ NE la_undef(),Nqsed)
  30 + indu=where(str_input_SED.STOKESU NE la_undef(),Nused)
  31 + qsed=dustem_initialize_internal_sed(Nqsed)
  32 + used=dustem_initialize_internal_sed(Nused)
  33 + indpolar=where(str_input_SED.STOKESU NE la_undef() AND str_input_SED.STOKESQ NE la_undef(),Npolar)
  34 +
  35 + psi_sed=dustem_initialize_internal_sed(Npolar)
  36 + smallp_sed=dustem_initialize_internal_sed(Npolar)
  37 + largep_sed=dustem_initialize_internal_sed(Npolar)
  38 +
  39 + syst_var.qsed=ptr_new(qsed)
  40 + syst_var.used=ptr_new(used)
  41 + (*syst_var.qsed).instru_names=dustem_filter2instru(str_input_SED[indq].filter)
  42 + (*syst_var.qsed).filt_names=strtrim(str_input_SED[indq].filter,2)
  43 + (*syst_var.qsed).wav=str_input_SED[indq].wavelength
  44 + (*syst_var.qsed).values=str_input_SED[indq].STOKESQ
  45 + (*syst_var.qsed).sigma=la_power(str_input_SED[indq].varianceQQ,0.5)
  46 + ;now U
  47 + (*syst_var.used).instru_names=dustem_filter2instru(str_input_SED[indu].filter)
  48 + (*syst_var.used).filt_names=strtrim(str_input_SED[indu].filter,2)
  49 + (*syst_var.used).wav=str_input_SED[indu].wavelength
  50 + (*syst_var.used).values=str_input_SED[indu].STOKESU
  51 + (*syst_var.used).sigma=la_power(str_input_SED[indu].varianceUU,0.5)
  52 + ;=== add up dependent fields
  53 + syst_var.psi_em=ptr_new(psi_sed)
  54 + syst_var.polfrac=ptr_new(smallp_sed)
  55 + syst_var.polsed=ptr_new(largep_sed)
  56 + polar_iqu2ippsi,str_input_SED[indpolar].STOKESI,str_input_SED[indpolar].STOKESQ,str_input_SED[indpolar].STOKESU,p,psi,largep=largep
  57 + (*syst_var.psi_em).values=psi
  58 + (*syst_var.polfrac).values=p
  59 + (*syst_var.polsed).values=largep
  60 + CIQ=str_input_SED[indpolar].VARIANCEII & CIQ[*]=0. ;because input SEDs don't include covariances, yet
  61 + CIU=str_input_SED[indpolar].VARIANCEII & CIU[*]=0. ;because input SEDs don't include covariances, yet
  62 + CQU=str_input_SED[indpolar].VARIANCEII & CQU[*]=0. ;because input SEDs don't include covariances, yet
  63 + polar_variance_iqu2ippsi,str_input_SED[indpolar].STOKESI $
  64 + ,str_input_SED[indpolar].STOKESQ $
  65 + ,str_input_SED[indpolar].STOKESU $
  66 + ,str_input_SED[indpolar].VARIANCEII $
  67 + ,str_input_SED[indpolar].VARIANCEQQ $
  68 + ,str_input_SED[indpolar].VARIANCEUU $
  69 + ,CIQ,CIU,CQU,variance_smallp,variance_psi,variance_largep,/lac,snr_largeP=snr_largeP
  70 +
  71 + (*syst_var.psi_em).instru_names=dustem_filter2instru(str_input_SED[indpolar].filter)
  72 + (*syst_var.psi_em).filt_names=strtrim(str_input_SED[indpolar].filter,2)
  73 + (*syst_var.psi_em).wav=str_input_SED[indpolar].wavelength
  74 +
  75 + (*syst_var.polfrac).instru_names = (*syst_var.psi_em).instru_names
  76 + (*syst_var.polfrac).filt_names = (*syst_var.psi_em).filt_names
  77 + (*syst_var.polfrac).wav = (*syst_var.psi_em).wav
  78 +
  79 + (*syst_var.polsed).instru_names = (*syst_var.psi_em).instru_names
  80 + (*syst_var.polsed).filt_names = (*syst_var.psi_em).filt_names
  81 + (*syst_var.polsed).wav = (*syst_var.psi_em).wav
  82 +
  83 + (*syst_var.psi_em).sigma=la_power(variance_psi,0.5)
  84 + (*syst_var.polfrac).sigma=la_power(variance_smallp,0.5)
  85 + (*syst_var.polsed).sigma=la_power(variance_largep,0.5)
  86 + ENDIF ELSE BEGIN
70 87  
71   - (*syst_var.psi_em).sigma=la_power(variance_psi,0.5)
72   - (*syst_var.polfrac).sigma=la_power(variance_smallp,0.5)
73   - (*syst_var.polsed).sigma=la_power(variance_largep,0.5)
  88 + ;===== note: those should include only points with polarization data, not all sed points
  89 + indq=where(str_input_SED.STOKESQ_EXT NE la_undef(),Nqsed)
  90 + indu=where(str_input_SED.STOKESU_EXT NE la_undef(),Nused)
  91 + qsed=dustem_initialize_internal_sed(Nqsed)
  92 + used=dustem_initialize_internal_sed(Nused)
  93 + indpolar=where(str_input_SED.STOKESU_EXT NE la_undef() AND str_input_SED.STOKESQ_EXT NE la_undef(),Npolar)
  94 +
  95 + psi_sed=dustem_initialize_internal_sed(Npolar)
  96 + smallp_sed=dustem_initialize_internal_sed(Npolar)
  97 + largep_sed=dustem_initialize_internal_sed(Npolar)
  98 +
  99 + syst_var.qext=ptr_new(qsed)
  100 + syst_var.uext=ptr_new(used)
  101 + (*syst_var.qext).instru_names=dustem_filter2instru(str_input_SED[indq].filter)
  102 + (*syst_var.qext).filt_names=strtrim(str_input_SED[indq].filter,2)
  103 + (*syst_var.qext).wav=str_input_SED[indq].wavelength
  104 + (*syst_var.qext).values=str_input_SED[indq].STOKESQ_EXT
  105 + (*syst_var.qext).sigma=la_power(str_input_SED[indq].varianceQQ_EXT,0.5)
  106 + ;now U
  107 + (*syst_var.uext).instru_names=dustem_filter2instru(str_input_SED[indu].filter)
  108 + (*syst_var.uext).filt_names=strtrim(str_input_SED[indu].filter,2)
  109 + (*syst_var.uext).wav=str_input_SED[indu].wavelength
  110 + (*syst_var.uext).values=str_input_SED[indu].STOKESU_EXT
  111 + (*syst_var.uext).sigma=la_power(str_input_SED[indu].varianceUU_EXT,0.5)
  112 + ;=== add up dependent fields
  113 + syst_var.psi_ext=ptr_new(psi_sed)
  114 + syst_var.fpolext=ptr_new(smallp_sed)
  115 + syst_var.polext=ptr_new(largep_sed)
  116 + polar_iqu2ippsi,str_input_SED[indpolar].STOKESI_EXT,str_input_SED[indpolar].STOKESQ_EXT,str_input_SED[indpolar].STOKESU_EXT,p,psi,largep=largep
  117 + (*syst_var.psi_ext).values=psi
  118 + (*syst_var.fpolext).values=p
  119 + (*syst_var.polext).values=largep
  120 + CIQ_EXT=str_input_SED[indpolar].VARIANCEII_EXT & CIQ_EXT[*]=0. ;because input SEDs don't include covariances, yet
  121 + CIU_EXT=str_input_SED[indpolar].VARIANCEII_EXT & CIU_EXT[*]=0. ;because input SEDs don't include covariances, yet
  122 + CQU_EXT=str_input_SED[indpolar].VARIANCEII_EXT & CQU_EXT[*]=0. ;because input SEDs don't include covariances, yet
  123 + polar_variance_iqu2ippsi,str_input_SED[indpolar].STOKESI_EXT $
  124 + ,str_input_SED[indpolar].STOKESQ_EXT $
  125 + ,str_input_SED[indpolar].STOKESU_EXT $
  126 + ,str_input_SED[indpolar].VARIANCEII_EXT $
  127 + ,str_input_SED[indpolar].VARIANCEQQ_EXT $
  128 + ,str_input_SED[indpolar].VARIANCEUU_EXT $
  129 + ,CIQ_EXT,CIU_EXT,CQU_EXT,variance_smallp_EXT,variance_psi_EXT,variance_largep_EXT,/lac,snr_largeP=snr_largeP_EXT
  130 +
  131 + (*syst_var.psi_ext).instru_names=dustem_filter2instru(str_input_SED[indpolar].filter)
  132 + (*syst_var.psi_ext).filt_names=strtrim(str_input_SED[indpolar].filter,2)
  133 + (*syst_var.psi_ext).wav=str_input_SED[indpolar].wavelength
  134 +
  135 + (*syst_var.fpolext).instru_names = (*syst_var.psi_ext).instru_names
  136 + (*syst_var.fpolext).filt_names = (*syst_var.psi_ext).filt_names
  137 + (*syst_var.fpolext).wav = (*syst_var.psi_ext).wav
  138 +
  139 + (*syst_var.polext).instru_names = (*syst_var.psi_ext).instru_names
  140 + (*syst_var.polext).filt_names = (*syst_var.psi_ext).filt_names
  141 + (*syst_var.polext).wav = (*syst_var.psi_ext).wav
  142 +
  143 + (*syst_var.psi_ext).sigma=la_power(variance_psi_EXT,0.5)
  144 + (*syst_var.fpolext).sigma=la_power(variance_smallp_EXT,0.5)
  145 + (*syst_var.polext).sigma=la_power(variance_largep_EXT,0.5)
  146 +
  147 + ENDELSE
  148 +
  149 +
74 150 ENDIF
75 151 END
76 152 \ No newline at end of file
... ...
src/idl/dustem_make_fits_predicted_ext.pro
  Diff comments   View file @3196a37
... ... @@ -51,16 +51,36 @@ ENDIF
51 51  
52 52  
53 53 ;===== define structures containing results
54   -one_str_input_EXT={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5),varianceII:la_undef(5),varianceQQ:la_undef(5),varianceUU:la_undef(5)}
55   -one_str_predicted_EXT={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5)}
  54 +one_str_input_EXT={FILTER:'',Wavelength:la_undef(4),StokesI_EXT:la_undef(5),StokesQ_EXT:la_undef(5),StokesU_EXT:la_undef(5),varianceII_EXT:la_undef(5),varianceQQ_EXT:la_undef(5),varianceUU_EXT:la_undef(5)}
  55 +one_str_predicted_EXT={FILTER:'',Wavelength:la_undef(4),StokesI_EXT:la_undef(5),StokesQ_EXT:la_undef(5),StokesU_EXT:la_undef(5)}
56 56  
57 57 Next=n_elements((*syst_var.ext).filt_names)
58 58 str_input_EXT=replicate(one_str_input_EXT,Next)
59 59 str_input_EXT.FILTER=(*syst_var.ext).filt_names
60 60 str_input_EXT.Wavelength=(*syst_var.ext).wav
61   -str_input_EXT.I=(*syst_var.ext).values
62   -str_input_EXT.varianceII=la_power((*syst_var.ext).sigma,2.)
63   -
  61 +str_input_EXT.StokesI_EXT=(*syst_var.ext).values
  62 +str_input_EXT.varianceII_EXT=la_power((*syst_var.ext).sigma,2.)
  63 +stop
  64 +;IC; if I didn't misunderstand the str_input_EXT output lacks this information
  65 +IF !run_pol THEN BEGIN
  66 + qext=aos2soa(*syst_var.qext)
  67 + FOR i=0L,n_elements(qext.wav)-1 DO BEGIN
  68 + ind=where(str_input_EXT.FILTER EQ qext[i].filt_names and str_input_EXT.wavelength EQ qext[i].wav,count)
  69 + IF count NE 0 THEN BEGIN
  70 + str_input_EXT[ind[0]].STOKESQ_EXT=qext[i].values
  71 + str_input_EXT[ind[0]].varianceQQ_EXT=la_power(qext[i].sigma,2.)
  72 + ENDIF
  73 + ENDFOR
  74 + uext=aos2soa(*syst_var.uext)
  75 + FOR i=0L,n_elements(uext.wav)-1 DO BEGIN
  76 + ind=where(str_input_EXT.FILTER EQ uext[i].filt_names and str_input_EXT.wavelength EQ uext[i].wav,count)
  77 + IF count NE 0 THEN BEGIN
  78 + str_input_EXT[ind[0]].STOKESU_EXT=uext[i].values
  79 + str_input_EXT[ind[0]].varianceUU_EXT=la_power(uext[i].sigma,2.)
  80 + ENDIF
  81 + ENDFOR
  82 +ENDIF
  83 +;IC: If we make an analogy with emission, why aren't there any predicted variances?
64 84 ;N_predicted_ext=n_elements(dustem_predicted_polext[0])
65 85 N_predicted_ext=n_elements(dustem_predicted_ext)
66 86 str_predicted_EXT=replicate(one_str_predicted_EXT,N_predicted_EXT)
... ... @@ -71,15 +91,15 @@ IF N_predicted_ext NE Next THEN BEGIN ;This is for cases when the extinction da
71 91 ENDIF ELSE BEGIN
72 92 str_predicted_EXT.FILTER=(*syst_var.ext).filt_names ;taken from input SED
73 93 str_predicted_EXT.Wavelength=(*syst_var.ext).wav ;taken from input SED
74   - str_predicted_EXT.I=dustem_predicted_ext
  94 + str_predicted_EXT.StokesI_EXT=dustem_predicted_ext
75 95 ENDELSE
76   -;stop
  96 +
77 97 IF !run_pol THEN BEGIN
78   - str_predicted_EXT.Q=interpol(dustem_predicted_Qext,(*syst_var.qext).wav,(*syst_var.ext).wav)
79   - str_predicted_EXT.U=interpol(dustem_predicted_Uext,(*syst_var.uext).wav,(*syst_var.ext).wav)
  98 + str_predicted_EXT.StokesQ_EXT=interpol(dustem_predicted_Qext,(*syst_var.qext).wav,(*syst_var.ext).wav)
  99 + str_predicted_EXT.StokesU_EXT=interpol(dustem_predicted_Uext,(*syst_var.uext).wav,(*syst_var.ext).wav)
80 100 ENDIF ELSE BEGIN
81   - str_predicted_EXT.Q=fully_undefined_ext.Q
82   - str_predicted_EXT.U=fully_undefined_ext.U
  101 + str_predicted_EXT.StokesQ_EXT=fully_undefined_ext.StokesQ_EXT
  102 + str_predicted_EXT.StokesU_EXT=fully_undefined_ext.StokesU_EXT
83 103 ENDELSE
84 104  
85 105 RETURN,str_predicted_EXT
... ...
src/idl/dustem_plot_dataset.pro
  Diff comments   View file @3196a37
... ... @@ -363,7 +363,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
363 363 if keyword_set(refresh) then begin ;The data points in the plot that are being refreshed
364 364  
365 365 ;making sure we're in the normalized plot so that overplotting happens here
366   - cgplot,1/vectw,vectx,/xlog,/nodata,/ys,xs=1,pos=position[1],noerase=1,xtickformat='(A1)',ytickformat='(A1)',color='Black',xr=xr,xtit='',yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,charsize=1.0
  366 + cgplot,1/vectw,vectx,/xlog,/nodata,ys=8,xs=8,pos=position[1],noerase=1,xtickformat='(A1)',ytickformat='(A1)',color='Black',xr=xr,xtit='',yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,charsize=1.0
367 367  
368 368 IF ct_spec NE 0 THEN BEGIN
369 369 xx=((*(*!dustem_data).ext).wav)[idx_spec]
... ... @@ -410,7 +410,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
410 410  
411 411  
412 412 ;making sure we're in the normal plot so that overplotting happens here
413   - cgplot,1/vectw,vectx,/xlog,/ylog,/nodata,/ys,xs=1,pos=position[0],noerase=1,xtickformat='(A1)',ytickformat='(A1)',color='Black',xr=xr,xtit='',yr=yr,yticks=2,ymino=2,xticklen=0.1,charsize=1.0
  413 + cgplot,1/vectw,vectx,/xlog,/ylog,/nodata,ys=8,xs=8,pos=position[0],noerase=1,xtickformat='(A1)',ytickformat='(A1)',color='Black',xr=xr,xtit='',yr=yr,yticks=2,ymino=2,xticklen=0.1,charsize=1.0
414 414  
415 415 ;Plotting of the spectra of the dust species
416 416  
... ... @@ -463,7 +463,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
463 463 if ct_spec ne 0 then begin
464 464  
465 465 ;Plotting of spectrum data points (to be fitted)
466   - cgplot,1/((*(*!dustem_data).ext).wav)(idx_spec),((*(*!dustem_data).ext).values)(idx_spec),/ylog,/xlog,/ys,xs=9,pos=position[0],noerase=1,xtit=' ',charsize=1.15,xtickformat='(A1)',color='Powder Blue',xr=xr,yr=yr,psym=16,syms=0.8,ytickformat='dstmwrp_exp'
  466 + cgplot,1/((*(*!dustem_data).ext).wav)(idx_spec),((*(*!dustem_data).ext).values)(idx_spec),/ylog,/xlog,ys=8,xs=8,pos=position[0],noerase=1,xtit=' ',charsize=1.15,xtickformat='(A1)',color='Powder Blue',xr=xr,yr=yr,psym=16,syms=0.8,ytickformat='dstmwrp_exp'
467 467  
468 468 ;Plotting of the spectrum error points
469 469 rms=2.*((*(*!dustem_data).ext).sigma)(idx_spec)/2.
... ... @@ -475,7 +475,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
475 475  
476 476 ;Plotting of filter data points (to be fitted) + testing for prior plotting
477 477 if ct_spec ne 0 then cgoplot,1/((*(*!dustem_data).ext).wav)(idx_filt),((*(*!dustem_data).ext).values)(idx_filt),charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,pos=position[0],/ys,xs=9,noerase=1,xtickformat='(A1)',xr=xr,yr=yr,ytickformat='dstmwrp_exp',/ylog,/xlog else $
478   - cgplot,1/((*(*!dustem_data).ext).wav)(idx_filt),((*(*!dustem_data).ext).values)(idx_filt),charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,pos=position[0],/ys,xs=9,noerase=1,xtickformat='(A1)',xtit=' ',xr=xr,yr=yr,ytickformat='dstmwrp_exp',/ylog,/xlog
  478 + cgplot,1/((*(*!dustem_data).ext).wav)(idx_filt),((*(*!dustem_data).ext).values)(idx_filt),charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,pos=position[0],ys=8,xs=8,noerase=1,xtickformat='(A1)',xtit=' ',xr=xr,yr=yr,ytickformat='dstmwrp_exp',/ylog,/xlog
479 479  
480 480  
481 481 ;Plotting of the filter error points
... ... @@ -502,7 +502,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
502 502 ;stop
503 503 if ct_spec_hdn ne 0 then begin
504 504 ;Plotting of hidden spectrum data points
505   - cgoplot,1/(((*(*!dustem_show).ext).wav)[idx_rmv_sed])(idx_spec_hdn),(((*(*!dustem_show).ext).values)[idx_rmv_sed])(idx_spec_hdn),pos=position[0],/ylog,/xlog,/ys,xs=9,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=yr,psym=16,syms=0.8
  505 + cgoplot,1/(((*(*!dustem_show).ext).wav)[idx_rmv_sed])(idx_spec_hdn),(((*(*!dustem_show).ext).values)[idx_rmv_sed])(idx_spec_hdn),pos=position[0],/ylog,/xlog,ys=8,xs=8,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=yr,psym=16,syms=0.8
506 506  
507 507 ;Plotting of hidden spectrum error points
508 508 rms=2.*(((*(*!dustem_show).ext).sigma)[idx_rmv_sed])(idx_spec_hdn)/2.
... ... @@ -511,7 +511,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
511 511  
512 512 if ct_filt_hdn then begin
513 513 ;Plotting of hidden filter data points
514   - cgoplot,1/(((*(*!dustem_show).ext).wav)[idx_rmv_sed])(idx_filt_hdn),(((*(*!dustem_show).ext).values)[idx_rmv_sed])(idx_filt_hdn),pos=position[0],/ylog,/xlog,/ys,xs=9,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=yr,psym=16,syms=0.8
  514 + cgoplot,1/(((*(*!dustem_show).ext).wav)[idx_rmv_sed])(idx_filt_hdn),(((*(*!dustem_show).ext).values)[idx_rmv_sed])(idx_filt_hdn),pos=position[0],/ylog,/xlog,ys=8,xs=8,noerase=1,charsize=1.15,xtickformat='(A1)',color='Black',xr=xr,yr=yr,psym=16,syms=0.8
515 515  
516 516 ;Plotting of hidden filter error point
517 517 rms=2.*(((*(*!dustem_show).ext).sigma)[idx_rmv_sed])(idx_filt_hdn)/2.
... ... @@ -1635,7 +1635,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
1635 1635 rms=2.*((*(*!dustem_data).qsed).sigma)(idx_filt)/2.
1636 1636  
1637 1637  
1638   - if ct_spec ne 0 then dustem_plot_mlog,((*(*!dustem_data).qsed).wav)(idx_filt),((*(*!dustem_data).qsed).values)(idx_filt),ppositions=position[0],charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,xs=1,xr=xr,noerase=1,/positive_only,xtickformat='(A1)',xtit='',/xlog,/ylog,yr=yr,ytickformat='(A1)';,ytickformat='(A1)'
  1638 + if ct_spec ne 0 then dustem_plot_mlog,((*(*!dustem_data).qsed).wav)(idx_filt),((*(*!dustem_data).qsed).values)(idx_filt),ppositions=position[0],charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,xs=1,xr=xr,noerase=1,/positive_only,xtickformat='(A1)',xtit='',/xlog,/ylog,yr=yr,/overplot,ytickformat='(A1)';,ytickformat='(A1)'
1639 1639 if ct_spec eq 0 then dustem_plot_mlog,((*(*!dustem_data).qsed).wav)(idx_filt),((*(*!dustem_data).qsed).values)(idx_filt),ppositions=position[0],charsize=1.15,color='Dodger Blue',psym=16,syms=0.8,xs=1,xr=xr,noerase=1,/positive_only,xtickformat='(A1)',xtit='',/xlog,/ylog,yr=yr;,ytickformat='(A1)'
1640 1640 dustem_plot_mlog,((*(*!dustem_data).qsed).wav)(idx_filt),((*(*!dustem_data).qsed).values)(idx_filt),ppositions=position[0],color='Dodger Blue',/positive_only,noerase=1,/overplot,rms=rms,xtit='',charsize=1.15,xtickformat='(A1)',xs=1,psym=16,syms=0.8,xr=xr,/xlog,/ylog,yr=yr,ytickformat='(A1)'
1641 1641  
... ... @@ -2218,7 +2218,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
2218 2218 if keyword_set(refresh) then begin ;The data points in the plot that are being refreshed
2219 2219  
2220 2220 ;Making sure we're placed in the right plot to refresh. Here normalized plot for QEXT.
2221   - cgplot,1/vectw,vectx,/nodata,/xlog,/ys,xs=1,pos=position[1],noerase=1,charsize=1.15,xtickformat='(A1)',ytickformat='(A1)',color='Powder Blue',xr=xr,xtit='',yr=[0.0,2.0],psym=8,syms=0.8
  2221 + cgplot,1/vectw,vectx,/nodata,/xlog,xs=1,pos=position[1],noerase=1,charsize=1.15,xtickformat='(A1)',ytickformat='(A1)',color='Powder Blue',xr=xr,xtit='',yr=[0.0,2.0],psym=8,syms=0.8
2222 2222  
2223 2223 IF ct_spec NE 0 THEN BEGIN
2224 2224  
... ... @@ -2255,7 +2255,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
2255 2255  
2256 2256 xtit=textoidl('1/\lambda (\mum^{-1})')
2257 2257  
2258   - cgplot,1/vectw,vectx,/xlog,/ys,xs=1,pos=position[1],noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0
  2258 + cgplot,1/vectw,vectx,/xlog,xs=1,pos=position[1],noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0
2259 2259 xyouts,pospltxt_n[0],pospltxt_n[1],textoidl('norm'),color=0,/normal,charsize=1.1
2260 2260  
2261 2261 endelse
... ... @@ -2574,7 +2574,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
2574 2574 if keyword_set(refresh) then begin ;The data points in the plot that are being refreshed
2575 2575  
2576 2576 ;Making sure we're placed in the right plot to refresh. Here normalized plot for QEXT.
2577   - cgplot,1/vectw,vectx,/nodata,/xlog,/ys,xs=1,pos=position[1],noerase=1,charsize=1.15,xtickformat='(A1)',ytickformat='(A1)',color='Powder Blue',xr=xr,xtit='',yr=[0.0,2.0],psym=8,syms=0.8
  2577 + cgplot,1/vectw,vectx,/nodata,/xlog,xs=1,pos=position[1],noerase=1,charsize=1.15,xtickformat='(A1)',ytickformat='(A1)',color='Powder Blue',xr=xr,xtit='',yr=[0.0,2.0],psym=8,syms=0.8
2578 2578  
2579 2579 IF ct_spec NE 0 THEN BEGIN
2580 2580  
... ... @@ -2607,7 +2607,7 @@ if keyword_set(dataset) then begin ;This shouldn't be a keyword. Should be chang
2607 2607  
2608 2608 xtit=textoidl('1/\lambda (\mum^{-1})')
2609 2609 ;if !run_pol then xtit = ''
2610   - cgplot,1/vectw,vectx,/xlog,/ys,xs=1,pos=position[1],noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0
  2610 + cgplot,1/vectw,vectx,/xlog,xs=1,pos=position[1],noerase=1,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0
2611 2611 xyouts,pospltxt_n[0],pospltxt_n[1],textoidl('norm'),color=0,/normal,charsize=1.1
2612 2612  
2613 2613 endelse
... ...
src/idl/dustem_read_fits_table.pro
  Diff comments   View file @3196a37
... ... @@ -224,7 +224,7 @@ IF unit_dustem_predicted_extinction_sed NE -1 THEN BEGIN
224 224 ENDIF ELSE BEGIN
225 225 ENDELSE
226 226 IF unit_dustem_shown_extinction_sed NE -1 THEN BEGIN
227   - str_shown_EXT=mrdfits(file,unit_shown_extinction_sed,header_shown_ext)
  227 + str_shown_EXT=mrdfits(file,unit_dustem_shown_extinction_sed,header_shown_ext)
228 228 ENDIF ELSE BEGIN
229 229 ENDELSE
230 230 ;=== read plugin informations and plugins headers
... ... @@ -252,7 +252,8 @@ ENDIF ELSE BEGIN
252 252 dustem_init,model=0
253 253 ENDELSE
254 254  
255   -defsysv,'!dustem_hcd',ptr_new(float(sxpar(header_input_sed,'NH')))
  255 +IF unit_input_emission_sed NE -1 THEN defsysv,'!dustem_hcd',ptr_new(float(sxpar(header_input_sed,'NH')))
  256 +IF unit_input_extinction_sed NE -1 THEN defsysv,'!dustem_hcd',ptr_new(float(sxpar(header_input_ext,'NH')))
256 257  
257 258 !dustem_model=used_model
258 259 IF !dustem_version.version NE used_version THEN BEGIN
... ... @@ -308,7 +309,6 @@ IF Nplugin NE 0 THEN BEGIN
308 309 (*!dustem_plugin).(i).PARAMTAG=ptr_new(paramtags)
309 310 ENDFOR
310 311 ENDIF
311   -;stop
312 312  
313 313 ;==== This is to create (*!dustem_data).sed in case it does not exist already
314 314 ; it doesn't exist at all times because dustem_init has been run
... ... @@ -319,7 +319,6 @@ IF unit_input_emission_sed NE -1 THEN BEGIN
319 319 ENDIF
320 320  
321 321 ;help,str_shown_SED
322   -;stop
323 322  
324 323 ;==== This is to create (*!dustem_show).sed in case it does not exist already
325 324 IF unit_dustem_shown_sed NE -1 THEN BEGIN
... ...
src/idl/dustem_write_fits_table.pro
  Diff comments   View file @3196a37
... ... @@ -187,7 +187,7 @@ IF Nplugins NE 0 THEN BEGIN
187 187 str_plugins[i]=plugin_tags[i]
188 188 unit=unit+1
189 189 ENDFOR
190   -ENDIF ELSE BEGIN
  190 +ENDIF ELSE BEGIN ;IC: This should be commented right? These quantities do not exist.
191 191 unit_plugins[i]=-1
192 192 str_plugins[i]=''
193 193 ENDELSE
... ... @@ -290,29 +290,51 @@ FOR i=0L,Nparams-1 DO BEGIN
290 290 ENDFOR
291 291 sxaddpar,hdr,'COMMENT','===================================================',after='PARF'+strtrim(i,2)
292 292  
293   -;stop
294 293  
295   -;===== complement extension header for input SED
296   -sxaddpar,header_input_sed,'TITLE','OBSERVED INPUT SED TO DUSTEMWRAP',before='COMMENT'
297   -sxaddpar,header_input_sed,'TTYPE1','FILTER','Dustemwrap Filter name'
298   -sxaddpar,header_input_sed,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
299   -sxaddpar,header_input_sed,'TTYPE3','I','Observed Total Stokes I Intensity [??]'
300   -sxaddpar,header_input_sed,'TTYPE4','Q','Observed Stokes Q Intensity [??]'
301   -sxaddpar,header_input_sed,'TTYPE5','U','Observed Stokes U Intensity [??]'
302   -sxaddpar,header_input_sed,'TTYPE6','VARIANCEII','Observed Total Stokes I Variance [??]'
303   -sxaddpar,header_input_sed,'TTYPE7','VARIANCEQQ','Observed Stokes Q Variance [??]'
304   -sxaddpar,header_input_sed,'TTYPE8','VARIANCEUU','Observed Stokes U Variance [??]'
305   -sxaddpar,header_input_sed,'NH',*!dustem_hcd,'Column density [H/cm2]'
  294 +IF unit_input_sed NE -1 THEN BEGIN
306 295  
307   -;stop
  296 + ;===== complement extension header for input SED
  297 + sxaddpar,header_input_sed,'TITLE','OBSERVED INPUT SED TO DUSTEMWRAP',before='COMMENT'
  298 + sxaddpar,header_input_sed,'TTYPE1','FILTER','Dustemwrap Filter name'
  299 + sxaddpar,header_input_sed,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
  300 + sxaddpar,header_input_sed,'TTYPE3','STOKESI','Observed Total Stokes I Intensity [??]'
  301 + sxaddpar,header_input_sed,'TTYPE4','STOKESQ','Observed Stokes Q Intensity [??]'
  302 + sxaddpar,header_input_sed,'TTYPE5','STOKESU','Observed Stokes U Intensity [??]'
  303 + sxaddpar,header_input_sed,'TTYPE6','VARIANCEII','Observed Total Stokes I Variance [??]'
  304 + sxaddpar,header_input_sed,'TTYPE7','VARIANCEQQ','Observed Stokes Q Variance [??]'
  305 + sxaddpar,header_input_sed,'TTYPE8','VARIANCEUU','Observed Stokes U Variance [??]'
  306 + sxaddpar,header_input_sed,'NH',*!dustem_hcd,'Column density [H/cm2]'
  307 +
  308 + ;===== complement extension header for predicted SED
  309 + sxaddpar,header_predicted_SED,'TITLE','OUTPUT BEST FIT SED BY DUSTEMWRAP'
  310 + sxaddpar,header_predicted_SED,'TTYPE1','FILTER','Dustemwrap Filter name'
  311 + sxaddpar,header_predicted_SED,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
  312 + sxaddpar,header_predicted_SED,'TTYPE3','STOKESI','Observed Total Stokes I Intensity [??]'
  313 + sxaddpar,header_predicted_SED,'TTYPE4','STOKESQ','Observed Stokes Q Intensity [??]'
  314 + sxaddpar,header_predicted_SED,'TTYPE5','STOKESU','Observed Stokes U Intensity [??]'
  315 +ENDIF ELSE BEGIN;IC: hesitated between this and 'IF unit_input_ext NE -1 THEN BEGIN'. This is under the assumtion that we always fit something.
  316 +
  317 + ;===== complement extension header for input extinction SED
  318 + sxaddpar,header_input_ext,'TITLE','OBSERVED INPUT EXTINCTION SED TO DUSTEMWRAP',before='COMMENT'
  319 + sxaddpar,header_input_ext,'TTYPE1','FILTER','Dustemwrap Filter name'
  320 + sxaddpar,header_input_ext,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
  321 + sxaddpar,header_input_ext,'TTYPE3','STOKESI_EXT','Observed Total Stokes I optical depth []'
  322 + sxaddpar,header_input_ext,'TTYPE4','STOKESQ_EXT','Observed Stokes Q optical depth []'
  323 + sxaddpar,header_input_ext,'TTYPE5','STOKESU_EXT','Observed Stokes U optical depth []'
  324 + sxaddpar,header_input_ext,'TTYPE6','VARIANCEII_EXT','Observed Total Stokes I optical depth Variance [??]'
  325 + sxaddpar,header_input_ext,'TTYPE7','VARIANCEQQ_EXT','Observed Stokes Q optical depth Variance [??]'
  326 + sxaddpar,header_input_ext,'TTYPE8','VARIANCEUU_EXT','Observed Stokes U optical depth Variance [??]'
  327 + sxaddpar,header_input_ext,'NH',*!dustem_hcd,'Column density [H/cm2]'
  328 +
  329 + ;===== complement extension header for predicted extinction SED
  330 + sxaddpar,header_predicted_EXT,'TITLE','OUTPUT BEST FIT EXTINCTION SED BY DUSTEMWRAP'
  331 + sxaddpar,header_predicted_EXT,'TTYPE1','FILTER','Dustemwrap Filter name'
  332 + sxaddpar,header_predicted_EXT,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
  333 + sxaddpar,header_predicted_EXT,'TTYPE3','STOKESI_EXT','Observed Total Stokes I optical depth []'
  334 + sxaddpar,header_predicted_EXT,'TTYPE4','STOKESQ_EXT','Observed Stokes Q optical depth []'
  335 + sxaddpar,header_predicted_EXT,'TTYPE5','STOKESU_EXT','Observed Stokes U optical depth []'
308 336  
309   -;===== complement extension header for predicted SED
310   -sxaddpar,header_predicted_SED,'TITLE','OUTPUT BEST FIT SED BY DUSTEMWRAP'
311   -sxaddpar,header_predicted_SED,'TTYPE1','FILTER','Dustemwrap Filter name'
312   -sxaddpar,header_predicted_SED,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
313   -sxaddpar,header_predicted_SED,'TTYPE3','I','Observed Total Stokes I Intensity [??]'
314   -sxaddpar,header_predicted_SED,'TTYPE4','Q','Observed Stokes Q Intensity [??]'
315   -sxaddpar,header_predicted_SED,'TTYPE5','U','Observed Stokes U Intensity [??]'
  337 +ENDELSE
316 338  
317 339 ;===== complement extension header for predicted spectrum total
318 340 sxaddpar,header_predicted_spectra,'TITLE','OUTPUT BEST FIT PER-GRAIN AND TOTAL EMISSION SPECTRA BY DUSTEMWRAP'
... ...