
import numpy as np
import matplotlib.pyplot as plt

def onda_plana(t, A, σ, ω_0, φ):
    """
    Genera una onda plana dada la frecuencia central y fase.
    Las unidades han de ser consistentes entre t, σ y ω_0.
    El pulso está normalizado.

    La onda plana se caracteriza por un término exponencial complejo que depende de la frecuencia central de la onda
    y la fase de la envolvente de la onda portadora.

    E(t) = E * exp(i * ω_0 * t) * exp(i * φ(t)) = A * exp(-t² / 2*σ) * exp(i * ( ω_0 * t + φ(t) ) )

    Argumentos:
        t (float): vector de tiempos
        A (float): amplitud del pulso
        σ (float): anchura del pulso
        ω_0 (float): frecuencia central (radianes / unidad de tiempo)
        φ (float): fase de la envolvente de la onda portadora (rad)

    Devuelve:
        E_pulso (float): forma del pulso gaussiano en el tiempo especificado
    """

    return A * np.exp(1j * ( ω_0 * t + φ ))

def pulso_gaussiano(t, A, σ, ω_0, φ):
    """
    Genera un pulso gaussiano dada su duración, frecuencia central y fase.
    Las unidades han de ser consistentes entre t, σ y ω_0.
    El pulso está normalizado.

    Un pulso gaussiano viene caracterizado por una envolvente en forma de gaussiana de expresión:

    E_envolvente = A * exp(-t² / 2*σ)

    Donde σ es la duración temporal del pulso, que está relacionada con el ancho de banda por la expresión:

    σ = FWHM / (2 * √log(2))

    FHWM es la anchura a media altura (full width half maximum).

    La envolvente viene modulada por un término exponencial complejo que depende de la frecuencia central de la onda,
    de manera que el pulso vendrá dado por el producto de la envolvente y esta exponencial, además del producto
    con la exponencial compleja que lleva la fase de la envolvente de la onda portadora.

    E(t) = E_envolvente * exp(i * ω_0 * t) * exp(i * φ(t)) = A * exp(-t² / 2*σ) * exp(i * ( ω_0 * t + φ(t) ) )

    Argumentos:
        t (float): vector de tiempos
        A (float): amplitud del pulso
        σ (float): anchura del pulso
        ω_0 (float): frecuencia central (radianes / unidad de tiempo)
        φ (float): fase de la envolvente de la onda portadora (rad)

    Devuelve:
        E_pulso (float): forma del pulso gaussiano en el tiempo especificado
    """

    return A * np.exp(-t*t / (2 * σ)) * np.exp(1j * ( ω_0 * t + φ ))

def pulso_sinc(t, A, σ, ω_0, φ):
    """
    Genera un pulso tipo sinc() gaussiano dada su duración, frecuencia central y fase.
    Las unidades han de ser consistentes entre t, σ y ω_0.
    El pulso está normalizado.

    Un pulso gaussiano viene caracterizado por una envolvente en forma de gaussiana de expresión:

    E_envolvente = A * sinc(t /*σ)


    La envolvente viene modulada por un término exponencial complejo que depende de la frecuencia central de la onda,
    de manera que el pulso vendrá dado por el producto de la envolvente y esta exponencial, además del producto
    con la exponencial compleja que lleva la fase de la envolvente de la onda portadora.

    E(t) = E_envolvente * exp(i * ω_0 * t) * exp(i * φ(t)) =A * sinc(t /σ) * exp(i * ( ω_0 * t + φ(t) ) )

    Los primeros mínimos se encuentran en +-pi*σ

    Argumentos:
        t (float): vector de tiempos
        A (float): amplitud del pulso
        σ (float): anchura del pulso
        ω_0 (float): frecuencia central (radianes / unidad de tiempo)
        φ (float): fase de la envolvente de la onda portadora (rad)

    Devuelve:
        E_pulso (float): forma del pulso gaussiano en el tiempo especificado
    """

    return  A*np.sinc(t / (np.pi* σ))*np.exp(1j * ( ω_0 * t + φ )) #np.sinc(x) = np.sin(pi*x)/(pi*x)) https://numpy.org/doc/stable/reference/generated/numpy.sinc.html


c=299792458.0     #m/s
pasos_temporales = 4000

