How to use presto - 10 common examples

Outputs:
                None

            *** Subbanding happens in-place ***
        """
        assert (self.numchans % nsub) == 0
        assert (subdm is None) or (subdm >= 0)
        nchan_per_sub = self.numchans // nsub
        sub_hifreqs = self.freqs[np.arange(nsub)*nchan_per_sub]
        sub_lofreqs = self.freqs[(1+np.arange(nsub))*nchan_per_sub-1]
        sub_ctrfreqs = 0.5*(sub_hifreqs+sub_lofreqs)
        
        if subdm is not None:
            # Compute delays
            ref_delays = psr_utils.delay_from_DM(subdm-self.dm, sub_ctrfreqs)
            delays = psr_utils.delay_from_DM(subdm-self.dm, self.freqs)
            rel_delays = delays-ref_delays.repeat(nchan_per_sub) # Relative delay
            rel_bindelays = np.round(rel_delays/self.dt).astype('int')
            # Shift channels
            self.shift_channels(rel_bindelays, padval)

        # Subband
        self.data = np.array([np.sum(sub, axis=0) for sub in \
                                np.vsplit(self.data, nsub)])
        self.freqs = sub_ctrfreqs
        self.numchans = nsub

Outputs:
                None

            *** Subbanding happens in-place ***
        """
        assert (self.numchans % nsub) == 0
        assert (subdm is None) or (subdm >= 0)
        nchan_per_sub = self.numchans // nsub
        sub_hifreqs = self.freqs[np.arange(nsub)*nchan_per_sub]
        sub_lofreqs = self.freqs[(1+np.arange(nsub))*nchan_per_sub-1]
        sub_ctrfreqs = 0.5*(sub_hifreqs+sub_lofreqs)
        
        if subdm is not None:
            # Compute delays
            ref_delays = psr_utils.delay_from_DM(subdm-self.dm, sub_ctrfreqs)
            delays = psr_utils.delay_from_DM(subdm-self.dm, self.freqs)
            rel_delays = delays-ref_delays.repeat(nchan_per_sub) # Relative delay
            rel_bindelays = np.round(rel_delays/self.dt).astype('int')
            # Shift channels
            self.shift_channels(rel_bindelays, padval)

        # Subband
        self.data = np.array([np.sum(sub, axis=0) for sub in \
                                np.vsplit(self.data, nsub)])
        self.freqs = sub_ctrfreqs
        self.numchans = nsub

short_nM = options.nM // 1000000

    # The basename of the data files
    if argv[1].endswith(".dat"):
        basename = "../"+argv[1][:-4]
    else:
        basename = "../"+argv[1]

    # Get the bird file (the first birdie file in the directory!)
    birdname = glob("../*.birds")
    if birdname:
        birdname = birdname[0]

    outnamebase = options.outdir+basename[3:]
    inf = read_inffile(basename)
    idata = infodata.infodata(basename+".inf")
    N = inf.N
    t0i = inf.mjd_i
    t0f = inf.mjd_f
    num = 0
    point = 0
    T = options.nM * inf.dt / 86400.0
    baryv = get_baryv(idata.RA, idata.DEC, idata.epoch, T, obs='GB')
    print("Baryv = ", baryv)
    inf.N = options.nM
    inf.numonoff = 0
    nM = options.nM // 1000000
    while point + options.nM < N:
        pM = point // 1000000
        outname = basename[3:]+'_%03dM'%nM+'_%02d'%num
        stdout.write('\n'+outname+'\n\n')
        inf.name = outname

else:
            self.bytes_per_spectra = (self.bits_per_sample * self.samples_per_spectra) // 8
        self.samples_per_subint = self.samples_per_spectra * self.spectra_per_subint
        self.bytes_per_subint = self.bytes_per_spectra * self.spectra_per_subint

        # Flip the band?
        if self.hi_freq < self.lo_freq:
            tmp = self.hi_freq
            self.hi_freq = self.lo_freq
            self.lo_freq = tmp
            self.df *= -1.0
            self.need_flipband = True
        # Compute the bandwidth
        self.BW = self.num_channels * self.df
        self.mjd = int(self.start_MJD[0])
        self.secs = (self.start_MJD[0] % 1)*pc.SECPERDAY

def gen_gaussians(params, N):
    """
    gen_gaussians(params, N):
        Return a model of a DC-component + M gaussians
            params is a sequence of 1+M*3 values
                the first value is the DC component.  Each remaining
                group of three represents the gaussians phase (0-1),
                FWHM (0-1), and amplitude (>0.0).
            N is the number of points in the model.
    """
    numgaussians = (len(params)-1) // 3
    model = Num.zeros(N, dtype='d') + params[0]
    for ii in range(numgaussians):
        phase, FWHM, amp = params[1+ii*3:4+ii*3]
        model += amp * gaussian_profile(N, phase, FWHM)
    return model

# Try to see how many freq derivs we have
        fs = [x.F0]
        for ii in range(1, 20):  # hopefully 20 is an upper limit!
            attrib = "F%d"%ii
            if hasattr(x, attrib):
                fs.append(getattr(x, attrib))
            else:
                break
        epoch = x.PEPOCH
        T = (x.FINISH - x.START) * 86400.0
    else:
        print("I don't recognize the file type for", filenm)
        sys.exit()
    newts = epoch + num.arange(int(T/10.0+0.5), dtype=num.float)/8640.0
    time = num.concatenate((time, newts))
    newps = 1.0 / pu.calc_freq(newts, epoch, *fs)
    period = num.concatenate((period, newps))
    print("%13.7f (%0.1f sec): " % (epoch, T), fs)

#!/usr/bin/env python
from __future__ import (print_function,division)
import presto.psr_utils as pu
import sys
from presto.infodata import infodata

