## Déviation d'un rayon lumineux dans un milieu d'indice optique variable
##JCourtin 06/2016


#importation des modules
from matplotlib import pyplot as plt
from math import pi,cos,sin,tan,acos,sqrt, atan
import time

############################ PARTIE A COMMENTER #################################

#Tracé de l'ensemble
def plotItALL(alpha, start, cadre, N):    #start--> (index, xPos, alpha, monte)
    Xmin, Xmax, Zmin, Zmax = cadre
    plt.close()
    plt.plot([Xmin,Xmax], [0,0], 'blue'  )#tracé de ??
    plt.plot([0,0],[Zmin,Zmax], 'blue')
    
    
    for i in range(N+1):  #tracé des ?????
        plt.plot([Xmin,Xmax], [z[i],z[i]], 'black')
    
    nP=[k*1e-3 for k in n]  #1e-3: facteur d'échelle
    plt.plot(nP,z)
    
    
    debut = time.clock()     #Chronomètre : que mesure-t-on ?
    xTraj , zTraj  = trajectoire(start)
    fin=time.clock()
    print("temps stratification : ",fin-debut)
    
    
    debut=time.clock()     #Chronomètre : que mesure-t-on ?
    xTrajE, zTrajE = trajectoireEuler(start, indiceQuad, gradIndiceQuad)
    fin=time.clock()
    print("temps Euler : ",fin-debut)

    
    #Que fait la commande suivante ?
    zTheo=[theorieQuad(x, alphaStart, deltaZ, Ninf, Nsup) for x in xTrajE]


    #tracé des courbes
    plt.plot(xTraj , zTraj , 'g')
    plt.plot(xTrajE, zTrajE, 'b')         #Que fait la ligne suivante ?
    plt.plot([xTrajE[k] for k in range(len(xTrajE)) if k%10==0], \
             [zTrajE[k] for k in range(len(xTrajE)) if k%10==0], 'bo')
    plt.plot(xTrajE, zTheo, 'r')    
    plt.show()


############################ PARTIE A MODIFIER #################################
#bien lire le MAIN avant de se lancer


##Trajectoire par la méthode en strate   --- Partie 1 ---
def trajectoire(start):
    
    
    #on pose alpha toujours positif [angle / horizontale]
    #Soit Descartes II : n1.cos(alpha1) = n2.cos(alpha2)
    #
    #on note "monte" un booléen True si on monte et False si on descend 
    #
    #Start -->  (index, xPos , alpha, monte)
    #Le but est simplement de construire le tuple start sur la strate suivante
    #On enregistre au passage les coordonnées du point start suivant.(xTraj, zTraj)
    #Enfin on boucle jusqu'à sortir du domaine d'étude.
    #
    #On utilise les listes z[k] et n[k] calculées dans le MAIN
    
    index=start[0]
    alpha=start[2]
    LOOP=True
    N=0
    
    xTraj=[start[1]]  #CI trajectoire par hypothèse
    zTraj=[0.]
    
    while(LOOP and xTraj[-1]<Xmax): #tracé du rayon et calul d'angle
    
        N+=0#corriger
        (index, xPos , alpha, monte) = stop#corriger
        nextX = xPos + (z[index+1]-z[index])/cos(alpha) #corriger 
                                                        #cf triangle rectangle
        n1=n[index]
    
        if (monte):
            n2=n[index+1];  nextIndex=index-1#monte    #!!! corriger
        else:
            n2=n[index-1];  nextIndex=index+1#descend  #!!! corriger
            
        xTraj+=[]    #compléter
        zTraj+=[]    #compléter
        
        if (n1/n2*cos(alpha)>1.): #nouveau start : Expliquer pourquoi ???
            monte=not(monte)                     # => change de sens
            nextAlpha=alpha
        else:
            nextAlpha=acos(n1/n2*cos(alpha)) #expliquer le cos (qui est correct)
        
        #point de départ suivant
        start=(nextIndex, nextX, nextAlpha, monte)
        
        #condition d'arret
        if (monte and z[index]>z[-2]): LOOP=False
        if (not(monte) and z[index]<z[2]):LOOP=False  
    
    #renvoie le graph de la trajectoire
    return (xTraj, zTraj)
 
 
 
 
    
