How to use the pint.toa function in Pint

def setUp(self):

        datadir = os.path.join(testdir, "datafile")
        # parfile = os.path.join(datadir, "J1744-1134.basic.par")
        # ngc = os.path.join(datadir, "NGC6440E")

        self.toas = toa.get_TOAs(os.path.join(datadir, "NGC6440E.tim"))
        self.model = pint.models.get_model(os.path.join(datadir, "NGC6440E.par"))
        self.modelWithTD = pint.models.get_model(os.path.join(datadir, "NGC6440E.par"))
        self.modelWithTD.CORRECT_TROPOSPHERE.value = True

        self.toasInvalid = toa.get_TOAs(os.path.join(datadir, "NGC6440E.tim"))

        for i in range(len(self.toasInvalid.table)):
            # adjust the timing by half a day to make them invalid
            self.toasInvalid.table["mjd"][i] += 0.5

        self.testAltitudes = np.arange(self.MIN_ALT, 90, 100) * u.deg
        self.testHeights = np.array([10, 100, 1000, 5000]) * u.m

        self.td = self.modelWithTD.components["TroposphereDelay"]

def setUpClass(cls):
        os.chdir(datadir)
        cls.parfileJ1853 = "J1853+1303_NANOGrav_11yv0.gls.par"
        cls.timfileJ1853 = "J1853+1303_NANOGrav_11yv0.tim"
        cls.toasJ1853 = toa.get_TOAs(cls.timfileJ1853, ephem="DE421", planets=False)
        cls.modelJ1853 = model.get_model(cls.parfileJ1853)
        cls.ltres, cls.ltbindelay = np.genfromtxt(
            cls.parfileJ1853 + ".tempo2_test", skip_header=1, unpack=True
        )

        cls.parfileJ0613 = "J0613-0200_NANOGrav_9yv1_ELL1H.gls.par"
        cls.timfileJ0613 = "J0613-0200_NANOGrav_9yv1.tim"
        cls.modelJ0613 = model.get_model(cls.parfileJ0613)
        cls.toasJ0613 = toa.get_TOAs(cls.timfileJ0613, ephem="DE421", planets=False)
        cls.parfileJ0613_STIG = "J0613-0200_NANOGrav_9yv1_ELL1H_STIG.gls.par"
        cls.modelJ0613_STIG = model.get_model(cls.parfileJ0613_STIG)

def setUpClass(self):
        os.chdir(datadir)
        self.parfileJ0023 = "J0023+0923_NANOGrav_11yv0.gls.par"
        self.timJ0023 = "J0023+0923_NANOGrav_11yv0.tim"
        self.toasJ0023 = toa.get_TOAs(self.timJ0023, ephem="DE436", planets=False)
        self.modelJ0023 = mb.get_model(self.parfileJ0023)
        # tempo result
        self.ltres, self.ltbindelay = np.genfromtxt(
            self.parfileJ0023 + ".tempo2_test", skip_header=1, unpack=True
        )

def setUpClass(cls):
        os.chdir(datadir)
        cls.par = "B1855+09_NANOGrav_9yv1.gls.par"
        cls.tim = "B1855+09_NANOGrav_9yv1.tim"
        cls.m = mb.get_model(cls.par)
        cls.t = toa.get_TOAs(cls.tim, ephem="DE436")
        cls.f = GLSFitter(cls.t, cls.m)
        # get tempo2 parameter dict
        with open("B1855+09_tempo2_gls_pars.json", "r") as fp:
            cls.t2d = json.load(fp)

def setUpClass(cls):
        os.chdir(datadir)
        cls.parf = "B1855+09_NANOGrav_dfg+12_TAI.par"
        cls.timf = "B1855+09_NANOGrav_dfg+12.tim"
        cls.JUMPm = mb.get_model(cls.parf)
        cls.toas = toa.get_TOAs(
            cls.timf, ephem="DE405", planets=False, include_bipm=False
        )
        # libstempo calculation
        cls.ltres = np.genfromtxt(
            cls.parf + ".tempo_test", unpack=True, names=True, dtype=np.longdouble
        )

def test_astropy_ephem(self):
        t_astropy = toa.get_TOAs(self.tim, ephem='DE436t')
        t_ephem = toa.get_TOAs(self.tim, ephem='DE436t', tdb_method="ephemeris")
        diff = (t_astropy.table['tdbld']-t_ephem.table['tdbld'])*86400.0
        assert np.all(np.abs(diff) < 5e-9), "Test TDB method, 'astropy' vs " \
                                            "'ephemeris' failed."