# -- Parámetros del pulso --
A = 1 # Amplitud del pulso
λ_0 = 100.55 # Longitud de onda de ejemplo (en micrómetros)
ω_0 = 2 * np.pi * c * 1e-12 / (λ_0 * 1e-6) # Frecuencia angular del pulso (rad / ps)
σ = 0.2 # Duración del pulso (ps)
t, Δt = np.linspace(0, 200*σ, num=pasos_temporales, retstep=True) # Vector de tiempos (centrado en cero, ps). Guardamos la separación entre datos
to=np.max(t)/2
print("Δt (ps)=",Δt)
fs=1/Δt
print("fs=1/Δt sampling rate (1/ps)=",fs)

Δt (ps)= 0.010002500625156289
fs=1/Δt sampling rate (1/ps)= 99.97500000000001


tipo_onda='gaussiana'
if(tipo_onda=='onda plana'):
  # Fase
  φ_0 =  0* np.ones(pasos_temporales) # Fase (constante en este caso)
  #φ_0 = (t-to)**2 # Variable con t
  onda =onda_plana(t-to, A, σ, ω_0, φ_0)


if(tipo_onda=='gaussiana'):
  φ_0 =  0* np.ones(pasos_temporales) # Fase (constante en este caso)
  onda=pulso_gaussiano(t-to, A, σ, ω_0, φ_0) # Vector con el campo complejo del pulso


if(tipo_onda=='sinc'):
  φ_0 =  0* np.ones(pasos_temporales) # Fase (constante en este caso)
  onda=pulso_sinc(t-to, A, σ, ω_0, φ_0) # Vector con el campo complejo del pulso


plt.figure(figsize=(16, 4))
plt.plot(t,np.real(onda),label='$Re(E(t))$')
plt.plot(t,np.abs(onda),label='$|E(t)|$')
delay=0.0
plt.plot(t+delay,np.real(onda),label=r'$Re(E(t+\tau))$',
                                color='red',linestyle='--')
plt.axhline(0.0,color='grey',linestyle='--')
plt.xlabel(r'$t$ (ps)')
plt.ylabel(r'$Re(E)$')
plt.axvline(to+np.pi*σ,color='grey',linestyle='--',label=r'$\pi \sigma$')
#plt.xlim(10,10+2*np.pi*1)
plt.legend()
plt.show()


def introduce_wave(x,N,pos):
    '''
    Genera un vector x "expandido" incrustando el original en un vector
    con más pasos temporales e iniciado con ceros
    '''
    n=x.shape[0]
    y=np.zeros((N*n),dtype=complex)
    if n > y.shape[0]:
        raise ValueError("n should be less than or equal to the number of columns in y.")
    # Copy the y array to avoid modifying it in-place
    x_expand = y.copy()
    # Determine the starting index for insertion
    # Range: 0, y.shape[0] - n + 1
    # Substitute the values from x into y
    x_expand[ pos:pos + n] = x
    return x_expand

def gamma_tau(E,tau,shape_ini,times_shape_ini,tau0):
  '''
  Calcula la función de correlación (promedio temporal de E*(t)|E(t+tau))
  salvo factores ctes (tiempo de simulación T y dt )
  '''
  N=times_shape_ini  # Los vectores de campo expandido tendrá N veces el tamaño original
  Et=introduce_wave(E,N,tau0) # Campo E(t)
  Ettau=introduce_wave(E,N,tau)     # Campo E(t+tau)
  promedio_temporal=np.sum(np.conjugate(Et)*Ettau) # "Integral (suma)" de <E*(t)|E(t+tau)>
  return promedio_temporal

'''
Se calcula la función de correlación, \Gamma, para muchos valores de tau
La idea es "barrer" E(t+tau) sobre el campo E(t)
'''
E=onda # Se ha definido previamente el campo en el array onda
n=E.shape[0]
N=10
tau0=int(N*n/2)
Gamma =[]
for tau in range(0, n*N - n + 1):
  Gamma.append(gamma_tau(E,tau,n,N,tau0))
Gamma=np.array(Gamma)
gamma=Gamma/Gamma[tau0] # Se calcula el grado de coherencia temporal \gamma


Et=introduce_wave(E,N,int(N*n/2)) #int(N*n/2))
t_rango_expand=Δt*(np.linspace(0,Et.shape[0],Et.shape[0]))
plt.plot(t_rango_expand,np.real(Et),label=r'Re($E(t))$',
                                color='red',linestyle='-')
plt.xlabel(r'$t$ (ps)')
#plt.xlim(0,0.3)
plt.legend()
plt.show()


t_max=np.argmax(Gamma)
texpand=Δt*(np.linspace(0,Gamma.shape[0],Gamma.shape[0])-t_max)
plt.plot(texpand,np.abs(gamma),label=r'|$\gamma(\tau)|$',
                              color='red',linestyle='-')
