| [1689] | 1 | import os
|
|---|
| 2 | import math
|
|---|
| 3 | from asap import scantable
|
|---|
| 4 | from asap import merge
|
|---|
| 5 | from asap import fitter
|
|---|
| 6 | from asap import selector
|
|---|
| 7 | from asap import rcParams
|
|---|
| [1722] | 8 | from asap import xyplotter
|
|---|
| [1689] | 9 |
|
|---|
| 10 | def _import_data(data):
|
|---|
| [1722] | 11 | if not isinstance(data, (list,tuple)):
|
|---|
| [1689] | 12 | if isinstance(data, scantable):
|
|---|
| 13 | return data
|
|---|
| 14 | elif isinstance(data, str):
|
|---|
| 15 | return scantable(data)
|
|---|
| 16 | tables = []
|
|---|
| 17 | for d in data:
|
|---|
| 18 | if isinstance(d, scantable):
|
|---|
| 19 | tables.append(d)
|
|---|
| 20 | elif isinstance(d, str):
|
|---|
| 21 | if os.path.exists(d):
|
|---|
| 22 | tables.append(scantable(d))
|
|---|
| 23 | else:
|
|---|
| 24 | raise IOError("Data file doesn't exists")
|
|---|
| 25 | else:
|
|---|
| 26 | raise TypeError("data is not a scantable or valid file")
|
|---|
| 27 | return merge(tables)
|
|---|
| 28 |
|
|---|
| 29 | def skydip(data, averagepol=True, tsky=300., plot=False):
|
|---|
| 30 | """Determine the opacity from a set of 'skydip' obervations.
|
|---|
| 31 | This can be any set of observations over a range of elevations,
|
|---|
| 32 | but will ususally be a dedicated (set of) scan(s).
|
|---|
| 33 | Return a list of 'n' opacities for 'n' IFs. In case of averagepol
|
|---|
| 34 | being 'False' a list of 'n*m' elements where 'm' is the number of
|
|---|
| 35 | polarisations, e.g.
|
|---|
| 36 | nIF = 3, nPol = 2 => [if0pol0, if0pol1, if1pol0, if1pol1, if2pol0, if2pol1]
|
|---|
| 37 |
|
|---|
| 38 | The opacity is determined by fitting a first order polynomial to:
|
|---|
| 39 |
|
|---|
| 40 |
|
|---|
| 41 | Tsys(airmass) = p0 + airmass*p1
|
|---|
| 42 |
|
|---|
| 43 | where
|
|---|
| 44 |
|
|---|
| 45 | airmass = 1/sin(elevation)
|
|---|
| 46 |
|
|---|
| 47 | tau = p1/Tsky
|
|---|
| 48 |
|
|---|
| 49 | Parameters:
|
|---|
| 50 | data: a list of file names or scantables or a single
|
|---|
| 51 | file name or scantable.
|
|---|
| 52 | averagepol: Return the average of the opacities for the polarisations
|
|---|
| 53 | This might be useful to set to 'False' if one polarisation
|
|---|
| 54 | is corrupted (Mopra). If set to 'False', an opacity value
|
|---|
| 55 | per polarisation is returned.
|
|---|
| 56 | tksy: The sky temperature (default 300.0K). This might
|
|---|
| 57 | be read from the data in the future.
|
|---|
| 58 | plot: Plot each fit (airmass vs. Tsys). Default is 'False'
|
|---|
| 59 | """
|
|---|
| [1722] | 60 | rcsave = rcParams['verbose']
|
|---|
| 61 | rcParams['verbose'] = False
|
|---|
| [1689] | 62 | scan = _import_data(data)
|
|---|
| 63 | f = fitter()
|
|---|
| 64 | f.set_function(poly=1)
|
|---|
| 65 | sel = selector()
|
|---|
| 66 | basesel = scan.get_selection()
|
|---|
| 67 | inos = scan.getifnos()
|
|---|
| 68 | pnos = scan.getpolnos()
|
|---|
| 69 | opacities = []
|
|---|
| 70 | for ino in inos:
|
|---|
| 71 | sel.set_ifs(ino)
|
|---|
| 72 | opacity = []
|
|---|
| [1722] | 73 | fits = []
|
|---|
| 74 | airms = []
|
|---|
| 75 | tsyss = []
|
|---|
| 76 |
|
|---|
| 77 | if plot:
|
|---|
| 78 | xyplotter.cla()
|
|---|
| 79 | xyplotter.ioff()
|
|---|
| 80 | xyplotter.clf()
|
|---|
| 81 | xyplotter.xlabel("Airmass")
|
|---|
| 82 | xyplotter.ylabel(r"$T_{sys}$")
|
|---|
| [1689] | 83 | for pno in pnos:
|
|---|
| 84 | sel.set_polarisations(pno)
|
|---|
| 85 | scan.set_selection(basesel+sel)
|
|---|
| [1722] | 86 | freq = scan.get_coordinate(0).get_reference_value()/1e9
|
|---|
| 87 | freqstr = "%0.4f GHz" % freq
|
|---|
| [1689] | 88 | tsys = scan.get_tsys()
|
|---|
| 89 | elev = scan.get_elevation()
|
|---|
| 90 | airmass = [ 1./math.sin(i) for i in elev ]
|
|---|
| [1722] | 91 | airms.append(airmass)
|
|---|
| 92 | tsyss.append(tsys)
|
|---|
| [1689] | 93 | f.set_data(airmass, tsys)
|
|---|
| 94 | f.fit()
|
|---|
| [1722] | 95 | fits.append(f.get_fit())
|
|---|
| [1689] | 96 | params = f.get_parameters()["params"]
|
|---|
| 97 | opacity.append(params[1]/tsky)
|
|---|
| 98 | if averagepol:
|
|---|
| 99 | opacities.append(sum(opacity)/len(opacity))
|
|---|
| 100 | else:
|
|---|
| 101 | opacities += opacity
|
|---|
| [1722] | 102 | if plot:
|
|---|
| 103 | colors = ['b','g','k']
|
|---|
| 104 | for i in range(len(airms)):
|
|---|
| 105 | xyplotter.plot(airms[i], tsyss[i], 'o', color=colors[i])
|
|---|
| 106 | xyplotter.plot(airms[i], fits[i], '-', color=colors[i])
|
|---|
| 107 | xyplotter.figtext(0.7,0.3-(i/30.0),
|
|---|
| 108 | r"$\tau_{fit}=%0.2f$" % opacity[i],
|
|---|
| 109 | color=colors[i])
|
|---|
| 110 | if averagepol:
|
|---|
| 111 | xyplotter.figtext(0.7,0.3-(len(airms)/30.0),
|
|---|
| 112 | r"$\tau=%0.2f$" % opacities[-1],
|
|---|
| 113 | color='r')
|
|---|
| 114 | xyplotter.title("IF%d : %s" % (ino, freqstr))
|
|---|
| 115 |
|
|---|
| 116 | xyplotter.ion()
|
|---|
| 117 | xyplotter.draw()
|
|---|
| 118 | raw_input("Hit <return> for next fit...")
|
|---|
| [1689] | 119 | sel.reset()
|
|---|
| [1722] | 120 |
|
|---|
| [1689] | 121 | scan.set_selection(basesel)
|
|---|
| [1722] | 122 | rcParams['verbose'] = rcsave
|
|---|
| 123 | if plot:
|
|---|
| 124 | xyplotter.close()
|
|---|
| [1689] | 125 | return opacities
|
|---|
