Page inachevée, encore en travaux...
xxxxxxxxxximport numpy as npimport numpy.ma as mafrom scipy.special import j1 # Fonction de Bessel d'ordre 1import matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Ddef IdiffCirc(r,d,l,f): """ Intensité I(r) diffractée par une pupille circulaire de diamètre d (longueur d'onde l, focale projection f) """ u = np.pi*d/(l*f)*10 # en cm-1 return 4*(j1(u*r)/(u*r))**2def IdiffRect(x,y,a,b,l,f): """ Intensité I(x,y) diffractée par une pupille rectangulaire de dimensions a, b (longueur d'onde l, focale projection f) """ k = 1/(l*f)*10 u = k*a # en cm-1 v = k*b # en cm-1 return (np.sinc(u*x))**2*(np.sinc(v*y))**2def nm2rgb(wl, gamma=0.8): ''' Conversion longueur d'onde -> RGB D'après Dan Bruton http://www.physics.sfasu.edu/astro/color/spectra.html and http://www.noah.org/wiki/Wavelength_to_RGB_in_Python ''' if wl >= 380 and wl <= 440: attenuation = 0.3 + 0.7 * (wl - 380) / (440 - 380) ; R = ((-(wl-440)/(440-380))*attenuation)**gamma ; G = 0 ; B = (1*attenuation)**gamma elif wl >= 440 and wl <= 490: R = 0 ; G = ((wl-440)/(490-440))**gamma ; B = 1 elif wl >= 490 and wl <= 510: R = 0 ; G = 1 ; B = (-(wl-510)/(510-490))**gamma elif wl >= 510 and wl <= 580: R = ((wl-510)/(580-510))**gamma ; G = 1 ; B = 0 elif wl >= 580 and wl <= 645: R = 1 ; G = (-(wl-645)/(645-580))**gamma ; B = 0 elif wl >= 645 and wl <= 750: attenuation = 0.3 + 0.7 * (750 - wl) / (750 - 645) ; R = (1*attenuation)**gamma ; G = 0 ; B = 0 else: R = 0 ; G = 0 ; B = 0 R *= 255 ; G *= 255 ; B *= 255 return np.array([int(R), int(G), int(B)])print('Initialisations ok')xxxxxxxxxxcourbes = ['Figure de diffraction', 'Intensité', 'Courbe I(r)'](layout=[['trace'],['d','lo'],['xm', 'Im', 'auto']])def _(trace = selector(courbes,label="Tracé"), d=slider(10,100,10,50,label="Diamètre d (µm)"), lo=slider(400,800,10,550,label="Longueur d'onde (nm)"), xm = input_box(np.round(1,1),label="xmax", type=str), Im = input_box(np.round(0.02,3),label="Imax", type=str), auto = checkbox(True,label="Echelle auto"), mesures = checkbox(False,label="Mesures")): #auto_update=false plt.close('all') ; plt.clf() d = int(d) ; lo = int(lo) ; xm = float(xm) ; Im = float(Im) f = 1 # Focale projection en m if auto: xm = 2*np.round(lo*f/d/10,1) # Demi largeur écran en cm p = xm/200 # Pas en cm x0 = np.arange(-xm,xm+p,p) ; x, y = np.meshgrid(x0,x0) r = 2*np.sqrt(x**2 + y**2) ; r = ma.masked_values(r, 0) Idiff = IdiffCirc(r,d,lo,f) ; Idiff = Idiff.filled(1) if trace in ['Figure de diffraction', 'Intensité']: Idiffmask = ma.masked_greater_equal(Idiff, Im) ; Idiffmask=Idiffmask.filled(Im)/Im rgb = nm2rgb(lo,0.8) l,h = Idiffmask.shape z = np.zeros((l,h,3)) for i in range(h): for j in range(l): z[i,j] = Idiffmask[i,j] * rgb z = ma.asarray(z,dtype='uint8') fig,(ax) = plt.subplots(1,1) if trace in['Figure de diffraction','Intensité']: if mesures: ax.grid() nx = 2*xm/p ; nticks = 10 xticks = np.linspace(0,nx,nticks+1) xticklabels = ["{:.1f}".format(i*xm*2) for i in range(-nticks//2,nticks//2+1)] ax.set_xticks(xticks) ; ax.set_yticks(xticks) ax.set_xticklabels(xticklabels) ; ax.set_yticklabels(xticklabels) ax.set_xlabel('x en mm') ; ax.set_ylabel('y en mm') else: ax.axis('off') ax.imshow(z,aspect='equal',origin='lower') elif trace == 'Courbe I(r)': x0m = ma.masked_values(x0, 0) Idiff0 = IdiffCirc(x0m**2,d,lo,f) ; Idiff0 = Idiff0.filled(1) Idiff0 = ma.masked_greater(Idiff0, Im) ax.plot(x0,Idiff0,color='red') plt.show() if trace == 'Intensité': r0, phi, z0 = var('r0 phi z0') trans=(r0*cos(phi),r0*sin(phi),z0) C=plot3d(lambda r0,phi : IdiffCirc(r0,d,lo,f),(r0,xm/1000,xm),(phi,0,2*pi),transformation=trans,opacity=0.87) show(C, aspect_ratio=(1,1,2))xxxxxxxxxxcourbes = ['Figure de diffraction', 'Intensité', 'Courbe I(x)'](layout=[['trace'],['a','b'],['lo'],['xm', 'ym', 'Im', 'auto']])def _(trace = selector(courbes,label="Tracé"), a=slider(10,300,10,50,label="Largeur a (µm)"), b=slider(10,300,10,0,label="Hauteur b (µm)"), lo=slider(400,800,10,550,label="Longueur d'onde (nm)"), xm = input_box(np.round(1,1),label="xmax", type=str), ym = input_box(np.round(1,1),label="ymax", type=str), Im = input_box(np.round(0.02,3),label="Imax", type=str), auto = checkbox(True,label="Echelle auto"), mesures = checkbox(False,label="Mesures")): #auto_update=false plt.close('all') ; plt.clf() a = int(a) ; b =int(b) ; lo = int(lo) ; xm = float(xm) ; ym = float(ym) ; Im = float(Im) f = 1 # Focale projection en m if auto: xm = 6*np.round(lo*f/a/10,1) # Demi largeur écran en cm ym = 6*np.round(lo*f/b/10,1) p = xm/200 # Pas en cm x0 = np.arange(-xm,xm+p,p) ; y0 = np.arange(-ym,ym+p,p) ; x, y = np.meshgrid(x0,y0) Idiff = IdiffRect(x,y,a,b,lo,f) if trace in ['Figure de diffraction', 'Intensité']: Idiffmask = ma.masked_greater_equal(Idiff, Im) ; Idiffmask=Idiffmask.filled(Im)/Im rgb = nm2rgb(lo,0.8) l,h = Idiffmask.shape z = np.zeros((l,h,3)) for i in range(h): for j in range(l): z[j,i] = Idiffmask[j,i] * rgb z = ma.asarray(z,dtype='uint8') fig,(ax) = plt.subplots(1,1) if trace in['Figure de diffraction','Intensité']: if mesures: ax.grid() nx = 2*xm/p ; ny = 2*ym/p ; nticks = 10 xticks = np.linspace(0,nx,nticks+1) ; yticks = np.linspace(0,ny,nticks+1) xticklabels = ["{:.1f}".format(i*xm*2) for i in range(-nticks//2,nticks//2+1)] yticklabels = ["{:.1f}".format(i*ym*2) for i in range(-nticks//2,nticks//2+1)] ax.set_xticks(xticks) ; ax.set_yticks(yticks) ax.set_xticklabels(xticklabels) ; ax.set_yticklabels(yticklabels) ax.set_xlabel('x en mm') ; ax.set_ylabel('y en mm') else: ax.axis('off') ax.imshow(z,aspect='equal',origin='lower') elif trace == 'Courbe I(x)': Idiff0 = IdiffRect(x0,0,a,b,lo,f) Idiff0 = ma.masked_greater(Idiff0, Im) ax.plot(x0,Idiff0,color='red') plt.show() if trace == 'Intensité': X, Y, z0 = var('X Y z0') C=plot3d(lambda X,Y : IdiffRect(X,Y,a,b,lo,f),(X,-xm,xm),(Y,-ym,ym),opacity=0.87) show(C, aspect_ratio=(1,1,2))xxxxxxxxxx