Changeset 1389 for branches

Timestamp:

07/27/07 02:08:12 (17 years ago)

Author:

TakTsutsumi

Message:

merged from NRAO version of ASAP2.1 with ALMA specific modifications

Location:

branches/alma/python

Files:

• asapfitter.py (modified) (6 diffs)
• asapmath.py (modified) (2 diffs)
• asapplotter.py (modified) (1 diff)
• scantable.py (modified) (4 diffs)

branches/alma/python/asapfitter.py

@@ -25 +25 @@
         self._p = None
         self._selection = None
+        self.uselinear = False
 
     def set_data(self, xdat, ydat, mask=None):
@@ -74 +75 @@
         Set the function to be fit.
         Parameters:
-            poly:    use a polynomial of the order given
+            poly:    use a polynomial of the order given with nonlinear least squares fit 
+            lpoly:   use polynomial of the order given with linear least squares fit
             gauss:   fit the number of gaussian specified
         Example:
             fitter.set_function(gauss=2) # will fit two gaussians
-            fitter.set_function(poly=3)  # will fit a 3rd order polynomial
+            fitter.set_function(poly=3)  # will fit a 3rd order polynomial via nonlinear method
+            fitter.set_function(lpoly=3)  # will fit a 3rd order polynomial via linear method
         """
         #default poly order 0
@@ -86 +89 @@
             n = kwargs.get('poly')
             self.components = [n]
+            self.uselinear = False 
+        elif kwargs.has_key('lpoly'):
+            self.fitfunc = 'poly'
+            n = kwargs.get('lpoly')
+            self.components = [n]
+            self.uselinear = True
         elif kwargs.has_key('gauss'):
             n = kwargs.get('gauss')
@@ -91 +100 @@
             self.fitfuncs = [ 'gauss' for i in range(n) ]
             self.components = [ 3 for i in range(n) ]
+            self.uselinear = False 
         else:
             msg = "Invalid function type."
@@ -146 +156 @@
             if len(fxdpar) and fxdpar.count(0) == 0:
                  raise RuntimeError,"No point fitting, if all parameters are fixed."
-            converged = self.fitter.fit()
+            if self.uselinear:
+                converged = self.fitter.lfit()
+            else:
+                converged = self.fitter.fit()
             if not converged:
                 raise RuntimeError,"Fit didn't converge."
@@ -501 +514 @@
                              array(self.data._getmask(self._fittedrow),
                                    copy=False))
-            
+                             
             ylab = self.data._get_ordinate_label()
 

branches/alma/python/asapmath.py

@@ -2 +2 @@
 from asap import rcParams
 from asap import print_log
+from asap import selector
 
 def average_time(*args, **kwargs):
@@ -113 +114 @@
     return s
 
+def dototalpower(calon, caloff, tcalval=0.0):
+    """
+    Do calibration for CAL on,off signals.
+    Adopted from GBTIDL dototalpower
+    Parameters:
+        calon:         the 'cal on' subintegration
+        caloff:        the 'cal off' subintegration
+        tcalval:       user supplied Tcal value
+    """
+    varlist = vars()
+    from asap._asap import stmath
+    stm = stmath()
+    stm._setinsitu(False)
+    s = scantable(stm._dototalpower(calon, caloff, tcalval))
+    s._add_history("dototalpower",varlist)
+    print_log()
+    return s
+
+def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
+    """
+    Calculate a quotient (sig-ref/ref * Tsys)
+    Adopted from GBTIDL dosigref
+    Parameters:
+        sig:         on source data
+        ref:         reference data
+        smooth:      width of box car smoothing for reference
+        tsysval:     user specified Tsys (scalar only)
+        tauval:      user specified Tau (required if tsysval is set)
+    """
+    varlist = vars()
+    from asap._asap import stmath
+    stm = stmath()
+    stm._setinsitu(False)
+    s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
+    s._add_history("dosigref",varlist)
+    print_log()
+    return s
+
+def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
+    """
+    Calibrate GBT position switched data
+    Adopted from GBTIDL getps
+    Currently calps identify the scans as position switched data if they
+    contain '_ps' in the source name. The data must contains 'CAL' signal
+    on/off in each integration. To identify 'CAL' on state, the word, 'calon'
+    need to be present in the source name field.
+    (GBT MS data reading process to scantable automatically append these
+    id names to the source names)
+
+    Parameters:
+        scantab:       scantable
+        scannos:       list of scan numbers
+        smooth:        optional box smoothing order for the reference
+                       (default is 1 = no smoothing)
+        tsysval:       optional user specified Tsys (default is 0.0,
+                       use Tsys in the data)
+        tauval:        optional user specified Tau
+        tcalval:       optional user specified Tcal (default is 0.0,
+                       use Tcal value in the data)
+    """
+    varlist = vars()
+    # check for the appropriate data
+    s = scantab.get_scan('*_ps*')
+    if s is None:
+        msg = "The input data appear to contain no position-switch mode data."
+        if rcParams['verbose']:
+            print msg
+            return
+        else:
+            raise TypeError(msg)
+    ssub = s.get_scan(scannos)
+    if ssub is None:
+        msg = "No data was found with given scan numbers!"
+        if rcParams['verbose']:
+            print msg
+            return
+        else:
+            raise TypeError(msg)
+    ssubon = ssub.get_scan('*calon')
+    ssuboff = ssub.get_scan('*[^calon]')
+    if ssubon.nrow() != ssuboff.nrow():
+        msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers."
+        if rcParams['verbose']:
+            print msg
+            return
+        else:
+            raise TypeError(msg)
+    cals = dototalpower(ssubon, ssuboff, tcalval)
+    sig = cals.get_scan('*ps')
+    ref = cals.get_scan('*psr')
+    if sig.nscan() != ref.nscan():
+        msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers."
+        if rcParams['verbose']:
+            print msg
+            return
+        else:
+            raise TypeError(msg)
+
+    #for user supplied Tsys
+    if tsysval>0.0:
+        if tauval<=0.0:
+            msg = "Need to supply a valid tau to use the supplied Tsys"
+            if rcParams['verbose']:
+                print msg
+                return
+            else:
+                raise TypeError(msg)
+        else:
+            sig.recalc_azel()
+            ref.recalc_azel()
+            #msg = "Use of user specified Tsys is not fully implemented yet."
+            #if rcParams['verbose']:
+            #    print msg
+            #    return
+            #else:
+            #    raise TypeError(msg)
+            # use get_elevation to get elevation and
+            # calculate a scaling factor using the formula
+            # -> tsys use to dosigref
+
+    #ress = dosigref(sig, ref, smooth, tsysval)
+    ress = dosigref(sig, ref, smooth, tsysval, tauval)
+    ress._add_history("calps", varlist)
+    print_log()
+    return ress
+
+def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
+    """
+    Do full (but a pair of scans at time) processing of GBT Nod data
+    calibration.
+    Adopted from  GBTIDL's getnod
+    Parameters:
+        scantab:     scantable
+        scannos:     a pair of scan numbers, or the first scan number of the pair
+        smooth:      box car smoothing order
+        tsysval:     optional user specified Tsys value
+        tauval:      optional user specified tau value (not implemented yet)
+        tcalval:     optional user specified Tcal value
+    """
+    varlist = vars()
+    from asap._asap import stmath
+    stm = stmath()
+    stm._setinsitu(False)
+
+    # check for the appropriate data
+    s = scantab.get_scan('*_nod*')
+    if s is None:
+        msg = "The input data appear to contain no Nod observing mode data."
+        if rcParams['verbose']:
+            print msg
+            return
+        else:
+            raise TypeError(msg)
+
+    # need check correspondance of each beam with sig-ref ...
+    # check for timestamps, scan numbers, subscan id (not available in
+    # ASAP data format...). Assume 1st scan of the pair (beam 0 - sig
+    # and beam 1 - ref...)
+    # First scan number of paired scans or list of pairs of
+    # scan numbers (has to have even number of pairs.)
+
+    #data splitting
+    scan1no = scan2no = 0
+
+    if len(scannos)==1:
+        scan1no = scannos[0]
+        scan2no = scannos[0]+1
+        pairScans = [scan1no, scan2no]
+    else:
+        #if len(scannos)>2:
+        #    msg = "calnod can only process a pair of nod scans at time."
+        #    if rcParams['verbose']:
+        #        print msg
+        #        return
+        #    else:
+        #        raise TypeError(msg)
+        #
+        #if len(scannos)==2:
+        #    scan1no = scannos[0]
+        #    scan2no = scannos[1]
+        pairScans = list(scannos)
+
+    if tsysval>0.0:
+        if tauval<=0.0:
+            msg = "Need to supply a valid tau to use the supplied Tsys"
+            if rcParams['verbose']:
+                print msg
+                return
+            else:
+                raise TypeError(msg)
+        else:
+            scantab.recalc_azel()
+    resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval))
+    resspec._add_history("calnod",varlist)
+    print_log()
+    return resspec
+
+def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
+    """
+    Calibrate GBT frequency switched data.
+    Adopted from GBTIDL getfs.
+    Currently calfs identify the scans as frequency switched data if they
+    contain '_fs' in the source name. The data must contains 'CAL' signal
+    on/off in each integration. To identify 'CAL' on state, the word, 'calon'
+    need to be present in the source name field.
+    (GBT MS data reading via scantable automatically append these
+    id names to the source names)
+
+    Parameters:
+        scantab:       scantable
+        scannos:       list of scan numbers
+        smooth:        optional box smoothing order for the reference
+                       (default is 1 = no smoothing)
+        tsysval:       optional user specified Tsys (default is 0.0,
+                       use Tsys in the data)
+        tauval:        optional user specified Tau
+    """
+    varlist = vars()
+    from asap._asap import stmath
+    stm = stmath()
+    stm._setinsitu(False)
+
+#    check = scantab.get_scan('*_fs*')
+#    if check is None:
+#        msg = "The input data appear to contain no Nod observing mode data."
+#        if rcParams['verbose']:
+#            print msg
+#            return
+#        else:
+#            raise TypeError(msg)
+    s = scantab.get_scan(scannos)
+    del scantab
+
+    resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
+    resspec._add_history("calfs",varlist)
+    print_log()
+    return resspec
+
 def simple_math(left, right, op='add', tsys=True):
     """

