Changeset 1725

Timestamp:

04/27/10 15:48:24 (15 years ago)

Author:

Malte Marquarding

Message:

Finishing touches to opacity calculations, docs, plotting and model

Location:

trunk

Files:

• doc/CHANGELOG (modified) (2 diffs)
• python/__init__.py (modified) (3 diffs)
• python/opacity.py (modified) (7 diffs)
• python/scantable.py (modified) (2 diffs)
• src/STFiller.cpp (modified) (1 diff)
• src/python_STAtmosphere.cpp (modified) (1 diff)
• test/test_scantable.py (modified) (3 diffs)

trunk/doc/CHANGELOG

@@ -2 +2 @@
 =========
 
-Release 3.0.0 []
+Release 3.0.0 [2010-04-??]
 
+* Ticket #177 Added function skydip to determine opacities.
+* Ticket #172 Fixed non-working scantable.resample
 * Ticket #155 Better output filenames. Ignore non-existsing beams/pols/ifs/scans
 * Ticket #157 numpy >= 1.1 support
 * Ticket #158 fixed plotter.set_font
-* Ticket #160 Aspect ration of pplotter customisable
-* Ticket #162 - fill in
+* Ticket #160 Aspect ratio of plotter is customisable now
 * Ticket #163 fixed for scantable.set_sourcetype
 * Ticket #164 Upgrade note in wiki FAQ
@@ -22 +23 @@
 * Interactive plotting annotations via optional argument interactive=True
 * Interactive creation of masks on the plotter - plotter.create_mask
+* Tidy up date range in asapplotter.plotazel/plotpointings
 
 Release 2.3.1 [2009-03-25]

trunk/python/__init__.py

@@ -393 +393 @@
 from linecatalog import linecatalog
 from opacity import skydip
+from opacity import model as opacity_model
 
 if rcParams['useplotter']:
@@ -412 +413 @@
 
 def is_ipython():
-    return '__IP' in dir(sys.modules["__main__"])
+    return 'IPython' in sys.modules.keys()
 
 if is_ipython():
@@ -630 +631 @@
                               scantable.create_mask
         skydip              - gain opacity values from a sky dip observation
+        opacity_model       - compute opacities fro given frequencies based on
+                              atmospheric model
 
     Note:

trunk/python/opacity.py

@@ -6 +6 @@
 from asap import selector
 from asap import rcParams
