基于折射率方向导数的自由梯度折射率介质的自适应光线追迹方法

本发明属于光学设计领域，特别是一种基于折射率方向导数的自由梯度折射率介质的自适应光线追迹方法。

背景技术：

1、近年来研究人员不断地拓展设计自由度，以在更复杂的光学系统中实现更高的光学质量。随着离子交换等工艺的发展，梯度折射率形式的材料越来越广泛地应用到光纤通信和成像光学等领域。另外，增材制造技术的发展使梯度折射率的形式可以变得更为复杂，向自由梯度折射率发展，为光学设计增加了折射率分布这一自由度。为了在设计过程中恰当地评价系统光学质量，自由梯度折射率的追迹方法是非常重要的。尽管可以使用解析的方式获得几种特殊形式梯度折射率的光线追迹结果，但对于自由梯度折射率这种复杂形式，解析法并不适用。现有光学设计软件使用数值计算方法逐步进行光线追迹，尽管可以通过减小步长来使结果精度提高，但过小的步长不仅极大地增加了计算耗时，并且精度上收益甚微。在光学设计过程中，像质评价和优化迭代都存在大量的光线追迹过程，因此需要光学设计者平衡计算效率与结果精度，通过反复调整获得较为合适的步长，但是这样会导致设计效率很低。

技术实现思路

1、本发明的目的在于针对上述现有技术存在的问题，提供一种用于自由梯度折射率介质的高效、高精度的光线追迹方法。

2、实现本发明目的的技术解决方案为：一种基于折射率方向导数的自由梯度折射率介质的自适应光线追迹方法，所述方法包括以下步骤：

3、步骤1，根据折射率分布n(r)在光线初始位置r0沿光线初始传播方向d0的折射率方向导数k，以及器件长度l，计算获得初始步长δh0，计算次数m＝1；

4、步骤2，使用步长δhm-1依据光线方程求解光线前进一步后的位置rm和方向dm；

5、步骤3，判断当前位置是否超出器件边界，若超出则执行步骤7，否则执行步骤4；

6、步骤4，计算折射率分布在步骤2获得的位置rm和方向dm的折射率变化系数cm，获得步长修正因子α，并求得新的步长δhm；

7、步骤5，判断δhm是否超出预设范围，若超出则进行修正，计算次数m＝m+1；

8、步骤6，循环步骤2至步骤5直至光线追迹超出器件范围；

9、步骤7，对光线最后一步的步长进行迭代修正，直至最终出射点的位置误差小于阈值。

10、进一步地，步骤1中计算获得初始步长δh0，具体过程包括：

11、步骤1-1，依据光线位置和方向，计算折射率在初始位置r0沿光线初始传播方向的方向导数：

12、

13、式中，n为三维分布形式的自由梯度折射率，r0＝(x0,y0,z0)为光线初始位置，d0＝(cosα0，cosβ0，cosγ0)为光线初始方向余弦，用于表示光线方向，其中α0，β0和γ0为光线与光轴正方向的夹角，k0为折射率在光线初始位置处沿传播方向的方向导数；

14、步骤1-2，根据所述方向导数获取折射率变化系数，计算公式为：

15、

16、式中，c0为折射率变化系数；

17、步骤1-3，根据所述折射率变化系数获取初始步长，计算公式为：

18、δh0＝(1.45c02-0.73c0+0.10)l-0.37

19、式中，δh0为初始步长，l为器件长度。

20、进一步地，步骤2所述使用步长δhm-1依据光线方程求解光线前进一步后的位置rm和方向dm，具体过程包括：

21、光线方程为：

22、

23、式中，r＝r(s)是光线上任一点p的位置矢量，并且为光线弧长s的函数；

24、定义变量dh＝ds/n，将上式重写为：

25、

26、引入光线矢量g＝ndr/ds＝(cosα,cosβ,cosγ)n，即g＝nd，将上式转化为一对与r和g相关的一阶常微分方程组：

27、

28、

29、采用微分方程数值解法求解所述一阶常微分方程组，则光线前进一步后的位置rm和光线矢量gm的计算公式为：

30、

31、gm＝gm-1+ψ(rm-1,gm-1,δhm-1,n)δhm-1

32、式中，和ψ均为增量函数；

33、由此获得：dm＝gm/n。

34、进一步地，步骤3中所述判断当前位置是否超出器件边界，具体过程包括：

35、步骤3-1，将步骤2中获得的光线位置rm＝(xm,ym,zm)的(xm,ym)代入后表面面型位置公式：

36、zsurface＝f(x,y)+l

37、式中，xm，ym，zm分别为光线位置的x，y，z坐标，f(x,y)为后表面的面型表达式，zsurface为后表面在(xm,ym)位置处的z坐标；

38、步骤3-2，比较zm与zsurface，若zm＞zsurface，则说明当前光线位置已超出器件边界。

39、进一步地，步骤4所述计算折射率分布在步骤2获得的位置rm和方向dm的折射率变化系数cm，获得步长修正因子α，并求得新的步长δhm，具体包括：

40、步骤4-1，计算折射率分布在光线当前位置rm和方向dm的折射率变化系数cm，计算公式为：

41、

42、式中，km为折射率在光线当前位置rm处沿传播方向的方向导数；

43、步骤4-2，使用cm计算步长修正因子α，计算公式为：

44、

45、式中，q＝|cm|-|cm-1|，由q构成的指数部分用于指示每一步δh的修正方向；cm-1为折射率分布在前一位置rm-1和前一方向dm-1的折射率变化系数；

46、步骤4-3，计算新的步长δhm，计算公式为：

47、δhm＝αδhm-1。

48、进一步地，步骤5中所述判断δhm是否超出预设范围，具体过程包括：

49、根据步骤4所得的步长δhm，判断其大小并进行修正，判断条件为：

50、0.001≤δhm≤0.1mm

51、若δhm符合上述判断条件则不进行修改；若δhm＞0.1mm则令δhm＝0.1mm，若δhm＜0.001mm则令δhm＝0.001mm。

52、进一步地，步骤7中所述对光线最后一步的步长进行迭代修正，直至最终出射点的位置误差小于阈值，具体过程包括：

53、步骤7-1，计算步长误差e：

54、e＝zm-zsurface

55、式中，zm为光线位置rm的z坐标，zsurface为后表面在(xm,ym)位置处的z坐标；

56、步骤7-2，利用所述步长误差e对光线最后一步的步长进行修正：

57、δh′last＝δhlast-e

58、式中，δhlast、δh′last分别是修正前和修正后的最后一步的步长；

59、步骤7-3，使用修正后的步长，利用步骤2的过程重新计算最后一步的光线位置rlast和光线矢量glast，迭代计算，直至满足下式：

60、|e|≤ε

61、式中，ε为阈值。

62、本发明与现有技术相比，其显著优点为：

63、1)本发明可以根据折射率的变化自适应的调整光线追迹过程中的步长，解决了传统方法中需要设计者手动设置并判断步长是否合理的问题。

64、2)本发明能够自适应修正每一步的步长，能够在逐步的光线追迹过程中判断沿当前光线的传输方向的折射率变化情况对步长做出合理修正。

65、3)本发明可以大幅提高计算效率，为自由梯度折射率光学系统优化设计带来便利。

66、4)本发明对于不同长度的器件计算精度都可以达到使用需求，并且能够在同样精度的情况下减少计算次数。

67、下面结合附图对本发明作进一步详细描述。