import scipy.integrate as integrate
import numpy as np
import matplotlib.pyplot as plt

## Maximas secondaires de RN

def RN(phi,N) :  # Fonction réseau
    if phi not in [0,2*np.pi,4*np.pi] :
        return ((np.sin(N*phi/2))**2)/((np.sin(phi/2))**2)/(N**2)
    else :
        return 1

def MaxSecondRN(N) :
    max = 0
    phimax = 0
    PhiListe = np.linspace(2*np.pi/N,4*np.pi/N,10000)
    for phi in PhiListe :
        val = RN(phi,N)
        if max < val :
            max, phimax = val,phi
    return N*phimax/np.pi,max

# x = w/w0   pulsation réduite
# b = tc*w0  produit durée de cohérence x pulsation centrale raie gaussienne
# lo         longueur d'onde centrale de la raie gaussienne
# On prend Jmax = 1

def JS(x) :  # densité spectrale
    b = 1700  # Bon ordre de grandeur pour raies usuelles
    res = np.exp(-b*b/4*(x-1)**2)
    return res

def f(x,delta,lo) :
    N = 1000
    res = JS(x)*RN(2*np.pi*x*delta/lo,N)  # Js(w)*R(w*delta/c)
    return res


## Utilisation de integrate

IS,er = integrate.quad(JS,0,np.inf)  # IS = intensité totale d'une raie

DeltaListe = np.linspace(-10,1300,10000)  # Liste des différences de marche en nm : de -10 nm à 1300 nm
L1,L2 = [],[]

lo = 546  # raie verte nm
for delta in DeltaListe :
    Ia,er = integrate.quad(f,0,np.inf,args=(delta,lo))
    L1.append(Ia/IS)  # Liste des intensités relatives.

lo = 600  # raie rouge en nm
for delta in DeltaListe :
    Ia,er = integrate.quad(f,0,np.inf,args=(delta,lo))
    L2.append(Ia/IS)

plt.close()
plt.figure()
plt.plot(DeltaListe,L1, "g.")
plt.plot(DeltaListe,L2, "r.")
plt.xlabel("Différence de marche $\delta$")
plt.ylabel("I($\delta$)/Imax")
plt.title("Intensité versus différence de marche")
plt.show()