前言

  初写这个程序时，用的是matlab。然而离开学校以后，价格高昂的matlab软件不再能使用。因此，用Python补充这些程序。如果有读者想要讨论电磁场理论相关的物理学或者编程以及数学，我的邮箱是707101557@qq.com。

程序

  下面两个程序用简单的方法实现仿真，算法的时间复杂度极高，且没有在编程技巧上进行改进，因此需要运行较长时间。在最后，我改进了变成技巧，加速了程序的运行

python程序

#---------------^_^-------------^_^------- #利用python仿真双缝实验 #程序再版与2019/1/5，初版写于2017年12月 #作者：cclplus #仅供学习交流使用 #如有疑问或者需求，可以联系作者707101557@qq.com import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import math n=int(50) #方镜上点的个数 rd=0.005#方镜的边长——单位m ld=1.0 #确定谐振腔的腔长 #假设起始时光强处处相等 I=np.ones((n,n),float) #这是一个需要长时间运行的程序 for t in range(30): In=np.zeros((n,n),float) for i in range(n): for j in range(n): for i1 in range(n): for j1 in range(n): In[i][j]=In[i][j]+1.0/(4.0*math.pi)*I[i1][j1]* (1.0+ld/math.sqrt((rd/float(n)*(i1-i))**2+ (rd/float(n)*(j1-j))**2+ld**2)) /float(n)/float(n) I=In*[1.0] S_I=sum(sum(I[i]) for i in range(len(I))) mul=float(n)*float(n)/S_I I*=[mul] I*=I fig = plt.figure() ax = Axes3D(fig) X=range(n) Y=range(n) ax.plot_surface(X, Y, I, rstride=1, cstride=1, cmap='rainbow') plt.show()

matlab程序

%作者：cclplus %时间:2017年12月 clc clear all; close all; %%%%%以下程序用于建立反射镜的模型top %%%方形镜 n=50;%取点,确定精度 rd=0.005;%确定方形镜的边长 ld=1;%确定谐振腔的长度 %%%%以上程序用于建立反射镜的模型bottom I=ones(n,n);%假设初始光强处处相等 for t=1:50 In=zeros(n,n); for i=1:n for j=1:n for i1=1:n for j1=1:n In(i,j)=In(i,j)+1/(4*pi)*I(i1,j1)*(1+ld/sqrt((rd/n*(i1-i))^2+(rd/n*(j1-j))^2+ld^2))*1/2500; end end end end I=In; end %%%bottom S_I=sum(sum(In)'); multipe=n^2/S_I; I=multipe.*In; %当找到有一组不符合的数据，即可跳出 %绘制图像 x=1:n; y=1:n; I=I.^2; surf(x,y,I); shading interp axis equal ##实验结果图
  受衍射的影响，激光的能量逐渐集中

对程序的改进

  我想有的人通过Python执行这段程序的时候，估计要怀疑计算机死机了，这确实是个问题，我在第一次写作时确实是等待了很长的运行时间，然后再把图片保存下来。
想要改进这个程序，需要采用动态链接库(dll)的方法 #---------------^_^-------------^_^------- #利用python仿真双缝实验 #程序再版与2019/1/5，初版写于2017年12月 #作者：cclplus #仅供学习交流使用 #如有疑问或者需求，可以联系作者707101557@qq.com import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import math from ctypes import * #dll = cdll.LoadLibrary('laser.dll'); dll=WinDLL('laser.dll') class StructPointer(Structure): _fields_ = [("arr", c_double * 2500)] def __init__( self ): for i in range(2500): self.arr[i]=0.0 dll.laser_mode.restype = POINTER(StructPointer) n=c_long(50) #方镜上点的个数 n_b=int(50) k=c_long(30) #ccl=np.ones((n,n),float) ccl=StructPointer() #dll.hello() ccl=dll.laser_mode(n,k) I=np.zeros((n_b,n_b),float) for i in range(n_b): for j in range(n_b): I[i][j]=ccl.contents.arr[i*n_b+j] X=range(n_b) Y=range(n_b) fig = plt.figure() ax = Axes3D(fig) ax.plot_surface(X,Y,I, rstride=1, cstride=1, cmap='rainbow') plt.show() 想要获取更详细的相关资料可以访问我的github主页https://github.com/YuruTu/laser

算法上的改进

  一力降十会，在绝对力量的面前，再多的技巧也只能黯然失 {MOD}。很多时候改进算法才是改良计算机程序的根本。个人认为真正的简化这个过程在于数学公式的推导。曹三松在《稳定腔激光模式理论的再研究》中推导了这个过程，并得出结论
共焦腔镜面上激光场振幅的表达式为 其中 C表示归一化系数，H表示Hermite多项式
光强I正比于光场振幅的平方
对推导过程持怀疑态度的读者可以联系本文作者cclplus，我会竭力为您解答
以下代码为小斐而写（如果有需要，我可以把这个程序打包成GUI模式） #写作时间：2019年2月26日 #为小斐而写 #作者cclplus import sympy#sympy的安装过程特殊，必须通过先在官网下载.gz文件 import matplotlib.pyplot as plt import math import numpy as np from mpl_toolkits.mplot3d import Axes3D eps = 0.001 #定义一个函数用于求Hermite多项式的平方值 #往后求导过程比较繁琐，暂不加以推广 def heimite(n,x): if n==0: return 1 elif n==1: return 4.0*x*x elif n==2: return (4.0*x*x-2.0)*(4.0*x*x-2.0) elif n==3: return (12.0*x-8.0*x**3.0)**2.0 else: print("暂时不支持计算比3更大的Hermite多项式的平方值 联系作者获得相关支持") if __name__ == "__main__": #没有用上全部参数,但也足以说明问题 num_plot = 200 m = 3 n = 3 R=np.zeros((num_plot,num_plot),float) for i in range( num_plot): for j in range( num_plot): x =float(i-num_plot/2)/float(num_plot/10) y =float(j-num_plot/2)/float(num_plot/10) R[i][j] =heimite(m,x)*heimite(n,y)*sympy.exp(-2.0*x*x-2.0*y*y) fig = plt.figure() X=range(num_plot) Y=range(num_plot) plt.contourf(X, Y, R, 8, alpha = 0.75, cmap = plt.cm.hot) plt.show() 结果