-from asap import xyplotter
+from asap._asap import atmosphere
+
+
+class model(object):
+    def _to_pascals(self, val):
+        if val > 2000:
+            return val
+        return val*100
+
+    def __init__(self, temperature=288, pressure=101325., humidity=0.5,
+                 elevation=700.):
+        """
+        This class implements opacity/atmospheric brightness temperature model
+        equivalent to the model available in MIRIAD. The actual math is a 
+        convertion of the Fortran code written by Bob Sault for MIRIAD.
+        It implements a simple model of the atmosphere and Liebe's model (1985) 
+        of the complex refractive index of air.
+
+        The model of the atmosphere is one with an exponential fall-off in
+        the water vapour content (scale height of 1540 m) and a temperature 
+        lapse rate of 6.5 mK/m. Otherwise the atmosphere obeys the ideal gas 
+        equation and hydrostatic equilibrium.
+
+        Note, the model includes atmospheric lines up to 800 GHz, but was not 
+        rigorously tested above 100 GHz and for instruments located at 
+        a significant elevation. For high-elevation sites it may be necessary to
+        adjust scale height and lapse rate.
+
+        Parameters:
+            temperature:    air temperature at the observatory (K)
+            pressure:       air pressure at the sea level if the observatory 
+                            elevation is set to non-zero value (note, by 
+                            default is set to 700m) or at the observatory
+                            ground level if the elevation is set to 0. (The 
+                            value is in Pascals or hPa, default 101325 Pa
+            humidity:       air humidity at the observatory (fractional), 
+                            default is 0.5
+            elevation:     observatory elevation about sea level (in meters)
+        """
+        self._atm = atmosphere(temp, self._to_pascals(pressure), humidity)
+        self.set_observatory_elevation(elevation):
+
+    def get_opacities(self, freq, elevation=None):
+        """Get the opacity value(s) for the fiven frequency(ies).
+        If no elevation is given the opacities for the zenith are returned.
+        If an elevation is specified refraction is also taken into account.
+        Parameters:
+            freq:       a frequency value in Hz, or a list of frequency values.
+                        One opacity value per frequency is returned as a scalar
+                        or list.
+            elevation:  the elevation at which to compute the opacity. If `None`
+                        is given (default) the zenith opacity is returned.
+
+
+        """
+        func = None
+        if isinstance(freq, (list, tuple)):
+            if elevation is None:
+                return self._atm.zenith_opacities(freq)
+            else:
+                elevation *= math.pi/180.
+                return self._atm.opacities(freq, elevation)
+        else:
+            if elevation is None:
+                return self._atm.zenith_opacity(freq)
+            else:
+                elevation *= math.pi/180.
+                return self._atm.opacity(freq, elevation)
+
+    def set_weather(self, temperature, pressure, humidity):
+        """Update the model using the given environmental parameters.
+        Parameters:
+            temperature:    air temperature at the observatory (K)
+            pressure:       air pressure at the sea level if the observatory 
+                            elevation is set to non-zero value (note, by 
+                            default is set to 700m) or at the observatory
+                            ground level if the elevation is set to 0. (The 
+                            value is in Pascals or hPa, default 101325 Pa
+            humidity:       air humidity at the observatory (fractional), 
+                            default is 0.5
+        """
+        pressure = self._to_pascals(pressure)
+        self._atm.set_weather(temperature, pressure, humidity)
+
+    def set_observatory_elevation(self, elevation:
+        """Update the model using the given the observatory elevation
+        Parameters:
+            elevation:  the elevation at which to compute the opacity. If `None`
+                        is given (default) the zenith opacity is returned.
+        """
+        self._atm.set_observatory_elevation(el)
+
 
 def _import_data(data):
@@ -27 +118 @@
     return merge(tables)
 
-def skydip(data, averagepol=True, tsky=300., plot=False):
+def skydip(data, averagepol=True, tsky=300., plot=False,
+           temperature=288, pressure=101325., humidity=0.5):
     """Determine the opacity from a set of 'skydip' obervations.
     This can be any set of observations over a range of elevations,
@@ -60 +152 @@
     rcsave = rcParams['verbose']
     rcParams['verbose'] = False
+    if plot:
+        from matplotlib import pylab
     scan = _import_data(data)
     f = fitter()
@@ -68 +162 @@
     pnos = scan.getpolnos()
     opacities = []
+    om = opacitymodel(temperature, pressure, humidity)
     for ino in inos:
         sel.set_ifs(ino)
@@ -74 +169 @@
         airms = []
         tsyss = []
-
         if plot:
-            xyplotter.cla()
-            xyplotter.ioff()
-            xyplotter.clf()
-            xyplotter.xlabel("Airmass")
-            xyplotter.ylabel(r"$T_{sys}$")
+            pylab.cla()
+            pylab.ioff()
+            pylab.clf()
+            pylab.xlabel("Airmass")
+            pylab.ylabel(r"$T_{sys}$")
         for pno in pnos:
             sel.set_polarisations(pno)
@@ -102 +196 @@
         if plot:
             colors = ['b','g','k']
-            for i in range(len(airms)):
-                xyplotter.plot(airms[i], tsyss[i], 'o', color=colors[i])
-                xyplotter.plot(airms[i], fits[i], '-', color=colors[i])
-                xyplotter.figtext(0.7,0.3-(i/30.0),
+            n = len(airms)
+            for i in range(n):
+                pylab.plot(airms[i], tsyss[i], 'o', color=colors[i])
+                pylab.plot(airms[i], fits[i], '-', color=colors[i])
+                pylab.figtext(0.7,0.3-(i/30.0),
                                   r"$\tau_{fit}=%0.2f$" % opacity[i],
                                   color=colors[i])
             if averagepol:
-                xyplotter.figtext(0.7,0.3-(len(airms)/30.0),
-                                  r"$\tau=%0.2f$" % opacities[-1],
+                pylab.figtext(0.7,0.3-(n/30.0),
+                                  r"$\tau_{avg}=%0.2f$" % opacities[-1],
                                   color='r')
-            xyplotter.title("IF%d : %s" % (ino, freqstr))
-
-            xyplotter.ion()
-            xyplotter.draw()
+                n +=1
+            pylab.figtext(0.7,0.3-(n/30.0),
+                          r"$\tau_{model}=%0.2f$" % om.get_opacities(freq*1e9),
+                          color='grey')
+            
+            pylab.title("IF%d : %s" % (ino, freqstr))
+
+            pylab.ion()
+            pylab.draw()
             raw_input("Hit <return> for next fit...")
         sel.reset()
@@ -122 +222 @@
     rcParams['verbose'] = rcsave
     if plot:
-        xyplotter.close()
+        pylab.close()
     return opacities

trunk/python/scantable.py

@@ -1115 +1115 @@
 
     @print_log_dec
-    def opacity(self, tau, insitu=None):
+    def opacity(self, tau=None, insitu=None):
         """
         Apply an opacity correction. The data
@@ -1126 +1126 @@
                          nIF*nPol or 1 and in order of IF/POL, e.g.
                          [opif0pol0, opif0pol1, opif1pol0 ...]
+                         if tau is `None` the opacities are determined from a
+                         model.
             insitu:      if False a new scantable is returned.
                          Otherwise, the scaling is done in-situ

trunk/src/STFiller.cpp

@@ -241 +241 @@
     if ( status != 0 ) break;
     n += 1;
-    
     Regex filterrx(".*[SL|PA]$");
     Regex obsrx("^AT.+");

trunk/src/python_STAtmosphere.cpp

@@ -44 +44 @@
         .def( init < double, double, double > () )
         .def("set_weather", &STAtmosphere::setWeather)
-        .def("set_observatory_elevation", &STAtmosphere::setObservatoryElevation)
+        .def("set_observatory_elevation", 
+             &STAtmosphere::setObservatoryElevation)
         .def("zenith_opacity", &STAtmosphere::zenithOpacity)
         .def("zenith_opacities", &STAtmosphere::zenithOpacities)

trunk/test/test_scantable.py

@@ -1 +1 @@
 import unittest
 import datetime
-from asap import scantable, selector, rcParams
+from asap import scantable, selector, rcParams, mask_not
 rcParams["verbose"] = False
 
@@ -8 +8 @@
 
         self.st = scantable("data/MOPS.rpf", average=True)
+        restfreqs = [86.243]     # 13CO-1/0, SiO the two IF
+        self.st.set_restfreqs(restfreqs,"GHz")
 
     def test_init(self):
@@ -135 +137 @@
         self.assertEqual(s1.getscannos(), (1,))
 
+    def test_flag(self):
+        q = self.st.auto_quotient()
+        q.set_unit('km/s')
+        q0 = q.copy()
+        q1 = q.copy()
+        msk = q0.create_mask([-10,20])
+        q0.flag(mask=mask_not(msk))
+        self.assertAlmostEqual(q0.stats(stat='max')[0], 95.62171936)
+        q1.flag(mask=msk)
+        self.assertAlmostEqual(q1.stats(stat='max')[0], 2.66563416)
 
 if __name__ == '__main__':