branches/alma/python/asapplotter.py

@@ -741 +741 @@
         return userlabel or d[mode]
 
+    def plotazel(self, scan=None):
+        """
+        plot azimuth and elevation  versus time of a scantable
+        """
+        import pylab as PL
+        from matplotlib.dates import DateFormatter, timezone, HourLocator, MinuteLocator, DayLocator
+        from matplotlib.ticker import MultipleLocator
+        from matplotlib.numerix import array, pi
+        self._data = scan
+        dates = self._data.get_time()
+        t = PL.date2num(dates)
+        tz = timezone('UTC')
+        PL.cla()
+        PL.ioff()
+        PL.clf()
+        tdel = max(t) - min(t)
+        ax = PL.subplot(2,1,1)
+        el = array(self._data.get_elevation())*180./pi
+        PL.ylabel('El [deg.]')
+        dstr = dates[0].strftime('%Y/%m/%d')
+        if tdel > 1.0:
+            dstr2 = dates[len(dates)-1].strftime('%Y/%m/%d')
+            dstr = dstr + " - " + dstr2
+            majloc = DayLocator()
+            minloc = HourLocator(range(0,23,12))
+            timefmt = DateFormatter("%b%d")
+        else:
+            timefmt = DateFormatter('%H')
+            majloc = HourLocator()
+            minloc = MinuteLocator(20)
+        PL.title(dstr)
+        PL.plot_date(t,el,'b,', tz=tz)
+        #ax.grid(True)
+        ax.yaxis.grid(True)
+        yloc = MultipleLocator(30)
+        ax.set_ylim(0,90)
+        ax.xaxis.set_major_formatter(timefmt)
+        ax.xaxis.set_major_locator(majloc)
+        ax.xaxis.set_minor_locator(minloc)
+        ax.yaxis.set_major_locator(yloc)
+        if tdel > 1.0:
+            labels = ax.get_xticklabels()
+        #    PL.setp(labels, fontsize=10, rotation=45)
+            PL.setp(labels, fontsize=10)
+        # Az plot
+        az = array(self._data.get_azimuth())*180./pi
+        if min(az) < 0:
+            for irow in range(len(az)):
+                if az[irow] < 0: az[irow] += 360.0
+
+        ax = PL.subplot(2,1,2)
+        PL.xlabel('Time (UT)')
+        PL.ylabel('Az [deg.]')
+        PL.plot_date(t,az,'b,', tz=tz)
+        ax.set_ylim(0,360)
+        #ax.grid(True)
+        ax.yaxis.grid(True)
+        #hfmt = DateFormatter('%H')
+        #hloc = HourLocator()
+        yloc = MultipleLocator(60)
+        ax.xaxis.set_major_formatter(timefmt)
+        ax.xaxis.set_major_locator(majloc)
+        ax.xaxis.set_minor_locator(minloc)
+        ax.yaxis.set_major_locator(yloc)
+        if tdel > 1.0:
+            labels = ax.get_xticklabels()
+            PL.setp(labels, fontsize=10)
+        PL.ion()
+        PL.draw()
+
+    def plotpointing(self, scan=None):
+        """
+        plot telescope pointings
+        """
+        import pylab as PL
+        from matplotlib.dates import DateFormatter, timezone
+        from matplotlib.ticker import MultipleLocator
+        from matplotlib.numerix import array, pi, zeros
+        self._data = scan
+        dir = array(self._data.get_directionval()).transpose()
+        ra = dir[0]*180./pi
+        dec = dir[1]*180./pi
+        PL.cla()
+        PL.ioff()
+        PL.clf()
+        ax = PL.axes([0.1,0.1,0.8,0.8])
+        ax = PL.axes([0.1,0.1,0.8,0.8])
+        ax.set_aspect('equal')
+        PL.plot(ra,dec, 'b,')
+        PL.xlabel('RA [deg.]')
+        PL.ylabel('Declination [deg.]')
+        PL.title('Telescope pointings')
+        [xmin,xmax,ymin,ymax] = PL.axis()
+        PL.axis([xmax,xmin,ymin,ymax])
+        PL.ion()
+        PL.draw()
+

