diff --git a/CHANGELOG.md b/CHANGELOG.md
index d78616a0f..a2b2da74d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project, at least loosely, adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Changed
+- Changed units of Phase to be u.dimensionless_unscaled instead of u.cycle, which was confusing
 
 ## [0.6.3] - 2020-05-04
 ### Added
diff --git a/src/pint/__init__.py b/src/pint/__init__.py
index 6748ce513..b77b6410f 100644
--- a/src/pint/__init__.py
+++ b/src/pint/__init__.py
@@ -17,7 +17,6 @@
     "ls",
     "dmu",
     "light_second_equivalency",
-    "dimensionless_cycles",
     "hourangle_second",
     "pulsar_mjd",
     "GMsun",
@@ -50,7 +49,6 @@
 
 # define equivalency for astropy units
 light_second_equivalency = [(ls, si.second, lambda x: x, lambda x: x)]
-dimensionless_cycles = [(u.cycle, None)]
 # hourangle_second unit
 hourangle_second = u.def_unit("hourangle_second", u.hourangle / np.longdouble(3600.0))
 
@@ -84,7 +82,7 @@
     "Tsun": Tsun,
     "GMsun": GMsun,
     "MJD": u.day,
-    "pulse phase": u.cycle,
+    "pulse phase": u.dimensionless_unscaled,
     "hourangle_second": hourangle_second,
 }
 
diff --git a/src/pint/mcmc_fitter.py b/src/pint/mcmc_fitter.py
index 951e10ed6..89af3e370 100644
--- a/src/pint/mcmc_fitter.py
+++ b/src/pint/mcmc_fitter.py
@@ -291,7 +291,6 @@ def get_event_phases(self):
         """
         phases = self.model.phase(self.toas)[1]
         # ensure all positive
-        phases = phases.to(u.cycle).value
         return np.where(phases < 0.0, phases + 1.0, phases)
 
     def lnposterior(self, theta):
@@ -642,7 +641,6 @@ def get_event_phases(self, index=None):
             print("Showing all %d phases" % len(phases))
         else:
             phases = self.model.phase(self.toas_list[index])[1]
-        phases = phases.to(u.cycle).value
         return np.where(phases < 0.0, phases + 1.0, phases)
 
     def get_template_vals(self, phases, index):
diff --git a/src/pint/models/glitch.py b/src/pint/models/glitch.py
index 87341e9ef..0facf1e62 100644
--- a/src/pint/models/glitch.py
+++ b/src/pint/models/glitch.py
@@ -5,8 +5,8 @@
 
 import astropy.units as u
 import numpy as np
+from astropy import log
 
-from pint import dimensionless_cycles
 from pint.models.parameter import prefixParameter
 from pint.models.timing_model import MissingParameter, PhaseComponent
 from pint.utils import split_prefixed_name
@@ -155,43 +155,39 @@ def glitch_phase(self, toas, delay):
         returns an array of phases in long double
         """
         tbl = toas.table
-        phs = np.zeros_like(tbl, dtype=np.longdouble) * u.cycle
+        phs = u.Quantity(np.zeros_like(tbl, dtype=np.longdouble))
         glepnames = [x for x in self.params if x.startswith("GLEP_")]
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            for glepnm in glepnames:
-                glep = getattr(self, glepnm)
-                eph = glep.value
-                idx = glep.index
-                dphs = getattr(self, "GLPH_%d" % idx).quantity
-                dF0 = getattr(self, "GLF0_%d" % idx).quantity
-                dF1 = getattr(self, "GLF1_%d" % idx).quantity
-                dF2 = getattr(self, "GLF2_%d" % idx).quantity
-                dt = (tbl["tdbld"] - eph) * u.day - delay
-                dt = dt.to(u.second)
-                affected = dt > 0.0  # TOAs affected by glitch
-                # decay term
-                dF0D = getattr(self, "GLF0D_%d" % idx).quantity
-                if dF0D != 0.0:
-                    tau = getattr(self, "GLTD_%d" % idx).quantity
-                    decayterm = (
-                        dF0D
-                        * tau
-                        * (1.0 - np.exp(-(dt[affected] / tau).to(u.Unit(""))))
-                    )
-                else:
-                    decayterm = 0.0
+        for glepnm in glepnames:
+            glep = getattr(self, glepnm)
+            eph = glep.value
+            idx = glep.index
+            dphs = getattr(self, "GLPH_%d" % idx).quantity
+            dF0 = getattr(self, "GLF0_%d" % idx).quantity
+            dF1 = getattr(self, "GLF1_%d" % idx).quantity
+            dF2 = getattr(self, "GLF2_%d" % idx).quantity
+            dt = (tbl["tdbld"] - eph) * u.day - delay
+            dt = dt.to(u.second)
+            affected = dt > 0.0  # TOAs affected by glitch
+            # decay term
+            dF0D = getattr(self, "GLF0D_%d" % idx).quantity
+            if dF0D != 0.0:
+                tau = getattr(self, "GLTD_%d" % idx).quantity
+                decayterm = dF0D * tau * (1.0 - np.exp(-(dt[affected] / tau)))
+            else:
+                decayterm = 0.0 * u.Unit("")
 
