code: def add_glitch(psr, epoch, amp):
      
      
      # Define the heaviside function
      heaviside = lambda x: 0.5 * (N.sign(x) + 1)
  
      # Glitches are spontaneous spin-up events.
      # Thus TOAs will be advanced, and resiudals will be negative.
      
      psr.stoas[:] -= amp * heaviside(psr.toas() - epoch) * \
        (psr.toas() - epoch)*86400.0
docs: Like pulsar term BWM event, but now differently parameterized: just an
      amplitude (not log-amp) parameter, and an epoch. [source: piccard]
      
      :param psr: pulsar object
      :param epoch: TOA time (MJD) the burst hits the earth
      :param amp: amplitude of the glitch