# import matplotlib
# matplotlib.use('TKAgg')
import matplotlib.pyplot as plt

import astropy.units as u
import os

datadir = os.path.dirname(os.path.abspath(str(__file__)))
parfile = os.path.join(datadir, "NGC6440E.par")
timfile = os.path.join(datadir, "NGC6440E.tim")

# Define the timing model
m = mb.get_model(parfile)

# Read in the TOAs
t = pint.toa.get_TOAs(timfile)

# Examples of how to select some subsets of TOAs
# These can be un-done using t.unselect()
#
# Use every other TOA
# t.select(np.where(np.arange(t.ntoas) % 2))

# Use only TOAs with errors < 30 us
# t.select(t.get_errors() < 30 * u.us)

# Use only TOAs from the GBT (although this is all of them for this example)
# t.select(t.get_obss() == 'gbt')

# Print a summary of the TOAs that we have
t.print_summary()

def get_toas(evtfile, flags, tcoords=None, minweight=0, minMJD=0, maxMJD=100000):
    if evtfile[:-3] == "tim":
        usepickle = False
        if "usepickle" in flags:
            usepickle = flags["usepickle"]
        ts = toa.get_TOAs(evtfile, usepickle=usepickle)
        # Prune out of range MJDs
        mask = np.logical_or(
            ts.get_mjds() < minMJD * u.day, ts.get_mjds() > maxMJD * u.day
        )
        ts.table.remove_rows(mask)
        ts.table = ts.table.group_by("obs")
    else:
        if "usepickle" in flags and flags["usepickle"]:
            try:
                picklefile = toa._check_pickle(evtfile)
                if not picklefile:
                    picklefile = evtfile
                ts = toa.TOAs(picklefile)
                return ts
            except:
                pass
        weightcol = flags["weightcol"] if "weightcol" in flags else None
        target = tcoords if weightcol == "CALC" else None
        tl = fermi.load_Fermi_TOAs(
            evtfile, weightcolumn=weightcol, targetcoord=target, minweight=minweight
        )
        tl = filter(lambda t: (t.mjd.value > minMJD) and (t.mjd.value < maxMJD), tl)
        ts = toa.TOAs(toalist=tl)
        ts.filename = evtfile
        ts.compute_TDBs()
        ts.compute_posvels(ephem="DE421", planets=False)

plt.close()


datadir = os.path.dirname(os.path.abspath(str(__file__)))
parfile = os.path.join(datadir, "NGC6440E.par.good")
timfile = os.path.join(datadir, "NGC6440E.tim")
print(parfile)
print(timfile)
nwalkers = 50
nsteps = 2000

# Define the timing model
m = mb.get_model(parfile)

# Read in the TOAs
t = pint.toa.get_TOAs(timfile)

# Print a summary of the TOAs that we have
t.print_summary()

# These are pre-fit residuals
rs = pint.residuals.Residuals(t, m).phase_resids
xt = t.get_mjds()
plt.plot(xt, rs, "x")
plt.title("%s Pre-Fit Timing Residuals" % m.PSR.value)
plt.xlabel("MJD")
plt.ylabel("Residual (phase)")
plt.grid()
plt.show()

# Now do the fit
print("Fitting...")

import astropy.units as u
import os

redge = 3
ledge = 0.5
datadir = os.path.dirname(os.path.abspath(str(__file__)))
parfile = os.path.join(datadir, 'NGC6440E.par.orig')
timfile = os.path.join(datadir, 'NGC6440E.tim')

# Define the timing model
m = mb.get_model(parfile)

# Read in the TOAs
t = pint.toa.get_TOAs(timfile)
t0 = pint.toa.get_TOAs(timfile)

print(t.table['groups'])
t.get_highest_density_range(70)
# Examples of how to select some subsets of TOAs
# These can be un-done using t.unselect()
#
# Use every other TOA
# t.select(np.where(np.arange(t.ntoas) % 2))

# Use only TOAs with errors < 30 us
# t.select(t.get_errors() < 30 * u.us)

# Use only TOAs from the GBT (although this is all of them for this example)
# t.select(t.get_obss() == 'gbt')

name = 'testing'