-                phs[affected] += (
-                    dphs
-                    + dt[affected]
-                    * (
-                        dF0
-                        + 0.5 * dt[affected] * dF1
-                        + 1.0 / 6.0 * dt[affected] * dt[affected] * dF2
-                    )
-                    + decayterm
+            log.info("{} {} ".format(dphs, dphs.unit))
+            phs[affected] += (
+                dphs
+                + dt[affected]
+                * (
+                    dF0
+                    + 0.5 * dt[affected] * dF1
+                    + 1.0 / 6.0 * dt[affected] * dt[affected] * dF2
                 )
-            return phs.to(u.cycle)
+                + decayterm
+            )
+        return phs
 
     def d_phase_d_GLPH(self, toas, param, delay):
         """Calculate the derivative wrt GLPH"""
@@ -206,8 +202,8 @@ def d_phase_d_GLPH(self, toas, param, delay):
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLPH = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLPH.units
-        dpdGLPH[affected] += 1.0 * u.cycle / par_GLPH.units
+        dpdGLPH = np.zeros(len(tbl), dtype=np.longdouble) / par_GLPH.units
+        dpdGLPH[affected] += 1.0 / par_GLPH.units
         return dpdGLPH
 
     def d_phase_d_GLF0(self, toas, param, delay):
@@ -223,9 +219,8 @@ def d_phase_d_GLF0(self, toas, param, delay):
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLF0 = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLF0.units
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            dpdGLF0[affected] = dt[affected]
+        dpdGLF0 = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF0.units
+        dpdGLF0[affected] = dt[affected]
         return dpdGLF0
 
     def d_phase_d_GLF1(self, toas, param, delay):
@@ -241,9 +236,8 @@ def d_phase_d_GLF1(self, toas, param, delay):
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLF1 = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLF1.units
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            dpdGLF1[affected] += np.longdouble(0.5) * dt[affected] * dt[affected]
+        dpdGLF1 = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF1.units
+        dpdGLF1[affected] += np.longdouble(0.5) * dt[affected] * dt[affected]
         return dpdGLF1
 
     def d_phase_d_GLF2(self, toas, param, delay):
@@ -259,11 +253,10 @@ def d_phase_d_GLF2(self, toas, param, delay):
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLF2 = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLF2.units
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            dpdGLF2[affected] += (
-                np.longdouble(1.0) / 6.0 * dt[affected] * dt[affected] * dt[affected]
-            )
+        dpdGLF2 = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF2.units
+        dpdGLF2[affected] += (
+            np.longdouble(1.0) / 6.0 * dt[affected] * dt[affected] * dt[affected]
+        )
         return dpdGLF2
 
     def d_phase_d_GLF0D(self, toas, param, delay):
@@ -281,11 +274,8 @@ def d_phase_d_GLF0D(self, toas, param, delay):
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLF0D = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLF0D.units
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            dpdGLF0D[affected] += tau * (
-                np.longdouble(1.0) - np.exp(-dt[affected] / tau)
-            )
+        dpdGLF0D = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF0D.units
+        dpdGLF0D[affected] += tau * (np.longdouble(1.0) - np.exp(-dt[affected] / tau))
         return dpdGLF0D
 
     def d_phase_d_GLTD(self, toas, param, delay):
@@ -300,17 +290,14 @@ def d_phase_d_GLTD(self, toas, param, delay):
         eph = np.longdouble(getattr(self, "GLEP_" + ids).value)
         par_GLTD = getattr(self, param)
         if par_GLTD.value == 0.0:
-            return np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLTD.units
+            return np.zeros(len(tbl), dtype=np.longdouble) / par_GLTD.units
         glf0d = getattr(self, "GLF0D_" + ids).quantity
         tau = par_GLTD.quantity
         dt = (tbl["tdbld"] - eph) * u.day - delay
         dt = dt.to(u.second)
         affected = np.where(dt > 0.0)[0]
-        dpdGLTD = np.zeros(len(tbl), dtype=np.longdouble) * u.cycle / par_GLTD.units
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            dpdGLTD[affected] += glf0d * (
-                np.longdouble(1.0) - np.exp(-dt[affected] / tau)
-            ) + glf0d * tau * (-np.exp(-dt[affected] / tau)) * dt[affected] / (
-                tau * tau
-            )
+        dpdGLTD = np.zeros(len(tbl), dtype=np.longdouble) / par_GLTD.units
+        dpdGLTD[affected] += glf0d * (
+            np.longdouble(1.0) - np.exp(-dt[affected] / tau)
+        ) + glf0d * tau * (-np.exp(-dt[affected] / tau)) * dt[affected] / (tau * tau)
         return dpdGLTD
diff --git a/src/pint/models/ifunc.py b/src/pint/models/ifunc.py
index d7a103695..0cfd58f34 100644
--- a/src/pint/models/ifunc.py
+++ b/src/pint/models/ifunc.py
@@ -132,5 +132,5 @@ def ifunc_phase(self, toas, delays):
         else:
             raise ValueError("Interpolation type %d not supported.".format(itype))
 
-        phase = (times * u.s) * self.F0.quantity * u.cycle
+        phase = ((times * u.s) * self.F0.quantity).to(u.dimensionless_unscaled)
         return phase
diff --git a/src/pint/models/jump.py b/src/pint/models/jump.py
index a04d24983..b8d24f276 100644
--- a/src/pint/models/jump.py
+++ b/src/pint/models/jump.py
@@ -6,7 +6,6 @@
 import astropy.units as u
 import numpy
 
-from pint import dimensionless_cycles
 from pint.models.parameter import maskParameter
 from pint.models.timing_model import DelayComponent, MissingParameter, PhaseComponent
 
@@ -120,8 +119,7 @@ def d_phase_d_jump(self, toas, jump_param, delay):
         d_phase_d_j = numpy.zeros(len(tbl))
         mask = jpar.select_toa_mask(toas)
         d_phase_d_j[mask] = self.F0.value
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            return (d_phase_d_j * self.F0.units).to(u.cycle / u.second)
+        return (d_phase_d_j * self.F0.units).to(1 / u.second)
 
     def print_par(self):
         result = ""
diff --git a/src/pint/models/spindown.py b/src/pint/models/spindown.py
index 557ac2f3c..e8f8f1098 100644
--- a/src/pint/models/spindown.py
+++ b/src/pint/models/spindown.py
@@ -7,7 +7,6 @@
 import numpy
 
 import pint.toa as toa
-from pint import dimensionless_cycles
 from pint.models.parameter import MJDParameter, floatParameter, prefixParameter
 from pint.models.timing_model import MissingParameter, PhaseComponent
 from pint.pulsar_mjd import Time
@@ -127,10 +126,9 @@ def spindown_phase(self, toas, delay):
         """
         dt = self.get_dt(toas, delay)
         # Add the [0.0] because that is the constant phase term
-        fterms = [0.0 * u.cycle] + self.get_spin_terms()
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            phs = taylor_horner(dt.to(u.second), fterms)
-            return phs.to(u.cycle)
+        fterms = [0.0 * u.dimensionless_unscaled] + self.get_spin_terms()
+        phs = taylor_horner(dt.to(u.second), fterms)
+        return phs.to(u.dimensionless_unscaled)
 
     def change_pepoch(self, new_epoch, toas=None, delay=None):
         """Move PEPOCH to a new time and change the related paramters.
@@ -197,13 +195,11 @@ def d_phase_d_F(self, toas, param, delay):
         fterms = [ft * numpy.longdouble(0.0) / unit for ft in fterms]
         fterms[order] += numpy.longdouble(1.0)
         dt = self.get_dt(toas, delay)
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            d_pphs_d_f = taylor_horner(dt.to(u.second), fterms)
-            return d_pphs_d_f.to(u.cycle / unit)
+        d_pphs_d_f = taylor_horner(dt.to(u.second), fterms)
+        return d_pphs_d_f.to(1 / unit)
 
     def d_spindown_phase_d_delay(self, toas, delay):
         dt = self.get_dt(toas, delay)
         fterms = [0.0] + self.get_spin_terms()
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            d_pphs_d_delay = taylor_horner_deriv(dt.to(u.second), fterms)
-            return -d_pphs_d_delay.to(u.cycle / u.second)
+        d_pphs_d_delay = taylor_horner_deriv(dt.to(u.second), fterms)
+        return -d_pphs_d_delay.to(1 / u.second)
diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py
index 12835b77b..de037847e 100644
--- a/src/pint/models/timing_model.py
+++ b/src/pint/models/timing_model.py
@@ -16,7 +16,6 @@
 from astropy import log
 
 import pint
-from pint import dimensionless_cycles
 from pint.models.parameter import strParameter, maskParameter
 from pint.phase import Phase
 from pint.utils import PrefixError, interesting_lines, lines_of, split_prefixed_name
@@ -906,8 +905,7 @@ def d_phase_d_toa(self, toas, sample_step=None):
         dp = sample_phase[1] - sample_phase[0]
         d_phase_d_toa = dp.int / (2 * sample_step) + dp.frac / (2 * sample_step)
         del copy_toas
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            return d_phase_d_toa.to(u.Hz)
+        return d_phase_d_toa.to(u.Hz)
 
     def d_phase_d_tpulsar(self, toas):
         """Return the derivative of phase wrt time at the pulsar.
@@ -922,7 +920,7 @@ def d_phase_d_param(self, toas, delay, param):
         # Is it safe to assume that any param affecting delay only affects
         # phase indirectly (and vice-versa)??
         par = getattr(self, param)
-        result = np.longdouble(np.zeros(toas.ntoas)) * u.cycle / par.units
+        result = np.longdouble(np.zeros(toas.ntoas)) / par.units
         param_phase_derivs = []
         phase_derivs = self.phase_deriv_funcs
         delay_derivs = self.delay_deriv_funcs
@@ -940,7 +938,7 @@ def d_phase_d_param(self, toas, delay, param):
             #                         d_delay_d_param
 
             d_delay_d_p = self.d_delay_d_param(toas, param)
-            dpdd_result = np.longdouble(np.zeros(toas.ntoas)) * u.cycle / u.second
+            dpdd_result = np.longdouble(np.zeros(toas.ntoas)) / u.second
             for dpddf in self.d_phase_d_delay_funcs:
                 dpdd_result += dpddf(toas, delay)
             result = dpdd_result * d_delay_d_p
@@ -978,8 +976,12 @@ def d_phase_d_param_num(self, toas, param, step=1e-2):
             h = ori_value * step
         parv = [par.value - h, par.value + h]
 
-        phase_i = np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.cycle
-        phase_f = np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.cycle
+        phase_i = (
+            np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled
+        )
+        phase_f = (
+            np.zeros((toas.ntoas, 2), dtype=np.longdouble) * u.dimensionless_unscaled
+        )
         for ii, val in enumerate(parv):
             par.value = val
             ph = self.phase(toas)
diff --git a/src/pint/models/wave.py b/src/pint/models/wave.py
index 0ba61df23..d8fb707b5 100644
--- a/src/pint/models/wave.py
+++ b/src/pint/models/wave.py
@@ -111,5 +111,5 @@ def wave_phase(self, toas, delays):
             times += wave_a * np.sin(wave_phase)
             times += wave_b * np.cos(wave_phase)
 
-        phase = (times) * self.F0.quantity * u.cycle
+        phase = ((times) * self.F0.quantity).to(u.dimensionless_unscaled)
         return phase
diff --git a/src/pint/phase.py b/src/pint/phase.py
index ee537a661..6fbcfde17 100644
--- a/src/pint/phase.py
+++ b/src/pint/phase.py
@@ -1,11 +1,3 @@
-# phase.py
-# Simple class representing pulse phase as integer and fractional
-# parts.
-# SUGGESTION(@paulray): How about adding some documentation here
-# describing why the fractional part is reduced to [-0.5,0.5] instead of [0,1].
-# I think I understand it, but it would be good to have it stated.
-# Also, probably one of the comparisons below should be <= or >=, so the
-# range is [-0.5,0.5) or (-0.5,0.5].
 from __future__ import absolute_import, division, print_function
 
 from collections import namedtuple
@@ -13,50 +5,72 @@
 import astropy.units as u
 import numpy
 
-from pint import dimensionless_cycles
-
 
 class Phase(namedtuple("Phase", "int frac")):
     """
-    Phase class array version
+    Class representing pulse phase as integer (.int) and fractional (.frac) parts.
+
+    The phase values are dimensionless Quantitys (u.dimensionless_unscaled == u.Unit("") == Unit(dimensionless))
 
     Ensures that the fractional part stays in [-0.5, 0.5)
+
+    SUGGESTION(@paulray): How about adding some documentation here
+    describing why the fractional part is reduced to [-0.5,0.5) instead of [0,1).
+
     """
 
     __slots__ = ()
 
     def __new__(cls, arg1, arg2=None):
-        # Assume inputs are numerical, could add an extra
-        # case to parse strings as input.
-        # if it is not a list, convert to a list
+        """Create new Phase object
+
+        Constructs a Phase object. 
+        Can be initialized with arrays or a scalar Quantity.
+        Can take inputs as plain numerical types, or dimensionaless Quantitys
+        Accepts either floating point argument (arg1) or pair of arguments with integer (arg1) and fractional (arg2) parts separate
+        Scalars are converted to length 1 arrays so `Phase.int` and `Phase.frac` are always arrays
+
+        Parameters
+        ----------
+        arg1 : array or dimensionless Quantity
+        arg2 : array or dimensionless Quantity 
+
+        Returns
+        -------
+        Phase : pulse phase object with arrays of dimensionless Quantitys as the int and frac parts
+        """
         if not hasattr(arg1, "unit"):
-            arg1 = arg1 * u.cycle
+            arg1 = u.Quantity(arg1)
+        else:
+            # This will raise an exception if the argument has any unit not convertable to Unit(dimensionless)
+            arg1 = arg1.to(u.dimensionless_unscaled)
+        #  If arg is scalar, convert to an array of length 1
         if arg1.shape == ():
             arg1 = arg1.reshape((1,))
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            arg1 = arg1.to(u.Unit(""))
-            # Since modf does not like dimensioned quantity
-            if arg2 is None:
-                ff, ii = numpy.modf(arg1)
+        if arg2 is None:
+            ff, ii = numpy.modf(arg1)
+        else:
+            if not hasattr(arg2, "unit"):
+                arg2 = u.Quantity(arg2)
             else:
-                if not hasattr(arg2, "unit"):
-                    arg2 = arg2 * u.cycle
-                if arg2.shape == ():
-                    arg2 = arg2.reshape((1,))
-                arg2 = arg2.to(u.Unit(""))
-                arg1S = numpy.modf(arg1)
-                arg2S = numpy.modf(arg2)
-                ii = arg1S[1] + arg2S[1]
-                ff = arg2S[0]
-            index = numpy.where(ff < -0.5)
-            ff[index] += 1.0
-            ii[index] -= 1
-            # The next line is >= so that the range is the interval [-0.5,0.5)
-            # Otherwise, the same phase could be represented 0,0.5 or 1,-0.5
-            index = numpy.where(ff >= 0.5)
-            ff[index] -= 1.0
-            ii[index] += 1
-            return super(Phase, cls).__new__(cls, ii.to(u.cycle), ff.to(u.cycle))
+                arg2 = arg2.to(u.dimensionless_unscaled)
+            if arg2.shape == ():
+                arg2 = arg2.reshape((1,))
+            arg1S = numpy.modf(arg1)
+            arg2S = numpy.modf(arg2)
+            # Prior code assumed that fractional part of arg1 was 0 if arg2 was present
+            # @paulray removed that assumption here
+            ff = arg1S[0] + arg2S[0]
+            ii = arg1S[1] + arg2S[1]
+        index = ff < -0.5
+        ff[index] += 1.0
+        ii[index] -= 1
+        # The next line is >= so that the range is the interval [-0.5,0.5)
+        # Otherwise, the same phase could be represented 0,0.5 or 1,-0.5
+        index = ff >= 0.5
+        ff[index] -= 1.0
+        ii[index] += 1
+        return super(Phase, cls).__new__(cls, ii, ff)
 
     def __neg__(self):
         # TODO: add type check for __neg__ and __add__
@@ -64,10 +78,8 @@ def __neg__(self):
 
     def __add__(self, other):
         ff = self.frac + other.frac
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            ii = numpy.modf(ff.to(u.Unit("")))[1]
-            ii = ii.to(u.cycle)
-            return Phase(self.int + other.int + ii, ff - ii)
+        ii = numpy.modf(ff)[1]
+        return Phase(self.int + other.int + ii, ff - ii)
 
     def __sub__(self, other):
         return self.__add__(other.__neg__())
diff --git a/src/pint/pintk/pulsar.py b/src/pint/pintk/pulsar.py
index 29d8ce9bf..a37a08113 100644
--- a/src/pint/pintk/pulsar.py
+++ b/src/pint/pintk/pulsar.py
@@ -168,7 +168,7 @@ def orbitalphase(self):
 
         phase = np.modf(tpb)[0]
         phase[phase < 0] += 1
-        return phase * u.cycle
+        return phase
 
     def dayofyear(self):
         """
diff --git a/src/pint/polycos.py b/src/pint/polycos.py
index 160dc70b4..eb5fd603a 100644
--- a/src/pint/polycos.py
+++ b/src/pint/polycos.py
@@ -615,9 +615,7 @@ def generate_polycos(
                 dt = (nodes * u.day - tmid).to("min")  # Use constant
                 rdcPhase = ph - refPhase
                 rdcPhase = (
-                    rdcPhase.int
-                    - (dt.value * model.F0.value * 60.0) * u.cycle
-                    + rdcPhase.frac
+                    rdcPhase.int - (dt.value * model.F0.value * 60.0) + rdcPhase.frac
                 )
                 dtd = dt.value.astype(float)  # Truncate to double
                 rdcPhased = rdcPhase.astype(float)
@@ -791,8 +789,8 @@ def eval_abs_phase(self, t):
             phaseInt += (absp.int,)
             phaseFrac += (absp.frac,)
             # Maybe add sort function here, since the time has been masked.
-        phaseInt = np.hstack(phaseInt).value * u.cycle
-        phaseFrac = np.hstack(phaseFrac).value * u.cycle
+        phaseInt = np.hstack(phaseInt).value
+        phaseFrac = np.hstack(phaseFrac).value
         absPhase = Phase(phaseInt, phaseFrac)
 
         return absPhase
diff --git a/src/pint/residuals.py b/src/pint/residuals.py
index 543101c8e..cac0222e9 100644
--- a/src/pint/residuals.py
+++ b/src/pint/residuals.py
@@ -5,7 +5,6 @@
 from scipy.linalg import LinAlgError
 from astropy import log
 
-from pint import dimensionless_cycles
 from pint.phase import Phase
 from pint.utils import weighted_mean
 
@@ -130,8 +129,7 @@ def calc_time_resids(self, weighted_mean=True):
         """Return timing model residuals in time (seconds)."""
         if self.phase_resids is None:
             self.phase_resids = self.calc_phase_resids(weighted_mean=weighted_mean)
-        with u.set_enabled_equivalencies(dimensionless_cycles):
-            return (self.phase_resids.to(u.Unit("")) / self.get_PSR_freq()).to(u.s)
+        return (self.phase_resids / self.get_PSR_freq()).to(u.s)
 
     def get_PSR_freq(self, modelF0=True):
         if modelF0:
diff --git a/src/pint/scripts/event_optimize.py b/src/pint/scripts/event_optimize.py
index 53635f4c7..e86439daa 100755
--- a/src/pint/scripts/event_optimize.py
+++ b/src/pint/scripts/event_optimize.py
@@ -279,7 +279,6 @@ def get_event_phases(self):
         """
         phss = self.model.phase(self.toas)[1]
         # ensure all postive
-        phss = phss.to(u.cycle).value
         return np.where(phss < 0.0, phss + 1.0, phss)
 
     def lnprior(self, theta):
diff --git a/src/pint/scripts/event_optimize_multiple.py b/src/pint/scripts/event_optimize_multiple.py
index 56519619d..740e23552 100755
--- a/src/pint/scripts/event_optimize_multiple.py
+++ b/src/pint/scripts/event_optimize_multiple.py
@@ -322,7 +322,7 @@ def main(argv=None):
                 gtemplate = cPickle.load(file(tname))
             except:
                 phases = (modelin.phase(ts)[1]).astype(np.float64)
-                phases[phases < 0] += 1 * u.cycle
+                phases[phases < 0] += 1 * u.dimensionless_unscaled
                 gtemplate = lctemplate.get_gauss2()
                 lcf = lcfitters.LCFitter(gtemplate, phases, weights=wlist[i])
                 lcf.fit(unbinned=False)
@@ -332,7 +332,7 @@ def main(argv=None):
                     protocol=2,
                 )
             phases = (modelin.phase(ts)[1]).astype(np.float64)
-            phases[phases < 0] += 1 * u.cycle
+            phases[phases < 0] += 1 * u.dimensionless_unscaled
             lcf = lcfitters.LCFitter(
                 gtemplate, phases.value, weights=wlist[i], binned_bins=200
             )
diff --git a/src/pint/scripts/fermiphase.py b/src/pint/scripts/fermiphase.py
index 798d4c754..87cbf1ad5 100755
--- a/src/pint/scripts/fermiphase.py
+++ b/src/pint/scripts/fermiphase.py
@@ -117,7 +117,6 @@ def main(argv=None):
     # Compute model phase for each TOA
     iphss, phss = modelin.phase(ts, abs_phase=True)
     # ensure all postive
-    phss = phss.to(u.cycle).value
     phases = np.where(phss < 0.0, phss + 1.0, phss)
     mjds = ts.get_mjds()
     weights = np.array([w["weight"] for w in ts.table["flags"]])
diff --git a/src/pint/scripts/photonphase.py b/src/pint/scripts/photonphase.py
index e34a644db..e0a922955 100755
--- a/src/pint/scripts/photonphase.py
+++ b/src/pint/scripts/photonphase.py
@@ -213,7 +213,6 @@ def main(argv=None):
     # Compute model phase for each TOA
     iphss, phss = modelin.phase(ts, abs_phase=True)
     # ensure all postive
-    phss = phss.to(u.cycle).value
     negmask = phss < 0.0
     phases = np.where(negmask, phss + 1.0, phss)
     h = float(hm(phases))
@@ -251,7 +250,7 @@ def main(argv=None):
             data_to_add["PULSE_PHASE"] = [phases, "D"]
 
         if args.absphase:
-            data_to_add["ABS_PHASE"] = [iphss - negmask * u.cycle, "K"]
+            data_to_add["ABS_PHASE"] = [iphss - negmask, "K"]
 
         if args.barytime:
             bats = modelin.get_barycentric_toas(ts)
diff --git a/src/pint/scripts/pintempo.py b/src/pint/scripts/pintempo.py
index 0d3212e52..101ec95cd 100755
--- a/src/pint/scripts/pintempo.py
+++ b/src/pint/scripts/pintempo.py
@@ -55,10 +55,11 @@ def main(argv=None):
     f = pint.fitter.WLSFitter(t, m)
     f.fit_toas()
 
-    # Print some basic params
-    print("Best fit has reduced chi^2 of", f.resids.chi2_reduced)
-    print("RMS in phase is", f.resids.phase_resids.std())
-    print("RMS in time is", f.resids.time_resids.std().to(u.us))
+    # Print fit summary
+    print(
+        "============================================================================"
+    )
+    f.print_summary()
 
     if args.plot:
         import matplotlib.pyplot as plt
diff --git a/src/pint/scripts/zima.py b/src/pint/scripts/zima.py
index 926cba74e..844776070 100755
--- a/src/pint/scripts/zima.py
+++ b/src/pint/scripts/zima.py
@@ -135,16 +135,13 @@ def main(argv=None):
         log.info("Computing observatory positions and velocities.")
         ts.compute_posvels(args.ephem, args.planets)
 
-    # F_local has units of Hz; discard cycles unit in phase to get a unit
-    # that TimeDelta understands
     log.info("Creating TOAs")
     F_local = m.d_phase_d_toa(ts)
-    rs = m.phase(ts).frac.value / F_local
+    rs = m.phase(ts).frac / F_local
 
     # Adjust the TOA times to put them where their residuals will be 0.0
-
     ts.adjust_TOAs(TimeDelta(-1.0 * rs))
-    rspost = m.phase(ts).frac.value / F_local
+    rspost = m.phase(ts).frac / F_local
 
     log.info("Second iteration")
     # Do a second iteration
@@ -160,24 +157,21 @@ def main(argv=None):
     if args.plot:
         # This should be a very boring plot with all residuals flat at 0.0!
         import matplotlib.pyplot as plt
-
-        with u.set_enabled_equivalencies(u.dimensionless_angles()):
-            rspost2 = m.phase(ts).frac / F_local
-            plt.errorbar(
-                ts.get_mjds().value,
-                rspost2.to(u.us).value,
-                yerr=ts.get_errors().to(u.us).value,
-            )
-            newts = pint.toa.get_TOAs(
-                args.timfile, ephem=args.ephem, planets=args.planets
-            )
-            rsnew = m.phase(newts).frac / F_local
-            plt.errorbar(
-                newts.get_mjds().value,
-                rsnew.to(u.us).value,
-                yerr=newts.get_errors().to(u.us).value,
-            )
-            # plt.plot(ts.get_mjds(),rspost.to(u.us),'x')
-            plt.xlabel("MJD")
-            plt.ylabel("Residual (us)")
-            plt.show()
+        from astropy.visualization import quantity_support
+
+        quantity_support()
+
+        rspost2 = m.phase(ts).frac / F_local
+        plt.errorbar(
+            ts.get_mjds(), rspost2.to(u.us), yerr=ts.get_errors().to(u.us), fmt="."
+        )
+        newts = pint.toa.get_TOAs(args.timfile, ephem=args.ephem, planets=args.planets)
+        rsnew = m.phase(newts).frac / F_local
+        plt.errorbar(
+            newts.get_mjds(), rsnew.to(u.us), yerr=newts.get_errors().to(u.us), fmt="."
+        )
+        # plt.plot(ts.get_mjds(),rspost.to(u.us),'x')
+        plt.xlabel("MJD")
+        plt.ylabel("Residual (us)")
+        plt.grid(True)
+        plt.show()
diff --git a/src/pint/toa.py b/src/pint/toa.py
index edc4b3e8e..68adfe475 100644
--- a/src/pint/toa.py
+++ b/src/pint/toa.py
@@ -749,7 +749,7 @@ def __init__(self, toafile=None, toalist=None):
                     self.get_freqs(),
                     self.get_obss(),
                     self.get_flags(),
-                    np.zeros(len(mjds)) * u.cycle,
+                    np.zeros(len(mjds)),
                     self.get_groups(),
                 ],
                 names=(
@@ -854,7 +854,7 @@ def get_pulse_numbers(self):
         # TODO: use a masked array?  Only some pulse numbers may be known
         if hasattr(self, "toas"):
             try:
-                return np.array([t.flags["pn"] for t in self.toas]) * u.cycle
+                return np.array([t.flags["pn"] for t in self.toas])
             except KeyError:
                 log.warning("Not all TOAs have pulse numbers, using none")
                 return None
@@ -864,7 +864,7 @@ def get_pulse_numbers(self):
                     raise ValueError(
                         "Pulse number cannot be both a column and a TOA flag"
                     )
-                return np.array(flags["pn"] for flags in self.table["flags"]) * u.cycle
+                return np.array(flags["pn"] for flags in self.table["flags"])
             elif "pulse_number" in self.table.colnames:
                 return self.table["pulse_number"]
             else:
@@ -1020,7 +1020,7 @@ def phase_columns_from_flags(self):
         try:
             pns = [flags["pn"] for flags in self.table["flags"]]
             self.table["pulse_number"] = pns
-            self.table["pulse_number"].unit = u.cycle
+            self.table["pulse_number"].unit = u.dimensionless_unscaled
 
             # Remove pn from dictionary to prevent redundancies
             for flags in self.table["flags"]:
@@ -1046,7 +1046,7 @@ def compute_pulse_numbers(self, model):
         # paulr: I think pulse numbers should be computed with abs_phase=True!
         phases = model.phase(self, abs_phase=True)
         self.table["pulse_number"] = phases.int
-        self.table["pulse_number"].unit = u.cycle
+        self.table["pulse_number"].unit = u.dimensionless_unscaled
 
     def adjust_TOAs(self, delta):
         """Apply a time delta to TOAs
diff --git a/tests/test_pintempo.py b/tests/test_pintempo.py
index 236b0ec5c..57473cb5b 100644
--- a/tests/test_pintempo.py
+++ b/tests/test_pintempo.py
@@ -19,12 +19,11 @@ def test_result():
     pintempo.main(cmd.split())
     lines = sys.stdout.getvalue()
     v = 999.0
+    # This line is in the output:
+    # Prefit residuals 1090.580262221985 us, Postfit residuals 21.182038051610704 us
     for l in lines.split("\n"):
-        if l.startswith("RMS in time is"):
-            v = float(l.split()[4])
-    # Check that RMS is less than 34 microseconds
-    from astropy import log
-
-    log.warning("%f" % v)
-    assert v < 34.0
+        if l.startswith("Prefit residuals"):
+            v = float(l.split()[6])
+    # Check that RMS is less than 30 microseconds
+    assert v < 30.0
     sys.stdout = saved_stdout