##First we load some packages.
from numpy import *
from matplotlib import *
from scipy import *
from pylab import *
from mpl_toolkits.mplot3d import *
import scipy as sp

#We define a function which is going to be the recursive function.
def tamari(x_t,y_t,z_t,h,a,b,c,d,e,f,g,u,i):
    dx_t=x_t+h*((x_t- a*y)*sp.cos(z_t) - b*y_t*sp.sin(z_t))
    dy_t=y_t+h*(sp.sin(z_t)*(c*y_t + x_t) + d*y_t*sp.cos(z_t))
    dz_t=z_t+h*(e + f*z_t + g*sp.arctan(((1-u)/(1-i))*(x_t*y_t)))
    return dx_t,dy_t,dz_t

#parametri
a=1.013
b=-0.011
c=0.02
d = 0.96
e = 0
f = 0.01
g = 1
u = 0.05
i = 0.05

#Definisco un rapporto incrementale lim h->0 di (t_finale -t_iniziale)/h
t_iniz=0
t_fin=100
h=0.0001
numsteps=int((t_fin-t_iniz)/h)

#using this parameters we build the time.
t=linspace(t_iniz,t_fin,numsteps)
#And the vectors for the solutions
x=zeros(numsteps)
y=zeros(numsteps)
z=zeros(numsteps)

#We set the initial conditions
x[0]=0
y[0]=0
z[0]=0

#Trange in size -1 perchè la prima condizione imposta è quella iniziale
for k in range(x.size-1):
    #We use the previous point to generate the new point using the recursion
    [x[k+1],y[k+1],z[k+1]]=tamari(x[k],y[k],z[k],t[k+1]-t[k],a,b,c,d,e,f,g,u,i)

#Now that we have the solution in vectors t,x,y,z is time to plot them.

#We create a figure and 4 axes on it. 3 of the axes are going to be 2D and the fourth one is a 3D plot.

fig = figure()
ax1 = fig.add_axes([0.1, 0.7, 0.4, 0.2])
ax2 = fig.add_axes([0.1, 0.4, 0.4, 0.2])
ax3 = fig.add_axes([0.1, 0.1, 0.4, 0.2])
ax4 = fig.add_axes([0.55, 0.25, 0.35, 0.5],projection='3d')

#And we add vectors to each plot

#grafico x,y,z, attrattore
ax1.plot(t, x,color='green',label='x(t)')
ax1.set_xlabel('t')
ax1.legend()
ax1.axis((t_iniz,t_fin,min(x),max(x)))

ax2.plot(t, y,color='orange',label='y(t)')
ax2.set_xlabel('t')
ax2.legend()
ax2.axis((t_iniz,t_fin,min(y),max(y)))

ax3.plot(t, z,color='blue',label='z(t)')
ax3.set_xlabel('t')
ax3.legend()
ax3.axis((t_iniz,t_fin,min(z),max(z)))

ax4.plot(x, y,z,color='black')
ax4.set_xlabel('x(t)')
ax4.set_ylabel('y(t)')
ax4.set_zlabel('z(t)')
ax4.set_title('Attrattore di Tamari')


show()