produce_orbitdiagram problem. about dynamicalsystems.jl HOT 7 CLOSED

using DataFrames, DifferentialEquations, DynamicalSystems, Plots

@inline @inbounds function vy0(C,p,v)
	"""
	Initial Velocity for cte jacobi conservation
		in:
			C = cte jacobi
			v = position and velocities of the particle
			cartesian coordinates
			p = mass ratio parameter == μ2
		out:
			ydot = velocity in y

	"""
		x, y, x̂= v
		μ = p
		#distance vectors
		r1 = sqrt((x+μ)^2 + y^2)
		r2 = sqrt((x - (1-μ))^2 + y^2)
		ŷ = sqrt(x^2 + y^2 - x̂^2 + 2*((1-μ)/r1 + μ/r2) - C)
		if x > 0
			return -ŷ
		end
		if x < 0
			return ŷ
		end
end


@inline @inbounds function eom!(du,u,p,t)
    """
    Planar Circular Restricted 3body problem
    Equations of motion in co-rotating frame

        in:
            u = position and velocities of the particle
            cartesian coordinates
                x = u[1]
                y = u[2]
                vx = u[3]
                vy = u[4]
            p = mass ratio parameter == μ2
            t = time span
        out:
            du = integration of the acel. and vel.

    """
        local μ = p

        #cartesian coordinates
        x = u[1]
        y = u[2]
        vx = u[3]
        vy = u[4]

        r1 = sqrt((x + μ)^2 + y^2)
        r2 = sqrt((x - (1-μ))^2 + y^2)

        du[1] = vx # F1
        du[2] = vy # F2
        du[3] =  2*vy  + x - (1 - μ)*(μ + x)/(r1^3) - μ*(x -(1-μ))/(r2^3) # F3
        du[4] = -2*vx  + y - ((1 - μ)/(r1^3) + μ/(r2^3))*y # F4
        return
    end

p=0.07

x0=-1.24
y0=0.19
vx0=0
vyu = vy0(-1.2,p,[x0,y0,vx0])

u0 = [x0,y0,vx0,vyu]
plane = (3, 0.0)

ds = ContinuousDynamicalSystem(eom!, u0, p)
psos = poincaresos(ds, plane, 6000.0; u0 = u0,direction = 0,rootkw = (xrtol = 1e-9, atol = 1e-9),
        diffeq = (alg = Vern9()))
scatter(psos[:, 1], psos[:, 2]; markersize = 2.0)

#code works fine until here

pvalues = collect(0.1:0.15:201)
i = Int(1)
plane = (3, 0.0)
tf = 200.0
p_index = 1

output = produce_orbitdiagram(ds, plane, i,pvalues;
                              tfinal = tf, Ttr = 200.0)

output = produce_orbitdiagram(ds, plane, i,p_index,pvalues;
                              tfinal = tf, Ttr = 200.0)

no method matching setindex!(::Float64, ::Float64, ::Int64)