docs: def calculate_splus_scross(nmax, mc, dl, F, e, t, l0, gamma, gammadot, inc):
      """
      Calculate splus and scross summed over all harmonics. 
      This waveform differs slightly from that in Taylor et al (2015) 
      in that it includes the time dependence of the advance of periastron.
      
      :param nmax: Total number of harmonics to use
      :param mc: Chirp mass of binary [Solar Mass]
      :param dl: Luminosity distance [Mpc]
      :param F: Orbital frequency of binary [Hz]
      :param e: Orbital Eccentricity
      :param t: TOAs [s]
      :param l0: Initial eccentric anomoly [rad]
      :param gamma: Angle of periastron advance [rad]
      :param gammadot: Time derivative of angle of periastron advance [rad/s]
      :param inc: Inclination angle [rad]
  
      """ 
      
      n = np.arange(1, nmax)
  
      # time dependent amplitudes
      an = get_an(n, mc, dl, F, e)
      bn = get_bn(n, mc, dl, F, e)
      cn = get_cn(n, mc, dl, F, e)
  
      # time dependent terms
      omega = 2*np.pi*F
      gt = gamma + gammadot * t
      lt = l0 + omega * t
  
      # tiled phase
      phase1 = n * np.tile(lt, (nmax-1,1)).T
      phase2 = np.tile(gt, (nmax-1,1)).T
      phasep = phase1 + 2*phase2
      phasem = phase1 - 2*phase2
  
      # intermediate terms
      sp = np.sin(phasem)/(n*omega-2*gammadot) + \
              np.sin(phasep)/(n*omega+2*gammadot)
      sm = np.sin(phasem)/(n*omega-2*gammadot) - \
              np.sin(phasep)/(n*omega+2*gammadot)
      cp = np.cos(phasem)/(n*omega-2*gammadot) + \
              np.cos(phasep)/(n*omega+2*gammadot)
      cm = np.cos(phasem)/(n*omega-2*gammadot) - \
              np.cos(phasep)/(n*omega+2*gammadot)
      
  
      splus_n = -0.5 * (1+np.cos(inc)**2) * (an*sp - bn*sm) + \
              (1-np.cos(inc)**2)*cn * np.sin(phase1)
      scross_n = np.cos(inc) * (an*cm - bn*cp)
          
  
      return np.sum(splus_n, axis=1), np.sum(scross_n, axis=1)