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
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