# -*- coding: utf-8 -*-
"""
Created on Wed Jan  1 14:30:33 2015
"""
# import basic modules
import numpy as np
import os
import matplotlib . pyplot as plt
import scipy.sparse as sp
import scipy.sparse.linalg  as la
import time
import math
import scipy.io as sio
from mpl_toolkits.mplot3d import Axes3D
import pylab as pyl

plt.close()
rhoe=1000 #mass volumique de l'eau
rhoa=1.2 #masse volumique de l'air
sa=900 #surface de référence des oiles fois le Cz des voile
se=.1 # surface de référence de la structure immergée
V=20 #vitesse du vent absolu
cx=0.1 # trainée des foils
#cx=.3 # trainée des flotteurss
#se=1 #surface des flotteurss
#=========
npvit=600 #nombre de points pour la vitesse du bateau
vitmin=0.1*V # %vitesse mini
vitmax=3.*V #%vitesse maxi
dx=(vitmax-vitmin)/npvit # pas de vitesse
#========
nptheta=3 # nombre de points dans la discretisation en allure
thetamin=.3 # %allure mini en radian
thetamax=1.4 # allure maxi en radian

dtheta=(thetamax-thetamin)/nptheta #pas en allure
#===========
npalpha=600 # nombre de points en bordage
alphamin=.1 # bordage mini
alphamax=thetamax # bordage maxi avant d'être à contre

dalpha=(alphamax-alphamin)/npalpha # pas en bordage

#=======
f=np.zeros((npalpha,npvit)) #fonction Fx-Tx
ord0=np.linspace(alphamin,alphamax,npalpha);
ab=np.linspace(vitmin,vitmax,npvit)
#figure
#============
for i in range (nptheta):
	theta=thetamin+dtheta*i
	for k in range (npvit):
		x=(vitmin+dx*k)/V
		ab[k]=x*V
		for j in range (npalpha):
			alpha=min(theta,alphamin+dalpha*j)
			ord0[j]=alpha
			cosmua=max(0.,(math.sin(theta-alpha)-x*math.sin(alpha))/(math.sqrt(1+2*x*math.cos(theta)+x**2)))
			cnumua=1.2*cosmua**2 # coefficient de propulsion vélique
			xi=rhoa*sa*cnumua*math.sin(alpha)/(rhoe*se*cx)#%poussée vélique
			f[j,k]=1*min(.00001,max(-.00001,(xi-1)*x**2+2*x*math.cos(theta)*xi));
			
		
        
	fig = plt.figure()
	ax = fig.gca()
	X, Y = pyl.meshgrid(ab, ord0)
	pyl.contourf(X, Y, f, 8, alpha=.75, cmap='summer')
	ax.set_xlabel('vitesse u du bateau')
	ax.set_ylabel(u'\u03B1<\u03B8 qui est le reglage du plan de voilure')
	plt.title('Vitesse du vent: '+str(V)+u' Angle amure \u03B8= '+str(theta))




plt.show()