Changeset 1389 for branches


Ignore:
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:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/alma/python/asapfitter.py

    r1306 r1389  
    2525        self._p = None
    2626        self._selection = None
     27        self.uselinear = False
    2728
    2829    def set_data(self, xdat, ydat, mask=None):
     
    7475        Set the function to be fit.
    7576        Parameters:
    76             poly:    use a polynomial of the order given
     77            poly:    use a polynomial of the order given with nonlinear least squares fit
     78            lpoly:   use polynomial of the order given with linear least squares fit
    7779            gauss:   fit the number of gaussian specified
    7880        Example:
    7981            fitter.set_function(gauss=2) # will fit two gaussians
    80             fitter.set_function(poly=3)  # will fit a 3rd order polynomial
     82            fitter.set_function(poly=3)  # will fit a 3rd order polynomial via nonlinear method
     83            fitter.set_function(lpoly=3)  # will fit a 3rd order polynomial via linear method
    8184        """
    8285        #default poly order 0
     
    8689            n = kwargs.get('poly')
    8790            self.components = [n]
     91            self.uselinear = False
     92        elif kwargs.has_key('lpoly'):
     93            self.fitfunc = 'poly'
     94            n = kwargs.get('lpoly')
     95            self.components = [n]
     96            self.uselinear = True
    8897        elif kwargs.has_key('gauss'):
    8998            n = kwargs.get('gauss')
     
    91100            self.fitfuncs = [ 'gauss' for i in range(n) ]
    92101            self.components = [ 3 for i in range(n) ]
     102            self.uselinear = False
    93103        else:
    94104            msg = "Invalid function type."
     
    146156            if len(fxdpar) and fxdpar.count(0) == 0:
    147157                 raise RuntimeError,"No point fitting, if all parameters are fixed."
    148             converged = self.fitter.fit()
     158            if self.uselinear:
     159                converged = self.fitter.lfit()
     160            else:
     161                converged = self.fitter.fit()
    149162            if not converged:
    150163                raise RuntimeError,"Fit didn't converge."
     
    501514                             array(self.data._getmask(self._fittedrow),
    502515                                   copy=False))
    503            
     516                             
    504517            ylab = self.data._get_ordinate_label()
    505518
  • branches/alma/python/asapmath.py

    r1362 r1389  
    22from asap import rcParams
    33from asap import print_log
     4from asap import selector
    45
    56def average_time(*args, **kwargs):
     
    113114    return s
    114115
     116def dototalpower(calon, caloff, tcalval=0.0):
     117    """
     118    Do calibration for CAL on,off signals.
     119    Adopted from GBTIDL dototalpower
     120    Parameters:
     121        calon:         the 'cal on' subintegration
     122        caloff:        the 'cal off' subintegration
     123        tcalval:       user supplied Tcal value
     124    """
     125    varlist = vars()
     126    from asap._asap import stmath
     127    stm = stmath()
     128    stm._setinsitu(False)
     129    s = scantable(stm._dototalpower(calon, caloff, tcalval))
     130    s._add_history("dototalpower",varlist)
     131    print_log()
     132    return s
     133
     134def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
     135    """
     136    Calculate a quotient (sig-ref/ref * Tsys)
     137    Adopted from GBTIDL dosigref
     138    Parameters:
     139        sig:         on source data
     140        ref:         reference data
     141        smooth:      width of box car smoothing for reference
     142        tsysval:     user specified Tsys (scalar only)
     143        tauval:      user specified Tau (required if tsysval is set)
     144    """
     145    varlist = vars()
     146    from asap._asap import stmath
     147    stm = stmath()
     148    stm._setinsitu(False)
     149    s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
     150    s._add_history("dosigref",varlist)
     151    print_log()
     152    return s
     153
     154def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
     155    """
     156    Calibrate GBT position switched data
     157    Adopted from GBTIDL getps
     158    Currently calps identify the scans as position switched data if they
     159    contain '_ps' in the source name. The data must contains 'CAL' signal
     160    on/off in each integration. To identify 'CAL' on state, the word, 'calon'
     161    need to be present in the source name field.
     162    (GBT MS data reading process to scantable automatically append these
     163    id names to the source names)
     164
     165    Parameters:
     166        scantab:       scantable
     167        scannos:       list of scan numbers
     168        smooth:        optional box smoothing order for the reference
     169                       (default is 1 = no smoothing)
     170        tsysval:       optional user specified Tsys (default is 0.0,
     171                       use Tsys in the data)
     172        tauval:        optional user specified Tau
     173        tcalval:       optional user specified Tcal (default is 0.0,
     174                       use Tcal value in the data)
     175    """
     176    varlist = vars()
     177    # check for the appropriate data
     178    s = scantab.get_scan('*_ps*')
     179    if s is None:
     180        msg = "The input data appear to contain no position-switch mode data."
     181        if rcParams['verbose']:
     182            print msg
     183            return
     184        else:
     185            raise TypeError(msg)
     186    ssub = s.get_scan(scannos)
     187    if ssub is None:
     188        msg = "No data was found with given scan numbers!"
     189        if rcParams['verbose']:
     190            print msg
     191            return
     192        else:
     193            raise TypeError(msg)
     194    ssubon = ssub.get_scan('*calon')
     195    ssuboff = ssub.get_scan('*[^calon]')
     196    if ssubon.nrow() != ssuboff.nrow():
     197        msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers."
     198        if rcParams['verbose']:
     199            print msg
     200            return
     201        else:
     202            raise TypeError(msg)
     203    cals = dototalpower(ssubon, ssuboff, tcalval)
     204    sig = cals.get_scan('*ps')
     205    ref = cals.get_scan('*psr')
     206    if sig.nscan() != ref.nscan():
     207        msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers."
     208        if rcParams['verbose']:
     209            print msg
     210            return
     211        else:
     212            raise TypeError(msg)
     213
     214    #for user supplied Tsys
     215    if tsysval>0.0:
     216        if tauval<=0.0:
     217            msg = "Need to supply a valid tau to use the supplied Tsys"
     218            if rcParams['verbose']:
     219                print msg
     220                return
     221            else:
     222                raise TypeError(msg)
     223        else:
     224            sig.recalc_azel()
     225            ref.recalc_azel()
     226            #msg = "Use of user specified Tsys is not fully implemented yet."
     227            #if rcParams['verbose']:
     228            #    print msg
     229            #    return
     230            #else:
     231            #    raise TypeError(msg)
     232            # use get_elevation to get elevation and
     233            # calculate a scaling factor using the formula
     234            # -> tsys use to dosigref
     235
     236    #ress = dosigref(sig, ref, smooth, tsysval)
     237    ress = dosigref(sig, ref, smooth, tsysval, tauval)
     238    ress._add_history("calps", varlist)
     239    print_log()
     240    return ress
     241
     242def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
     243    """
     244    Do full (but a pair of scans at time) processing of GBT Nod data
     245    calibration.
     246    Adopted from  GBTIDL's getnod
     247    Parameters:
     248        scantab:     scantable
     249        scannos:     a pair of scan numbers, or the first scan number of the pair
     250        smooth:      box car smoothing order
     251        tsysval:     optional user specified Tsys value
     252        tauval:      optional user specified tau value (not implemented yet)
     253        tcalval:     optional user specified Tcal value
     254    """
     255    varlist = vars()
     256    from asap._asap import stmath
     257    stm = stmath()
     258    stm._setinsitu(False)
     259
     260    # check for the appropriate data
     261    s = scantab.get_scan('*_nod*')
     262    if s is None:
     263        msg = "The input data appear to contain no Nod observing mode data."
     264        if rcParams['verbose']:
     265            print msg
     266            return
     267        else:
     268            raise TypeError(msg)
     269
     270    # need check correspondance of each beam with sig-ref ...
     271    # check for timestamps, scan numbers, subscan id (not available in
     272    # ASAP data format...). Assume 1st scan of the pair (beam 0 - sig
     273    # and beam 1 - ref...)
     274    # First scan number of paired scans or list of pairs of
     275    # scan numbers (has to have even number of pairs.)
     276
     277    #data splitting
     278    scan1no = scan2no = 0
     279
     280    if len(scannos)==1:
     281        scan1no = scannos[0]
     282        scan2no = scannos[0]+1
     283        pairScans = [scan1no, scan2no]
     284    else:
     285        #if len(scannos)>2:
     286        #    msg = "calnod can only process a pair of nod scans at time."
     287        #    if rcParams['verbose']:
     288        #        print msg
     289        #        return
     290        #    else:
     291        #        raise TypeError(msg)
     292        #
     293        #if len(scannos)==2:
     294        #    scan1no = scannos[0]
     295        #    scan2no = scannos[1]
     296        pairScans = list(scannos)
     297
     298    if tsysval>0.0:
     299        if tauval<=0.0:
     300            msg = "Need to supply a valid tau to use the supplied Tsys"
     301            if rcParams['verbose']:
     302                print msg
     303                return
     304            else:
     305                raise TypeError(msg)
     306        else:
     307            scantab.recalc_azel()
     308    resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval))
     309    resspec._add_history("calnod",varlist)
     310    print_log()
     311    return resspec
     312
     313def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
     314    """
     315    Calibrate GBT frequency switched data.
     316    Adopted from GBTIDL getfs.
     317    Currently calfs identify the scans as frequency switched data if they
     318    contain '_fs' in the source name. The data must contains 'CAL' signal
     319    on/off in each integration. To identify 'CAL' on state, the word, 'calon'
     320    need to be present in the source name field.
     321    (GBT MS data reading via scantable automatically append these
     322    id names to the source names)
     323
     324    Parameters:
     325        scantab:       scantable
     326        scannos:       list of scan numbers
     327        smooth:        optional box smoothing order for the reference
     328                       (default is 1 = no smoothing)
     329        tsysval:       optional user specified Tsys (default is 0.0,
     330                       use Tsys in the data)
     331        tauval:        optional user specified Tau
     332    """
     333    varlist = vars()
     334    from asap._asap import stmath
     335    stm = stmath()
     336    stm._setinsitu(False)
     337
     338#    check = scantab.get_scan('*_fs*')
     339#    if check is None:
     340#        msg = "The input data appear to contain no Nod observing mode data."
     341#        if rcParams['verbose']:
     342#            print msg
     343#            return
     344#        else:
     345#            raise TypeError(msg)
     346    s = scantab.get_scan(scannos)
     347    del scantab
     348
     349    resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
     350    resspec._add_history("calfs",varlist)
     351    print_log()
     352    return resspec
     353
    115354def simple_math(left, right, op='add', tsys=True):
    116355    """
  • branches/alma/python/asapplotter.py

    r1358 r1389  
    741741        return userlabel or d[mode]
    742742
     743    def plotazel(self, scan=None):
     744        """
     745        plot azimuth and elevation  versus time of a scantable
     746        """
     747        import pylab as PL
     748        from matplotlib.dates import DateFormatter, timezone, HourLocator, MinuteLocator, DayLocator
     749        from matplotlib.ticker import MultipleLocator
     750        from matplotlib.numerix import array, pi
     751        self._data = scan
     752        dates = self._data.get_time()
     753        t = PL.date2num(dates)
     754        tz = timezone('UTC')
     755        PL.cla()
     756        PL.ioff()
     757        PL.clf()
     758        tdel = max(t) - min(t)
     759        ax = PL.subplot(2,1,1)
     760        el = array(self._data.get_elevation())*180./pi
     761        PL.ylabel('El [deg.]')
     762        dstr = dates[0].strftime('%Y/%m/%d')
     763        if tdel > 1.0:
     764            dstr2 = dates[len(dates)-1].strftime('%Y/%m/%d')
     765            dstr = dstr + " - " + dstr2
     766            majloc = DayLocator()
     767            minloc = HourLocator(range(0,23,12))
     768            timefmt = DateFormatter("%b%d")
     769        else:
     770            timefmt = DateFormatter('%H')
     771            majloc = HourLocator()
     772            minloc = MinuteLocator(20)
     773        PL.title(dstr)
     774        PL.plot_date(t,el,'b,', tz=tz)
     775        #ax.grid(True)
     776        ax.yaxis.grid(True)
     777        yloc = MultipleLocator(30)
     778        ax.set_ylim(0,90)
     779        ax.xaxis.set_major_formatter(timefmt)
     780        ax.xaxis.set_major_locator(majloc)
     781        ax.xaxis.set_minor_locator(minloc)
     782        ax.yaxis.set_major_locator(yloc)
     783        if tdel > 1.0:
     784            labels = ax.get_xticklabels()
     785        #    PL.setp(labels, fontsize=10, rotation=45)
     786            PL.setp(labels, fontsize=10)
     787        # Az plot
     788        az = array(self._data.get_azimuth())*180./pi
     789        if min(az) < 0:
     790            for irow in range(len(az)):
     791                if az[irow] < 0: az[irow] += 360.0
     792
     793        ax = PL.subplot(2,1,2)
     794        PL.xlabel('Time (UT)')
     795        PL.ylabel('Az [deg.]')
     796        PL.plot_date(t,az,'b,', tz=tz)
     797        ax.set_ylim(0,360)
     798        #ax.grid(True)
     799        ax.yaxis.grid(True)
     800        #hfmt = DateFormatter('%H')
     801        #hloc = HourLocator()
     802        yloc = MultipleLocator(60)
     803        ax.xaxis.set_major_formatter(timefmt)
     804        ax.xaxis.set_major_locator(majloc)
     805        ax.xaxis.set_minor_locator(minloc)
     806        ax.yaxis.set_major_locator(yloc)
     807        if tdel > 1.0:
     808            labels = ax.get_xticklabels()
     809            PL.setp(labels, fontsize=10)
     810        PL.ion()
     811        PL.draw()
     812
     813    def plotpointing(self, scan=None):
     814        """
     815        plot telescope pointings
     816        """
     817        import pylab as PL
     818        from matplotlib.dates import DateFormatter, timezone
     819        from matplotlib.ticker import MultipleLocator
     820        from matplotlib.numerix import array, pi, zeros
     821        self._data = scan
     822        dir = array(self._data.get_directionval()).transpose()
     823        ra = dir[0]*180./pi
     824        dec = dir[1]*180./pi
     825        PL.cla()
     826        PL.ioff()
     827        PL.clf()
     828        ax = PL.axes([0.1,0.1,0.8,0.8])
     829        ax = PL.axes([0.1,0.1,0.8,0.8])
     830        ax.set_aspect('equal')
     831        PL.plot(ra,dec, 'b,')
     832        PL.xlabel('RA [deg.]')
     833        PL.ylabel('Declination [deg.]')
     834        PL.title('Telescope pointings')
     835        [xmin,xmax,ymin,ymax] = PL.axis()
     836        PL.axis([xmax,xmin,ymin,ymax])
     837        PL.ion()
     838        PL.draw()
     839
  • branches/alma/python/scantable.py

    r1373 r1389  
    516516        """
    517517        return self._get_column(self._getdirection, row)
     518
     519    def get_directionval(self, row=-1):
     520        """
     521        Get a list of Positions on the sky (direction) for the observations.
     522        Return a float for each integration in the scantable.
     523        Parameters:
     524            row:    row no of integration. Default -1 return all rows
     525        Example:
     526            none
     527        """
     528        return self._get_column(self._getdirectionvec, row)
     529
    518530
    519531    def set_unit(self, unit='channel'):
     
    12151227
    12161228
    1217     def poly_baseline(self, mask=None, order=0, plot=False, insitu=None):
     1229    def poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None):
    12181230        """
    12191231        Return a scan which has been baselined (all rows) by a polynomial.
     
    12241236                        indivual fit has to be approved, by typing 'y'
    12251237                        or 'n'
     1238            uselin:     use linear polynomial fit
    12261239            insitu:     if False a new scantable is returned.
    12271240                        Otherwise, the scaling is done in-situ
     
    12401253            f = fitter()
    12411254            f.set_scan(self, mask)
    1242             f.set_function(poly=order)
     1255            #f.set_function(poly=order)
     1256            if uselin:
     1257                f.set_function(lpoly=order)
     1258            else:
     1259                f.set_function(poly=order)
    12431260            s = f.auto_fit(insitu, plot=plot)
    12441261            s._add_history("poly_baseline", varlist)
Note: See TracChangeset for help on using the changeset viewer.