# -*- coding: utf-8 -*-
"""
Created on Tue Mar 14 09:47:22 2017

@author: valerie
"""
from scipy.special import jn
from scipy.special import kn
import math
import numpy

from scipy.optimize import bisect
import pylab as plt;
import numpy as np
from math import sqrt


#On donne les arguments i et m comme definie dans le rapport => LP[i,m] plus besoin de toucher le reste du programme
m=1
i=1
v=72.6












#//////////////////////////////////////////////////////DEBUT PROGRAMME////////////////////////////////////////////////////////////////



def isNaN(x):
    return (x == x) == False

#Initialisation de la fonction
def f(x):
    if i!=0:
        t=v*sqrt(1-x)*(jn(n-1,v*sqrt(1-x)))/jn(n,v*sqrt(1-x))+v*sqrt(x)*(kn(n,v*sqrt(x))/kn(n-1,v*sqrt(x)))
    else:t=v*sqrt(1-x)*(jn(1,v*sqrt(1-x)))/jn(0,v*sqrt(1-x))+v*sqrt(x)*(kn(1,v*sqrt(x))/kn(0,v*sqrt(x))) # j'ai rajouté une condition sur i pour le cas ou i=0.
    return t

#Fonction qui incrémente a et b de dx, b=a+dx jusqu'à ce que le produit f(a)*f(b) soit négatif avec une pente entre a et b petite. La fonction donne alors l'intervalle a, b.
def recherche(f,a,b,dx):
    x1=a
    x2=a+dx
    f1=f(a)
    f2=f(x2)
    delta=10000000
    
    #Ajout d'un test supplémentaire: on avance aussi tant que la pente entre x2 et x1 reste grande (branche infinie)
    while (f1*f2>0.0) or (abs((f2-f1)/dx)>delta):
        if x1 >= b:
            return x1,x2
            print ('Il n y a plus de racines sur cet intervalle [%f,%f]'%(a,b))
            break
           
        x1 = x2; f1 = f2
        x2 = x1 + dx; f2 = f(x2)
        
        if isNaN(x2)== True:
            print isNaN(x2)
            return x1,x2
            break
       # raw_input() 
    return x1,x2
    
#Fonction qui calcule la racine de la fonction f sur l'intervalle [a,b]    
def racine(f,a,b):
    root = bisect(f,a,b)
    return root

#Fonction qui applique recherche et racine en boucle tant que la borne sup de l'intervalle reste inférieure à b.
def recherche_racines_multiples(f,a,b,dx):
    
    x1temp=a
    x2temp=b
    while x1temp<b:
        x1,x2=recherche(f,x1temp,x2temp,dx)
       
        if x1>=b:
            break
        if isNaN(x2)== True:
            break
        x=round(racine(f,x1,x2),6)
        resultat.append(x)
       
       
#        print ('%f est une racine' %x) #j'ai commenté afin de pas encombrer affichage
        x1temp=x+dx
        x2temp=b+dx
        xs=x
        
    return xs

  
#
Listedesvaleurs = {}
for n in range(i,i+1):
    resultat=[]
    recherche_racines_multiples(f,0.0001,0.9999,0.00001)
    #print [xs for xs in resultat] #j'ai commenté afin de pas encombrer affichage

#    print n, resultat #j'ai commenté afin de pas encombrer affichage
    Listedesvaleurs[n]= resultat
#print Listedesvaleurs #j'ai commenté afin de pas encombrer affichage
    
print 'les differentes racine possible sont',resultat

if m>len(resultat): # au cas ou il n'y ait pas de racine correspondant au m choisi
    print 'il n y a pas de racine correspondant a cette valeur de m (',m,') veuillez baisser la valeur de m'
    print 'il y a',len(resultat),'racine(s) possible(s)'
#//////////////////////////////////////////////////TRACE////////////////////////////////////////////////////////////////////



A=E0=1
a=62.5
n1=1.458
n2=1.453
lda=0.6328
K0=2*math.pi/lda


b=resultat[len(resultat)-m]
U=v*(1-b)**0.5
W=v*b**0.5

x = np.linspace(-70, 70,1000)
y = np.linspace(-70, 70,1000)
X, Y = np.meshgrid(x, y)



R=np.sqrt((X**2)+(Y**2))
phi=np.arctan(Y/X)


Z=(A/jn(i,U))*jn(i,(R*U/a))*np.cos(i*phi)
P=(A/kn(i,W))*kn(i,(R*W/a))*np.cos(i*phi)

Z[R>=a]=P[R>=a]
z1=abs(Z*Z)
z2=np.max(z1)
z3=255*(z1/z2)


 # calcul du tableau des valeurs de Z
plt.imshow(z3, interpolation="bicubic",origin="lower", extent=[-63,63,-63,63],cmap=plt.cm.gray)

#plt.plot(x,Z)
#plt.axis([60,70,0,1])



plt.colorbar()
 

plt.title('Modes')

plt.show()