• 几何光学与 Mathematica

    写在前面

            光路是很重要的,用mma可以很好得模拟它.但是很遗憾,mma跑得不快,如果是因为sanction而用不上正版的话...

            嗯,所以会python是很重要的

            如果有能力自己去解一个方程的话,也可以...

            但是题主选用了Scipy  

            至于可视化吗,题主选用了matplotlib

            有些地方为了算的快一点,也用了multiprocessing  和 numpy

    几何光学:在连续折射率介质中进行光线追踪

    • 光线方程\frac{d}{ds}[n(\mathbf{r})\frac{d\mathbf{r}}{ds}]=\bigtriangledown n(\mathbf{r})
      • 模拟光在大气中的折射\left\{\begin{matrix} n(x,y,z)=n(y)=n_0-ay\\ \frac{d}{ds}[n(y(s))*y'(s)]=\frac{\partial n}{\partial y}=-a\\ \frac{d}{ds}[n(y(s))*z'(s)]=\frac{\partial n}{\partial z}=0 \end{matrix}\right.
      • 模拟光在光纤中的传播
        • 纤芯的折射率模型n(r)=n_1[1-\frac{n_1-n_2}{n_1}(\frac{r}{a})^2]
        • 纤芯中心折射率n1  包层区域折射率n2
          • 对于单模光纤  n2/n1  通常在99.4%~99.7%
          • 对于多模光纤  n2/n1  通常在98%~99%
        • 光线围绕纤芯中心周期性地向前传播,表现出自聚焦的特性
        • 光线路径的周期大小与初始入射角有关,不同的光线将不严格聚集在相同的点上,一个入射的光脉冲将逐渐弥漫散开

    几何光学:在折射率跃变介质中进行光线追踪

    • Snell 定律
    • 球面凸透镜
      • 球面方程
        • f_1(x,y)=R^2-[y^2+(x-R+a)^2]=0
      • 透镜半径
        • r=\sqrt{a(2R-a)}
      • 曲面的法线方向
        • \Omega =\frac{\bigtriangledown f}{\begin{vmatrix} \bigtriangledown f \end{vmatrix}}
      • 平行光束正入射

        • 光束的粗细相对于凸透镜半径r来说  y方向上宽度为 -0.1r~0.1r
      • 平行光斜入射

      • 轴上光源

      • 轴外光源

      • 凸透镜的近轴焦点

      • 厚透镜的形式焦距的倒数

        • \frac{1}{f}=(n-1)[\frac{1}{r_1}-\frac{1}{r_2}+\delta \frac{n_2-1}{n_2r_1r_2}]

      • 物点经过凸透镜之后的像距

        • 待写

    • 三棱镜

    •  白光的色散和色光的合成

      • 单个三棱镜的色散

      • 两个三棱镜的色散  

    • 消色差透镜(凹凸组合透镜)

    波动光学:光衍射的计算

    • 根据夫琅禾费衍射理论,衍射场的振幅分布为:
      • U(\alpha,\beta)=\iint_{\sum}u(x,y)e^{2\pi i (xsin\alpha + ysin\beta)/\lambda}dxdy
      • u(x,y)衍射屏的透射函数
      • 衍射强度的分布I(\alpha,\beta)=\begin{vmatrix} U(\alpha,\beta) \end{vmatrix}^2
      • U=U_1+iU_2=\iint_{\sum}ucosfdxdy+i\iint_{\sum}sinfdxdy
    • 单个圆孔的衍射

      • u(x,y)=UnitStep[r^2-(x-x_0)^2-(y-y_0)^2]
      • 硬核的来了,这个程序等他跑完也就该下班了,几个小时也不一定跑的完,我们还是换Python取编写.当然,如果读者会C或者Fortran啥的那更好的,不过,写个几百行累死人了,不是吗?
      • 热力图是表征强度的一个很好的工具,
      • 下面给出单个圆孔的夫琅禾费衍射计算程序（以方向角形式展示的，所以他不是一个标砖的圆，这个演示完会有用X，Y坐标表示的程序）
      • 当然,以Python语言的一般特性,这个n取到100的时候,,,,
    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. import seaborn as sns
    4. from scipy.integrate import dblquad
    5. import math
    6.  
    7. #lambda 用k表示
    8. r = 3.0
    9. k = 0.6
    10. q = 2*math.pi / k
    11. #alpha 用a表示 beta用b表示
    12. a0 = math.pi
    13. b0 = math.pi
    14. n = 10
    15.  
    16. def fun(D,k):
    17.     a = D[0]
    18.     b = D[1]
    19.     u = "1"
    20.     x1 = 3
    21.     x2 = -3
    22.     y1 = lambda x: (9-x**2)**0.5
    23.     y2 = lambda x: -(9-x**2)**0.5
    24.     U1 = dblquad(lambda x,y:eval(u)*math.cos(2*math.pi*(x*math.sin(a)+y*math.sin(b))/k)
    25.                  ,x1,x2,y1,y2)
    26.     U2 = dblquad(lambda x,y:eval(u)*math.sin(2*math.pi*(x*math.sin(a)+y*math.sin(b))/k)
    27.                  ,x1,x2,y1,y2)
    28.     U = U1[0]**2 + U2[0]**2
    29.     return U
    30.  
    31. data0 = np.array([[[i,j] for j in np.linspace(-a0,a0,n)] for i in np.linspace(-b0,b0,n)])
    32. data = [[fun(j,k) for j in i] for i in data0]
    33. #注意!
    34. data =  np.array(data)#设置二维矩阵
    35. f,ax = plt.subplots(figsize=(11,12))
    36. #数据
    37.  
    38. f = open("diffraction.txt","w")
    39. f.write(str(data))
    40. f.close()
    41.  
    42. #注意!
    43. sns.heatmap(data, ax=ax,vmin=0,vmax=6,cmap='YlOrRd',linewidths=2,cbar=True)
    44. #热力图绘制代码
    45.  
    46.  
    47. plt.title('heatmap') 
    48. plt.ylabel('y_label')  
    49. plt.xlabel('x_label')
    50. plt.savefig("diffraction.jpg")
    51. plt.show()
    • 这时候必须multiprocessing辅助一下了
      1. import matplotlib.pyplot as plt
      2. import numpy as np
      3. import seaborn as sns
      4. from scipy.integrate import dblquad
      5. import math
      6. import datetime
      7. import multiprocessing as mp
      8. import time
      9.  
      10. a0 = math.pi/10
      11. b0 = math.pi/10
      12. n = 100
      13. k = 0.6
      14.  
      15.  
      16. def final_fun(name, param):
      17.     result = []
      18.     global a0
      19.     global b0
      20.     global n
      21.     global k
      22.     A = np.linspace(-a0,a0,n)
      23.     B = np.linspace(-b0,b0,n)
      24.     for params in param:
      25.         a = A[params[0]]
      26.         b = B[params[1]]
      27.         u = "1"
      28.         x1 = 3
      29.         x2 = -3
      30.         y1 = lambda x: (9-x**2)**0.5
      31.         y2 = lambda x: -(9-x**2)**0.5
      32.         U1 = dblquad(lambda x,y:eval(u)*math.cos(2*math.pi*(x*math.sin(a)+y*math.sin(b))/k)
      33.                      ,x1,x2,y1,y2)
      34.         U2 = dblquad(lambda x,y:eval(u)*math.sin(2*math.pi*(x*math.sin(a)+y*math.sin(b))/k)
      35.                      ,x1,x2,y1,y2)
      36.         U = U1[0]**2 + U2[0]**2
      37.         result.append(U)
      38.     return {name:result}
      39.  
      40.  
      41. if __name__ == '__main__': 
      42.     start_time = time.time()
      43.     num_cores = int(mp.cpu_count())
      44.     pool = mp.Pool(num_cores)
      45.     param_dict = {"task1":[[0,i] for i in range(n)],
      46. "task2":[[1,i] for i in range(n)],
      47. "task3":[[2,i] for i in range(n)],
      48. "task4":[[3,i] for i in range(n)],
      49. "task5":[[4,i] for i in range(n)],
      50. "task6":[[5,i] for i in range(n)],
      51. "task7":[[6,i] for i in range(n)],
      52. "task8":[[7,i] for i in range(n)],
      53. "task9":[[8,i] for i in range(n)],
      54. "task10":[[9,i] for i in range(n)],
      55. "task11":[[10,i] for i in range(n)],
      56. "task12":[[11,i] for i in range(n)],
      57. "task13":[[12,i] for i in range(n)],
      58. "task14":[[13,i] for i in range(n)],
      59. "task15":[[14,i] for i in range(n)],
      60. "task16":[[15,i] for i in range(n)],
      61. "task17":[[16,i] for i in range(n)],
      62. "task18":[[17,i] for i in range(n)],
      63. "task19":[[18,i] for i in range(n)],
      64. "task20":[[19,i] for i in range(n)],
      65. "task21":[[20,i] for i in range(n)],
      66. "task22":[[21,i] for i in range(n)],
      67. "task23":[[22,i] for i in range(n)],
      68. "task24":[[23,i] for i in range(n)],
      69. "task25":[[24,i] for i in range(n)],
      70. "task26":[[25,i] for i in range(n)],
      71. "task27":[[26,i] for i in range(n)],
      72. "task28":[[27,i] for i in range(n)],
      73. "task29":[[28,i] for i in range(n)],
      74. "task30":[[29,i] for i in range(n)],
      75. "task31":[[30,i] for i in range(n)],
      76. "task32":[[31,i] for i in range(n)],
      77. "task33":[[32,i] for i in range(n)],
      78. "task34":[[33,i] for i in range(n)],
      79. "task35":[[34,i] for i in range(n)],
      80. "task36":[[35,i] for i in range(n)],
      81. "task37":[[36,i] for i in range(n)],
      82. "task38":[[37,i] for i in range(n)],
      83. "task39":[[38,i] for i in range(n)],
      84. "task40":[[39,i] for i in range(n)],
      85. "task41":[[40,i] for i in range(n)],
      86. "task42":[[41,i] for i in range(n)],
      87. "task43":[[42,i] for i in range(n)],
      88. "task44":[[43,i] for i in range(n)],
      89. "task45":[[44,i] for i in range(n)],
      90. "task46":[[45,i] for i in range(n)],
      91. "task47":[[46,i] for i in range(n)],
      92. "task48":[[47,i] for i in range(n)],
      93. "task49":[[48,i] for i in range(n)],
      94. "task50":[[49,i] for i in range(n)],
      95. "task51":[[50,i] for i in range(n)],
      96. "task52":[[51,i] for i in range(n)],
      97. "task53":[[52,i] for i in range(n)],
      98. "task54":[[53,i] for i in range(n)],
      99. "task55":[[54,i] for i in range(n)],
      100. "task56":[[55,i] for i in range(n)],
      101. "task57":[[56,i] for i in range(n)],
      102. "task58":[[57,i] for i in range(n)],
      103. "task59":[[58,i] for i in range(n)],
      104. "task60":[[59,i] for i in range(n)],
      105. "task61":[[60,i] for i in range(n)],
      106. "task62":[[61,i] for i in range(n)],
      107. "task63":[[62,i] for i in range(n)],
      108. "task64":[[63,i] for i in range(n)],
      109. "task65":[[64,i] for i in range(n)],
      110. "task66":[[65,i] for i in range(n)],
      111. "task67":[[66,i] for i in range(n)],
      112. "task68":[[67,i] for i in range(n)],
      113. "task69":[[68,i] for i in range(n)],
      114. "task70":[[69,i] for i in range(n)],
      115. "task71":[[70,i] for i in range(n)],
      116. "task72":[[71,i] for i in range(n)],
      117. "task73":[[72,i] for i in range(n)],
      118. "task74":[[73,i] for i in range(n)],
      119. "task75":[[74,i] for i in range(n)],
      120. "task76":[[75,i] for i in range(n)],
      121. "task77":[[76,i] for i in range(n)],
      122. "task78":[[77,i] for i in range(n)],
      123. "task79":[[78,i] for i in range(n)],
      124. "task80":[[79,i] for i in range(n)],
      125. "task81":[[80,i] for i in range(n)],
      126. "task82":[[81,i] for i in range(n)],
      127. "task83":[[82,i] for i in range(n)],
      128. "task84":[[83,i] for i in range(n)],
      129. "task85":[[84,i] for i in range(n)],
      130. "task86":[[85,i] for i in range(n)],
      131. "task87":[[86,i] for i in range(n)],
      132. "task88":[[87,i] for i in range(n)],
      133. "task89":[[88,i] for i in range(n)],
      134. "task90":[[89,i] for i in range(n)],
      135. "task91":[[90,i] for i in range(n)],
      136. "task92":[[91,i] for i in range(n)],
      137. "task93":[[92,i] for i in range(n)],
      138. "task94":[[93,i] for i in range(n)],
      139. "task95":[[94,i] for i in range(n)],
      140. "task96":[[95,i] for i in range(n)],
      141. "task97":[[96,i] for i in range(n)],
      142. "task98":[[97,i] for i in range(n)],
      143. "task99":[[98,i] for i in range(n)],
      144. "task100":[[99,i] for i in range(n)]}
      145.     
      146.     results = [pool.apply_async(final_fun, args=(name, param)) for name, param in param_dict.items()]
      147.     results = [p.get() for p in results]
      148.     end_time = time.time()
      149.     use_time = end_time - start_time
      150.     print("多进程计算 共消耗: " + "{:.2f}".format(use_time) + " 秒")
      151.     data = [i["task"+str(results.index(i)+1)] for i in results]
      152.     data = np.array(data)
      153.     f = open("diffraction.txt","w")
      154.     f.write(str(data))
      155.     f.close()
      156.     fig,ax = plt.subplots(figsize=(11,12))
      157.     sns.heatmap(data, ax=ax,vmin=0,vmax=6,cmap='YlOrRd',\
      158.         linewidths=2,cbar=True)
      159.  
      160.  
      161.     plt.title('diffraction') 
      162.     plt.ylabel('y_label')  
      163.     plt.xlabel('x_label')   
      164.     plt.savefig("diffraction.jpg")
      165.     plt.show()

    用x，y坐标表示，这个更符合人的观察

    • 单个矩形孔衍射

     

    • 随机孔衍射

          

  • 相关阅读:
    AntBlazor Theme in ABP Framework
    C练题笔记之：Leetcode-654. 最大二叉树
    消息队列RabbitMQ核心：简单（Hello World）模式、队列（Work Queues）模式、发布订阅模式
    第一章 Google软件测试介绍
    首发Yolov8涨点神器：华为诺亚2023极简的神经网络模型 VanillaNet---VanillaBlock助力检测，实现暴力涨点
    CSDN 编程竞赛第七期题解
    3D建模基础教程：石墨工具介绍
    AI无处不在，科技改变生活：开放原子全球开源峰会参会感悟
    【软件测试】毕业失眠的夜晚，望着天花板呆住了 到底该何去何从......
    @ConditionalOnProperty注解作用
  • 原文地址：https://blog.csdn.net/Chandler_river/article/details/126107360