Open In Colab

.Licencia Creative Commons
<span xmlns:dct="http://purl.org/dc/terms/" property="dct:title"></span> The following notes written by <span xmlns:cc="http://creativecommons.org/ns#" property="cc:attributionName">Sergio Gutiérrez Rodrigo (sergut@unizar.es) </span>. Distributed under License Creative Commons Atribución-NoComercial-CompartirIgual 4.0 Internacional

Departamento de Física Aplicada
Universidad de Zaragoza
Instituto de Nanociencia y Materiales de Aragón (INMA)
C/ Pedro Cerbuna, 12, 50009, Zaragoza, España

Óptica - Tema 1- Fermat's principle, Lagrangian optics and the ray equation

La ecuación de la Eikonal, formulada por primera vez por William Rowan Hamilton en 1831, se empleada en numerosos campos dentro de la ciencia y la ingeniería. Su origen se encuentra tanto en la teoría de propagación de ondas como en la óptica geométrica. Obteniendo sus soluciones seremos capaces de conocer las trayectorias de los rayos en cualquier tipo de medio, independientemente de lo complejo que sea.

La trayectoria de un rayo en un medio óptico es una geodésica, cuyas características dependerán del índice de refracción ($n$). En un medio homogéneo, es decir donde $n=\text{cte}$, los rayos serán líneas rectas. Sin embargo, nos centraremos en sistemas ópticos donde el índice de refracción cambia de forma continua con las direcciones del espacio provocando desviaciones en los rayos.

En un ámbito puramente geométrico, las geodésicas se definen como el camino de mínima distancia entre un número finito de puntos pertenecientes a una superficie. Sin embargo, Fermat dedujo que en óptica las curvas que definían los rayos eran caminos que aseguraban el mínimo tiempo transcurrido para ir de un punto a otro. Este último resultado puede ser reformulado como un problema de cálculo variacional, donde el rayo entre dos puntos A y B será aquel para el cual el tiempo transcurrido ($\tau$) sea un extremo:

El camino óptico se define como:

$$\mathcal{L} = \int_A^B n \, ds$$

Principio de Fermat:

$$\delta \mathcal{L} = \delta \left( \int_A^B n \, ds \right)=0$$

Ya que $n=c/v$, el tiempo $\tau$ es un funcional del camino óptico $n \, ds$, \begin{equation} \tau=\int_A^B \frac{ds}{v} \rightarrow \delta\tau=0 \end{equation}

que se conoce como principio de mínima acción donde $ds$ representa la longitud de arco. Esta última relación implica que el camino óptico es una variable estacionaria dentro de nuestro sistema, lo cual significa que es un máximo o un mínimo tal y como buscábamos.

Podemos entonces deducir la ecuación de rayos empleando las ecuaciones de Euler-Lagrange.

