spectralindex

Return to directory of Todd's CASA extensions

This function is designed to fit a power-law spectral index to the results output from casa's fluxscale task, OR from a simple list of frequencies, flux densities and uncertainties. In the former use case, currently it also requires the output from the listobs task to determine the center frequencies of each spectral window. It has two options for determining the uncertainty of the fitted parameters: using the covariance matrix (obtained by setting trials=0) or running a Monte Carlo simulation (default, with trials=2000). The latter appears to be more accurate in most cases. It produces a plot summarizing the results. Finally, if you have only one flux density measurement and a spectral index, you can use this function to ask what the flux density would be at a different frequency.

Usage:

au.spectralindex(self, filename='', yfilename='', source='', verbose=False, maxpoints=0, 
                       trials=2000, spw='', plotdir='', labelspw=False, referenceFrame='TOPO', 
                       plaintext=False, lineNumbers=None,columns=None, yscale=1.0,
                       plotunits='linear', freqs=[], fluxes=[], errors=[], plotfile='', 
                       showplot=True, xaxis=[], spectralIndex=None, silent=False,
                       axisEqual=True) 
  • filename: contains a listobs output file
  • yfilename: contains a fluxscale output file
  • source: sourcename to choose from the (possibly) multi-source fluxscale file
  • maxpoints: the maximum number of spws to select for the fit (0=no max)
  • trials: if > 0, use a Monte-Carlo technique estimate the fit uncertainties, otherwise, use the sqrt(covar_matrix) from scipy.optimize.leastsq (default) There is a minimum number of 100 trials, and ~1000 is recommended.
  • spw: the spws to use, e.g. ''=all, '1~3,5,6~8'=[1,2,3,5,6,7,8]
  • plotdir: the directory in which to write the plotfile
  • labelspw: draw spw numeric labels adjacent to each point
  • referenceFrame: the frequency reference frame, if not 'TOPO'
  • plaintext: if True, then read the freqs, flux densities and uncertainties from a single plain text file, specified by filename
  • lineNumbers: which lines to read, where 1 is the first line in the file
  • columns: which columns to read for freq, flux, error (starting at 1)
  • yscale: factor to scale the y values after reading from the file
  • plotunits: 'linear' or 'log' (only used if trials==0)
  • freqs,fluxes,errors: alternative approach to using filename, just pass 3 lists
  • xaxis: a list of points for which to compute model y values
  • spectralindex: If you have only one flux measurement, you can specify it along with a spectral index (via this parameter) and it will compute expected flux densities at the frequencies specified by the xaxis parameter.
  • axisEqual: if True, then when plotting, set pl.axis('equal')

Input file formats

The listobs file should look like this. It is keying off the string TOPO, so there can be any number of leading comment characters on each line, and any number of unrelated lines of text at the top. It can also have an additional column (the spw name) between the spw and the #Chans, as this is the new format of listobs output.
 SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)  TotBW(kHz)  Corrs
  0         128 TOPO  42738       1000          128000      RR  LL
  1         128 TOPO  42866       1000          128000      RR  LL
  2         128 TOPO  42994       1000          128000      RR  LL
  3         128 TOPO  43122       1000          128000      RR  LL
  4         128 TOPO  43250       1000          128000      RR  LL
  5         128 TOPO  43378       1000          128000      RR  LL
  6         128 TOPO  43506       1000          128000      RR  LL
  7         128 TOPO  43634       1000          128000      RR  LL
  8         128 TOPO  44488       1000          128000      RR  LL
  9         128 TOPO  44616       1000          128000      RR  LL
  10        128 TOPO  44744       1000          128000      RR  LL
  11        128 TOPO  44872       1000          128000      RR  LL
  12        128 TOPO  45000       1000          128000      RR  LL
  13        128 TOPO  45128       1000          128000      RR  LL
  14        128 TOPO  45256       1000          128000      RR  LL
  15        128 TOPO  45384       1000          128000      RR  LL

