def f_multiple_sin(x, *params):
    y = np.zeros_like(x)
    for i in range(0, len(params), 4): #4=amplitude, freq, phase, offset
        amplitude = params[i]
        freq = params[i+1]
        phase = params[i+2]
        offset = params[i+3]
        y = y + amplitude*np.sin(np.multiply(freq, x)+phase)+offset
    return y

def fit_multiple(t, flux, flux_err, guess):
    popt, pcov = optimize.curve_fit(
        f_multiple_sin, t, flux, sigma=flux_err, p0=guess,
        bounds=(guess-np.multiply(guess,0.1),guess+np.multiply(guess,0.1)),
        absolute_sigma=True
        )

    perr = np.sqrt(np.diag(pcov))

    return popt, perr

guess = [4.50148944e-03, 2.40000040e+02, 3.01766641e-03, 8.99996136e-01, 3.14546648e-03, 2.61818207e+02, 2.94282247e-03, 5.56770657e-06]
popt, perr = fit_multiple(time, flux, flux_err, guess)