• 几何光学与 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坐标表示，这个更符合人的观察

    • 单个矩形孔衍射

     

    • 随机孔衍射

          

  • 相关阅读:
    我的NVIDIA开发者之旅-Jetson Nano 2gb教你怎么训练模型（完整的模型训练套路）
    Excel - 选择性粘贴和单元格引用规则
    SQL for Data Science Quiz Week1
    遥感和GEE不正式告别
    Lazarus网络编程
    NCP81239MNTXG 开关降压/升压控制器，USB 功率传递和 Type-C 应用
    Metabase学习教程：仪表盘-5
    【Python实战】零基础实战教程(一) Hello World!
    力扣刷题记录62.1-----538. 把二叉搜索树转换为累加树
    基于Python的性能优化
  • 原文地址：https://blog.csdn.net/Chandler_river/article/details/126107360