#!/usr/bin/env python
"""
solve and plot a system of coupled partial differential equations which result in a Lorenz attractor

\dot{x} = a(y-x)
\dot{y} = x(b-z) -y
\dot{x} = xy - cz

(see http://de.wikipedia.org/wiki/Lorenz-Attraktor)

author: Bob Wimmer
wimmer@physik.uni-kiel.de

date: June 3, 2012

"""



from numpy import *
from math import *
from scipy.integrate import odeint
import Gnuplot, Gnuplot.funcutils

a= 10.
b= 28.
c= 8./3.

def f(v,t):
    x,y,z = v[0], v[1], v[2]
    f0 = a*(y - x)
    f1 = x*(b - z) - y
    f2 = x*y - c*z
    return [f0, f1, f2]



def main():
    #initial conditions:
    x0 = 1.
    y0 = 1.
    z0 = 1.
    v0 = [x0,y0,z0]
    t = linspace(0.,50., 10000)
    sol = odeint(f,v0,t)
    X = sol[:,0]
    Y = sol[:,1]
    Z = sol[:,2]
    g = Gnuplot.Gnuplot(persist=True)
    g('set xlabel "x"')
    g('set ylabel "y"')
    g('set zlabel "z"')
    g('set style data lines')
    g.splot(sol,v0)

main()