# -*- coding: utf-8 -*-
from scipy.integrate import odeint
from scipy.optimize import fsolve
from pylab import figure, plot, show, xlabel, ylabel, quiver ,legend
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])




## całkowanie


# prametry
I=1
time = np.linspace(0.0,200.0,1000) 
# warunki początkowe
yinit = np.array([0.,.0]) 
y = odeint(FH, yinit, time, args =(I,))




figure(1) 
plot(time,y[:,0]) # y[:,0] pierwsza kolumna y - zawiera zmienną v
#plot(time,y[:,1]) # y[:,1] druga kolumna y - zawiera zmienną w
xlabel('t') 
ylabel('y') 
legend(('v','w'))

## obraz fazowy
figure(2)
wezly_v = np.arange(-2.5,2.5, 0.1)
wezly_w = 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')

## 2) pole wektorowe
#v,w = np.meshgrid(wezly_v , wezly_w )
#dv = v -v**3/3. -w +I
#dw = 0.08 *(v+0.7-0.8*w)
#quiver(v,w,dv,dw)

                
# 3) przykładowa krzywa całkowa
plot(y[:,0],y[:,1]) # y[:,0] is the first column of y 
xlabel('v') 
ylabel('w') 
show()