# -*- coding: utf-8 -*-
from scipy.integrate import odeint
from pylab import figure, plot, show, xlabel, ylabel ,legend,draw,ion,xlim, subplot,clf
import numpy as np

def w_pionowa(v,I):
    #izoklina stycznych pionowych (dv/dt=0)
    return v-v**3/3+I
 
def w_pozioma(v):
    #izoklina stycznych poziomych (dw/dt =0)
    return (v+0.7)/0.8   
    
def FH(y,t,I):
    # równania FH
    v = y[0]
    w = y[1]
    dv = v -v**3/3. -w +I
    dw = 0.08*(v + 0.7 - 0.8*w)  
    return np.array([dv, dw])

def rysuj(y,t,I):
        
    figure(1) 
    clf()
    subplot(2,1,2)
    plot(t,y[:,0],'b') # y[:,0] pierwsza kolumna y - zawiera zmienną v
    xlabel('t') 
    ylabel('y') 
    legend(('v',))
    xlim([0,600])
    
    ## obraz fazowy
    subplot(2,1,1)
    wezly_v = np.arange(-2.5,2.5, 0.1)
    # 1) izokliny
    plot(wezly_v, w_pionowa(wezly_v, I),'r')
    plot(wezly_v, w_pozioma(wezly_v),'g')
    
    # 3) przykładowa krzywa całkowa
    plot(y[:,0],y[:,1],'b') # y[:,0] is the first column of y 
    xlabel('v') 
    ylabel('w') 
    draw()
    raw_input()

## całkowanie

ion()
# prametry
I=0
time1 = np.linspace(0.0,200.0,1000) 
# warunki początkowe
yinit = np.array([0.,.0]) 
y1 = odeint(FH, yinit, time1, args =(I,))
rysuj(y1,time1,I)


yinit=y1[-1,:]
time2 = time1[-1]+np.linspace(0.0,200.0,1000)
I =  1 #6.125  #2.222
y2 = odeint(FH, yinit, time2, (I,))
rysuj(np.concatenate((y1,y2)),np.concatenate((time1,time2)),I)

yinit=y2[-1,:]
time3 = time2[-1]+np.linspace(0.0,200.0,1000)
I =  0 #6.125  #2.222
y3 = odeint(FH, yinit, time3, (I,))
rysuj(np.concatenate((y1,y2,y3)),np.concatenate((time1,time2,time3)),I)
show()