[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 |
---|
| 8 | |
---|
| 9 | def _import_data(data): |
---|
| 10 | if not hasattr(data, "__len__"): |
---|
| 11 | if isinstance(data, scantable): |
---|
| 12 | return data |
---|
| 13 | elif isinstance(data, str): |
---|
| 14 | return scantable(data) |
---|
| 15 | tables = [] |
---|
| 16 | for d in data: |
---|
| 17 | if isinstance(d, scantable): |
---|
| 18 | tables.append(d) |
---|
| 19 | elif isinstance(d, str): |
---|
| 20 | if os.path.exists(d): |
---|
| 21 | tables.append(scantable(d)) |
---|
| 22 | else: |
---|
| 23 | raise IOError("Data file doesn't exists") |
---|
| 24 | else: |
---|
| 25 | raise TypeError("data is not a scantable or valid file") |
---|
| 26 | return merge(tables) |
---|
| 27 | |
---|
| 28 | def skydip(data, averagepol=True, tsky=300., plot=False): |
---|
| 29 | """Determine the opacity from a set of 'skydip' obervations. |
---|
| 30 | This can be any set of observations over a range of elevations, |
---|
| 31 | but will ususally be a dedicated (set of) scan(s). |
---|
| 32 | Return a list of 'n' opacities for 'n' IFs. In case of averagepol |
---|
| 33 | being 'False' a list of 'n*m' elements where 'm' is the number of |
---|
| 34 | polarisations, e.g. |
---|
| 35 | nIF = 3, nPol = 2 => [if0pol0, if0pol1, if1pol0, if1pol1, if2pol0, if2pol1] |
---|
| 36 | |
---|
| 37 | The opacity is determined by fitting a first order polynomial to: |
---|
| 38 | |
---|
| 39 | |
---|
| 40 | Tsys(airmass) = p0 + airmass*p1 |
---|
| 41 | |
---|
| 42 | where |
---|
| 43 | |
---|
| 44 | airmass = 1/sin(elevation) |
---|
| 45 | |
---|
| 46 | tau = p1/Tsky |
---|
| 47 | |
---|
| 48 | Parameters: |
---|
| 49 | data: a list of file names or scantables or a single |
---|
| 50 | file name or scantable. |
---|
| 51 | averagepol: Return the average of the opacities for the polarisations |
---|
| 52 | This might be useful to set to 'False' if one polarisation |
---|
| 53 | is corrupted (Mopra). If set to 'False', an opacity value |
---|
| 54 | per polarisation is returned. |
---|
| 55 | tksy: The sky temperature (default 300.0K). This might |
---|
| 56 | be read from the data in the future. |
---|
| 57 | plot: Plot each fit (airmass vs. Tsys). Default is 'False' |
---|
| 58 | """ |
---|
| 59 | scan = _import_data(data) |
---|
| 60 | f = fitter() |
---|
| 61 | f.set_function(poly=1) |
---|
| 62 | sel = selector() |
---|
| 63 | basesel = scan.get_selection() |
---|
| 64 | inos = scan.getifnos() |
---|
| 65 | pnos = scan.getpolnos() |
---|
| 66 | opacities = [] |
---|
| 67 | for ino in inos: |
---|
| 68 | sel.set_ifs(ino) |
---|
| 69 | opacity = [] |
---|
| 70 | for pno in pnos: |
---|
| 71 | sel.set_polarisations(pno) |
---|
| 72 | scan.set_selection(basesel+sel) |
---|
| 73 | tsys = scan.get_tsys() |
---|
| 74 | elev = scan.get_elevation() |
---|
| 75 | airmass = [ 1./math.sin(i) for i in elev ] |
---|
| 76 | f.set_data(airmass, tsys) |
---|
| 77 | f.fit() |
---|
| 78 | if plot: |
---|
| 79 | f.plot(residual=True, plotparms=True) |
---|
| 80 | raw_input("Hit <return> for next fit...") |
---|
| 81 | params = f.get_parameters()["params"] |
---|
| 82 | opacity.append(params[1]/tsky) |
---|
| 83 | if averagepol: |
---|
| 84 | opacities.append(sum(opacity)/len(opacity)) |
---|
| 85 | else: |
---|
| 86 | opacities += opacity |
---|
| 87 | sel.reset() |
---|
| 88 | scan.set_selection(basesel) |
---|
| 89 | return opacities |
---|