Similarly, the fluxscale file should look like this. It can contain the results for multiple sources. It is keying off of the string 'SpW', so there can be any number of leading comment characters on each line, and any number of unrelated lines of text at the top.

# Opening MS: 10C186_B_Q_d1_3s.ms for calibration.
# Found reference field(s): J1331+3030
# Found transfer field(s):  J1717-3342 J1924-2914
# Flux density for J1717-3342 in SpW=0 is: 2.33891 ± 0.0620232 (SNR = 37.7102, N= 44)
# Flux density for J1717-3342 in SpW=1 is: 2.33168 ± 0.0608224 (SNR = 38.3358, N= 44)
# Flux density for J1717-3342 in SpW=2 is: 2.33139 ± 0.0642111 (SNR = 36.3082, N= 44)
# Flux density for J1717-3342 in SpW=3 is: 2.33639 ± 0.0697579 (SNR = 33.4929, N= 44)
# Flux density for J1717-3342 in SpW=4 is: 2.31025 ± 0.0631668 (SNR = 36.5738, N= 44)
# Flux density for J1717-3342 in SpW=5 is: 2.32929 ± 0.0667446 (SNR = 34.8986, N= 44)
# Flux density for J1717-3342 in SpW=6 is: 2.35698 ± 0.0698231 (SNR = 33.7564, N= 44)
# Flux density for J1717-3342 in SpW=7 is: 2.32728 ± 0.0668059 (SNR = 34.8364, N= 44)
# Flux density for J1717-3342 in SpW=8 is: 2.23395 ± 0.0394774 (SNR = 56.588, N= 46)
# Flux density for J1717-3342 in SpW=9 is: 2.23898 ± 0.0404523 (SNR = 55.3487, N= 46)
# Flux density for J1717-3342 in SpW=10 is: 2.20212 ± 0.0408321 (SNR = 53.9311, N= 46)
# Flux density for J1717-3342 in SpW=11 is: 2.18275 ± 0.0392774 (SNR = 55.5727, N= 46)
# Flux density for J1717-3342 in SpW=12 is: 2.17859 ± 0.0388576 (SNR = 56.0659, N= 46)
# Flux density for J1717-3342 in SpW=13 is: 2.16182 ± 0.0405544 (SNR = 53.3066, N= 46)
# Flux density for J1717-3342 in SpW=14 is: 2.16303 ± 0.0413542 (SNR = 52.305, N= 46)
# Flux density for J1717-3342 in SpW=15 is: 2.19681 ± 0.0400311 (SNR = 54.8777, N= 46)
# Flux density for J1924-2914 in SpW=0 is: 14.6731 ± 0.138529 (SNR = 105.921, N= 44)
# Flux density for J1924-2914 in SpW=1 is: 14.5996 ± 0.136761 (SNR = 106.753, N= 44)
# Flux density for J1924-2914 in SpW=2 is: 14.6242 ± 0.146673 (SNR = 99.7061, N= 44)
# Flux density for J1924-2914 in SpW=3 is: 14.6802 ± 0.158394 (SNR = 92.6816, N= 44)
# Flux density for J1924-2914 in SpW=4 is: 14.5259 ± 0.141069 (SNR = 102.97, N= 44)
# Flux density for J1924-2914 in SpW=5 is: 14.6075 ± 0.151812 (SNR = 96.2206, N= 44)
# Flux density for J1924-2914 in SpW=6 is: 14.7261 ± 0.157892 (SNR = 93.2666, N= 44)
# Flux density for J1924-2914 in SpW=7 is: 14.5783 ± 0.14938 (SNR = 97.5919, N= 44)
# Flux density for J1924-2914 in SpW=8 is: 13.9264 ± 0.0748849 (SNR = 185.971, N= 46)
# Flux density for J1924-2914 in SpW=9 is: 13.9368 ± 0.0751502 (SNR = 185.452, N= 46)
# Flux density for J1924-2914 in SpW=10 is: 13.7826 ± 0.0785357 (SNR = 175.494, N= 46)
# Flux density for J1924-2914 in SpW=11 is: 13.7178 ± 0.0778214 (SNR = 176.273, N= 46)
# Flux density for J1924-2914 in SpW=12 is: 13.6833 ± 0.0748109 (SNR = 182.905, N= 46)
# Flux density for J1924-2914 in SpW=13 is: 13.6308 ± 0.0772038 (SNR = 176.555, N= 46)
# Flux density for J1924-2914 in SpW=14 is: 13.6372 ± 0.0792322 (SNR = 172.117, N= 46)
# Flux density for J1924-2914 in SpW=15 is: 13.7463 ± 0.0735882 (SNR = 186.8, N= 46)