plt.xlabel(r'$\tau-\tau_0$ (ps)')
plt.xlim(-100*σ,100*σ)
plt.legend()
plt.show()


def gamma_sinc_analitica(tau,ω_0):
  '''
  Podemos obtener gamma analítica para la distribución espectral de tipo salto
  Sabiendo que su transformada de Fourier, el pulso E(t), es una función proporcional
  a sinc(t).
  Definida la onda tipo sinc(t), la ventana de la distribución escalar es:
  delta_nu=1/(pi*sigma), donde el mínimo de la onda sinc() ocurre en σ*np.pi
  '''
  delta_nu=1.0/(σ*np.pi) # en t=σ*np.pi ocurre el mínimo de la onda sinc()
  return np.sinc(np.pi*delta_nu*tau*Δt/np.pi)*np.exp(1j*ω_0*tau*Δt) #np.sin(np.pi*delta_nu*tau*Δt)*np.exp(1j*ω_0*tau*Δt)/(np.pi*delta_nu*tau*Δt)

Gamma_sinc =[]
for tau in range(0, n*N - n + 1):
  Gamma_sinc.append(gamma_sinc_analitica(tau,ω_0))

Gamma_sinc=np.array(Gamma_sinc)


if (tipo_onda=='sinc'):
  t_max=np.argmax(Gamma)
  texpand=Δt*(np.linspace(0,Gamma.shape[0],Gamma.shape[0])-t_max)
  plt.plot(texpand,np.abs(gamma),label=r'|$\gamma(\tau)|$ numérica' ,
                                color='red',linestyle='-')
  texpand=Δt*(np.linspace(0,Gamma.shape[0],Gamma.shape[0]))
  plt.plot(texpand,np.abs(Gamma_sinc),label=r'|$\gamma(\tau)|$ analítica',
                                color='blue',linestyle='--')

  plt.axvline(np.pi*σ,color='grey',linestyle='--',label=r'$1/\Delta \nu=\pi \sigma$')
  plt.xlabel(r'$\tau-\tau_0$ (ps)')
  #plt.ylabel(r'$\gamma(\tau)$')
  plt.xlim(0,5*np.pi*σ)
  plt.legend()
  plt.show()


import numpy as np

def intensidad_young(xp,I1,I2,lambda_nm,n,d,D1,D2,xs,gamma_tipo='onda plana'):
  # Longitud de onda, k y frecuencia angular
  c=299792458.0     #m/s
  lambda0=lambda_nm*1e-9
  k = 2.0 * np.pi / lambda0
  omega0=c*k

  # Cálculo de delta
  delta= (2.0*np.pi*n/lambda0)*(xs*d/D1+xp*d/D2)

  # Se selecciona el tipo de onda (gamma_c)
  # Por defecto 'onda plana'
  # Se pueden añadir tantas como se quiera
  if(gamma_tipo=='onda plana'):
    gamma_c=1.0

  if(gamma_tipo=='perfil cuadrado'):
    tau=delta/omega0
    delta_lambda=20.0 #nm # delta_lambda=0.0 --> onda plana
    lambda_i=lambda0-delta_lambda*1e-9
    lambda_f=lambda0+delta_lambda*1e-9
    delta_nu =c*(lambda_i-lambda_f)/(lambda_f*lambda_i)
    gamma_c= np.sinc(np.pi*delta_nu*tau/np.pi) #np.sinc(x)=sin(pi*x)/(pi*x)

  # Cálculo de la intensidad
  return I1+I2+2.0*np.sqrt(I1*I2)*gamma_c*np.cos(delta)


# Intensidades
I1,I2=1.0,1.0 # Unidades arbitrarias

# Índice de refracción
n=1.0

# Distancias
D1=1.0  # m
D2=1.0  # m
d=0.001 # m
xs=0.0 # m

# Parámetros electromagnéticos
c=299792458.0     #m/s
lambda_nm = 600.0 # nm
lambda0=lambda_nm*1e-9
k = 2.0 * np.pi / lambda0
omega0=c*k

print("Frecuencia central=",np.round(omega0*1e-12/(2.0*np.pi),2),"THz")

# gamma_c
gamma_tipo='perfil cuadrado'


