docs: def add_ecc_cgw(psr, gwtheta, gwphi, mc, dist, F, inc, psi, gamma0,
                  e0, l0, q, nmax=100, nset=None, pd=None, periEv=True,
                  psrTerm=True, tref=0, check=True, useFile=True):
      """
      Simulate GW from eccentric SMBHB. Waveform models from
      Taylor et al. (2015) and Barack and Cutler (2004).
  
      WARNING: This residual waveform is only accurate if the
      GW frequency is not significantly evolving over the 
      observation time of the pulsar.
  
      :param psr: pulsar object
      :param gwtheta: Polar angle of GW source in celestial coords [radians]
      :param gwphi: Azimuthal angle of GW source in celestial coords [radians]
      :param mc: Chirp mass of SMBMB [solar masses]
      :param dist: Luminosity distance to SMBMB [Mpc]
      :param F: Orbital frequency of SMBHB [Hz]
      :param inc: Inclination of GW source [radians]
      :param psi: Polarization of GW source [radians]
      :param gamma0: Initial angle of periastron [radians]
      :param e0: Initial eccentricity of SMBHB
      :param l0: Initial mean anomaly [radians]
      :param q: Mass ratio of SMBHB
      :param nmax: Number of harmonics to use in waveform decomposition
      :param nset: Fix the number of harmonics to be injected
      :param pd: Pulsar distance [kpc]
      :param periEv: Evolve the position of periapsis [boolean] 
      :param psrTerm: Option to include pulsar term [boolean] 
      :param tref: Fiducial time at which initial parameters are referenced [s]
      :param check: Check if frequency evolves significantly over obs. time
      :param useFile: Use pre-computed table of number of harmonics vs eccentricity
  
      :returns: Vector of induced residuals
      """
      
      # define variable for later use
      cosgwtheta, cosgwphi = N.cos(gwtheta), N.cos(gwphi)
      singwtheta, singwphi = N.sin(gwtheta), N.sin(gwphi)
      sin2psi, cos2psi = N.sin(2*psi), N.cos(2*psi)
  
      # unit vectors to GW source
      m = N.array([singwphi, -cosgwphi, 0.0])
      n = N.array([-cosgwtheta*cosgwphi, -cosgwtheta*singwphi, singwtheta])
      omhat = N.array([-singwtheta*cosgwphi, -singwtheta*singwphi, -cosgwtheta])
      
      # pulsar location
      if 'RAJ' and 'DECJ' in psr.pars():
          ptheta = N.pi/2 - psr['DECJ'].val
          pphi = psr['RAJ'].val
      elif 'ELONG' and 'ELAT' in psr.pars():
          fac = 180./N.pi
          coords = ephem.Equatorial(ephem.Ecliptic(str(psr['ELONG'].val*fac), str(psr['ELAT'].val*fac)))
  
          ptheta = N.pi/2 - float(repr(coords.dec))
          pphi = float(repr(coords.ra))
  
  
      # use definition from Sesana et al 2010 and Ellis et al 2012
      phat = N.array([N.sin(ptheta)*N.cos(pphi), N.sin(ptheta)*N.sin(pphi),\
              N.cos(ptheta)])
  
      fplus = 0.5 * (N.dot(m, phat)**2 - N.dot(n, phat)**2) / (1+N.dot(omhat, phat))
      fcross = (N.dot(m, phat)*N.dot(n, phat)) / (1 + N.dot(omhat, phat))
      cosMu = -N.dot(omhat, phat)
  
      # get values from pulsar object
      toas = N.double(psr.toas())*86400 - tref
      
      if check:
          # check that frequency is not evolving significantly over obs. time
          y = eu.solve_coupled_ecc_solution(F, e0, gamma0, l0, mc, q,
                                            N.array([0.0,toas.max()]))
          
          # initial and final values over observation time
          Fc0, ec0, gc0, phic0 = y[0,:]
          Fc1, ec1, gc1, phic1 = y[-1,:]
  
          # observation time
          Tobs = 1/(toas.max()-toas.min())
  
          if N.abs(Fc0-Fc1) > 1/Tobs:
              print('WARNING: Frequency is evolving over more than one frequency bin.')
              print('F0 = {0}, F1 = {1}, delta f = {2}'.format(Fc0, Fc1, 1/Tobs))
      
      # get gammadot for earth term
      if periEv==False:
          gammadot = 0.0
      else:
          gammadot = eu.get_gammadot(F, mc, q, e0)
  
      if nset is not None:
          nharm = nset
      elif useFile:
          if e0 > 0.001 and e0 < 0.999:
              nharm = min(int(ecc_interp(e0)), nmax) + 1
          elif e0 < 0.001:
              nharm = 3
          else:
              nharm = nmax
      else:
          nharm = nmax
      
      ##### earth term #####
      splus, scross = eu.calculate_splus_scross(nharm, mc, dist, F, e0, toas,
                                                l0, gamma0, gammadot, inc)
      
      ##### pulsar term #####
      if psrTerm:
  
          # convert units
          pd *= eu.KPC2S   # convert from kpc to seconds
      
          # get pulsar time
          tp = toas - pd * (1-cosMu)
  
          # solve coupled system of equations to get pulsar term values
          y = eu.solve_coupled_ecc_solution(F, e0, gamma0, l0, mc, q,
                                            N.array([0.0, tp.min()]))
          
          # get pulsar term values
          if N.any(y):
              Fp, ep, gp, lp = y[-1,:]
  
              # get gammadot at pulsar term
              gammadotp = eu.get_gammadot(Fp, mc, q, ep)
  
              if useFile:
                  if ep > 0.001 and ep < 0.999:
                      nharm = min(int(ecc_interp(ep)), nmax)
                  elif ep < 0.001:
                      nharm = 3
                  else:
                      nharm = nmax
              else:
                  nharm = nmax
  
  
              splusp, scrossp = eu.calculate_splus_scross(nharm, mc, dist, Fp,
                                                          ep, toas, lp, gp,
                                                          gammadotp, inc)
  
              rr = (fplus*cos2psi - fcross*sin2psi) * (splusp - splus) + \
                      (fplus*sin2psi + fcross*cos2psi) * (scrossp - scross)
  
          else:
              rr = N.zeros(len(p.toas))
              
      else:
          rr = - (fplus*cos2psi - fcross*sin2psi) * splus - \
                  (fplus*sin2psi + fcross*cos2psi) * scross
           
      psr.stoas[:] += rr/86400