##calcul de la trajectoire par la méthode d'Euler     --- Partie 2 ---
def trajectoireEuler(start, indice, gradIndice):


    #Dans cet méthode on fixe un déplacement élémentaire ds dont la direction
    #est donnée par le vecteur unitaire u.
    #On doit donc simplement trouver le prochain vecteur u et avancer de ds 
    #dans sa direction. On part du meme point Start ["monte" ne sert à rien ici]
    #
    #On doit donc fournir à la méthode la fonction indice ainsi que son gradient


    #ds=1.0e-5#m
    ds=(Zmax - Zmin)/N#m
    
    alpha=0#corriger
    xTraj=[42]  #Corriger #CI trajectoire par hypothèse
    zTraj=["Toto"] #corriger
    
    n=indice(zTraj[-1], deltaZ, Ninf, Nsup) #expliquer
    
    u=(sin(alpha), tan(alpha))  #calcul de u    #!!! corriger
   
    while(xTraj[-1]<Zmin):  #!!! corriger et expliquer 
    
    
        xTraj+=[xTraj[1]-ds*arcos(alpha)]  #avance de ds sur l'arc #!!! corriger
        zTraj+=[zTraj[1]+ds*arctan(alpha)] #!!! corriger
        
        
        nNew=indice(zTraj[-1], deltaZ, Ninf, Nsup)#nouvel indice
        dn=nNew-n                                 #dn variation d'indice
        n=nNew
        
        grad_R=gradIndiceQuad(zTraj[-1], deltaZ, Ninf, Nsup)#gradient selon er
        
        u0=0#formule            #re-calcul la direction u
        u1=0#formule  #cf methode Euler 
        u=(u0,u1)
        
        alpha=atan(u[0]*u[1]) #!!! corriger #re-calcul alpha
        
    return (x, z) #corriger 




##Trajectoire théorique du rayon
def theorieQuad(x, alphaStart, deltaZ, Nmin, Nmax): #cf cours
    kappa=(Nmax**2-Nmin**2)/(2*Nmax**2)
    zTheo=deltaZ*sin(alphaStart)/sqrt(2*kappa)
    zTheo=zTheo*sin(sqrt(2*kappa)/(deltaZ*cos(alphaStart))*x)
    return zTheo
    
    

##différentes lois d'indice
def indiceQuad(z, deltaZ, Nmin, Nmax):      #loi d'indice "quadratique"
    kappa=(Nmax**2-Nmin**2)/(2*Nmax**2)
    if (abs(z)>deltaZ):
        return #???
    else:
        return Nmax*sqrt(1-2*kappa*(z/deltaZ)**2)

def gradIndiceQuad(z, deltaZ, Nmin, Nmax):  #gradient associé
    kappa=(Nmax**2-Nmin**2)/(2*Nmax**2)
    if (abs(z)>deltaZ):
        return #?????
    else:
        return -2*kappa*Nmax/deltaZ**2/sqrt(1-2*kappa*(z/deltaZ)**2)*z

def indiceLin(z, deltaZ, Nmin, Nmax):        #loi d'indice "linéaire"
    if (abs(z)>deltaZ):
        return #????
    else:
        return #Trouver la loi d'indice



    
############ MAIN ##############  NE PAS MODIFIER ##############################

if (__name__=="__main__"):
        
##definition de l'indice en fonction de z
    
    
    N=50    #Nombre de strates
    
    z=[]     #altitudes des strates
    n=[]     #indice de la strate
    
    Zmin=-3.0e-4  #m
    Zmax=+3.0e-4  #m
    Xmin=0.       #m
    Xmax=9.9e-3  #m
    
    cadre =(Xmin, Xmax, Zmin, Zmax) #domaine du tracé
    
    Nsup=2.4
    Ninf=1.
    deltaZ=2.0e-4 #demi-largeur de la fibre
    
##Calcul de la loi d'indice
    for i in range(N+1):
        z+=[Zmin+(Zmax-Zmin)/N*i]  #position des strates
        n+=[indiceQuad(z[i], deltaZ, Ninf, Nsup)] #indice de la strate

        
##conditions initiales
    alphaStart=pi/4.
    start = (N//2, Xmin, alphaStart, True) #--> (index, xPos , alpha, monte)
    
    plotItALL(alphaStart, start, cadre, N) #trace les courbes