import Numeric.GSL.ODE
import Numeric.LinearAlgebra
import Graphics.Plot

vanderpol mu = do
    let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)]
        ts = linspace 1000 (0,50)
        sol = toColumns $ odeSolve (xdot mu) [1,0] ts
    mplot (ts : sol)
    mplot sol


harmonic w d = do
    let xdot w d t [x,v] = [v, a*x + b*v] where a = -w^2; b = -2*d*w
        ts = linspace 100 (0,20)
        sol = odeSolve (xdot w d) [1,0] ts
    mplot (ts : toColumns sol)


kepler v a = mplot (take 2 $ toColumns sol) where
    xdot t [x,y,vx,vy] = [vx,vy,x*k,y*k]
        where g=1
              k=(-g)*(x*x+y*y)**(-1.5)
    ts = linspace 100 (0,30)
    sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts
    degree = pi/180


main = do
    vanderpol 2
    harmonic 1 0
    harmonic 1 0.1
    kepler 0.3 60
    kepler 0.4 70