plot(time, m.resid, seriestype=:scatter)
plot!([0], seriestype = "hline")

using LsqFit

# a two-parameter exponential model
# x: array of independent variables
# p: array of model parameters
# model(x, p) will accept the full data set as the first argument `x`.
# This means that we need to write our model function so it applies
# the model to the full dataset. We use `@.` to apply the calculations
# across all rows.
@. model(x, p) = p[1]*exp(-x*p[2])

# some example data
# xdata: independent variables
# ydata: dependent variable
xdata = range(0, stop=10, length=20)
ydata = model(xdata, [1.0 2.0]) + 0.01*randn(length(xdata))
p0 = [0.5, 0.5]

fit = curve_fit(model, xdata, ydata, p0)

prediction = model.(xdata, Ref(fit.param))

using Plots 
scatter(xdata, ydata)
plot!(xdata, prediction)

function fourier!(y::AbstractArray,
                  x::AbstractArray,
                  h::Int64,
                  p0::Vector{Float64})

    if (h == 1)
        @. model1(t, p)= p[1] + p[2]*cos(2*π*t/24) + p[3]*sin(2*π*t/24)
        fit = curve_fit(model1,x, y, p0)
        return(fit)
    elseif (h==2)
        @. model2(t, p) = p[1] + p[2]*cos(2*π*t/24) + p[3]*sin(2*π*t/24) +
                          p[4]*cos(2*π*t/12) + p[5]*sin(2*π*t/12)
                fit = curve_fit(model2,x, y, p0)
            return(fit)
            return(model1)
    elseif (h==3)
        @. model3(t, p) = p[1] + p[2]*cos(2*π*t/24) + p[3]*sin(2*π*t/24) +
                          p[4]*cos(2*π*t/12) + p[5]*sin(2*π*t/12) +
                          p[6]*cos(2*π*t/8) + p[7]*sin(2*π*t/8)
                fit = curve_fit(model3,x, y, p0)
            return(fit)
        else
            error("Only  three harmonics allowed")
    end

end