if len(sys.argv) != 2:
    print("chooseN ")
    print("    Prints a good value for fast FFTs to be used for -numout in prepdata/prepsubband")
    sys.exit(1)

if sys.argv[1].endswith('.inf'):
    inf = infodata(sys.argv[1])
    n = inf.N
else:
    try:
        n = int(sys.argv[1])
    except:
        print("chooseN ")
        print("    Prints a good value for fast FFTs to be used for -numout in prepdata/prepsubband")
        sys.exit(2)

print(pu.choose_N(n))

opts.is_no_detrend = True  # we don't do detrending for events
		if opts.mjd == '' and not opts.is_chandra:
			print("Error: for events' file start MJD _must_ be given with --mjd option or --chandra option!")
			sys.exit(1)
		if opts.nbins == -1:
			print("Error: number of bins _must_ be given with --nbins option!")
			sys.exit(1)
		if opts.rebin != 1:
			print("Event data can not be re-binned")
			opts.rebin = 1
	else:
		if (opts.mjd == '' and not opts.is_chandra) or opts.tsamp == '' or opts.psrname == '':
			# reading inf-file to get corresponding info
			inffile = datfile.split(".dat")[0] + ".inf"
			try:
				id = inf.infodata(inffile)
				tsamp = id.dt        # sampling time
				startmjd = id.epoch  # start MJD
                		source = id.object   # pulsar name
			except:
				print("Error: Can't read the inf-file '%s'!" % (inffile,))
				sys.exit(1)

	# overwriting MJD, tsamp, pulsarname from the command line if given
	if opts.is_chandra:
		opts.mjd = "50814.0"
		print("Chandra event file. Reference start MJD is %s" % (opts.mjd))
	if opts.mjd != '':
		startmjd = float(opts.mjd)
	if opts.tsamp != '':
		tsamp = float(opts.tsamp)
	if opts.psrname != '':

def parse_eph(filenm):
    global period, time
    suffix = filenm.split(".")[-1]
    if suffix=="bestprof":
        x = bestprof.bestprof(filenm)
        fs = pu.p_to_f(x.p0_bary, x.p1_bary, x.p2_bary)
        epoch = x.epochi_bary + x.epochf_bary
        T = x.T
    elif suffix=="par":
        x = parfile.psr_par(filenm)
        # Try to see how many freq derivs we have
        fs = [x.F0]
        for ii in range(1, 20):  # hopefully 20 is an upper limit!
            attrib = "F%d"%ii
            if hasattr(x, attrib):
                fs.append(getattr(x, attrib))
            else:
                break
        epoch = x.PEPOCH
        T = (x.FINISH - x.START) * 86400.0
    else:
        print("I don't recognize the file type for", filenm)

sp_pgplot.ppgplot.pgslw(5)
            sweepstart = sweeped_start - 0.2 * sweep_duration
            sp_pgplot.ppgplot.pgsci(0)
            sp_pgplot.ppgplot.pgline(delays + sweepstart, freqs)
            sp_pgplot.ppgplot.pgsci(1)

        #### Figure texts 
        if integrate_spec:
            sp_pgplot.ppgplot.pgsvp(0.81, 0.97, 0.64, 0.909)
            sp_pgplot.ppgplot.pgsch(0.62)
        else:
            sp_pgplot.ppgplot.pgsvp(0.745, 0.97, 0.64, 0.909)
            sp_pgplot.ppgplot.pgsch(0.7)
        sp_pgplot.ppgplot.pgslw(3)
        sp_pgplot.ppgplot.pgmtxt('T', -1.1, 0.01, 0.0, "RA: %s" % RA)
        sp_pgplot.ppgplot.pgmtxt('T', -2.6, 0.01, 0.0, "DEC: %s" % dec)
        sp_pgplot.ppgplot.pgmtxt('T', -4.1, 0.01, 0.0, "MJD: %f" % MJD)
        sp_pgplot.ppgplot.pgmtxt('T', -5.6, 0.01, 0.0, "Obs date: %s %s %s" % (date[0], date[1], date[2]))
        sp_pgplot.ppgplot.pgmtxt('T', -7.1, 0.01, 0.0, "Telescope: %s" % telescope)
        sp_pgplot.ppgplot.pgmtxt('T', -8.6, 0.01, 0.0, "DM: %.2f pc cm\\u-3\\d" % dm)
        if sigma:
            sp_pgplot.ppgplot.pgmtxt('T', -10.1, 0.01, 0.0, "S/N\\dMAX\\u: %.2f" % sigma)
        else:
            sp_pgplot.ppgplot.pgmtxt('T', -10.1, 0.01, 0.0, "S/N\\dMAX\\u: N/A")
        sp_pgplot.ppgplot.pgmtxt('T', -11.6, 0.01, 0.0, "Number of samples: %i" % nbins)
        sp_pgplot.ppgplot.pgmtxt('T', -13.1, 0.01, 0.0, "Number of subbands: %i" % nsub)
        sp_pgplot.ppgplot.pgmtxt('T', -14.6, 0.01, 0.0, "Pulse width: %.2f ms" % (pulse_width * 1e3))
        sp_pgplot.ppgplot.pgmtxt('T', -16.1, 0.01, 0.0, "Sampling time: %.3f \\gms" % (tsamp * 1e6))
        sp_pgplot.ppgplot.pgmtxt('T', -17.6, 0.0, 0.0, "Bary pulse peak time: %.2f s" % (bary_start))
        sp_pgplot.ppgplot.pgsvp(0.07, 0.7, 0.01, 0.05)
        sp_pgplot.ppgplot.pgmtxt('T', -2.1, 0.01, 0.0, "%s" % fn)