(57)【要約】【課題】屈折率異方性がある媒質内でも、gradient  index光線方程式の微分の部分を差分近似せずに、高速度で計算する。【解決手段】  屈折率が一様でなく屈折率異方性のある媒質内での光線追跡方法であって、ｎ’を光線屈折率とするとき、屈折率異方性のある媒質内でのgradient index光線方程式n'd2x/dz2=[1+(dx/dz)2+(dy/dz)2][(∂n'/∂x)-(∂n'/∂z)(dx/dz)n'd2y/dz2=[1+(dx/dz)2+(dy/dz)2][(∂n'/∂x)-(∂n'/∂z)(dy/dz)について、ｘ、ｙを光軸方向ｚの関数として冪級数展開することにより解く。ただし、ｘ、ｙは互いに直交し且つ前記光軸方向ｚと直交する方向である。【選択図】図１
2006003287.tif 000002
【特許請求の範囲】
【請求項１】
屈折率が一様でなく屈折率異方性のある媒質内での光線追跡方法であって、ｎ’を光線屈折率とするとき、屈折率異方性のある媒質内でのgradient index光線方程式
【数１】
2006003287.tif 000094
【数２】
2006003287.tif 000095
について、ｘ、ｙを光軸方向ｚの関数として冪級数展開することにより解くことを特徴とする光線追跡方法。
ただし、ｘ、ｙは互いに直交し且つ前記光軸方向ｚと直交する方向である。
【請求項２】
前記光線追跡方法は、
ｘ、ｙの冪級数展開の係数Ａi，Ｂi（ｉ＝０，．．．Ｋ＋１）から、
「光線方向に直交する初期ベクトルの成分に関する冪級数展開の係数」を求め、
次に、「初期ベクトルの回転に関する冪級数展開の係数」を求め、
次に、「電気ベクトルの成分に関する冪級数展開の係数」を求め、
次に、「電気ベクトルの成分の積に関する冪級数展開の係数」を求め、
次に、「光線屈折率およびその微分に関する冪級数展開の係数」を求め、
次に、「gradient index 光線方程式を解く前段階の冪級数展開の係数」を求め、
その係数と、係数Ａi，Ｂi（ｉ＝０，...，ｋ＋１）より、高次の係数Ａk+2，Ｂk+2を計算することを特徴とする請求項１に記載の光線追跡方法。
【請求項３】
請求項１または２に記載の光線追跡方法を記憶したことを特徴とするコンピュータ読み取り可能な記憶媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、偏光を利用する光学系または、半導体素子製造用の投影光学系や、天体観測用の光学系等において、特に光線通過部材の、加熱や外力による内部応力分布の変動による光線の偏光分離が問題となるような光学系の解析や、電離層のように異方性不均質な磁気プラズマ中を伝搬する光線の解析に関するものである。
【背景技術】
【０００２】
従来、複屈折がある硝材を光線が透過する場合の光線追跡は、例えば、１軸結晶を光線が入射する場合の屈折光線の屈折率と波面放線方向を求める方法については、M．C．Simonによって、“Ray tracing formulas for monoaxial optical components“という題名の文献で紹介されている。（Appl．Opt．(1983)pp354-360）
また、一般的に２軸結晶に光線が入射する場合の屈折光線の屈折率と波面法線ベクトルを求める方法の原理については、ボルン・ウオルフ著：「光学の原理III」（邦訳：東海大学出版会）において、Fresnelの公式（同文献、P980参照）を利用して求める方法が示されている。（同文献、P999〜P1000参照）
しかし、ガラスの温度上昇分布に差がある場合や、ガラスの一部に外力を加えた場合、ガラス内で歪みが発生し、内部応力が発生するが、 発生する内部応力は、通常方向異方性があるので、屈折率異方性があり、光線がガラスを透過すると、偏光になる。一方、内部応力は、ガラス内の位置により異なっていて、それに従って、屈折率も異なるので、gradient  index光線追跡をする必要がある。
【０００３】
内部応力のある媒質内でのgradient index光線追跡では、計算誤差が蓄積しないように、屈折率の計算は高精度に行う必要があるが、前記のM．C．Simonの文献等で述べられているような屈折率を計算する方法では、gradient index光線追跡に絶えられるような精度が保証できないことがわかった。そこで、本願発明者は、波面法線楕円体によって、屈折率を計算する方法を開発した。（特開平１１−１０１７１４号公報、および橋本純夫：「内部応力分布があるガラス内の光線追跡」光学26(1997) 677を参照）
【特許文献１】特開平１１−１０１７１４号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
しかし、この方法では、内部応力のある媒質内での屈折率は正確に計算できるが、媒質内では、波面法線方向と光線方向が異なるために、光線の微小進行毎に、光線の進行方向を電場ベクトルと磁場ベクトルの外積すなわち、Poynting ベクトルを求めることによって光線方向を修正する必要があった。このような計算方法であると、微小進行毎の修正によって、計算誤差が蓄積する恐れがあった。
【０００５】
そこで、Maxwell方程式から、屈折率異方性がある媒質内でもgradient  index光線追跡が可能な光線方程式（光線ベクトルの方向を直接求める微分方程式）を導出することが考えられる。ここで、その光線方程式は、屈折率異方性が無い場合の光線方程式における屈折率の代わりに、光線屈折率を用いれば良いものである。
【０００６】
しかしながら、光路を有限個に分割して、gradient  index光線方程式の微分の部分を差分近似して求める場合、分割数を多くしないと高精度な計算値が得られないが、分割数を多くすることにより、多大な計算時間が必要であるという問題が生じる。
【０００７】
そこで、本発明は、屈折率異方性がある媒質内でも、gradient  index光線方程式の微分の部分を差分近似せずに、高速度で計算する方法を提供するものである。
【課題を解決するための手段】
【０００８】
上述の目的を達成するため、本発明にかかる光線追跡方法は、屈折率が一様でなく屈折率異方性のある媒質内での光線追跡方法であって、ｎ’を光線屈折率とするとき、屈折率異方性のある媒質内でのgradient index光線方程式
【０００９】
【数３】
2006003287.tif 000003
【００１０】
【数４】
2006003287.tif 000004
【００１１】
について、ｘ、ｙを光軸方向ｚの関数として冪級数展開することにより解くことを特徴とする。
【発明の効果】
【００１２】
本発明によれば、屈折率異方性がある媒質内でも、たとえば数値解析よりも、高速かつ高精度に光線追跡を行うことができる。
【発明を実施するための最良の形態】
【００１３】
まず、屈折率異方性が無い場合においては、光軸をｚ軸とすると、それに垂直なｘ,ｙの光線到達位置をｚの関数としてgradient  index光線方程式を冪級数展開することにより、数値解析よりも高速で高精度に計算できる。（Duncan  T.  Mooreの方法、J．Opt．Soc．Am．６５，４５１（１９７５）を参照）
しかし、この方法を、屈折率異方性がある媒質内でのgradient  index光線方程式に対して適用するのは容易ではない。特に、光線が進行するに従って偏光方向が回転する場合は、かなり困難である。
【００１４】
そこでまず、屈折率異方性の分布が子午面に関して対称で、かつ光線が子午面にのみ存在する場合において、gradient  index光線方程式を冪級数展開する方法を検討する。
子午面をｘｚ面とすると、子午面と光線楕円体の交線の楕円は、
【００１５】
【数５】
2006003287.tif 000005
で表わすことができる。光線の傾きがｍ＝ｄｘ／ｄｚの場合、xz面内で偏光する光線の単位電気ベクトルは、光線の方向に垂直なので、
【００１６】
【数６】
2006003287.tif 000006
【００１７】
と表わすことができる。原点から電気ベクトルと楕円の交点の距離をγとすると、
【００１８】
【数７】
2006003287.tif 000007
【００１９】
となり、この(3)を(1)に代入することにより、光線屈折率は、
【００２０】
【数８】
2006003287.tif 000008
【００２１】
となり、その距離による微分は、
【００２２】
【数９】
2006003287.tif 000009
【００２３】
となる。ｘｚ面に関して対称な場合、屈折率異方性の場合のgradient  index光線方程式
【００２４】
【数１０】
2006003287.tif 000010
【００２５】
に代入して、
【００２６】
【数１１】
2006003287.tif 000011
【００２７】
特に、γxz＝０の場合  γxx＝ｎx2，  γzz＝ｎz2  …(9)
【００２８】
【数１２】
2006003287.tif 000012
【００２９】
さらに、
【００３０】
【数１３】
2006003287.tif 000013
【００３１】
の場合、
【００３２】
【数１４】
2006003287.tif 000014
【００３３】
となる。このように、gradient  index光線方程式を冪級数展開するためには、微分方程式を冪級数展開しやすい形にすることが必要である。
【００３４】
先ずｘをｚで
【００３５】
【数１５】
2006003287.tif 000015
【００３６】
というように冪級数展開した場合、
【００３７】
【数１６】
2006003287.tif 000016
【００３８】
【数１７】
2006003287.tif 000017
【００３９】
となる。
【００４０】
【数１８】
2006003287.tif 000018
【００４１】
とおくと、(15)より、
【００４２】
【数１９】
2006003287.tif 000019
【００４３】
となる。
【００４４】
また、(12)式の場合のように、ｎxがｘのみの関数の場合、
【００４５】
【数２０】
2006003287.tif 000020
【００４６】
と表わすことができる。この場合、
【００４７】
【数２１】
2006003287.tif 000021
【００４８】
とおくと、
【００４９】
【数２２】
2006003287.tif 000022
【００５０】
より、
【００５１】
【数２３】
2006003287.tif 000023
【００５２】
となる。
【００５３】
【数２４】
2006003287.tif 000024
【００５４】
とおくと、
【００５５】
【数２５】
2006003287.tif 000025
【００５６】
となる。従って、
【００５７】
【数２６】
2006003287.tif 000026
【００５８】
とおくと、
【００５９】
【数２７】
2006003287.tif 000027
【００６０】
となる。
また、(12)式の場合のように、ｎzが一定の場合、
【００６１】
【数２８】
2006003287.tif 000028
【００６２】
とおくと、
【００６３】
【数２９】
2006003287.tif 000029
【００６４】
となる。(10)式に、(11)式、(26)式、(28)式を代入し、ｚの係数を比較することにより、
【００６５】
【数３０】
2006003287.tif 000030
【００６６】
というように、(13)式の(k＋2)次の係数を、0次から(k＋1)次までの係数により表せるので、低次の係数から順に高次の係数を求めることができる。
【００６７】
次に、屈折率異方性の分布が子午面に関して非対称な場合において、gradient  index光線方程式を冪級数展開する方法を述べる。その場合、「波面法線楕円体との双対性により、光線方向のベクトルｔと直交する平面と、光線楕円体との交線である楕円の2つの軸の方向が電気ベクトルＥの方向であり、軸の長さが、光線屈折率の逆数に比例する」[ボルン・ウオルフ著：「光学の原理III,P988」（邦訳：東海大学出版会）]ことを利用する。
【００６８】
まず、光線方向に垂直で、互いに垂直な2ベクトルを選定する。光線方向をベクトルｔ＝（ｔx，ｔy，ｔz）とすると、それに垂直で、かつ互いに垂直な初期ベクトルｐ，ｑの成分は、例えば、
【００６９】
【数３１】
2006003287.tif 000031
【００７０】
【数３２】
2006003287.tif 000032
【００７１】
となる。光線ベクトルの成分を微分で表わすと、
【００７２】
【数３３】
2006003287.tif 000033
【００７３】
となるので、
【００７４】
【数３４】
2006003287.tif 000034
【００７５】
と置くと、初期ベクトルｐ，ｑの成分は、
【００７６】
【数３５】
2006003287.tif 000035
【００７７】
となる。そこで、
【００７８】
【数３６】
2006003287.tif 000036
【００７９】
と置くと、ベクトルｐ，ｑの微分形式には平方根は含まれないので、後の計算に便利である。
【００８０】
次に、ベクトルｐ，ｑをθ回転させることにより電気ベクトルの方向になる場合、その方向は、光線方向のベクトルと直交する平面と、光線楕円体との交線である楕円の2つの軸の方向であるので、
光線楕円体が、
【００８１】
【数３７】
2006003287.tif 000037
【００８２】
と表わせる場合、
【００８３】
【数３８】
2006003287.tif 000038
【００８４】
となる。ただし、(41)式におけるａ，ｂ，ｃは、
【００８５】
【数３９】
2006003287.tif 000039
【００８６】
【数４０】
2006003287.tif 000040
【００８７】
【数４１】
2006003287.tif 000041
【００８８】
である。ただし、これらの式で、ｌ，ｍが１，２，３の値をとる時、それぞれｘ，ｙ，ｚの添え字を表わしている。以後の場合もこの規則に準じるとする。このとき、ａ’，ｂ’，ｃ’を、
【００８９】
【数４２】
2006003287.tif 000042
【００９０】
【数４３】
2006003287.tif 000043
【００９１】
【数４４】
2006003287.tif 000044
【００９２】
により定義する。
【００９３】
【数４５】
2006003287.tif 000045
【００９４】
より、
(41)、(45)、(46)、(47)を代入して、
【００９５】
【数４６】
2006003287.tif 000046
【００９６】
とすると、
【００９７】
【数４７】
2006003287.tif 000047
【００９８】
となる。ベクトルｐ，ｑをθ回転させることにより電気ベクトルの方向になる場合、そのベクトルの成分は２つの互いに垂直な方向が存在し、
【００９９】
【数４８】
2006003287.tif 000048
【０１００】
または、
【０１０１】
【数４９】
2006003287.tif 000049
【０１０２】
となる。
(51)式の場合、
【０１０３】
【数５０】
2006003287.tif 000050
【０１０４】
となる。ただし、
【０１０５】
【数５１】
2006003287.tif 000051
【０１０６】
である。
(52)式の場合、
【０１０７】
【数５２】
2006003287.tif 000052
【０１０８】
となる。ただし、
【０１０９】
【数５３】
2006003287.tif 000053
【０１１０】
である。なお、(54)式と(56)式において、Ｈだけ比が異なっているが、これは、(39)式と(40)式において、後で計算しやすいように、ベクトルｐ，ｑの微分形式に平方根を含まないように、√Ｈだけ比が異なって定義したことに起因している。以後の過程においても、２つの異なる電気ベクトルの方向に対して、微分形式の比がＨまたは、その2乗だけ異なる場合がある。
【０１１１】
光線楕円体が、(40)式で表わせる場合、光線屈折率は
【０１１２】
【数５４】
2006003287.tif 000054
【０１１３】
となるので、その微分は、
【０１１４】
【数５５】
2006003287.tif 000055
【０１１５】
等となる。(57)式，(58)式を屈折率異方性の場合のgradient  index  光線方程式
【０１１６】
【数５６】
2006003287.tif 000056
【０１１７】
【数５７】
2006003287.tif 000057
【０１１８】
に代入して整理し、さらに、(53)式、または(55)式を代入して整理することにより、
【０１１９】
【数５８】
2006003287.tif 000058
【０１２０】
【数５９】
2006003287.tif 000059
【０１２１】
が得られる。
【０１２２】
(61)式，(62)式を冪級数展開により求める場合、
【０１２３】
【数６０】
2006003287.tif 000060
【０１２４】
とする。最初、係数Ａ0，Ｂ0，Ａ1，Ｂ1は入射した光線の位置および方向により与えられる。ｉ≧２の場合の係数Ａi，Ｂiは以下のようにして求められる。まず、
【０１２５】
【数６１】
2006003287.tif 000061
【０１２６】
と冪級数展開したとき、係数Ｈkは
【０１２７】
【数６２】
2006003287.tif 000062
【０１２８】
により求められる。(36)〜(39)式より、光線方向に直交する初期ベクトルｐ，ｑの成分に関する冪級数展開の係数は、
【０１２９】
【数６３】
2006003287.tif 000063
【０１３０】
【数６４】
2006003287.tif 000064
【０１３１】
となる。次に(42)〜(47)式より、初期ベクトルの回転に関するａ’，ｂ’，ｃ’の冪級数展開の係数は、
【０１３２】
【数６５】
2006003287.tif 000065
【０１３３】
【数６６】
2006003287.tif 000066
【０１３４】
【数６７】
2006003287.tif 000067
【０１３５】
とすることにより、γlmの冪級数展開の係数を、γlm,kとすると、
【０１３６】
【数６８】
2006003287.tif 000068
【０１３７】
【数６９】
2006003287.tif 000069
【０１３８】
【数７０】
2006003287.tif 000070
【０１３９】
となる。(49)式より、
【０１４０】
【数７１】
2006003287.tif 000071
【０１４１】
の冪級数展開の係数は、
【０１４２】
【数７２】
2006003287.tif 000072
【０１４３】
となる。次に(54)式、または、(56)式の冪級数展開の係数を求める必要があるが、微分方程式に平方根が含まれている。この場合、冪級数展開で求めると、計算精度が劣化する恐れがある。そこで、
【０１４４】
【数７３】
2006003287.tif 000073
【０１４５】
【数７４】
2006003287.tif 000074
【０１４６】
とし、
【０１４７】
【数７５】
2006003287.tif 000075
【０１４８】
とすると、(77)式の両辺を2乗して、ｚの同じ次数の項を等しいとおくことにより、
【０１４９】
【数７６】
2006003287.tif 000076
【０１５０】
というように、低次の係数から順に高次の係数を求める。従って、(54)式、または、(56)式における電気ベクトルの成分に関する冪級数展開の係数は、
【０１５１】
【数７７】
2006003287.tif 000077
【０１５２】
または、
【０１５３】
【数７８】
2006003287.tif 000078
【０１５４】
【数７９】
2006003287.tif 000079
【０１５５】
となる。従って、ベクトルe'の成分どうしの積を冪級数展開した場合の係数は、
【０１５６】
【数８０】
2006003287.tif 000080
【０１５７】
となり、γlmとそのｘ，ｙ，ｚで微分した場合の冪級数展開の係数をそれぞれ、
γlm,k，γlm,x,k，γlm,y,k，γlm,z,kとすると、光線屈折率およびその微分に関する冪級数展開の係数は、
【０１５８】
【数８１】
2006003287.tif 000081
【０１５９】
となり、gradient index 光線方程式を解く前段階の冪級数展開の係数は、
【０１６０】
【数８２】
2006003287.tif 000082
【０１６１】
となり、(61)式,(62)式より, (63)式の高次の係数Ａk+2，Ｂk+2は、
【０１６２】
【数８３】
2006003287.tif 000083
【０１６３】
【数８４】
2006003287.tif 000084
【０１６４】
となる。
【０１６５】
次に発明の実施の形態を、図１により説明する。図１において、冪級数展開の係数Ａi，Ｂi（ｉ＝０，．．．，Ｋ＋１）から、係数Ｈkを(65)式により求め、また、「光線方向に直交する初期ベクトルの成分に関する冪級数展開の係数」ｐ’l,k，ｑ’l,kを(66)式,(67)式により求め、次に、「初期ベクトルの回転に関する冪級数展開の係数」ａ’k，ｂ’k，ｃ’kおよび
【０１６６】
【数８５】
2006003287.tif 000085
【０１６７】
を(68)式〜(74)式より求める。それにより、「電気ベクトルの成分に関する冪級数展開の係数」ｅ’l,kを(75)式〜(81)式により求め、「電気ベクトルの成分の積に関する冪級数展開の係数」Ｅlm,kを(82)式により求める。さらに、「光線屈折率およびその微分に関する冪級数展開の係数」Ｆk，Ｆx,k，Ｆy,k，Ｆz,kを (83)式により求め、それにより、「gradient index 光線方程式を解く前段階の冪級数展開の係数」Ｔx,k，Ｔy,kを(84)式により求める。それと、係数Ａi，Ｂi（ｉ＝０，．．．ｋ＋１），Ｈｉ（ｉ＝０，．．．ｋ）より、高次の係数Ａk+2，Ｂk+2を(85)式,(86)式により求める。この計算が終わると、図１の最初の段階に戻り、より高次の冪級数展開の係数を求める。
【実施例】
【０１６８】
次に、本発明による効果を実施例により説明する。先ず、屈折率異方性の分布が子午面に関して対称で、かつ光線が子午面にのみ存在する場合において、gradient  index光線方程式を冪級数展開により求め、数値解析の結果と比較した。実施例１として、間隔が20ｍｍで、平面の材質の場合で、屈折率がｎx＝１．６，ｎy＝ｎz＝１．５で、楕円体の主軸と座標軸が一致し、
【０１６９】
【数８６】
2006003287.tif 000086
【０１７０】
光線出発点が（ｘ，ｙ，ｚ）=(0,0,0)で、最初の光線方向が（ｔx，ｔy，ｔz）=(0.3,0.3,0.953939)  の場合、
(31)式のように冪級数展開によりｘの到達位置を求め、数値解析による計算結果と比較した。その結果を表1に示す。
【０１７１】
【表１】
2006003287.tif 000087
【０１７２】
表１において、分割数∞での到達位置は、分割数20,  200,  2000での到達位置から、到達位置を分割数Nの関数ｆ（Ｎ）として、推定式
【０１７３】
【数８７】
2006003287.tif 000088
【０１７４】
より、定数を求め、
【０１７５】
【数８８】
2006003287.tif 000089
【０１７６】
より、分割数∞での到達位置を推定した値である。数値解析による分割数∞での推定値と、冪級数展開の値は、10の-10乗まで一致していることが分かる。
次に、実施例２として、間隔が5ｍｍで、平面の材質の場合で、光線楕円体の係数γlmが
【０１７７】
【数８９】
2006003287.tif 000090
【０１７８】
のように冪級数展開できる場合、
【０１７９】
【数９０】
2006003287.tif 000091
【０１８０】
でそれ以外の係数は全てゼロの場合で、光線出発点が（ｘ，ｙ，ｚ）=(10,0,0)で、最初の光線方向が（ｔx，ｔy，ｔz）=(-0.3,0.3,0.905539)  の場合、(86)式，(87)式により、冪級数展開によりｘ、ｙの到達位置を求め、数値解析による計算結果と比較した。その結果を表２およびに表３に示す。表２は電気ベクトルの方向が(51)式で表わされる場合の光線追跡であり、表３は電気ベクトルの方向が(52)式で表わされる場合の光線追跡である。
【０１８１】
【表２】
2006003287.tif 000092
【０１８２】
【表３】
2006003287.tif 000093
【０１８３】
表２および表３において、分割数∞での到達位置は、分割数20,  200,  2000での到達位置から、到達位置を分割数Nの関数ｆ（Ｎ）として、(87)式の推定式により、定数Ｃ0，Ｃ1，αを求め、(88)式より、分割数∞での到達位置を推定した値である。数値解析による分割数∞での推定値と、冪級数展開の値は、10の-8乗まで一致していることが分かる。
【０１８４】
上述の各実施例によれば、屈折率異方性がある媒質内でも、gradient  index光線方程式を冪級数展開することにより、数値解析よりも高速で高精度に計算することが可能である。
【図面の簡単な説明】
【０１８５】
【図１】本発明の計算の順序を示す図である。
【図１】
2006003287.tif 000096