docs: def add_cgw(psr, gwtheta, gwphi, mc, dist, fgw, phase0, psi, inc, pdist=1.0, \
                          pphase=None, psrTerm=True, evolve=True, \
                          phase_approx=False, tref=0):
      """
      Function to create GW-induced residuals from a SMBMB as 
      defined in Ellis et. al 2012,2013. Tries to be smart about it...
  
      :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 fgw: Frequency of GW (twice the orbital frequency) [Hz]
      :param phase0: Initial Phase of GW source [radians]
      :param psi: Polarization of GW source [radians]
      :param inc: Inclination of GW source [radians]
      :param pdist: Pulsar distance to use other than those in psr [kpc]
      :param pphase: Use pulsar phase to determine distance [radian]
      :param psrTerm: Option to include pulsar term [boolean] 
      :param evolve: Option to exclude evolution [boolean]
      :param tref: Fidicuial time at which initial parameters are referenced
  
      :returns: Vector of induced residuals
      """
  
      # convert units
      mc *= eu.SOLAR2S         # convert from solar masses to seconds
      dist *= eu.MPC2S    # convert from Mpc to seconds
  
      # define initial orbital frequency 
      w0 = N.pi * fgw
      phase0 /= 2 # orbital phase
      w053 = w0**(-5/3)
  
      # 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)
      incfac1, incfac2 = 0.5*(3+N.cos(2*inc)), 2*N.cos(inc)
  
      # 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])
  
      # various factors invloving GW parameters
      fac1 = 256/5 * mc**(5/3) * w0**(8/3) 
      fac2 = 1/32/mc**(5/3)
      fac3 = mc**(5/3)/dist
  
      # 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 = psr.toas()*86400 - tref
      if pphase is not None:
          pd = pphase/(2*N.pi*fgw*(1-cosMu)) / eu.KPC2S
      else:
          pd = pdist
  
      # convert units
      pd *= eu.KPC2S   # convert from kpc to seconds
      
      # get pulsar time
      tp = toas-pd*(1-cosMu)
  
      # evolution
      if evolve:
  
          # calculate time dependent frequency at earth and pulsar
          omega = w0 * (1 - fac1 * toas)**(-3/8)
          omega_p = w0 * (1 - fac1 * tp)**(-3/8)
  
          # calculate time dependent phase
          phase = phase0 + fac2 * (w053 - omega**(-5/3))
          phase_p = phase0 + fac2 * (w053 - omega_p**(-5/3))
      
      # use approximation that frequency does not evlolve over observation time
      elif phase_approx:
          
          # frequencies
          omega = w0
          omega_p = w0 * (1 + fac1 * pd*(1-cosMu))**(-3/8)
          
          # phases
          phase = phase0 + omega * toas
          phase_p = phase0 + fac2 * (w053 - omega_p**(-5/3)) + omega_p*toas
            
      # no evolution
      else: 
          
          # monochromatic
          omega = w0
          omega_p = omega
          
          # phases
          phase = phase0 + omega * toas
          phase_p = phase0 + omega * tp
          
  
      # define time dependent coefficients
      At = N.sin(2*phase) * incfac1
      Bt = N.cos(2*phase) * incfac2
      At_p = N.sin(2*phase_p) * incfac1
      Bt_p = N.cos(2*phase_p) * incfac2
  
      # now define time dependent amplitudes
      alpha = fac3 / omega**(1/3)
      alpha_p = fac3 / omega_p**(1/3)
  
      # define rplus and rcross
      rplus = alpha * (At*cos2psi + Bt*sin2psi)
      rcross = alpha * (-At*sin2psi + Bt*cos2psi)
      rplus_p = alpha_p * (At_p*cos2psi + Bt_p*sin2psi)
      rcross_p = alpha_p * (-At_p*sin2psi + Bt_p*cos2psi)
  
      # residuals
      if psrTerm:
          res = fplus*(rplus_p-rplus)+fcross*(rcross_p-rcross)
      else:
          res = -fplus*rplus - fcross*rcross
  
      psr.stoas[:] += res/86400