{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[{"file_id":"1JrJRijkDwZ11GZ0Mg68TmlwUZKdbbDjS","timestamp":1693486578147}],"toc_visible":true,"authorship_tag":"ABX9TyMDABHkoelRZeb9GYo6LeSg"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\")
"],"metadata":{"id":"znzSo4wE5fes"}},{"cell_type":"markdown","source":["![\"Licencia](\"https://i.creativecommons.org/l/by-nc-sa/4.0/88x31.png\")
The following notes written by Sergio Gutiérrez Rodrigo () . Distributed under License Creative Commons Atribución-NoComercial-CompartirIgual 4.0 Internacional"],"metadata":{"id":"tYZjdh8TvE7-"}},{"cell_type":"markdown","source":["```\n","Departamento de Física Aplicada\n","Universidad de Zaragoza\n","Instituto de Nanociencia y Materiales de Aragón (INMA)\n","C/ Pedro Cerbuna, 12, 50009, Zaragoza, España\n","```\n","\n","\n","\n","\n"],"metadata":{"id":"j2H5UHGe1Uf8"}},{"cell_type":"markdown","source":["---\n","# **Óptica - Tema 1 - Fraunhofer Diffraction: circular aperture**\n","\n","---"],"metadata":{"id":"D7vYGeB21ZQU"}},{"cell_type":"markdown","source":[""],"metadata":{"id":"Xl4pJdL7UOLK"}},{"cell_type":"markdown","source":["Ecuación general:\n","$$E_{cP'}=cte \\int_\\sigma e^{-\\imath k [(\\alpha_s-\\alpha_p)\\xi+(\\beta_s-\\beta_p)\\eta]}d\\sigma$$\n","\n","\n"],"metadata":{"id":"vTa6SHwTrAjp"}},{"cell_type":"markdown","source":["Planteamiento del problema:\n","\n","1. Fuente en el eje de simetría ($z$): $\\alpha_s=\\beta_s=0$\n","2. Diferencial en coordenadas cilíndricas:\n","$$d\\sigma=\\rho d\\rho d\\psi$$\n","\n","3. Además:\n","$$\\xi=\\rho \\cos(\\psi)$$\n","$$\\eta=\\rho \\sin(\\psi)$$\n","\n","4. Por simetría se puede tomar $\\psi'=0 \\rightarrow \\beta_p=0$ (coordenadas en la pantalla)\n","\n","5. $\\alpha_p=\\sin(\\theta)$\n"],"metadata":{"id":"jfkF1Nwhs0mq"}},{"cell_type":"markdown","source":["Ecuación final difracción por abertura circular de radio $R$:\n","$$E_{cP'}=cte \\int_0^R \\rho d\\rho \\int_0^{2\\pi}e^{\\imath k \\rho \\sin(\\theta)\\cos(\\psi)}d\\psi$$\n","\n","\n"],"metadata":{"id":"Lg2i9Bj5ty6y"}},{"cell_type":"markdown","source":[],"metadata":{"id":"U4pFjRvDUj2y"}},{"cell_type":"markdown","source":["En la situación habitual donde el objeto difractor se encuentra entre dos lentes convergentes:\n","$$\\sin(\\theta)\\approx \\dfrac{\\rho'}{f'}$$ donde $f'$ es la focal de la lente colectora y $\\rho'$ la distancia al eje desde las coordenadas ${X'Y'}$, en el plano focal imagen de la lente colectora."],"metadata":{"id":"O-MLNcvLu9LR"}},{"cell_type":"markdown","source":["La integral resulta ser:\n","$$E_{cP'}=cte \\int_0^R \\rho d\\rho \\int_0^{2\\pi}e^{\\imath k \\dfrac{\\rho \\rho'}{f'}\\cos(\\psi)}d\\psi$$\n","\n"],"metadata":{"id":"pS0KzsRvvet6"}},{"cell_type":"markdown","source":["Utilizando la definición de la función de Bessel de orden $0$:\n","$$J_0(u)=\\dfrac{1}{2\\pi}\\int_0^{2\\pi}e^{\\imath u \\cos(\\psi)}d\\psi$$\n","\n","Tomando $u=k\\dfrac{\\rho \\rho'}{f'}$\n","\n","Resulta que:\n","\n","$$E_{cP'}=cte \\int_0^R \\rho d\\rho [2\\pi J_0(k\\dfrac{\\rho \\rho'}{f'})]$$\n","\n","Utilizando la definición de la función de Bessel de orden $1$:\n","$$\\dfrac{R}{a}J_1(aR)=\\int_0^R J_0(a\\rho)\\rho d\\rho$$\n","\n","Encontramos que:\n","\n","$$E_{cP'}=cte \\times 2\\pi\\dfrac{J_1(aR)R}{a}$$\n","\n","donde $a=k\\dfrac{\\rho'}{f'}$ \n","\n","\n"],"metadata":{"id":"mck0dby6v6IS"}},{"cell_type":"markdown","source":["Es habitual definir $$Z=k\\dfrac{\\rho' R}{f'}=\\dfrac{\\pi}{\\lambda}\\dfrac{\\rho' D}{f'}$$ donde $D=2R$\n","\n","El resultado final, salvo constante de proporcionalidad ($\\in \\mathbb{C}$) es:\n","\n","$$E_{cP'} =cte R^2 [\\dfrac{2 J_1(Z)}{Z}]$$ y por lo tanto la intensidad será\n","\n","$$I_{P'}=E_{cP'} E^*_{cP'} =cte R^4(\\dfrac{2 J_1(Z)}{Z})^2$$\n","\n","*Note: $J_1(Z)/Z \\rightarrow 1/2 \\text{ for } Z \\rightarrow 0$*"],"metadata":{"id":"bF-7vNSSx9xj"}},{"cell_type":"markdown","source":["# Bessel functions (Bessel function of the first kind of real order and complex argument)"],"metadata":{"id":"bJevGc571jIc"}},{"cell_type":"code","source":["import numpy as np\n","import matplotlib.pyplot as plt\n","from scipy.special import jv,jn_zeros\n"," #https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.jv.html#scipy.special.jv\n","fig, ax = plt.subplots(dpi=150)\n","x = np.linspace(-15., 15., 1000)\n","for i in range(0,2):\n"," ax.plot(x, jv(i, x), label=f'$J_{i!r}$')\n","#ax.axvline(jn_zeros(1,1),color='green',linestyle='--') # Z1\n","ax.axvline(0,color='black',linestyle='--')\n","ax.axhline(0.0,color='black',linestyle='--')\n","plt.xlabel('$u$')\n","ax.legend()\n","plt.show()"],"metadata":{"id":"1DGR2OI_oNnI"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# Diffraction by an aperture (numerical)"],"metadata":{"id":"j-M41YaeL6pm"}},{"cell_type":"markdown","source":["For $r_p,r_s \\gg \\lambda \\Longrightarrow$\n","$E_{cP}=\\dfrac{\\imath k A}{4\\pi} \\int_\\sigma \\dfrac{e^{-\\imath k (r_p+r_s)}}{r_p r_s} \\left( \\cos(\\theta)+\\sin(\\chi)\\right ) d\\sigma$"],"metadata":{"id":"Ppg5Zso5jqu6"}},{"cell_type":"markdown","source":[""],"metadata":{"id":"QvpjvrMm_GBy"}},{"cell_type":"markdown","source":["# Fraunhofer diffraction by an aperture (numerical)"],"metadata":{"id":"x_TfJJ2IME4j"}},{"cell_type":"markdown","source":["For $l_p,l_s \\gg \\text{aperture dimensions} \\Longrightarrow$\n","$E_{cP}=cte \\int_\\sigma e^{-\\imath k [(\\alpha_s-\\alpha_p)\\xi+(\\beta_s-\\beta_p)\\eta]} d\\sigma$"],"metadata":{"id":"K_jROSgtME4q"}},{"cell_type":"markdown","source":[""],"metadata":{"id":"MAY7nJdYMDOJ"}},{"cell_type":"code","execution_count":null,"metadata":{"id":"H6mQ5eVnKq8j"},"outputs":[],"source":["import math\n","import numpy as np\n","\n","def diff(theta, phi, k, rs, rp):\n"," coss = (np.cos(theta) + np.cos(phi)) / (rs * rp)\n"," E = coss * np.exp(-1j*k* (rp + rs))\n"," return E\n","\n","def diffraction_intensity_circle_numerical(xs, ys, xp, yp, S, P, A, R,lambda_nm):\n"," k = 2.0 * math.pi / (lambda_nm*1e-9)\n"," dx = R/100 # cartesian mesh (X and Y directions)\n"," # Circle defined inside a square, sides (2R,2R)\n"," Nx = round(2 * R / dx)\n"," Ny = round(2 * R / dx)\n"," E = 0.0+0.0j\n"," for i in range(Nx):\n"," xa = R * (Nx - 2 * i) / Nx\n"," for j in range(Ny):\n"," ya = R * (Ny - 2 * j) / Ny\n"," if(xa**2+ya**2<=R**2):\n"," # Definition of the shape: points inside a circle\n"," rs = np.sqrt(S**2 + (xa - xs)**2 +(ya - ys)**2)\n"," rp = np.sqrt(P**2 + (xa - xp)**2 + (ya - yp)**2)\n"," # S=|zs| y P=|zp| cos(theta)=cos(N,-rs)=/|rs| & cos(phi)=/|rp|\n"," theta = math.acos(S / rs)\n"," phi = math.acos(P / rp)\n"," E+= diff(theta, phi, k, rs, rp)\n"," E = 1j*E * k * A * dx * dx / (4 * np.pi)\n"," return np.real(E*np.conjugate(E))"]},{"cell_type":"markdown","source":["# Fraunhofer diffraction (analytic)"],"metadata":{"id":"msMRElpxqPG9"}},{"cell_type":"code","source":["def diffraction_intensity_circle_analytic(xp,yp,fp, R, lambda_nm):\n"," k = 2.0 * math.pi / (lambda_nm*1e-9)\n"," rho=np.sqrt(xp**2+yp**2)\n"," Z=k*rho*R/fp\n"," Ecp=(R**2)*2*jv(1,Z)/Z\n"," Ip=Ecp*np.conjugate(Ecp)\n"," return Ip"],"metadata":{"id":"3HIsZOLO1p7i"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# Parameters (geometry)"],"metadata":{"id":"31WL-PawSmg-"}},{"cell_type":"markdown","source":["Infrarrojo lejano con abertura rectangular\n","+ 𝑆 = 100.0 𝑚 (distancia al emisor)\n","+ $P \\approx f'$ = 1.0 𝑚 (distancia a la pantalla)\n","+ $R$ = 1 𝑚𝑚 (anchura abertura horizontal)\n","+ λ= 600.0 n𝑚"],"metadata":{"id":"_r612Tx2Q3eY"}},{"cell_type":"code","source":["xs, ys = 0.0, 0.0 # meters\n","S, P = 20.0, 10.0 # meters\n","A = 1 # amplitude (arbitrary units)\n","R=1e-3 # meters\n","lambda_nm = 600.0 # nm\n","k = 2.0 * math.pi / (lambda_nm*1e-9)"],"metadata":{"id":"zeQjDEdwSplf"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# Diffraction pattern for $\\varphi=0$ (both analytic and numerical)"],"metadata":{"id":"Hs2F_KgvDsaS"}},{"cell_type":"code","source":["def calc_diff_pattern_phi0():\n"," Xp=100\n"," m=5\n"," xpm = m*1.22*P*lambda_nm*1e-9/(2.0*R)\n"," ypm = xpm\n"," diff_pattern = []\n"," for i in range(Xp):\n"," print(\"iteration\",i,' of ',Xp)\n"," xp=xpm * (Xp - 2 * i) / Xp - xs +1e-5\n"," yp=1.0e-9\n"," I_numerical= diffraction_intensity_circle_numerical(xs, ys, xp, yp, S, P, A,R, lambda_nm)\n"," I_analytic=diffraction_intensity_circle_analytic(xp-xs, yp-ys, P, R, lambda_nm)\n"," diff_pattern.append([xp,I_analytic,I_numerical])\n"," return np.array(diff_pattern)"],"metadata":{"id":"wPzGpwenEAp2"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["diff_pattern=calc_diff_pattern_phi0()"],"metadata":{"id":"xMTYICw6NYUh"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["## Plot static figure (both analytic and numerical)"],"metadata":{"id":"mgEJ9uPPr5y-"}},{"cell_type":"code","source":["import matplotlib.pyplot as plt\n","#Z1=1.22*P*lambda_nm*1e-9/(2.0*R)\n","Z1=jn_zeros(1,1)*P*lambda_nm*1e-9/(2.0*np.pi*R)\n","print(\"Z0 bessel=\",Z1)\n","x=diff_pattern[:,0]\n","y_analytic=diff_pattern[:,1]/np.max(diff_pattern[:,1]) # Normalized to 1\n","y_numerical=diff_pattern[:,2]/np.max(diff_pattern[:,2]) # Normalized to 1\n","#plt.plot(x,np.log(y_analytic),label='Analytic')\n","plt.plot(x,y_analytic,label='Analytic')\n","plt.axvline(Z1,color='green',linestyle='--')\n","plt.axvline(-Z1,color='green',linestyle='--')\n","plt.plot(x,y_numerical,label='Numerical',linestyle='--')\n","\n","plt.xlabel(r'$\\rho$ (m)')\n","plt.ylabel('Intensity (normalized to unity)')\n","plt.legend()\n","plt.show()"],"metadata":{"id":"T6Py4-VcDspb"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["## Plot dynamic figure (analytic)"],"metadata":{"id":"qGrtluxkT2Wf"}},{"cell_type":"code","source":["import matplotlib.pyplot as plt\n","import numpy as np\n","import ipywidgets as widgets\n","from ipywidgets import interactive\n","from ipywidgets import SelectionSlider\n","\n","def plot_function_with_sliders(f, labels, x_min, y_min,ref_x,max_y,value_scale):\n"," # Extract the parameter names from the function signature, excluding 'x'\n"," parameters = list(f.__code__.co_varnames)[1:f.__code__.co_argcount]\n","\n"," # Create sliders for each parameter\n"," # Define the base, min exponent, and max exponent for the slider\n"," sliders = {}\n"," i=0\n"," for param in parameters:\n"," #print(\"parameter=\",param)\n"," #print(\"[initial value, min, max, step]=\",value_scale[i])\n"," sliders[param] = widgets.FloatSlider(value=value_scale[i][0],\n"," min=value_scale[i][1],\n"," max=value_scale[i][2],\n"," step=value_scale[i][3],\n"," description=param,\n"," readout_format='.3f')\n"," i+=1\n","\n"," # Create sliders for x_max and y_max\n"," #print(\"r max\")\n"," #print(\"[initial value, min, max, step]=\",value_scale[i])\n"," x_max_slider = widgets.FloatSlider(value=value_scale[i][0],\n"," min=value_scale[i][1],\n"," max=value_scale[i][2],\n"," step=value_scale[i][3],\n"," description='radial scale',\n"," readout_format='.3f')\n"," #print(\"I max\")\n"," #print(\"[initial value, min, max, step]=\",value_scale[i+1])\n"," y_max_slider = widgets.FloatSlider(value=value_scale[i+1][0],\n"," min=value_scale[i+1][1],\n"," max=value_scale[i+1][2],\n"," step=value_scale[i+1][3],\n"," description='I scale',\n"," readout_format='.3f')\n","\n"," # Define a function to update the plot\n"," def update_plot(**kwargs):\n"," plt.figure(figsize=(8, 6))\n"," plt.xlabel(labels[0])\n"," plt.ylabel(labels[1])\n","\n"," x = np.linspace(x_min-x_max_slider.value, abs(x_min)+x_max_slider.value, 200)\n","\n"," # Pass slider values as keyword arguments to the input function\n"," params = {param: slider.value for param, slider in sliders.items()}\n"," y = f(x, **params)\n","\n"," plt.plot(x/ref_x, y/max_y)\n"," plt.scatter(x/ref_x, y/max_y,color='red')\n"," plt.title(\"Intensity normalized at $(x_p,y_p)=(0,0)$ value \"+str(max_y)+'$m^2$')\n"," plt.grid(True)\n"," plt.xlim((x_min-x_max_slider.value)/ref_x, (abs(x_min)+x_max_slider.value)/ref_x)\n"," plt.ylim(y_min, y_max_slider.value)\n"," plt.show()\n","\n"," # Create an interactive plot with the sliders\n"," interactive_plot = interactive(update_plot, **sliders, x_max=x_max_slider, y_max=y_max_slider)\n"," return interactive_plot"],"metadata":{"id":"HXhPvAk_T2cu"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Example of an input function\n","function=diffraction_intensity_circle_analytic #analytic_rectangle(xp, yp, P, a,b, k)\n","\n","\n","# Define labels, axis limits, and initial parameter values\n","labels = [r'$\\rho/R$', \"Intensity (norm.)\"]\n","m=5\n","Z1=jn_zeros(1,1)*P*lambda_nm*1e-9/(2.0*np.pi*R)\n","x_min = -m*Z1\n","y_min = 0.0\n","max_y = R**4\n","ref_x = R\n","# value_scale (each): [initial value, min,max,step]\n","value_scale=[[1e-9,-1e-2,1e-2,1e-3], # yp\n"," [P,0.1*P,10*P,0.1*P], # P\n"," [R,0.1*R,10*R,0.1*R], # R\n"," [lambda_nm,0.1*lambda_nm,3.0*lambda_nm,0.1*lambda_nm], # lambda (wavelength) in nm\n"," [-R,0.1*R,10*R,0.05*R], # x scale\n"," [1.0,0.01,50.0,0.05]] # y scale\n","\n","# Create the interactive plot\n","interactive_plot = plot_function_with_sliders(function,\n"," labels,\n"," x_min, y_min,\n"," ref_x,max_y,\n"," value_scale)\n","\n","# Display the interactive plot\n","interactive_plot"],"metadata":{"id":"-yGvE2xSXuHI"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# 2D diffraction pattern (analytic)"],"metadata":{"id":"40xc5XwuN2dy"}},{"cell_type":"code","source":["def diff_pattern_2D():\n"," Xp=200\n"," Yp=200\n"," m=5\n"," xpm = m*1.22*P*lambda_nm*1e-9/(2.0*R)\n"," ypm = xpm\n","\n"," diff_pattern = []\n"," for i in range(Xp):\n"," xp=xpm * (Xp - 2 * i) / Xp - xs +1e-5\n"," for j in range(Yp):\n"," yp=ypm * (Yp - 2 * j) / Yp - ys +1e-5\n"," I_analytic=diffraction_intensity_circle_analytic(xp, yp, P, R, lambda_nm)\n"," diff_pattern.append([xp,yp,I_analytic])\n"," return np.array(diff_pattern)"],"metadata":{"id":"qm2CKK5dro-J"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["diff_pattern_2D=diff_pattern_2D()\n","print(diff_pattern_2D.shape)"],"metadata":{"id":"VRPmse2otTpQ"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["import numpy as np\n","import matplotlib.pyplot as plt\n","\n","# Generate some sample data (replace this with your numpy array)\n","data =diff_pattern_2D\n","\n","# Extract x, y, and the value from the data\n","x = data[:, 0]\n","y = data[:, 1]\n","value = data[:, 2]/np.max(data[:, 2])\n","\n","# Reshape the value to a 2D grid for contour plotting\n","x_unique = np.unique(x)\n","y_unique = np.unique(y)\n","X, Y = np.meshgrid(x_unique, y_unique)\n","Z = value.reshape(len(y_unique), len(x_unique))\n","\n","# Create a contour plot\n","plt.figure(figsize=(6, 6))\n","max_value=np.max(Z)/100\n","contour = plt.contourf(X, Y, Z, levels=500, cmap='viridis',\n"," vmax=max_value)\n","plt.colorbar(contour) # Use scientific notation for colorbar\n","\n","\n","plt.xlabel('X')\n","plt.ylabel('Y')\n","plt.title('Intensity (normalized to unity)')\n","plt.gca().set_aspect('equal') # Set equal aspect ratio\n","plt.show()\n"],"metadata":{"id":"i1N1KpD4NaGo"},"execution_count":null,"outputs":[]}]}