code: def get_an(n, mc, dl, F, e):
      
      
      # convert to seconds
      mc *= SOLAR2S
      dl *= MPC2S
      
      omega = 2 * np.pi * F
      
      amp = n * mc**(5/3) * omega**(2/3) / dl
      
      ret = -amp * (ss.jn(n-2,n*e) - 2*e*ss.jn(n-1,n*e) +
                    (2/n)*ss.jn(n,n*e) + 2*e*ss.jn(n+1,n*e) -
                    ss.jn(n+2,n*e))
  
      return ret
docs: Compute a_n from Eq. 22 of Taylor et al. (2015).
      
      :param n: Harmonic number
      :param mc: Chirp mass of binary [Solar Mass]
      :param dl: Luminosity distance [Mpc]
      :param F: Orbital frequency of binary [Hz]
      :param e: Orbital Eccentricity
      
      :returns: a_n