def lomb_scargle_model(time, signal, error, sys_err=0.05, nharm=8, nfreq=3, tone_control=5.0):

    dy0 = np.sqrt(error**2 + sys_err**2)

    wt = 1. / dy0**2
    time = time.copy() - min(time) # speeds up lomb_scargle code to have min(time)==0
    signal = signal.copy()

    chi0 = np.dot(signal**2, wt)

# TODO parametrize?
    f0 = 1. / max(time)
    df = 0.8 / max(time) # 20120202 :    0.1/Xmax
    fmax = 33. #pre 20120126: 10. # 25
    numf = int((fmax - f0) / df) # TODO !!! this is off by 1 point, fix?

    model_dict = {'freq_fits' : []}
    lambda0_range = [-np.log10(len(time)), 8] # these numbers "fix" the strange-amplitude effect
    for i in range(nfreq):
        if i == 0:
            fit = fit_lomb_scargle(time, signal, dy0, f0, df, numf,
                    tone_control=tone_control, lambda0_range=lambda0_range,
                    nharm=nharm, detrend_order=1)
            model_dict['trend'] = fit['trend_coef'][1]
        else:
            fit = fit_lomb_scargle(time, signal, dy0, f0, df, numf,
                    tone_control=tone_control, lambda0_range=lambda0_range,
                    nharm=nharm, detrend_order=0)
        model_dict['freq_fits'].append(fit)
        signal -= fit['model']
        model_dict['freq_fits'][-1]['resid'] = signal.copy()
        if i == 0:
            model_dict['varrat'] = np.dot(signal**2, wt) / chi0

    model_dict['nfreq'] = nfreq
    model_dict['nharm'] = nharm
    model_dict['chi2'] = fit['chi2']
    model_dict['f0'] = f0
    model_dict['df'] = df
    model_dict['numf'] = numf

    return model_dict
Simultaneous fit of a sum of sinusoids by weighted least squares:
dataset/ETHPy150Open cesium-ml/cesium/cesium/science_features/lomb_scargle.py/lomb_scargle_model