# -*- coding: utf-8 -*-
"""
Created on Wed Jan  1 14:30:33 2015
ok
@author: joseorellana
"""
# import basic modules
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import pylab as pyl
import math
import cmath
import os  

def cm0(alpha):
	return alpha*(.4-alpha)

def dcm0(alpha):
	return .4-2*alpha
	
	
NPA=1000
V=200.
L=1.
S=.1
ro=1.2
x=np.linspace(0,1.1,NPA)
xdot=np.linspace(-8,8,NPA)
pb=np.zeros((NPA,NPA))
ce=np.zeros((NPA,NPA))

coef1=(ro*S*L*V**2)/2

for i in range (NPA):
	for j in range (NPA):
		xa=math.atan(math.tan(x[i])-L*xdot[j]/(V*math.cos(x[i])))
		dxa=-L/(V*math.cos(x[i]))*math.cos(xa)**2
		t1=(xdot[j]*L/V-math.sin(x[i]))*cm0(xa)
		t2=(1-2*(xdot[j]*L/V)*math.sin(x[i])+(xdot[j]*L/V)**2)*dcm0(xa)*dxa/2
		pb[j,i]=(t1+t2)*coef1
		pb[j,i]=min(.01,max(pb[j,i],0))
		
		

figure=plt.subplot(1, 3,1)	
fig = plt.figure(1)
ax = fig.gca()
Y, X = pyl.meshgrid(x, xdot)
pyl.contourf(X, Y, pb, 8, alpha=.75, cmap='summer')
ax.set_ylabel('vitesse angulaire')
ax.set_xlabel('angle d incidence')
plt.title('critere de Poincare-Bendixson')



for i in range (NPA):
	for j in range (NPA):
		xa=math.atan(math.tan(x[i]-L*xdot[j]/(V*math.cos(x[i]))))
		t3=(1-(2*xdot[j]*L*math.sin(x[i]))/V+(xdot[j]*L/V)**2)*cm0(xa)-cm0(x[i])
		if np.abs(xdot[j])>0:
			t4=coef1*t3/xdot[j]
		else:
			t4=coef1*t3/.000001
		ce[j,i]=t4
		ce[j,i]=min(.01,max(ce[j,i],0))

figure=plt.subplot(1, 3,2)	
fig = plt.figure(1)
ax = fig.gca()
Y, X = pyl.meshgrid(x, xdot)
pyl.contourf(X, Y,ce, 8, alpha=.75, cmap='summer')
ax.set_ylabel('vitesse angulaire')
ax.set_xlabel('angle d incidence')
plt.title('critere de l energie')

figure=plt.subplot(1, 3,3)	
fig = plt.figure(1)
ax = fig.gca()
Y, X = pyl.meshgrid(x, xdot)
pyl.contourf(X, Y, pb-ce, 8, alpha=.75, cmap='summer')
ax.set_ylabel('vitesse angulaire')
ax.set_xlabel('angle d incidence')
plt.title('Difference des criteres')
plt.show()