from pyndamics import Simulation

sim=Simulation()   # get a simulation object

sim.add("mice'=b*mice - d*mice",    # the equations
    100,                            # initial value
    plot=True)                      # display a plot

sim.params(b=1.1,d=0.08)

sim.run(0,4)

sim=Simulation()   # get a simulation object

sim.stock("mice",100,plot=False)
sim.inflow("mice","b*mice")
sim.outflow("mice","d*mice")

sim.params(b=1.1,d=0.08)

sim.run(0,4)

x,y=sim.t,sim.mice
plot(x,y,'r--')
xlabel('Days')
ylabel('Number of Mice')

sim=Simulation()

sim.add("deer' = r*deer*(1-deer/K)-c*deer*wolf",
                initial_value=350,
                plot=True)

sim.add("wolf' = -Wd*wolf+D*deer*wolf",
                initial_value=50,
                plot=True)

sim.params(r=0.25,D=0.001,c=0.005,Wd=0.3,K=1e500)

print sim.equations()
sim.run(0,500)

from pyndamics import phase_plot
phase_plot(sim,'deer','wolf')

sim=Simulation()

# logistic growth
sim.add("pop' = r*pop*(1-pop/K)",
                initial_value=350,
                plot=1)

# exponential growth
sim.add("pop2' = r*pop2",
                initial_value=350,
                plot=1)


sim.params(r=0.25,K=3000)

sim.run(0,5)


sim=Simulation()
sim.add("x''=-k*x/m -b*x'",[10,0],plot=True)
sim.params(k=1.0,m=1.0,b=0.5)
print sim.equations()
sim.run(0,20)

sim=Simulation()
sim.add("p'=p*(1-p)",0.1,plot=True)
sim.run(0,10)

from pyndamics import vector_field
vector_field(sim,p=linspace(-1,2,20))

vector_field(sim,rescale=True,p=linspace(-1,2,20))

sim=Simulation()
sim.add("x'=sigma*(y-x)",14,plot=True)
sim.add("y'=x*(rho-z)-y",8.1,plot=True)
sim.add("z'=x*y-beta*z",45,plot=True)
sim.params(sigma=10,beta=8.0/3,rho=15)
sim.run(0,50,num_iterations=10000)  # increase the resolution

phase_plot(sim,'x','z')

phase_plot(sim,'x','y','z')

sim=Simulation()
sim.add("x1''=(-b1 * x1' - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1",[0.5,0.0],plot=[1,2])
sim.add("x2''=(-b2 * x2' - k2 * (x2 - x1 - L2)) / m2",[2.25,0.0],plot=[1,2])
sim.params( m1 = 1.0, # masses
            m2 = 1.5, # masses
            k1 = 8.0, # Spring constants
            k2 = 40.0,
            L1 = 0.5, # Natural lengths
            L2 = 1.0,
            b1 = 0.8, # Friction coefficients
            b2 = 0.5,
        )


sim.run(0,10) 

N=2
sim=Simulation()
sim.add_system(N,"x[i]''=(-b[i] * x[i]' - k[i] * (x[i] -x[i-1] -L[i]) + k[i+1] * (x[i+1] - x[i] - L[i+1])) / m[i]",
           [[0.5,0.0],[2.25,0.0]],  # initial conditions
           [[1,2],[1,2]],  # plot arguments
           )

sim.system_params(N,
            k=[8.0,40.0],
            m=[1,1.5],
            L=[0.5,1],
            b=[0.8,0.5],
            
            x_1=0,   # these are the boundary conditions on each side
            x_N=0,
            k_N=0,
            L_N=0,
        )

print "Here are the equations:"
print sim.equations()
sim.run(0,10) 

N=3
sim=Simulation()
sim.add_system(N,"x[i]''=(-b[i] * x[i]' - k[i] * (x[i] -x[i-1] -L[i]) + k[i+1] * (x[i+1] - x[i] - L[i+1])) / m[i]",
           [[0.5,0.0],[2.25,0.0],[2.5,0.0]],  # initial conditions
           [[1,2],[1,2],[1,2]],  # plot arguments
           )

sim.system_params(N,
            k=[8.0,40.0,30.0],
            m=[1,1.5,0.5],
            L=[0.5,1,2],
            b=[0.8,0.5,.1],
            
            x_1=0,   # these are the boundary conditions on each side
            x_N=0,
            k_N=0,
            L_N=0,
        )

print "Here are the equations:"
print sim.equations()
sim.run(0,10)