# -*- coding: utf-8 -*-
"""

Source : Vincent Grenard, merci à  lui !

XXX brouillon XXX
pointillé noir = rayons correspondant à  un objet sur l'axe à  la position 
vérifiant les relations de conjugaison

rayons rouges : correspondent à  l'objet le plus loin net
rayons bleus : correspondent à  l'objet le plus proche net


Les calculs sont fait dans le cadre de l'approximation des petits angles, ne 
tient pas compte du caractére non stigmatique des lentilles, en particulier aux 
grandes ouvertures. Calculs non garantis car fait vite fait sur un brouillon.
(Basé sur le formalisme de l'optique matricielle, plus facile qu'avec les 
relations de conjugaison et de la trigo)

Ne tiens pas compte non plus du phénomène de diffraction qui est le phénomène
(à  ma connaissance) qui limite la résolution lorsque les diaphragmes sont très
fermés.

"""
import matplotlib.path as mpath
import matplotlib.patches as mpatches
import numpy as np #testé avec np.__version__ 1.16.1
import matplotlib.pyplot as plt #testé avec matplotlib.__version__ 2.0.0
# import PIL.Image as im#attention à  la majuscule ; testé avec PIL.__version__ '5.4.1'
from matplotlib.widgets import Slider, Button#, RadioButtons
# from scipy.signal import fftconvolve #testé avec scipy.__version__ 1.2.1
# plt.ion()

####### définition des constante et valeurs initiales

fp = 1 #distance focale de la lentille
xAp = 1.5 #position algébrique du capteur
epsi = 0.5 #taille(rayon) du pixel
R = 2.3 #rayon du diaphragme

inf = float("inf")
#####   limite de l'espace  pour matplotlib  #####
xA = 1/(-1/fp+1/xAp) #position de l'image
lim = [2*xA,xAp*2,-4,4]   # xdeb xfin ydeb yfin
marge = 0.05
x_deb = lim[0]
x_fin = lim[1]
y_deb = lim[2]
y_fin = lim[3]
delta_x = x_fin - x_deb
delta_y = y_fin - y_deb
###### fin de déf des constantes

def traverse_vide(x_old,y_old,theta_old,x_new):
    """ optique matricielle 
        traversé du vide : theta inchangé, y change : ynew = yold + theta*dist
    """
    return y_old + (x_new-x_old)*theta_old #theta_new = theta_old
def traverse_lentille(y_old,theta_old,fp):
    """ optique matricielle 
        traversé d'une lentille : y inchangé, theta_new = theta_old - y/f'
    """
    theta_new = -y_old/fp + theta_old
    return theta_new #y_new = y_old

def traverse_lentille_arriere(y_new, theta_new, fp):
    """ optique matricielle 
        traversé d'une lentille : y inchangé, theta_new = theta_old - y/f'
        ici, on connait theta_new et on veut theta_old
    """
    theta_old = theta_new + y_new/fp
    return theta_old





fig = plt.figure(figsize = [10,6])

#axe optique
plt.plot([x_deb,x_fin],[0,0],'-k')#axe optique
deb_fleche = x_fin - 0.025*delta_x
altitude_fleche = delta_y * 0.025
plt.plot( [deb_fleche,x_fin,deb_fleche],[altitude_fleche,0,-altitude_fleche],'-k') # flèche

#lentille
x_O, f_p = 0, fp #une seule lentille ici
plt.plot([x_O,x_O],[0.995*y_deb,0.995*y_fin],'-k',lw=2)
if f_p<0:
    sg = -1
else:
    sg = 1
lx_fl = 0.05*delta_x
ly_fl = sg * 0.05*delta_y
plt.plot([x_O-lx_fl,x_O,x_O+lx_fl],[0.995*y_fin-ly_fl,0.995*y_fin,0.995*y_fin-ly_fl],'-k',lw=2)
plt.plot([x_O-lx_fl,x_O,x_O+lx_fl],[0.995*y_deb+ly_fl,0.995*y_deb,0.995*y_deb+ly_fl],'-k',lw=2)

#foyers
plt.plot([x_O+f_p]*2,[delta_y*0.01,-delta_y*0.01])
plt.plot([x_O-f_p]*2,[delta_y*0.01,-delta_y*0.01])
plt.axis('off')
plt.axis(lim)

#cercle de confusion ; position xAp
plot_cercleConf1, = plt.plot([0.5*x_fin]*2,[-y_fin, y_fin],'k--')
plot_cercleConf, = plt.plot([0.5*x_fin]*2,[-epsi, epsi],lw = 3,color = "green")


#tracé du diaphragme
global diaphragme_bas, diaphragme_haut
diaphragme_bas, = plt.plot([-0.075,-0.075],[0.9*y_deb, -R],'-',lw = 5, color = [0.5,0.5,0.5])
diaphragme_haut, =plt.plot([-0.075,-0.075],[0.9*y_fin, R],'-',lw = 5, color = [0.5,0.5,0.5])


#rayons limites = rayon qui passe par le haut du diaphragme + haut et bas du cercle de confusion (et symétrique)
theta_lim_new1 = (epsi-R)/xAp
theta_lim_new2 = (-epsi-R)/xAp
theta_lim_old1 = traverse_lentille_arriere(R,theta_lim_new1,fp)
theta_lim_old2 = traverse_lentille_arriere(R,theta_lim_new2,fp)
x_old1 = -R/theta_lim_old1#point sur l'axe de départ pour le rayon
x_old2 = -R/theta_lim_old2
#XXX remarque, si x_old1>0, c'est qu'on a atteint l'hyperfocale => x_old1 réel = -inf