Examples:

1. Using a single plain text file:
CASA <4>: au.spectralindex('plain.txt', plaintext=True)
Error-weighted fit: Slope: 2.250+-0.109  Flux D. @ 115.000GHz: 0.233+-0.005
   Un-weighted fit: Slope: 2.234         Flux D. @ 115.000GHz: 0.239
Plot saved in ./plain.txt.png

2. Using a listobs output file and a fluxscale output file:
CASA <5>: au.spectralIndex('listobs.txt','fluxscale_d2_new.txt','J17')
Read 16 x values for J1717-3342 (avg=44125.000)
Read 16 y values for J1717-3342 (avg=2.264)
Error-weighted fit: Slope: -1.454+-0.640  Flux D. @ 42.802GHz: 2.362+-0.014
   Un-weighted fit: Slope: -1.453         Flux D. @ 42.802GHz: 2.364
Plot saved in ./fluxscale_d2_new.txt.J1717-3342.png
fluxscale d2 new.txt.J1717-3342.png

3. Using a list of flux density measurements, and finding the value at another frequency:
CASA <9>: au.spectralindex(freqs=[22.3,43.3,222], fluxes=[.0003,.001,.044], errors=[.00001,.0001,.001], xaxis=[115])
Value at 22.300000 = 0.000295 +- 0.000022 (scaled MAD = 0.000022)
Error-weighted fit: Slope: 2.177+-0.222  Flux D. @ 22.300GHz: 0.295+-0.022 mJy
   Un-weighted fit: Slope: 2.198         Flux D. @ 22.300GHz: 0.270 mJy
Value at 115.000000 = 0.010475 +- 0.000451 (scaled MAD = 0.000448)
  Out[9]:
{'covar': array([[ 0.00043185, -0.00038152],
       [-0.00038152,  0.0015951 ]]),
 'intercept': -2.5972443931713598,
 'interceptUncertainty': 0.50968131460806199,
 'meanOfLogX': 1.7770485779507217,
 'spectralIndex': 2.1767855043259088,
 'spectralIndexUncertainty': 0.22226760634713091,
 'yaxis': [0.010476089728640246],
 'yaxisUncertainty': [0.00045094643771614152]}

-- ToddHunter - 2012-03-07
Topic attachments
I Attachment Action Size Date Who Comment
fluxscale_d2_new.txttxt fluxscale_d2_new.txt manage 2 K 2013-08-08 - 08:39 ToddHunter  
fluxscale_d2_new.txt.J1717-3342.pngpng fluxscale_d2_new.txt.J1717-3342.png manage 40 K 2013-08-08 - 12:03 ToddHunter  
fluxscale_d2_new.txt.J1924-2914.pngpng fluxscale_d2_new.txt.J1924-2914.png manage 42 K 2012-03-07 - 10:37 ToddHunter  
listobs.txttxt listobs.txt manage 1 K 2013-08-08 - 08:13 ToddHunter  
plain.txttxt plain.txt manage 58 bytes 2013-06-22 - 16:24 ToddHunter example file for the plaintext=True option
Topic revision: r10 - 2015-05-14, ToddHunter
This site is powered by FoswikiCopyright © by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding NRAO Public Wiki? Send feedback