(57)【要約】【課題】屈折率異方性がある媒質内でもｇｒａｄｉｅｔ  ｉｎｄｅｘ光線追跡が可能な光線方程式を導出し、その光線方程式により、高精度に光線追跡を行う方法を提供するする。【解決手段】光線位置から内部応力を求め、内部応力から光線楕円体の係数を求める。次に求めた光線楕円体の係数と、与えられた光線方向から、光線楕円体の幾何学的な性質から単位電気ベクトルｅを求めるか、残差ベクトルを最小にするように多元のＮｅｗｔｏｎ法で単位電気ベクトルｅを求める。さらにその距離による微分を求める。次に、単位電気ベクトルｅから、光線屈折率ｎ’を求める。さらにその距離による微分を求める。求めた光線屈折率と、その距離による微分、および光線方向の単位ベクトルを、光線方程式に代入して、次の微小距離進行後の光線方向と光線位置を求める。光線方程式は、その微分方程式を差分化するか、Ｒｕｎｇｅ−Ｋｕｔｔａ法により数値計算する。【選択図】  図１
2004361134.tif 000002
【特許請求の範囲】
【請求項１】
屈折率が一様でなく屈折率異方性のある媒質内での光線追跡方であって、光線追跡に用いる光線方程式
【数１】
2004361134.tif 000096
において、ｎ’を光線屈折率とすることを特徴とする光線追跡方法。
【請求項２】
請求項１に記載の光線追跡方法であって、前記光線屈折率を、「光線方向のベクトルｔと直交する平面と、光線楕円体との交線である楕円の２つの軸の方向が電気ベクトルＥの方向であり、軸の長さが、光線屈折率の逆数に比例する」という関係を利用して、幾何学的に計算することを特徴とする光線追跡方法。
【請求項３】
請求項１に記載の光線追跡方法であって、光線屈折率を、残差ベクトル：
【数２】
2004361134.tif 000097
（ただし、ｎ’は光線屈折率、ｅは単位電気ベクトル、ｄは単位電気変位ベクトル、ｔは光線方向の単位ベクトル）を最小にすることにより、単位電気ベクトルｅを求め、その成分により、式：
【数３】
2004361134.tif 000098
により求めることを特徴とする光線追跡方法。
但し、γｉｊは光線楕円体の係数であり、ｉ、ｊについては、１をｘ−ｙ−ｚ直交座標系のｘ、２をｙ、３をｚとするものとする。
【請求項４】
請求項３に記載の光線追跡方法であって、前記残差ベクトルｆ（ｅ）を最小にする方法が、単位電気ベクトルｅの成分を未知数とするニュートン法であることを特徴とする光線追跡方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、偏光を利用する光学系または、半導体素子製造用の投影光学系や、天体観測用の光学系等において、特に光線通過部材の、加熱や外力による内部応力分布の変動による光線の偏光分離が問題となるような光学系中における光線の追跡方法や、電離層のように異方性不均質な磁気プラズマ中を伝搬する光線の追跡方法に関するものである。
【０００２】
【従来の技術】
従来、複屈折がある硝材を光線が透過する場合の光線追跡は、例えば、１軸結晶を光線が入射する場合の屈折光線の屈折率と波面放線方向を求める方法については、Ｍ．Ｃ．Ｓｉｍｏｎによって、“Ｒａｙ  ｔｒａｃｉｎｇ  ｆｏｒｍｕｌａｓ  ｆｏｒ  ｍｏｎｏａｘｉａｌ  ｏｐｔｉｃａｌ  ｃｏｍｐｏｎｅｎｔｓ“という題名の文献で紹介されている（Ａｐｐｌ．Ｏｐｔ．（１９８３）ｐｐ３５４−３６０）。
【０００３】
また、一般的に２軸結晶に光線が入射する場合の屈折光線の屈折率と波面法線ベクトルを求める方法の原理については、ボルン・ウオルフ著：「光学の原理ＩＩＩ」（邦訳：東海大学出版会）において、Ｆｒｅｓｎｅｌの公式（同文献、Ｐ９８０参照）を利用して求める方法が示されている。（同文献、Ｐ９９９〜Ｐ１０００参照）
しかし、ガラスの温度上昇分布に差がある場合や、ガラスの一部に外力を加えた場合、ガラス内で歪みが発生し、内部応力が発生するが、発生する内部応力は、通常方向異方性があるので、屈折率異方性があり、光線がガラスを透過すると、偏光になる。一方、内部応力は、ガラス内の位置により異なっていて、それに従って、屈折率も異なるので、ｇｒａｄｉｅｔ  ｉｎｄｅｘ光線追跡をする必要がある。
【０００４】
内部応力のある媒質内でのｇｒａｄｉｅｎｔ  ｉｎｄｅｘ光線追跡では、計算誤差が蓄積しないように、屈折率の計算は高精度に行う必要があるが、前記のＭ．Ｃ．Ｓｉｍｏｎの文献等で述べられているような屈折率を計算する方法では、ｇｒａｄｉｅｎｔ  ｉｎｄｅｘ光線追跡に絶えられるような精度が保証できないことがわかった。そこで、発明者は、波面法線楕円体によって、屈折率を計算する方法を開発した（橋本純夫：「内部応力分布があるガラス内の光線追跡」光学２６（１９９７）  ６７７）。
【０００５】
【非特許文献１】Ａｐｐｌ．Ｏｐｔ．（１９８３）ｐｐ３５４−３６０
【非特許文献２】「光学の原理ＩＩＩ」
【非特許文献３】光学２６（１９９７）  ６７７
【０００６】
【発明が解決しようとする課題】
しかし、前記「内部応力分布があるガラス内の光線追跡」に示される方法では、内部応力のある媒質内での屈折率は正確に計算できるが、媒質内では、波面法線方向と光線方向が異なるために、光線の微小進行毎に、光線の進行方向を電場ベクトルと磁場ベクトルの外積すなわち、Ｐｏｙｎｔｉｎｇ  ベクトルを求めることによって光線方向を修正する必要があった。このような計算方法であると、微小進行毎の修正によって、計算誤差が蓄積する恐れがあった。
【０００７】
本発明はこのような事情に鑑みてなされたもので、屈折率異方性がある媒質内でもｇｒａｄｉｅｔ  ｉｎｄｅｘ光線追跡が可能な光線方程式（光線ベクトルの方向を直接求める微分方程式）を導出し、その光線方程式により、高精度に光線追跡を行う方法を提供するすることを課題とする。
【０００８】
【課題を解決するための手段】
前記課題を解決するための第１の手段は、屈折率が一様でなく屈折率異方性のある媒質内での光線追跡方であって、光線追跡に用いる光線方程式
【０００９】
【数４】
2004361134.tif 000003
【００１０】
において、ｎ’を光線屈折率とすることを特徴とする光線追跡方法（請求項１）である。
【００１１】
なお、本明細書においては、ベクトルを示すのに式中では通常のベクトル記号を使用するが、文章中でイメージデータを用いると煩雑になるので、例えば上記式における
【００１２】
【数５】
2004361134.tif 000004
【００１３】
を、ベクトルｒ、ｒ（ベクトル）のようにして表す場合があり、又、前後の関係からベクトルであることが明らかである場合は、これらの記載を省略して単にｒと表す場合がある。
【００１４】
まず、屈折率異方性がある媒質内でもｇｒａｄｉｅｔ  ｉｎｄｅｘ光線追跡が可能な前記光線方程式を導いた方法について述べる。
【００１５】
光線の磁場、および電場、電気変位ベクトルはそれぞれ次のように表わされる。
【００１６】
【数６】
2004361134.tif 000005
【００１７】
ここで、ｋ０は真空中の波数であり、φはｅｉｋｏｎａｌである。
【００１８】
【数７】
2004361134.tif 000006
【００１９】
であり、ｒは光線の位置ベクトル、Ｓは波面法線方向の単位ベクトル、ｎは屈折率である。屈折率ｎは光速ｃと位相速度ｖｐとにより次のように表される。
【００２０】
【数８】
2004361134.tif 000007
【００２１】
（１）〜（３）をＭａｘｗｅｌｌ方程式に代入することにより、
【００２２】
【数９】
2004361134.tif 000008
【００２３】
となる。そして、（５）に（６）を代入することにより、
【００２４】
【数１０】
2004361134.tif 000009
【００２５】
であるから、
【００２６】
【数１１】
2004361134.tif 000010
【００２７】
となる。光線方向はポインテイングベクトル
【００２８】
【数１２】
2004361134.tif 000011
【００２９】
の方向である。ここで、
【００３０】
【数１３】
2004361134.tif 000012
【００３１】
はベクトルＨの複素共役である。しかし、複屈折のある媒質内では、２本の直線偏光に分かれ、直線偏光ではベクトルＨが実数なので、光線方向は
【００３２】
【数１４】
2004361134.tif 000013
【００３３】
とすることができる。従って、（６）より、
【００３４】
【数１５】
2004361134.tif 000014
【００３５】
となる。等方性媒質の場合、
【００３６】
【数１６】
2004361134.tif 000015
【００３７】
であるので、（７）より、
【００３８】
【数１７】
2004361134.tif 000016
【００３９】
となることにより、前式は簡単になる。しかし、異方性媒質の場合、ベクトルＤとベクトルＥは同じ方向のベクトルでないので、
【００４０】
【数１８】
2004361134.tif 000017
【００４１】
となり、（１２）式の第２項を無視することはできない。従って、この式からｅｉｋｏｎａｌ方程式
【００４２】
【数１９】
2004361134.tif 000018
【００４３】
を導出することはできない。しかし、幾何光学的な波面の考察から、
【００４４】
【数２０】
2004361134.tif 000019
【００４５】
が成立するので、この式の両辺を２乗すると、ｅｉｋｏｎａｌ方程式（１６）は成立する。そして、（１２）式より、
【数２１】
2004361134.tif 000020
従って、光線方向の単位ベクトルは、
【００４６】
【数２２】
2004361134.tif 000021
【００４７】
（１９）式で
【００４８】
【数２３】
2004361134.tif 000022
【００４９】
とおくと、
【００５０】
【数２４】
2004361134.tif 000023
【００５１】
となる。ｒを光線の位置ベクトル、ｓを光線に沿って測った長さとすると、
【００５２】
【数２５】
2004361134.tif 000024
【００５３】
が成立する。例えばベクトルｔのｘ成分において、
【００５４】
【数２６】
2004361134.tif 000025
【００５５】
であるので、
【００５６】
【数２７】
2004361134.tif 000026
【００５７】
従って、（２３）式に（２１）式を代入して（２５）式を用いることにより、
【００５８】
【数２８】
2004361134.tif 000027
【００５９】
となる。
【００６０】
【数２９】
2004361134.tif 000028
【００６１】
に（２３）式および（２１）式を代入して、
【００６２】
【数３０】
2004361134.tif 000029
【００６３】
つまり、異方性媒質の場合、周知の波面法線方向の屈折率の代わりに上記ｎ’を用いることにより、Ｇｒａｄｉｅｎｔ  Ｉｎｄｅｘ  光線追跡が可能であるということがわかる。
【００６４】
次に、このｎ’が光線屈折率に等しいことを照明する。（２０）式において、
【００６５】
【数３１】
2004361134.tif 000030
【００６６】
は、波面の幾何学的な考察から、（１６）式により与えられる。そして、（２０）式における残りの未知量は
【００６７】
【数３２】
2004361134.tif 000031
【００６８】
であるが、（５）、（６）より、
【００６９】
【数３３】
2004361134.tif 000032
【００７０】
であるから、
【００７１】
【数３４】
2004361134.tif 000033
【００７２】
となる。両辺を２乗して整理することにより、
【００７３】
【数３５】
2004361134.tif 000034
【００７４】
（１６）式を代入して整理することにより、
【００７５】
【数３６】
2004361134.tif 000035
従って、（２２）式より、
【００７６】
【数３７】
2004361134.tif 000036
【００７７】
「光学の原理ＩＩＩ」Ｐ９７８より、
【００７８】
【数３８】
2004361134.tif 000037
【００７９】
だから、
【００８０】
【数３９】
2004361134.tif 000038
【００８１】
となる。これは、「光学の原理ＩＩＩ」Ｐ９７８によると、光線屈折率（エネルギー屈折率）になる。なおベクトルＤとベクトルＥのなす角をｃｏｓαとすると、
【００８２】
【数４０】
2004361134.tif 000039
【００８３】
となる。
【００８４】
前記課題を解決するための第２の手段は、前記第１の手段であって、前記光線屈折率を、「光線方向のベクトルｔと直交する平面と、光線楕円体との交線である楕円の２つの軸の方向が電気ベクトルＥの方向であり、軸の長さが、光線屈折率の逆数に比例する」という関係を利用して、幾何学的に計算することを特徴とするものである。
【００８５】
前記課題を解決するための第３の手段は、前記第１の手段であって、光線屈折率を、残差ベクトル：
【００８６】
【数４１】
2004361134.tif 000040
【００８７】
（ただし、ｎ’は光線屈折率、ｅは単位電気ベクトル、ｄは単位電気変位ベクトル、ｔは光線方向の単位ベクトル）を最小にすることにより、単位電気ベクトルｅを求め、その成分により、式：
【００８８】
【数４２】
2004361134.tif 000041
【００８９】
により求めることを特徴とするもの（請求項３）である。但し、γｉｊは光線楕円体の係数であり、ｉ、ｊについては、１をｘ−ｙ−ｚ直交座標系のｘ、２をｙ、３をｚとするものとする。
【００９０】
前記課題を解決するための第４の手段は、前記第３の手段であって、前記残差ベクトルｆ（ｅ）を最小にする方法が、単位電気ベクトルｅの成分を未知数とするニュートン法であることを特徴とするもの（請求項４）である。
【００９１】
【発明の実施の形態】
以下、前記第２の手段から第４の手段の実施の形態について説明する。
【００９２】
まず、光線屈折率を幾何学的に求める方法の例を説明する。光線方向の屈折率は、光線楕円体（原点からの距離が光線屈折率の逆数に比例する楕円体）により計算することができる。すなわち、「波面法線楕円体との双対性により、光線方向のベクトルｔと直交する平面と、光線楕円体との交線である楕円の２つの軸の方向が電気ベクトルＥの方向であり、軸の長さが、光線屈折率の逆数に比例する」
（「光学の原理ＩＩＩ」，Ｐ９８８）ことを利用する。そのために応力分布から光線楕円体を求める。すなわち、光線楕円体を、
【００９３】
【数４３】
2004361134.tif 000042
【００９４】
と置く。但し、γｉｊについては、ｉとｊについては、１のときｘ、２のときｙ、３のときｚとし、
ｘ１＝ｘ、ｘ２＝ｙ、ｘ３＝ｚ    …（３８）
とする
光線方向の単位ベクトルｔ＝（ｔｘ，ｔｙ，ｔｚ）に垂直で、互いに垂直な任意の２つの単位ベクトルをベクトルｕ＝（ｕｘ，ｕｙ，ｕｚ）、ベクトルｖ＝（ｖｘ，ｖｙ，ｖｚ）とすると、（ｘ１，ｘ２，ｘ３）を、ベクトルｔ、ベクトルｕ、ベクトルｖを直交座標軸とする座標上の点（ｔ，ｕ，ｖ）に変換する場合、１のときｘ、２のときｙ、３のときｚとし、（３８）式のように表わすと、
【００９５】
【数４４】
2004361134.tif 000043
【００９６】
となる。従って、（ｘ１，ｘ２，ｘ３）が光線方向の単位ベクトルｔに垂直な平面ｔ＝０上にある場合、
【００９７】
【数４５】
2004361134.tif 000044
【００９８】
となる。（４０）式を（３７）式に代入して、
【００９９】
【数４６】
2004361134.tif 000045
【０１００】
となる。これは、光線方向のベクトルと直交する平面と、光線楕円体の光線である楕円を表わしているので、この楕円を
【０１０１】
【数４７】
2004361134.tif 000046
【０１０２】
と置くと、
【０１０３】
【数４８】
2004361134.tif 000047
【０１０４】
となる。ベクトルｕ、ベクトルｖの座標軸を角度α回転させて、ベクトルｐ、ベクトルｑになったとすると、ベクトルｐ、ベクトルｑの座標軸による座標値Ｐ，Ｑは
【０１０５】
【数４９】
2004361134.tif 000048
【０１０６】
となるので、
【０１０７】
【数５０】
2004361134.tif 000049
【０１０８】
となる。（４７）式を（４２）式に代入して整理すると、
【０１０９】
【数５１】
2004361134.tif 000050
【０１１０】
となる。ベクトルｐ、ベクトルｑが楕円の主軸となるためには、（４８）式にあるＰＱの係数がゼロになる必要があるので、角度αの満たすべき条件は、
【０１１１】
【数５２】
2004361134.tif 000051
【０１１２】
となる。（４９）式に（４３）式から（４５）式を代入することにより、角度αを求める。（４６）式より、
【０１１３】
【数５３】
2004361134.tif 000052
【０１１４】
となり、（５０）式より、電気ベクトルの単位ベクトルｅがベクトルｐ、ベクトルｑのうちどちらかになる。
【０１１５】
ベクトルｅがわかると、原点からベクトルｅと光線楕円体の交点との距離は光線屈折率の逆数になるので、（３７）式  にｘｉ＝ｅｉ／ｎ’、ｘｊ＝ｅｉｊ／ｎ’を代入して、
【０１１６】
【数５４】
2004361134.tif 000053
【０１１７】
となる。
また、光線屈折率の距離による微分は、
【０１１８】
【数５５】
2004361134.tif 000054
【０１１９】
となる。
【０１２０】
このようにして、（５２）式から光線屈折率ｎ’が求まる。なお、（５２）式におけるγｉｊの求め方については後述する。
【０１２１】
以下、残差ベクトルを最小にすることにより、単位電気ベクトルｅを求め、その成分から光線屈折率を求める方法の例について説明する。この方法においては、電気ベクトルＥの他に、電気変位ベクトルＤも必要とするので、それに先立ち、電気ベクトルＥから電気変位ベクトルＤを求める方法を示す。
【０１２２】
（３７）式の光線楕円体が座標変換
【０１２３】
【数５６】
2004361134.tif 000055
【０１２４】
によって楕円体の主軸と座標軸とが一致するように直交変換した場合、すなわち、光線楕円体を
【０１２５】
【数５７】
2004361134.tif 000056
【０１２６】
という形に直交変換した場合、（５５）式に（５４）式を代入して、
【０１２７】
【数５８】
2004361134.tif 000057
【０１２８】
となる。（３７）式と比較して、
【０１２９】
【数５９】
2004361134.tif 000058
【０１３０】
（５５）式のように、光線楕円体の主軸と座標軸が一致している場合、
【０１３１】
【数６０】
2004361134.tif 000059
【０１３２】
であり、また、
【０１３３】
【数６１】
2004361134.tif 000060
【０１３４】
であるので、
【０１３５】
【数６２】
2004361134.tif 000061
【０１３６】
となる。Ｃｉｊが単位行列であることから、（５４）式の逆変換において、
【０１３７】
【数６３】
2004361134.tif 000062
【０１３８】
であり、（６１）式に（６０）式を代入し、また、（５４）式の変換により、
【０１３９】
【数６４】
2004361134.tif 000063
【０１４０】
であるので、（６１）式に（６２）式を代入し、（５７）式により、ベクトルＥからベクトルＤを求める関係式は、
【０１４１】
【数６５】
2004361134.tif 000064
【０１４２】
となる。このようにして、電気ベクトルＥから電気変位ベクトルＤが求まる。なお、γｉｊの求め方については後述する。
【０１４３】
光線楕円体による計算方法の精度を確認する為に、残差ベクトルを導入する。（１０）式に、（２１）式、（２２）式の関係を代入すると、
【０１４４】
【数６６】
2004361134.tif 000065
となり、従って、
【０１４５】
【数６７】
2004361134.tif 000066
【０１４６】
となる。また、（７）式に、（２１）式を代入して、
【０１４７】
【数６８】
2004361134.tif 000067
【０１４８】
（６５）式に（６７）式を代入して、
【０１４９】
【数６９】
2004361134.tif 000068
電気ベクトルＥを単位電気ベクトルとし、（６３）式の代わりに
【０１５０】
【数７０】
2004361134.tif 000069
【０１５１】
とすると、
【０１５２】
【数７１】
2004361134.tif 000070
【０１５３】
となるので、（３５）式より、
【０１５４】
【数７２】
2004361134.tif 000071
【０１５５】
また、ｅは単位ベクトルなので、
【０１５６】
【数７３】
2004361134.tif 000072
これらを（６９）式に代入すると、
【０１５７】
【数７４】
2004361134.tif 000073
【０１５８】
となる。そこで、
【０１５９】
【数７５】
2004361134.tif 000074
【０１６０】
を残差ベクトルとして、この計算方法の精度を見積もる。そして、この残差ベクトルを最小にするような電気単位ベクトルｅを求めれば、（５２）式を使用することにより、光線屈折率を求めることができる。以下、この方法について説明する。
【０１６１】
残差ベクトルｆ（ｅ）を最小にするような電気単位ベクトルｅは、波面法線方向の屈折率を多元のＮｅｗｔｏｎ法で求めたのと同様な方法で、残差ベクトルを用いて求めることができる。すなわち、ある電気単位ベクトルｅの近似値に対して、
【０１６２】
【数７６】
2004361134.tif 000075
【０１６３】
のような未知数が３個の連立方程式を満たす解Δｅｊ（ｊ＝１，３）を計算し、ベクトルｅのさらによい近似をｅ＋Δｅ（Δｅはベクトル）とし、この計算を繰り返して、高精度なベクトルｅの値を求める。ベクトルｅを求めた後、（５２）式により、光線屈折率ｎ’を求める。（７６）式において、ベクトルｔが与えられた場合における残差ベクトルｆ（ｅ）の成分ｆｉを、ベクトルｅの成分ｅｊで微分すると、（５２）式、（７０）式、（７５）式より、
【０１６４】
【数７７】
2004361134.tif 000076
【０１６５】
となる。ここで、δｉｊはＫｒｏｎｅｃｋｅｒのδである。
【０１６６】
光線屈折率ｎ’の距離による微分を、残差ベクトルから求めることもできる。多元のＮｅｗｔｏｎ法で、電気単位ベクトルｅとｎ’を求めて後、
【０１６７】
【数７８】
2004361134.tif 000077
【０１６８】
の両辺を距離ｘで微分して、先ず
【０１６９】
【数７９】
2004361134.tif 000078
【０１７０】
を求めた後、
【０１７１】
【数８０】
2004361134.tif 000079
を求める。（７５）式において、（５２）式、（７０）式より、
【０１７２】
【数８１】
2004361134.tif 000080
【０１７３】
（５２）式より、  また、
【０１７４】
【数８２】
2004361134.tif 000081
【０１７５】
より、
【０１７６】
【数８３】
2004361134.tif 000082
【０１７７】
と置くと、
【０１７８】
【数８４】
2004361134.tif 000083
【０１７９】
となり、∂ｅｊ／∂ｘを未知数とする連立方程式になるので、（８２）式を解けば、∂ｅｊ／∂ｘが求められ、
（５３）式により、∂ｎ’／∂ｘが求められる。
【０１８０】
以上のようにして、幾何学的方法又は逐次近似法により光線屈折率ｎとその微分を求めることができるが、光線楕円体の係数γｉｊが求まっていることを前提としてきた。以下、光線楕円体の係数γｉｊを求める方法について説明する。
【０１８１】
光線楕円体の主軸が全て、ｘ’ｙ’ｚ’軸と一致する場合、それぞれの方向の屈折率は、
【０１８２】
【数８５】
2004361134.tif 000084
【０１８３】
となる。ｘｙ軸をθ回転させて、ｘ’ｙ’軸になったとすると、ｘ’ｙ’座標系での主応力σｘ’ｘ’、σｙ’ｙ’は、ｘｙ座標系での主応力σｘｘ、σｙｙおよびせん断応力σｘｙにより、
【０１８４】
【数８６】
2004361134.tif 000085
【０１８５】
と表わされる。この場合、光線楕円体の係数は、原点と、ｘ’軸と光線楕円体の交点の距離をｒ’とすると、ｘ，ｙ，ｚをｘ’，ｙ’，ｚ’に置き換えた式に、（ｘ’，ｙ’，ｚ’）＝（ｒ’，０，０）を代入して、
【０１８６】
【数８７】
2004361134.tif 000086
【０１８７】
となる。一方、ｘ，ｙ，ｚ’座標軸において、
（ｘ，ｙ，ｚ）＝（ｒｃｏｓθ，ｒｓｉｎθ，０）を光線楕円体に代入すると、
【０１８８】
【数８８】
2004361134.tif 000087
【０１８９】
となる。ｒとｒ’は同じ距離なので、（８８）式と（９０）式が任意の角度θに対して成り立つように三角関数の係数を比較すると、
【０１９０】
【数８９】
2004361134.tif 000088
【０１９１】
となる。
【０１９２】
又、光線楕円体の係数は、以下のようにして精密に求めることができる。ｘｙｚ軸を回転させて、ｘ’ｙ’ｚ’軸になる場合に光線楕円体の主軸と座標軸が一致した場合、（５５）式のようになる。光線楕円体の主軸と座標軸が一致した場合、この座標系では、せん断応力は存在しない。従って、せん断応力が存在しないように座標系を回転させると、光線楕円体の主軸と座標軸が一致することになる。先ず、ｘｙｚ座標系を、ｚ軸を中心にθ回転させて、ｘ’ｙ’ｚ’座標系になったとすると、ｘ’ｙ’ｚ’座標系でのせん断応力σｘ’ｙ’はｘｙ座標系での主応力σｘｘ、σｙｙおよびせん断応力σｘｙにより、
【０１９３】
【数９０】
2004361134.tif 000089
【０１９４】
と表わされるが、σｘ’ｙ’となるようにするには、
【０１９５】
【数９１】
2004361134.tif 000090
【０１９６】
となり、
【０１９７】
【数９２】
2004361134.tif 000091
【０１９８】
回転させることが必要である。同様にして、他の方向の軸を中心に回転させて、その方向に垂直な面でのせん断応力をゼロにして、これらの操作を２，３度繰り返すことにより、ほぼ完全に光線楕円体の主軸と座標軸が一致することになる。この場合の、ｘｙｚ軸からｘ’ｙ’ｚ’軸への直交変換が（５４）式のような場合、光線楕円体は（５５）式のようになり、（５５）式に（５４）式を代入して、（５６）式のようになる。
【０１９９】
従って、ｘｊｘｋの項の係数は、
【０２００】
【数９３】
2004361134.tif 000092
【０２０１】
となる。ここで、
【０２０２】
【数９４】
2004361134.tif 000093
【０２０３】
である。従って、光線楕円体の係数の距離による微分は、
【０２０４】
【数９５】
2004361134.tif 000094
【０２０５】
となる。
【０２０６】
以上のように述べた計算方法をアルゴリズム図で示すと、図１のようになる。図１で、矢印のある実線は、同じ微小距離のステップ内でのアルゴリズムの流れを示す。そして、矢印のある破線は、次の微小距離のステップに移行するときのアルゴリズムの流れを示す。
【０２０７】
先ず、光線位置から内部応力を求め、内部応力から光線楕円体の係数を（９１）式〜（９３）式、または、（９６）式〜（９８）式により求める。次に求めた光線楕円体の係数と、与えられた光線方向から、（４３）式〜（４５）式、（４９）式、（５０）式のように光線楕円体の幾何学的な性質から単位電気ベクトルｅを求めるか、（７５）式の残差ベクトルを最小にするように（７６）式のような多元のＮｅｗｔｏｎ法で単位電気ベクトルｅを求める。さらにその距離による微分を、（４３）式〜（４５）式、（４９）式、（５０）式を距離で微分するか、（８１）式、（８２）式により求める。
【０２０８】
次に、（５２）式により、単位電気ベクトルｅから、光線屈折率ｎ’を求める。さらにその距離による微分を、（５３）式により求める。求めた光線屈折率と、その距離による微分、および光線方向の単位ベクトルを、（２８）式の光線方程式に代入して、次の微小距離進行後の光線方向と光線位置を求める。光線方程式は、その微分方程式を差分化するか、Ｒｕｎｇｅ−Ｋｕｔｔａ法により数値計算する。
【０２０９】
なお以前開発した波面法線を追跡する方法（被特許文献３）に比べて、微小距離進行毎に波面法線方向から光線方向を計算する過程が無くなり、やや簡単な構成となっている。
【０２１０】
なお、光線屈折率を求める方法を２つ示したが、（７５）式の残差ベクトルによって計算の精度を見積もると、光線楕円体の幾何学的な性質から単位電気ベクトルｅを求める方法も、（７６）式のような多元のＮｅｗｔｏｎ法で単位電気ベクトルｅを求める方法も、倍精度の限界の程度まで残差を小さくすることができ、精度的にはどちらも大差無いことがわかった。また、それぞれの方法により、光線屈折率の微分を求めることができるが、光線楕円体による方法であると、（４９）式の距離による微分を求める必要があり、この数式計算はかなり煩雑である。むしろ、多元のＮｅｗｔｏｎ法から、（８０）式、（８１）式を用いてから（５３）式を用いる方法の方が、簡単である。
【０２１１】
以上述べた計算方法は、ガラス内での光線追跡であるが、ガラス外部の空気などの異なった媒質から光線が入射する場合は、光線方向に対してはスネルの法則が適用できないので、以前のように、先ず波面法線楕円体により波面法線方向を決定してから、光線方向を求める。そして、ガラスから出射する際には、（６３）式により電気ベクトルから電気変位ベクトルを求め、電気変位ベクトルと磁気ベクトルの双方に垂直な方向が波面法線方向となる。その波面法線方向から、スネルの法則によりガラス出射後の波面法線方向を求める。
【０２１２】
【実施例】
以上に述べた計算方法で、実施例として、外部応力によるガラス内の屈折率変動による収差変動を計算した。計算例として、図２のようなレンズＧ１〜Ｇ２０からなる投影光学系（１／４倍、像側のＮＡ：０．５６，最大像高１２ｍｍ）に、波長２４８ｎｍの光線を透過させた。この投影レンズのレンズデータは、表１のようになっている。表１で、第０面は物体面（レチクル）であり、第４１面は像面（ウェハ）である。また、媒質でＡＩＲとなっているのは、空気であり、ＳＩＯ２（ＳｉＯ２）となっているのは石英ガラスである。計算に用いたガラスのヤング率は、７２×１０５Ｎ／ｍｍ２，ポアソン比は０．１６４，光弾性定数は、Ｃ１＝６．３９５×１０−７ｍｍ２／Ｎ、
Ｃ２＝４．１７６×１０−６ｍｍ２／Ｎである。
【０２１３】
【表１】
2004361134.tif 000095
【０２１４】
外半径２７．７ｍｍのレンズＧ２０の子午線方向（ｙ軸とする）の両端に、１０Ｎ／ｃｍ２の外部応力を加えた場合のガラス内の応力分布は、有限要素法により、図３のようになっている。図３（Ａ）はｙ軸方向の応力分布で、図３（Ｂ）はｚ軸方向（光軸およびｙ軸の両方に垂直な方向）の応力分布である。図３の等高線の間隔は、１Ｎ／ｃｍ２である。
【０２１５】
このような外力を加えた場合において、レンズＧ２０に光束を透過させ、前記のように、光線方程式に光線屈折率を用いてｇｒａｄｉｅｎｔ  ｉｎｄｅｘ光線追跡をして（この方法を略して光線屈折率法ということにする）波面収差を計算し、外力を加えない場合の波面収差との差をとり、波面収差変動を計算した。その結果、波面収差変動は図４のようになる。電気変位ベクトルがｙ軸方向の偏光を、前記のような応力分布のガラスに透過させて、像面上でｙ軸方向の偏光成分の波面収差変動を計算したのが図４（Ａ）であり、電気変位ベクトルがｚ軸方向の偏光を、同様な応力分布のガラスに透過させて、像面上でｚ軸方向の偏光成分の波面収差変動を計算したのが図４（Ｂ）である。
【０２１６】
図４（Ａ）の差の最大値は、０．０００３８λであり、図４（Ｂ）の差の最大値は、０．００３２λである。また、図４（Ａ）の等高線の間隔は、０．００００３８λであり、図４（Ｂ）の等高線の間隔は、０．０００３２λである。但し、応力分布を計算するときの有限要素法の分割数や、応力分布のｆｉｔｔｉｎｇ法の違いにより、波面収差変動分布に若干の相違がある。
【０２１７】
なお、ｙ軸方向の偏光を同じ状態のガラスに透過させて、像面上でｚ軸方向の偏光成分を計算したものと、ｚ軸方向の偏光を同じ状態のガラスに透過させて、像面上でｙ軸方向の偏光成分を計算したものとは、図４（Ａ）および図４（Ｂ）の偏光よりも、かなり、強度が小さく無視できる程である。
【０２１８】
また、このような応力状態の場合、今回新たに発明した光線屈折率法に対して、従来のように波面法線方向の屈折率を求め、光線の微小進行毎に光線の進行方向を電場ベクトルと磁場ベクトルの外積を求めることによって光線方向を修正する方法（波面法線法と略する）とを比較した。従来の計算方法による波面収差変動は図５のようになる。電気変位ベクトルがｙ軸方向の偏光を、前記のような応力分布のガラスに透過させて、像面上でｙ軸方向の偏光成分の波面収差変動を計算したのが図５（Ａ）であり、電気変位ベクトルがｚ軸方向の偏光を、同様な応力分布のガラスに透過させて、像面上でｚ軸方向の偏光成分の波面収差変動を計算したのが図５（Ｂ）である。図５（Ａ）の差の最大値は、０．００１２λであり、図５（Ｂ）の差の最大値は、０．００３２λである。また、図５（Ａ）の等高線の間隔は、０．０００１２λであり、図５（Ｂ）の等高線の間隔は、０．０００３２λである。
【０２１９】
図４と図５を比較してみると、波面収差変動分布がかなり違っていることが分かる。さらに図３を参照してみると、今回新たに開発した光線屈折率法の分布の方が、応力の分布に近く、理想的であることが分かる。なお、投影レンズのＧ２０のように光線が集光しているような場合は、光線屈折率法と波面法線法によって、分布がかなり違っているが、ガラスに平行な光線を通すような場合は、光線屈折率法と波面法線法による計算結果に、それ程違いは無いことが分かった。
【０２２０】
【発明の効果】
以上説明したように、本発明によれば、レンズによって集光する場合のように、波面法線方向と光線方向が異なる場合においても、内部応力による収差変動を、高精度に計算することが可能となる。
【図面の簡単な説明】
【図１】本発明の実施の形態である光線追跡方法の手順を示す図である。
【図２】本発明の実施例に使用した投影レンズの光路図である。
【図３】本発明の実施例に使用した、レンズに加えた応力分布を示す図である。
【図４】本発明の実施例である光線屈折率法による波面収差変動分布計算結果を示す図である。
【図５】波面法線法による波面収差変動分布計算結果を示す図である。
【符号の説明】
Ｇ１〜Ｇ２０…レンズ
【図１】
2004361134.tif 000099
【図２】
2004361134.tif 000100
【図３】
2004361134.tif 000101
【図４】
2004361134.tif 000102
【図５】
2004361134.tif 000103