# plotlim1_haut, = plt.plot([x_old1,0,x_fin],[0,R,traverse_vide(0,R,theta_lim_new1,x_fin)],"-b")
# plotlim1_bas, = plt.plot([x_old1,0,x_fin],[0,-R,-traverse_vide(0,R,theta_lim_new1,x_fin)],"-b")

# plotlim2_haut, = plt.plot([x_old2,0,x_fin],[0,R,traverse_vide(0,R,theta_lim_new2,x_fin)],"-r")
# plotlim2_bas, = plt.plot([x_old2,0,x_fin],[0,-R,-traverse_vide(0,R,theta_lim_new2,x_fin)],"-r")

#netteté «parfaite»
plotnet_haut, = plt.plot([xA,0,x_fin],[0,R,traverse_vide(0,R,traverse_lentille(R,-R/xA,fp),x_fin)],"-k")
plotnet_bas, =plt.plot([xA,0,x_fin],[0,-R,-traverse_vide(0,R,traverse_lentille(R,-R/xA,fp),x_fin)],"-k")

pp3 = plt.Polygon([[xA, 0],
                   [0, R],
                   [traverse_lentille(R,-R/xA,fp), 0]])


####### création des sliders
axcolor = 'lightgoldenrodyellow'
# ax_R = plt.axes([0.05, 0.05, 0.2, 0.1], facecolor=axcolor)
ax_epsi = plt.axes([0.05, 0.05, 0.2, 0.1], facecolor=axcolor)
ax_xAp = plt.axes([0.05, 0.2, 0.2, 0.1], facecolor=axcolor)

# slider_R = Slider(ax_R, 'R', 0.001,0.8*y_fin, valinit=0.4*y_fin)#, valstep=1)
slider_epsi = Slider(ax_epsi, u'E', 0.001,0.8*x_fin, valinit=0.5*x_fin)
slider_xAp = Slider(ax_xAp, "xA", 1.01*fp,x_fin, valinit=(2*fp+x_fin)/3)



def update(osef):
    # R = slider_R.val
    x_epsi = slider_epsi.val
    xAp = slider_xAp.val
    xA = 1/(-1/fp+1/xAp) #position de l'image
    plot_cercleConf.set_data([x_epsi]*2,[-epsi, epsi])
    plot_cercleConf1.set_data([x_epsi]*2,[-y_fin, y_fin])
    diaphragme_bas.set_data([-0.075,-0.075],[0.9*y_deb, -R])
    diaphragme_haut.set_data([-0.075,-0.075],[0.9*y_fin, R])
    #rayons limites = rayon qui passe par le haut du diaphragme + haut et bas du cercle de confusion (et symétrique)
    theta_lim_new1 = (epsi-R)/xAp
    theta_lim_new2 = (-epsi-R)/xAp
    theta_lim_old1 = traverse_lentille_arriere(R,theta_lim_new1,fp)
    theta_lim_old2 = traverse_lentille_arriere(R,theta_lim_new2,fp)
    if theta_lim_old2 != 0:
        x_old2 = -R/theta_lim_old2#point sur l'axe de départ pour le rayon
    else:
        x_old2 = -inf
    x_old1 = -R/theta_lim_old1
    #XXX remarque, si x_old1>0, c'est qu'on a atteint l'hyperfocale => x_old1 réel = -inf
    
    # if x_old2 < 0:
    #     plotlim2_haut.set_data([x_old2,0,x_fin],[0,R,traverse_vide(0,R,theta_lim_new2,x_fin)])
    #     plotlim2_bas.set_data([x_old2,0,x_fin],[0,-R,-traverse_vide(0,R,theta_lim_new2,x_fin)])
    # else:
    #     plotlim2_haut.set_data([x_deb,0,x_fin],[R,R,traverse_vide(0,R,traverse_lentille(R,0,fp),x_fin)])
    #     plotlim2_bas.set_data([x_deb,0,x_fin],[-R,-R,-traverse_vide(0,R,traverse_lentille(R,0,fp),x_fin)])
    # plotlim1_haut.set_data([x_old1,0,x_fin],[0,R,traverse_vide(0,R,theta_lim_new1,x_fin)])
    # plotlim1_bas.set_data([x_old1,0,x_fin],[0,-R,-traverse_vide(0,R,theta_lim_new1,x_fin)])
    
    #netteté «parfaite»
    plotnet_haut.set_data([xA,0,x_fin],[0,R,traverse_vide(0,R,traverse_lentille(R,-R/xA,fp),x_fin)])
    plotnet_bas.set_data([xA,0,x_fin],[0,-R,-traverse_vide(0,R,traverse_lentille(R,-R/xA,fp),x_fin)])
    profondeur = x_old1-x_old2
    if profondeur < 0:
        print("profondeur de champ infini")
        
    else:
            
        print("profondeur de champ =",profondeur)
    if x_old2 < 0:
        print("dont",int(round(100*(x_old1-xA)/profondeur)),"% proche")
        print("et",int(round(-100*(x_old2-xA)/profondeur)),"% loin")
# slider_R.on_changed(update) #connection signal/slot
slider_epsi.on_changed(update)
slider_xAp.on_changed(update)




# --- Affichage de la fenêtre ---
plt.show()