branches/alma/python/scantable.py

@@ -516 +516 @@
         """
         return self._get_column(self._getdirection, row)
+
+    def get_directionval(self, row=-1):
+        """
+        Get a list of Positions on the sky (direction) for the observations.
+        Return a float for each integration in the scantable.
+        Parameters:
+            row:    row no of integration. Default -1 return all rows
+        Example:
+            none
+        """
+        return self._get_column(self._getdirectionvec, row)
+
 
     def set_unit(self, unit='channel'):
@@ -1215 +1227 @@
 
 
-    def poly_baseline(self, mask=None, order=0, plot=False, insitu=None):
+    def poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None):
         """
         Return a scan which has been baselined (all rows) by a polynomial.
@@ -1224 +1236 @@
                         indivual fit has to be approved, by typing 'y'
                         or 'n'
+            uselin:     use linear polynomial fit
             insitu:     if False a new scantable is returned.
                         Otherwise, the scaling is done in-situ
@@ -1240 +1253 @@
             f = fitter()
             f.set_scan(self, mask)
-            f.set_function(poly=order)
+            #f.set_function(poly=order)
+            if uselin:
+                f.set_function(lpoly=order)
+            else:
+                f.set_function(poly=order)
             s = f.auto_fit(insitu, plot=plot)
             s._add_history("poly_baseline", varlist)