def calc_intensity_pattern():
    Xp=400
    m=10
    xpm = (2*m+1)*lambda0*D2/(2.0*d)-xs*D2/D1
    diff_pattern = []
    for i in range(Xp):
      print("x iteration",i,' of ',Xp)
      xp=xpm * (Xp - 2 * i) / Xp
      I_analytic=intensidad_young(xp,I1,I2,lambda_nm,n,d,D1,D2,xs,
                                  gamma_tipo=gamma_tipo)
      diff_pattern.append([xp,I_analytic])
    return np.array(diff_pattern)

# Calculo de la intensidad
intensity_pattern=calc_intensity_pattern()


import matplotlib.pyplot as plt
x=intensity_pattern[:,0]
y_analytic=intensity_pattern[:,1]
plt.plot(x,y_analytic,label='Analytic',color='red')

plt.axvline((2*0+1)*lambda0*D2/(2.0*d)-xs*D2/D1,color='grey',linestyle='--',label=r'$x_{min}(m=0)$')
plt.axvline((1)*lambda0*D2/d-xs*D2/D1,color='orange',linestyle='--',label=r'$x_{max}(m=1)$')
plt.xlabel(r'$x$ (m)')
plt.ylabel('Intensity')
plt.legend()
plt.show()


def calc_intensity_pattern_2D():
    Xp,Yp=400,400
    m=10
    xpm = (2*m+1)*lambda0*D2/(2.0*d)-xs*D2/D1
    ypm = xpm
    diff_pattern = []
    for i in range(Xp):
      print("x iteration",i,' of ',Xp)
      xp=xpm * (Xp - 2 * i) / Xp
      for j in range(Yp):
        yp=ypm * (Yp - 2 * j) / Yp
        I_analytic=intensidad_young(xp,I1,I2,lambda_nm,n,d,D1,D2,xs,
                                  gamma_tipo=gamma_tipo)
        diff_pattern.append([xp,yp,I_analytic])
    return np.array(diff_pattern)

intensity_pattern_2D=calc_intensity_pattern_2D()


import numpy as np
import matplotlib.pyplot as plt

# Generate some sample data (replace this with your numpy array)
data =intensity_pattern_2D

# Extract x, y, and the value from the data
x = data[:, 0]
y = data[:, 1]
value = data[:, 2]

# Reshape the value to a 2D grid for contour plotting
x_unique = np.unique(x)
y_unique = np.unique(y)
X, Y = np.meshgrid(x_unique, y_unique)
Z = value.reshape(len(y_unique), len(x_unique))


# Create a contour plot
plt.figure(figsize=(12, 2))
max_value=np.max(Z)/1.0
contour = plt.contourf(X, Y,  np.rot90(Z), levels=500, cmap='viridis',
                       vmax=max_value)
plt.colorbar(contour)  # Use scientific notation for colorbar


plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('Intensity')
plt.show()


import matplotlib.pyplot as plt
import numpy as np
import ipywidgets as widgets
from ipywidgets import interactive
from ipywidgets import SelectionSlider

