function parametrized_kuramoto(du,u,p,t)
    # NumberOfOscillators, CouplingConstant, and NaturalFreqs
    N,K,w = p
    # Centroid
    u[1] = 0
    for i in 1:N
        u[1] += cos(u[i+1]) + sin(u[i+1])*im
    end
    u[1] = u[1]/N
    r = abs(u[1])
    ψ = angle(u[1])
    # Oscillators
    for i in 1:N
    	du[i+1] = w[i] + K*r*sin(ψ-u[i+1])
    end
end

function getparameters_kuramoto(N,K,γ,distribution)
    distribution == "Uniform" ? nothing : error("$(distribution) is not yet supported")
    j = collect(1:N)
    w = broadcast(j->((2j-N-1)/(N-1))*((pi*γ)/2),j)
    p = (N=N, K=K, w=w)
    return p
end


function getstartpoints_kuramoto(N)
    d = Uniform(-pi,pi)
    u0::Array{ComplexF64} = rand(d,N)
    pushfirst!(u0,mean(u0))
    return u0
end

function kuramoto(N,K,T,fs,γ,distribution)
    tspan = (0.0,T)
    u0 = getstartpoints_kuramoto(N)
    p = getparameters_kuramoto(N,K,γ,distribution)
    prob = ODEProblem(parametrized_kuramoto,u0,tspan,p)
    sol = solve(prob,RK4(),saveat=collect(0:1/fs:T),progress=true)
end