La longitud de arco $ds$ y la velocidad $v$ vienen definidas como: \begin{equation} ds=v dt=\sqrt{x'^2+y'^2+z'^2} dt \hspace{0.75cm} v=\frac{c}{n} \end{equation} donde se utiliza la notación $'$ para las derivadas temporales $d/dt$.

Teniendo en cuenta que las variables espaciales en este caso dependen del tiempo. Introduciendo estas consideraciones podemos definir un lagrangiano mediante en cual derivaremos la ecuación de la Eikonal. \begin{equation} \delta \mathcal{L}=\delta \left( \int_{t_1}^{t_2}n(x,y,z)\sqrt{x'^2+y'^2+z'^2}dt\right)=0 \end{equation}

$$\implies L(x,x',y,y',z,z',t)=n(x,y,z)\sqrt{x'^2+y'^2+z'^2}$$

Ecuaciones de Euler-Lagrange aplicadas al lagrangiano óptico $L(x,x',y,y',z,z',t)=n(x,y,z)\sqrt{x'^2+y'^2+z'^2}$:

$$\dfrac{\partial L}{\partial q_i}-\dfrac{d}{dt}\left( \dfrac{\partial L}{\partial q_i'} \right)=0$$

Aquí $q_i=x,y,z$ , y $q_i'=x',y',z'$

Entonces:

$ \begin{equation} \dfrac{\partial L}{\partial x}-\dfrac{d}{dt}\left( \dfrac{\partial L}{\partial x'} \right)=0 \\ \dfrac{\partial L}{\partial y}-\dfrac{d}{dt}\left( \dfrac{\partial L}{\partial y'} \right)=0 \\ \dfrac{\partial L}{\partial z}-\dfrac{d}{dt}\left( \dfrac{\partial L}{\partial z'} \right)=0 \end{equation} $

Ya que $ds=\sqrt{x'^2+y'^2+z'^2} dt$

$$ \rightarrow \dfrac{\partial L}{\partial x}=v\left( \dfrac{\partial n}{\partial x} \right)$$
$$ \rightarrow \dfrac{\partial L}{\partial x'}=n\left( \dfrac{x'}{\sqrt{x'^2+y'^2+z'^2}} \right)=\dfrac{n}{v}\dfrac{dx}{dt}=\dfrac{n}{v}\dfrac{dx}{ds}\dfrac{ds}{dt}=n\dfrac{dx}{ds}$$

Finalmente:

$$\rightarrow \dfrac{d}{dt}\left( n\dfrac{d x}{ds}\right)=\dfrac{ds}{dt}\dfrac{d}{ds}\left( n\dfrac{d x}{ds}\right)=v\dfrac{d}{ds}\left( n\dfrac{d x}{ds}\right)$$

Y por lo tanto: $$\dfrac{\partial n}{\partial x}=\dfrac{d}{ds}\left( n\dfrac{d x}{ds}\right)$$

Empleando las ecuaciones de Euler-Lagrange al resto de componentes se llega a la ecuación de las trayectores de la ótica geométrica: \begin{equation} \frac{d}{ds}\left(n\frac{d\textbf{r}}{ds}\right)=\nabla n \end{equation} donde $\textbf{r}=(x,y,z)$.

La parte de la izquierda de la igualdad representa la curvatura del rayo, mientras que en el otro lado encontramos el gradiente del índice de refracción. En consecuencia, el rayo siempre se doblará hacia la región de máximo gradiente. Si consideramos un medio isótropo, el lado derecho sería igual a cero ya que $n$ no dependería de las direcciones del espacio, es por eso que las trayectorias son líneas rectas como se comentó anteriormente.

About waves and rays

El momento $\vec p$ se define a partir del lagrangiano $L$ como

  1. $p_i=\dfrac{\partial L}{\partial x_i'}$ where $x_1 = x, x_2= y, x_3= z$
  2. Además $p_i=\dfrac{\partial L}{\partial x_i'}=n\dfrac{d x_i}{ds}$ para el lagrangiano óptico.

  1. $L=n\dfrac{x_1'^2}{\sqrt{x_1'^2+x_2'^2+x_3'^2}}+n\dfrac{x_2'^2}{\sqrt{x_1'^2+x_2'^2+x_3'^2}}+n\dfrac{x_3'^2}{\sqrt{x_1'^2+x_2'^2+x_3'^2}} \implies L=\vec v \vec p$

  2. Usando (3) en $$\mathcal{L} = \int_A^B n \, ds=\int_A^B L(x_i,x_i',t) \, dt$$ Se llega a que $$\mathcal{L} = \int_A^B \vec p \,. d\vec l$$

    • Se ha utilizado que $\int_{t_1}^{t_2} p_i v_i dt=\int_{t_1}^{t_2} p_i dx_i$*
    • $\delta \mathcal{L} = 0$ usando el momento se conoce como princpio de mínima acción (Maupertius).

  3. Otro resultado notable es que: $$\dfrac{\partial \mathcal{L}}{\partial x_i}=\int \dfrac{\partial L}{\partial x_i}\, dt=\int \dfrac{d p_i}{d t}\, dt=\int dp_i=p_i$$ Y por lo tanto $$\vec p = \nabla \mathcal{L}$$

    • $\vec p$ es tangente a la trayectoria de la luz (los rayos) en cada punto
    • $\vec p$ es ortogonal a cualquier superficie $\ni \mathcal{L} =cte$ (respecto de una referencia), ya que el gradiente define superficies equipotenciales por ser conservativo $\int_A^B \nabla \mathcal{L} \, d\vec l =\mathcal{L}(B)-\mathcal{L}(A)$, $\implies \int_C \nabla \mathcal{L} \, d\vec l =0 $ si $C \in S \, \ni \mathcal{L}(S)=cte$
    • $\vec p$ es un campo conservativo, ya que $\int_A^B \nabla \mathcal{L} \, d\vec l =\mathcal{L}(B)-\mathcal{L}(A) \, \forall A,B$ +El camino óptico entre dos frentes de onda es el mismo para cualquier trayectoria seguida. Si $A$ y $B$ pertenecen a un frente de onda $S_1$ y, $C$ y $D$ pertenecen a un frente de onda $S_2 \implies \int_B^C n\,ds=\int_A^D n\,ds$

GRadient-INdex medium (GRIN)

En esta sección mostraremos las ecuaciones a resolver mediante el uso de PINN para el caso de un medio GRIN, en particular una fibra óptica.

Consideramos un índice de refracción cuadrático \begin{equation} n^2(x,y)=n_o^2\left(1-g^2(x^2+y^2)\right) \end{equation}

donde g será una constante real característica del medio y $n_o$ el índice de refracción en el eje de simetría ($z$).

Se puede observar cómo $n(x,y)$ decrece conforme nos alejamos en dirección radial de este. Principalmente, resolveremos dos casos: uno dará lugar a una ecuación diferencial ordinaria y el último a dos ecuaciones diferenciales acopladas. Previamente motivaremos el interés de conocer las soluciones de esta ecuación con este particular índice de refracción.

Los medios no homogéneos han sido siempre de gran interés para aquellos que trabajan en campos de óptica puesto que son sistemas que se encuentran en la propia naturaleza. Las lentes del ojo humano y la atmósfera son claros ejemplos de estos. Además ofrecen grandes ventajas en cuanto a la construcción de instrumentos ópticos: reducción de tamaños, costes, peso... A la hora de crear un sistema óptico, el primer paso es construir un método adecuado para el trazado de rayos. El principio teórico de esta metodología es la ecuación de la Eikonal, puesto que permite la simulación de todo tipo de medios. En general, se suelen emplear aproximaciones de primer orden para estudiar las ecuaciones, es decir, rayos paraxiales (aquellos cuyas trayectorias forman ángulos pequeños con el eje óptico).

Aunque la ecuación de rayos apareciese hace bastante tiempo, gracias a los nuevos desarrollos y estudios podemos considerar nuevas simetrías de forma eficiente e incluso en algunos casos podemos integrar la Eikonal casi de forma completa. Actualmente, uno de los grandes usos de este tipo de medios es la creación de fibras ópticas ya que transmiten la información de manera más eficiente. La curvatura de los rayos provoca que estos no se reflejen en las superficies de la fibra y por tanto las pérdidas disminuyen

Forma aproximada

Considerando $g^2(x^2+y^2)\ll 1$ podemos aproximar el índice de refracción empleando el desarrollo en serie de Taylor

Desarrollo serie de Taylor: \begin{equation} \sqrt{1+x}\approx 1+\frac{x}{2}-\frac{x^2}{8}\dots \end{equation}

con el fin de obtener una expresión del gradiente más sencilla.

Usando coordenadas cilíndricas con $r=\sqrt{x^2+y^2}$ tenemos

\begin{equation} n(r)=n_o\sqrt{1-(gr)^2}\approx n_o\left(1-\frac{(gr)^2}{2}\right) \end{equation}

Para el caso de rayos paraxiales (casi paralelos a la dirección de propagación): $$d/ds \approx d/dz$$ y por lo tanto podemos escribir \begin{equation} \frac{d}{dz}\left(n\frac{d\textbf{r}}{dz}\right)=-n_og^2\textbf{r} \hspace{0.5cm} \textbf{r}=(x,y)\,. \end{equation}

En la parte izquierda de la última igualdad aplicamos la regla de la cadena.

Parte izquierda:

  • $n\neq n(z)$
  • $1-g^2(x^2+y^2)\approx 1 \rightarrow n\approx n_0$.

Obtenemos:

\begin{equation} \frac{d^2\textbf{r}}{d^2z}+g^2\textbf{r}=0 \hspace{0.5cm} \textbf{r}=(x,y) \,. \end{equation}

Rayos meridionales (contenidos en un plano que contiene al eje óptico $z$)

\begin{equation} \frac{d^2r}{dz^2}+g^2r=0 \hspace{0.5cm} r=\sqrt{x^2+y^2} \end{equation}

donde $r$ representa la distancia radial al eje $z$.

Esta ecuación puede ser resuelta analíticamente de forma sencilla ya que es de la forma del oscilador armónico, por lo tanto nos servirá de gran ayuda para contrastar nuestros resultados.

La solución es:

\begin{equation} r(z)=r_o\cos(gz)+r_o'\sin(gz)/g \end{equation}

con $r_o$ y $r'_o$ condiciones iniciales que determinan respectivamente la posición y la tangente del ángulo con el eje $z$.

Vemos como en este caso nuestra solución estará en un plano que contiene al eje $z$.

  • La trayectoria del rayo es periódica con período $\Lambda=2\pi/g$.
  • Si consideramos varios rayos todos paralelos al eje $z$ es decir con $r'_o=0$ tras una distancia $d$ tal que $gd=\pi/2$ todos ellos colapsarán en un punto del eje $z$ que se denomina punto focal.

Forma no aproximada

En esta situación, no consideramos un gradiente aproximado. Además, no particularizaremos la ecuación obtenida para rayos meridionales. La parte izquierda de la expresión \eqref{eq: paso intermedio aproximacion meridional} no variará, puesto que seguimos manteniendo la aproximación $1-g^2(x^2+y^2)\approx 1$. El gradiente vendrá definido por: \begin{equation} \nabla n(x,y)=\left(\frac{\partial n}{\partial x}, \frac{\partial n}{\partial y}\right) \end{equation} donde por la forma cuadrática de $n(x,y)$ tendremos que las derivadas son simétricas. Finalmente obtenemos dos ecuaciones diferenciales ordinarias acopladas \begin{equation} \left\lbrace\begin{array}{cc} \displaystyle\frac{d^2x}{dz^2}=& \displaystyle\frac{-g^2x}{\sqrt{1-g^2(x^2+y^2)}} \\ \displaystyle\frac{d^2y}{dz^2}=& \displaystyle\frac{-g^2y}{\sqrt{1-g^2(x^2+y^2)}} \end{array}\right. \end{equation} Notar que en este caso debemos de tener especial cuidado con los valores de $(x,y)$ para los cuales la raíz se anule, ya que provocarán que nuestras ecuaciones diverjan. Otro aspecto importante a tener en cuenta es que el índice de refracción siempre sea real, ya que en una fibra óptica no tiene sentido tener medios absorbentes, que son aquellos donde $n$ posee una parte imaginaria. Al igual que en el caso anterior podremos definir soluciones en el plano XZ o YZ, pero también más generales. Podremos estudiar las trayectorias de los rayos que no se encuentran confinadas en un solo plano. En las siguientes secciones fijaremos uno de los dos planos anteriores mediante la imposición de las condiciones iniciales, con el fin de poder comprar resultados con la solución aproximada.

Important libraries

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import cmath
import math
import time

Analytical ray path ($g^2(x^2+y^2)\ll 1$)

In [ ]:
def PathRay(z,r_o=0.5,r_1=-0.8,g=0.25):
    # Period = 2*pi/g
    return r_o*np.cos(g*z)+r_1*np.sin(g*z)/g

z=np.linspace(0,L,n_steps)
z=np.reshape(z,(n_steps,1))

print('g=',g)
r_o=np.sqrt(x_inicial**2+y_inicial**2)
r=PathRay(z,r_o=r_o,r_1=x_prime_inicial,g=g)

g= 0.5

Runge Kutta

In [ ]:
def rungekutta4(f, x0, t):
    n = len(t)
    x = np.zeros((n, len(x0)))
    x[0] = x0
    for i in range(n - 1):
        h = t[i+1] - t[i]
        k1 = f(x[i])
        k2 = f(x[i] + k1 * h / 2.)
        k3 = f(x[i] + k2 * h / 2.)
        k4 = f(x[i] + k3 * h)
        x[i+1] = x[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)
    return x
In [ ]:
#Function that returns the equations of eikonal fiber
def eikonal_fiber(x_sol):
  dx_dz = x_sol[2]
  dy_dz = x_sol[3]
  dx2_dz2 = -(g**2)*x_sol[0]/np.real(cmath.sqrt(1-(g**2*(x_sol[0]**2+x_sol[1]**2))))
  dy2_dz2 = -(g**2)*x_sol[1]/np.real(cmath.sqrt(1-(g**2*(x_sol[0]**2+x_sol[1]**2))))
  return np.array([dx_dz, dy_dz,dx2_dz2, dy2_dz2])

def generate_data(n_steps,tmax,init_cond,batch_size):
  data=[]
  targets=[]
  t = np.linspace(0, tmax, n_steps)
  for i in range (1,batch_size+1):
    x0 = np.random.uniform(-1.0, 1.0)
    y0 = np.random.uniform(-1.0, 1.0)
    z0 = np.random.uniform(-1.0, 1.0)
    results =  rungekutta4(f=eikonal_fiber, x0=init_cond,t=t)
    data.append(results[:n_steps])
  return np.array(data)

Parameters

In [ ]:
g=0.5

# Period = 2.0*np.math.pi/g
P=2.0*np.math.pi/g
L =3*P

n_steps=500 # Number of points in the domain

#Initial conditions of space and solution
# Initial conditions
z0=0.0
x0=[1.0,0.0]
dx_dz0=[0.0,0.2]
x_inicial,y_inicial=x0
x_prime_inicial,y_prime_inicial=dx_dz0

Generate data with Runge-Kutta

In [ ]:
batch_size=1
init_cond = [x_inicial, y_inicial, x_prime_inicial, y_prime_inicial]
inicial=time.time()
data=np.squeeze(generate_data(n_steps,L,init_cond,batch_size))
final=time.time()
print(final-inicial)
print(data.shape)
x_final=data[:,0]
y_final=data[:,1]

0.022270679473876953
(500, 4)

Plot the results

In [ ]:
plt.plot(z, r, color="green",marker='o',markersize=0.0, linestyle='-.', linewidth=1, label="x(z)=y(z) aprox.")
plt.plot(z,x_final,color="red",markersize=0.0, linestyle='--', linewidth=1,  label="x(z) Runge Kutta")
plt.plot(z,y_final,color="blue",markersize=0.0, linestyle='--', linewidth=1,  label="y(z) Runge Kutta")

plt.legend(fontsize=12)
plt.xlabel("z",fontsize=16)
plt.ylabel("Position at plane XY",fontsize=16)
plt.show()
In [ ]:
ax = plt.axes(projection='3d')
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

ax.scatter(x_final,y_final,z, c='g',s=1,label="(x(z),y(z)) RK")
ax.legend()
plt.show()