def plot_function_with_sliders(f, labels, x_min, y_min,ref_x,max_y,value_scale):
    # Extract the parameter names from the function signature, excluding 'x'
    # and the last
    parameters = list(f.__code__.co_varnames)[1:f.__code__.co_argcount-1]

    # Create sliders for each parameter
    # Define the base, min exponent, and max exponent for the slider
    sliders = {}
    i=0
    for param in parameters:
        #print("parameter=",param)
        #print("[initial value, min, max, step]=",value_scale[i])
        sliders[param] = widgets.FloatSlider(value=value_scale[i][0],
                                             min=value_scale[i][1],
                                             max=value_scale[i][2],
                                             step=value_scale[i][3],
                                             description=param,
                                             readout_format='.4f')
        i+=1
    # Create sliders for x_max and y_max
    #print("x max")
    #print("[initial value, min, max, step]=",value_scale[i])
    x_max_slider = widgets.FloatSlider(value=value_scale[i][0],
                                       min=value_scale[i][1],
                                       max=value_scale[i][2],
                                       step=value_scale[i][3],
                                       description='x scale',
                                       readout_format='.4f')
    i+=1
    #print("y max")
    #print("[initial value, min, max, step]=",value_scale[i])
    y_max_slider = widgets.FloatSlider(value=value_scale[i][0],
                                       min=value_scale[i][1],
                                       max=value_scale[i][2],
                                       step=value_scale[i][3],
                                       description='y scale',
                                       readout_format='.4f')

    # Define a function to update the plot
    def update_plot(**kwargs):
        N=200
        # Extract x, y, and the value from the data
        x = np.linspace(x_min-x_max_slider.value, abs(x_min)+x_max_slider.value, N)
        y = np.linspace(y_min-y_max_slider.value, abs(y_min)+y_max_slider.value, N)

        # Pass slider values as keyword arguments to the input function
        params = {param: slider.value for param, slider in sliders.items()}
        params['gamma_tipo'] = value_scale[-1][0]

        # Cálculo de la función
        Z =f(x, **params)

        # Reshape the value to a 2D grid for contour plotting
        x_unique = np.unique(x)
        y_unique = np.unique(y)
        X, Y = np.meshgrid(x_unique, y_unique)

        # Repeat the row N times to create a (N,N) array
        Z = np.tile(Z, (N, 1))

        # Result as a (N, N) array
        Z = Z.T  # Transpose the array
        Z = Z.reshape(len(y_unique), len(x_unique))

        # Create a contour plot
        plt.figure(figsize=(12, 2))
        max_value=np.max(Z)
        print("Valor mínimo=",np.min(Z))
        print("Valor máximo=",max_value)
        contour = plt.contourf(X, Y,  np.rot90(Z), levels=50, cmap='viridis')

        # Set the minimum and maximum values for the colorbar
        contour.set_clim(0.0,max_value)  # Set the desired min and max values here

        plt.colorbar(contour)

        plt.xlim((x_min-x_max_slider.value)/ref_x, (abs(x_min)+x_max_slider.value)/ref_x)
        plt.ylim((y_min-y_max_slider.value)/ref_x, (abs(y_min)+y_max_slider.value)/ref_x)
        plt.xlabel(labels[0])
        plt.ylabel(labels[1])
        plt.title('Intensity')
        plt.show()

    # Create an interactive plot with the sliders
    interactive_plot = interactive(update_plot, **sliders, x_max=x_max_slider, y_max=y_max_slider)
    return interactive_plot


# Example of an input function
function=intensidad_young


# Define labels, axis limits, and initial parameter values
labels = [r'$x(m)$', r'$y(m)$']
m=10
xpm = (2*m+1)*lambda0*D2/(2.0*d)-xs*D2/D1
x_min = -xpm
y_min = -xpm
max_y = 1.0
ref_x = 1.0 #lambda0*D2/d
# value_scale (each): [initial value, min,max,step]
value_scale=[[I1,1e-4*I1,10*I1,0.01*I1], # I1
             [I2,1e-4*I2,10*I2,0.01*I2], # I2
             [lambda_nm,0.1*lambda_nm,3.0*lambda_nm,0.1*lambda_nm],      # lambda (wavelength) in nm
             [n,1.0,3.5,0.1],           # n
             [d,0.1*d,10*d,0.05*d],     # d
             [D1,0.1*D1,10*D1,0.05*D1], # D1
             [D2,0.1*D2,10*D2,0.05*D2], # D2
             [xs,0.0,1e-2,1e-3],        # xs
             [-xpm,0.1*xpm,10*xpm,0.05*xpm],    # x scale
             [-xpm,0.1*xpm,10*xpm,0.05*xpm],    # y scale
             [gamma_tipo]]       # tipo gamma

# Create the interactive plot
interactive_plot = plot_function_with_sliders(function,
                                              labels,
                                              x_min, y_min,
                                              ref_x,max_y,
                                              value_scale)

# Display the interactive plot
interactive_plot

Óptica - Tema 1 -Coherence and interference phenomena¶

Coherencia¶

Definición de ondas planas y pulsos (gaussiano, sinc)¶

Parámetros¶

Cálculo del pulso¶

Representación gráfica¶

Cálculo numérico del grado de coherencia $\gamma(\tau)$ a partir de $E(t)$¶

Ejemplo de campo "expandido"¶

Representación de $\gamma (\tau)$¶

Caso de onda tipo $sinc()$: comparación cálculo numérico vs analítico¶

Función analítica $\gamma (\tau)$ para una distribución espectral tipo salto¶

Representación de $\gamma (\tau)$¶

Interferencias: franjas de Young¶

Observación en función de las coordenadas $x,y$ de la pantalla.¶

Ejemplo $g(\nu)\rightarrow$ Tipo "perfil cuadrado" (simétrico)¶

Función de Python para calcular la intensidad de las franjas de Young (aproximación)¶

Parámetros geométricos¶

Experimento de Young en la dirección $x$ ($y=0$)¶

Cálculo de la intensidad¶

Figura estática¶

Experimento de Young plano $XY$¶

Cálculo de la intensidad¶

Figura estática¶

Figura dinámica¶