def METHOD_NAME(self, toas, delays):
    times = 0
    wave_names = ["WAVE%d" % ii for ii in range(1, self.num_wave_terms + 1)]
    wave_terms = [getattr(self, name) for name in wave_names]
    wave_om = self.WAVE_OM.quantity
    base_phase = (
        wave_om
        * (
            toas.table["tdbld"] * u.day
            - self.WAVEEPOCH.value * u.day
            - delays.to(u.day)
        )
    ).value
    for k, wave_term in enumerate(wave_terms):
        wave_a, wave_b = wave_term.quantity
        METHOD_NAME = (k + 1) * base_phase
        times += wave_a * np.sin(METHOD_NAME)
        times += wave_b * np.cos(METHOD_NAME)
    return ((times) * self._parent.F0.quantity).to(u.dimensionless_unscaled)