def period_folding(x, y, dy, lomb_model, sys_err=0.05):
    out_dict = {}
    model_vals = np.zeros(len(y))
    freq_2p = lomb_model['freq_fits'][0]['freq'] * 0.5
    ytest_2p = y.copy()
    dy0 = np.sqrt(dy**2 + sys_err**2)

    # Here we force the freq to just 2*freq1_Period; we also do not use linear
    # detrending since we are not searching for freqs, and we want the
    # resulting model to be smooth when in phase-space. Detrending would result
    # in non-smooth model when period folded
    lambda0_range = [-np.log10(len(x)), 8.]
    fit = ls.fit_lomb_scargle(x, ytest_2p, dy0, freq_2p, lomb_model['df'], 1,
            lambda0_range=lambda0_range, nharm=lomb_model['nharm'], detrend_order=0)
    model_vals += fit['model']

    ytest_2p -= fit['model']
    for i in range(1, lomb_model['nfreq']):
        fit = ls.fit_lomb_scargle(x, ytest_2p, dy0, lomb_model['f0'],
                lomb_model['df'], lomb_model['numf'], 
                lambda0_range=lambda0_range, nharm=lomb_model['nharm'],
                detrend_order=0)
        ytest_2p -= fit['model']

    out_dict['1p_resid'] = lomb_model['freq_fits'][-1]['resid']
    out_dict['2p_resid'] = ytest_2p

    # So the following uses the 2*Period model, and gets a time-sorted, folded t and m:
    # NOTE: if this is succesful, I think a lot of other features could characterize the
    # shapes of the 2P folded data (not P or 2P dependent).
    # the reason we choose 2P is that occasionally for eclipsing
    # sources the LS code chooses 0.5 of true period (but never 2x
    # the true period).  slopes are not dependent upon the actual
    # period so 2P is fine if it gives a high chance of correct fitting.
    # NOTE: we only use the model from freq1 because this with its harmonics seems to
    # adequately model shapes such as RRLyr skewed sawtooth, multi minima of rvtau
    # without getting the scatter from using additional LS found frequencies.
    t_2per_fold = np.array(x % (1. / freq_2p))
    t_2per_sort_inds = np.argsort(t_2per_fold)
    t_2per_fold = t_2per_fold[t_2per_sort_inds]
    y_2per_fold = np.array(model_vals)[t_2per_sort_inds]
    out_dict['folded_slopes'] = slopes = np.diff(y_2per_fold) / np.diff(t_2per_fold)

    return out_dict
This section is used to calculate Dubath (10. Percentile90:2P/P),
dataset/ETHPy150Open cesium-ml/cesium/cesium/science_features/period_folding.py/period_folding