(57)【要約】【課題】高精度で、またＭａｘｗｅｌｌ方程式を満足するような形で、光線追跡を行う方法を提供する。【解決手段】屈折率異方性のある媒質内での光線追跡であって、残差ベクトル【数１】2004028843.tif 000066（ただし、ｎは屈折率、ｅ（ベクトル）は電場ベクトル、ｄ（ベクトル）は電気変位ベクトル、ｓ（ベクトル）は屈折光の波面法線方向の単位ベクトル）を最小にすることにより、電気変位ベクトルｄ（ベクトル）を求める工程を有することを特徴とする光線追跡方法。【選択図】　なし【特許請求の範囲】
【請求項１】
屈折率異方性のある媒質内での光線追跡であって、残差ベクトル
【数１】
2004028843.tif 000064
（ただし、ｎは屈折率、ｅ（ベクトル）は電場ベクトル、ｄ（ベクトル）は電気変位ベクトル、ｓ（ベクトル）は屈折光の波面法線方向の単位ベクトル）を最小にすることにより、電気変位ベクトルｄ（ベクトル）を求める工程を有することを特徴とする光線追跡方法。
【請求項２】
請求項１に記載の光線追跡方法であって、残差ベクトルを最小にする方法が、電気変位ベクトルｄ（ベクトル）の成分を未知数とするニュートン法であることを特徴とする光線追跡方法。
【請求項３】
屈折率異方性のある媒質に入射して屈折する光線の波面法線ベクトルと、屈折率を求める工程を有する光線追跡方法であって、請求項１に記載の残差ベクトルおよび、残差ベクトル
【数２】
2004028843.tif 000065
（ただし、ｃ（ベクトル）およびｓ（ベクトル）は、それぞれ入射光および屈折光の波面法線方向の単位ベクトル、ｑ（ベクトル）は境界面の法線ベクトル、ｎｃは屈折媒質に入射する前の媒質（空気等）の屈折率、ｎは入射後の媒質の屈折率）
の両者を最小にすることにより、ｄ（ベクトル）とｓ（ベクトル）を求める工程を有することを特徴とする光線追跡方法。
【請求項４】
請求項３に記載の光線追跡方法であって、残差ベクトルを最小にする方法が、電気変位ベクトルｄ（ベクトル）の成分と、波面法線方向の単位ベクトルｓ（ベクトル）の成分を未知数とするニュートン法であることを特徴とする光線追跡方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、偏光を利用する光学系または、半導体素子製造用の投影光学系や、天体観測用の光学系等であって、特に光線通過部材の、加熱や外力による内部応力分布の変動による光線の偏光分離が問題となるような光学系における光線追跡方法に関するものである。
【０００２】
なお、本明細書においてはベクトルが多用されるが、電子出願の制約上、文章中においてベクトル記号を使用することができないので、数式中においてベクトル記号を用いて例えば
【０００３】
【数３】
2004028843.tif 000002
【０００４】
と表されるものを、文章中においてはｘ（ベクトル）と表すことにする。
【０００５】
【従来の技術】
従来、複屈折がある硝材を光線が透過する場合の光線追跡については、例えば、１軸結晶を光線が入射する場合の屈折光線の屈折率と波面放線方向を求める方法について、Ｍ．Ｃ．Ｓｉｍｏｎによって、“Ｒａｙ　ｔｒａｃｉｎｇ　ｆｏｒｍｕｌａｓ　ｆｏｒ　ｍｏｎｏａｘｉａｌ　ｏｐｔｉｃａｌ　ｃｏｍｐｏｎｅｎｔｓ“という題名の文献で紹介された。（Ａｐｐｌ．Ｏｐｔ．（１９８３）ｐｐ３５４−３６０）に記載されている。
【０００６】
また、一般的に２軸結晶に光線が入射する場合の屈折光線の屈折率と波面法線ベクトルを求める方法の原理については、ボルン・ウオルフ著：「光学の原理ＩＩＩ」（邦訳：東海大学出版会）において、Ｆｒｅｓｎｅｌの公式（同文献、Ｐ９８０参照）を利用して求める方法が示されている。（同文献、Ｐ９９９〜Ｐ１０００参照）
【０００７】
【発明が解決しようとする課題】
しかし、ガラスの温度上昇分布に差がある場合や、ガラスの一部に外力を加えた場合、ガラス内で歪みが発生し、内部応力が発生するが、発生する内部応力は、通常方向異方性があるので、屈折率異方性があり、光線がガラスを透過すると、偏光になる。一方、内部応力は、ガラス内の位置により異なっていて、それに従って、屈折率も異なるので、ｇｒａｄｉｅｔ　ｉｎｄｅｘ光線追跡をする必要がある。
【０００８】
内部応力のある媒質内でのｇｒａｄｉｅｎｔ　ｉｎｄｅｘ光線追跡では、計算誤差が蓄積しないように、屈折率の計算は高精度に行う必要があるが、前記のＭ．Ｃ．Ｓｉｍｏｎの文献等で述べられているような屈折率を計算する方法では、ｇｒａｄｉｅｎｔ　ｉｎｄｅｘ光線追跡に絶えられるような精度が保証できないことがわかった。そこで、本発明者は、波面法線楕円体によって、屈折率を計算する方法を開発した（橋本純夫：「内部応力分布があるガラス内の光線追跡」光学２６（１９９７）　６７７：以下この論文を「先行論文」という）。
【０００９】
この方法は、光学パラメータ及び光弾性パラメータを用いて屈折光の波面法線ベクトルを求め、これから「波面法線楕円体を考えて、波面法線楕円体の原点を通り、波面法線ベクトルと直交する平面と楕円体の光線は楕円になる。このとき、この楕円の主軸の長さは、位相速度の逆数に比例し、その主軸の方向は、電気変位ベクトルの振動方向に一致している。」という原理に基づいて、電気変位ベクトルを求め、さらに、この電気変位ベクトルと誘電テンソルの関係式により電気ベクトルを求めると共に、２つの電気変位ベクトルに直交する方向を磁気ベクトルの方向とし、このようにして求められた電気変位ベクトルと磁気ベクトルの外積からポインティングベクトル、すなわち光線の進行方法を求めるものである。
【００１０】
しかしながら、この方法での計算精度がどの程度であるかは、明確にされてはいなかった。また、ｇｒａｄｉｅｎｔ　ｉｎｄｅｘ光線追跡によって、Ｍａｘｗｅｌｌ方程式を満足しているかどうかも明白でなかった。
【００１１】
本発明はこのような事情に鑑みてなされたもので、前記本発明者が開発した先行論文の方法よりも高精度で、またＭａｘｗｅｌｌ方程式を満足するような形で、光線追跡を行う方法を提供することを課題とする。
【００１２】
【課題を解決するための手段】
前記課題を解決するための第１の手段は、屈折率異方性のある媒質内での光線追跡であって、残差ベクトル
【００１３】
【数４】
2004028843.tif 000003
【００１４】
（ただし、ｎは屈折率、ｅ（ベクトル）は電場ベクトル、ｄ（ベクトル）は電気変位ベクトル、ｓ（ベクトル）は屈折光の波面法線方向の単位ベクトル）を最小にすることにより、電気変位ベクトルｄ（ベクトル）を求める工程を有することを特徴とする光線追跡方法（請求項１）である。
【００１５】
前記課題を解決するための第２の手段は、前記第１の手段であって、残差ベクトルを最小にする方法が、電気変位ベクトルｄ（ベクトル）の成分を未知数とするニュートン法であることを特徴とするもの（請求項２）である。
【００１６】
前記課題を解決するための第３の手段は、屈折率異方性のある媒質に入射して屈折する光線の波面法線ベクトルと、屈折率を求める工程を有する光線追跡方法であって、請求項１に記載の残差ベクトルおよび、残差ベクトル
【００１７】
【数５】
2004028843.tif 000004
【００１８】
（ただし、ｃ（ベクトル）およびｓ（ベクトル）は、それぞれ入射光および屈折光の波面法線方向の単位ベクトル、ｑ（ベクトル）は境界面の法線ベクトル、ｎｃは屈折媒質に入射する前の媒質（空気等）の屈折率、ｎは入射後の媒質の屈折率）
の両者を最小にすることにより、ｄ（ベクトル）とｓ（ベクトル）を求める工程を有することを特徴とする光線追跡方法（請求項３）である。
【００１９】
前記課題を解決するための第４の手段は、前記第３の手段であって、残差ベクトルを最小にする方法が、電気変位ベクトルｄ（ベクトル）の成分と、波面法線方向の単位ベクトルｓ（ベクトル）の成分を未知数とするニュートン法であることを特徴とするもの（請求項４）である。
【００２０】
本発明者は、先ず、屈折率異方性のある媒質内で、Ｍａｘｗｅｌｌ方程式を満足しているかどうかを検証するための、残差ベクトルを定義した。そして、その残差ベクトルを最小にするように、電気変位ベクトルｄ（ベクトル）の成分の値を求めることにより、光線追跡を行うことにした。ｄ（ベクトル）の成分を求めるには、これを未知数とするニュートン法により求めることが好ましい。
【００２１】
同様に、空気、または真空などの等方性媒質から屈折率異方性のある媒質内へ光線が入射する場合において、前記の残差ベクトルとともに、入射光の光線ベクトルと屈折光の波面法線ベクトルｓ（ベクトル）に関する残差ベクトルを定義して、両方の残差ベクトルを共に最小にするような、電気変位ベクトルｄ（ベクトル）の成分と、波面法線ベクトルｓ（ベクトル）の成分の値を求めることにより光線追跡を行うことにした。電気変位ベクトルｄ（ベクトル）の成分と、波面法線ベクトルｓ（ベクトル）の成分の値を求めるには、これらを未知数とするニュートン法により求めることが好ましい。
【００２２】
以下、これらの事項について詳細に説明する。
先ず、Ｍａｘｗｅｌｌ方程式から始めて、屈折率異方性のある媒質内での波面法線方向の単位ベクトルｓ（ベクトル）と屈折率ｎを計算する方法の精度を検証する方法について説明する。
【００２３】
光線の磁場、および電場、電気変位ベクトルをそれぞれ、
【００２４】
【数６】
2004028843.tif 000005
【００２５】
とする。ただし、ｋ０は真空中の波数であり、φはｅｉｋｏｎａｌ
【００２６】
【数７】
2004028843.tif 000006
【００２７】
である。ここで、ｒ（ベクトル）は光線の位置ベクトル、ｓ（ベクトル）は波面法線方向の単位ベクトル、ｎは光速ｃと位相速度ｖｐとによりｎ＝ｃ／ｖｐ　と表される屈折率である。ｓ（ベクトル）は、例えば先行論文に詳述されているような方法により求めることができるが、前記第１の手段、第２の手段においては、これを求める方法については特に規定しておらず、他の公知の方法によって求めてもよい。これらをＭａｘｗｅｌｌ方程式に代入すると、
【００２８】
【数８】
2004028843.tif 000007
となり、（４）式に（５）式を代入することにより、
【００２９】
【数９】
2004028843.tif 000008
【００３０】
であるから、
【００３１】
【数１０】
2004028843.tif 000009
【００３２】
となる。非等方媒質の場合、等方媒質の場合とは違って、ｄ（ベクトル）とｅ（ベクトル）は同じ方向のベクトルでないため、一般に、
【００３３】
【数１１】
2004028843.tif 000010
【００３４】
である。よって、（６）式から、ｅｉｋｏｎａｌ方程式
【００３５】
【数１２】
2004028843.tif 000011
【００３６】
を導出することはできない。しかし、ｎをｃ／ｖｐと表されるような屈折率とする限り、幾何光学的な波面の考察から、
【００３７】
【数１３】
2004028843.tif 000012
【００３８】
が成立するので、この式の両辺を２乗すると、ｅｉｋｏｎａｌ方程式（１１）は成立する。（９）式に（１２）式を代入すると、
【００３９】
【数１４】
2004028843.tif 000013
【００４０】
となる。そこで、
【００４１】
【数１５】
2004028843.tif 000014
【００４２】
と置いて、残差ベクトル
【００４３】
【数１６】
2004028843.tif 000015
【００４４】
を計算すると、先行論文において波面法線楕円体によりｄ（ベクトル）とｅ（ベクトル）を計算した方法の、計算精度がどの程度であるかがわかる。波面法線楕円体から計算された残差ベクトルは、内部応力が大きく成る程、大きくなる。
【００４５】
前記第１の手段においては、先行論文と異なり、（１４）式で定義される残差ベクトルが最小となるようにｄ（ベクトル）を定める。この方法は、先行文献と異なり、Ｍａｘｗｅｌｌの方程式を前提としているので、より正確にｄ（ベクトル）を求めることができる。
【００４６】
後述するように、ｎはｄ（ベクトル）の成分に関する２次式になるので、残差ベクトル
【００４７】
【数１７】
2004028843.tif 000016
【００４８】
は、ｄ（ベクトル）の成分に関する非線形な式になる。そこで、
【００４９】
【数１８】
2004028843.tif 000017
【００５０】
をゼロにするようなｄ（ベクトル）を求めるためには、Ｎｅｗｔｏｎ法を用いることが好ましい。［Ｔ．Ｒ．マッカーラ：計算機のための数値計算法概論（サイエンス社、１９８４）　ｐ．７３参照］、以下、Ｎｅｗｔｏｎ法を用いた例について説明する。
すなわち、ｉ＝１〜３を、直交座標系のｘ、ｙ、ｚに対応させた場合、
【００５１】
【数１９】
2004028843.tif 000018
【００５２】
の未知数が３個の連立方程式を満たす解Δｄｊ（ｊ＝１，３）を計算し、ｄ（ベクトル）のさらに良い近似値を
【００５３】
【数２０】
2004028843.tif 000019
【００５４】
とし、この計算を繰り返して、高精度なｄ（ベクトル）の値を求める。
（１０）式において、ｓ（ベクトル）が与えられた場合の
【００５５】
【数２１】
2004028843.tif 000020
【００５６】
の成分ｆｉを、ｄ（ベクトル）の成分ｄｊで微分すると、
【００５７】
【数２２】
2004028843.tif 000021
【００５８】
となる。ここで、δｉｊはＫｒｏｎｅｃｋｅｒのδである。
【００５９】
以下、内部応力がある場合の屈折率の計算方法を示す。先ず、応力光学定数ｑ１１、ｑ１２と、光弾性定数Ｃ１、Ｃ２との関係は、
【００６０】
【数２３】
2004028843.tif 000022
【００６１】
である。そして、楕円の主軸と座標軸とが一致しない一般の屈折率楕円体の式を、
【００６２】
【数２４】
2004028843.tif 000023
【００６３】
とすると、βｉｊとｑ１１、ｑ１２の関係式は、応力が加わる前の屈折率ｎ０と、境界面と入射光線の交点上の応力分布σｉｊ（ｉ，ｊ＝ｘ，ｙ，ｚ）により、
【００６４】
【数２５】
2004028843.tif 000024
【００６５】
となる。（２０）式が成り立つためには、応力光学定数ｑ１１，ｑ１２とｑ４４の関係において、
ｑ４４＝２（ｑ１１−ｑ１２）　…（２１）
の関係式が成立する必要があるが、これは以下のようにして成立することが証明される。すなわち、応力によって変化した波面法線楕円体
【００６６】
【数２６】
2004028843.tif 000025
【００６７】
において、立方晶系以上の対称性がある場合、
【００６８】
【数２７】
2004028843.tif 000026
【００６９】
等となるが　［ボルン・ウオルフ：光学の原理ＩＩＩ（東海大学出版会、１９７６）　　ｐ．１０２６〜１０２９参照］、ガラスのように等方媒質の場合、応力のない場合の屈折率をｎ０とすると、
【００７０】
【数２８】
2004028843.tif 000027
【００７１】
より、
【００７２】
【数２９】
2004028843.tif 000028
【００７３】
となる。等方媒質は任意の座標変換に対して不変であるので、ｘｙ座標系をθ回転させて、ｘ’ｙ’座標系になったとすると、ｘ’ｙ’座標系での主応力σｘ’ｘ’，σｙ’ｙ’は、ｘｙ座標系での主応力σｘｘ，σｙｙおよびせん断応力σｘｙにより、
【００７４】
【数３０】
2004028843.tif 000029
【００７５】
と表される。原点と、ｘ’軸と波面法線楕円体の交点の距離をｒ’とすると、（２２）式のｘ，ｙ，ｚをｘ’，ｙ’，ｚ’に置き換えた式に、
（ｘ’，ｙ’，ｚ’）＝（ｒ’，０，０）
を代入して、（２７）式、（２８）式より、
【００７６】
【数３１】
2004028843.tif 000030
【００７７】
一方、（２２）式において、
（ｘ，ｙ，ｚ）＝（ｒｃｏｓθ，ｒｓｉｎθ，０）
を代入すると、（２４）式、（２５）式により、
【００７８】
【数３２】
2004028843.tif 000031
【００７９】
ｘ軸をθ回転させたものがｘ’軸になることにより、（２９）式のｒ’と（３０）式のｒは同じなので、σｘｙの係数が等しいことにより、
ｑ４４＝２（ｑ１１−ｑ１２）　…（３１）
となり、（２１）式が成立する。
【００８０】
また、波面法線楕円体を考慮することにより、
【００８１】
【数３３】
2004028843.tif 000032
【００８２】
となる。また、
【００８３】
【数３４】
2004028843.tif 000033
【００８４】
である。これを証明すると、以下のようになる。すなわち、一般の、主軸が座標軸と一致しない屈折率楕円体
【００８５】
【数３５】
2004028843.tif 000034
【００８６】
が座標変換
【００８７】
【数３６】
2004028843.tif 000035
【００８８】
によって楕円体の主軸と座標軸が一致するように直交変換した場合、すなわち、屈折率楕円体を
【００８９】
【数３７】
2004028843.tif 000036
【００９０】
という形に直交変換した場合を考える。（３６）式に（３５）式を代入して、ｉとｊを入れ替えて、
【００９１】
【数３８】
2004028843.tif 000037
【００９２】
（３４）式と比較して、
【００９３】
【数３９】
2004028843.tif 000038
【００９４】
また、（３４）式の座標系に相当する電気ベクトル、および電気変位ベクトルの成分をｅｉ、ｄｉ、（３６）式の座標系に相当する電気ベクトルおよび電気変位ベクトルの成分を、ｅ’ｉ、ｄ’ｉとし、また、直交変換マトリックスはＣｉｊ−１＝Ｃｊｉであるので、
【００９５】
【数４０】
2004028843.tif 000039
【００９６】
（３６）式のように、屈折率楕円体の主軸と座標軸が一致している場合、
β’ｉｉ＝１／ｎ’ｉ２　…（４１）
であり、また、
ｎ’ｉ２＝ε’ｉμ　…（４２）
ｄ’ｉ＝ε’ｉｅ’ｉ　…（４３）
ｅ’ｉ＝ｄ’ｉ／ε’ｉ＝μβ’ｉｉｄ’ｉ　…（４４）
（４０）式に、（４４）式、（３９）式を代入して、
【００９７】
【数４１】
2004028843.tif 000040
【００９８】
（３８）式より、
【００９９】
【数４２】
2004028843.tif 000041
【０１００】
が求められ、これは（３３）式と同等となる。
【０１０１】
（１９）式、（２０）式、（３２）式、（３３）式等をｄｊで微分して後、（１６）式に代入して
∂ｆｉ／∂ｄｊを求め、（１５）式に代入することにより、この問題を解くことができる。
【０１０２】
そして、ｇｒａｄｉｅｎｔ　ｉｎｄｅｘ光線追跡には、∂ｎ／∂ｘｊの値が必要であるが、それは以下のような方法で求められる。先ず、（１４）式に（３２）式を代入すると、その成分は、
【０１０３】
【数４３】
2004028843.tif 000042
【０１０４】
となり、（１３）式の両辺をｘｊで微分すると、
∂ｆｉ／∂ｘｊ＝０…（４８）
となることから、
【０１０５】
【数４４】
2004028843.tif 000043
【０１０６】
また、（３２）式より、
【０１０７】
【数４５】
2004028843.tif 000044
【０１０８】
となり、これを（４９）式に代入することにより、∂ｄｋ／∂ｘｊに関する連立一次方程式が成立し、求めた∂ｄｋ／∂ｘｊの値を、（５０）式に代入することにより∂ｎ／∂ｘｊが求められる。
【０１０９】
以上のように、Ｎｅｗｔｏｎ法でｄ（ベクトル）を定めた後、（３３）式により、ｅ（ベクトル）を定める。光線の進行方向のベクトル（ポインティングベクトル）ｒ（ベクトル）は、ｐｏｉｎｔｉｎｇ　ｖｅｃｔｏｒの方向なので、ｒ（ベクトル）は、ｄ（ベクトル）とｓ（ベクトル）の双方に垂直な磁気ベクトルｈ（ベクトル）により、
【０１１０】
【数４６】
2004028843.tif 000045
【０１１１】
により求められる。以上の説明に置いては、ｄ（ベクトル）を定めるのにニュートン法を用いたが、（１４）式で定義される残差ベクトルを最小にするｄ（ベクトル）を求めるためには、他の公知の最適化方法を用いてもよい。
【０１１２】
以上を要約すると、前記第１の手段においては、公知の方法によりｓ（ベクトル）を求め、これから（１４）式で定義される残差ベクトルを最小にするｄ（ベクトル）を求め（（１４）式で必要とされるｅ（ベクトル）は（３３）式や、その他の公知の関係式を用いて求める）、求まったｄ（ベクトル）とｅ（ベクトル）に垂直なｈ（ベクトル）を求め、光線の進行方向を示すポインティングベクトルｒ（ベクトル）をｅ（ベクトル）とｈ（ベクトル）の外積として求める。
【０１１３】
微小距離進行後のｓ（ベクトル）の変化は、光線の微分方程式を差分近似するか、Ｒｕｎｇｅ−Ｋｕｔｔａ法　［Ｔ．Ｒ．マッカーラ：計算機のための数値計算法概論（サイエンス社、１９８４）　ｐ．２７７　参照］　により求める。
【０１１４】
前記先行論文では、「ガラスを層状に分割した場合、各層に入射するごとに、電気変位ベクトルが互いに垂直な２つの偏光に分離するために、ガラス透過後の光線は光束となると考えられる」としたが、直角に近い角度へ電気変位ベクトルが推移するような偏光の成分はごく僅かであると考えられるので、無限に近い本数の光線になるとは考えにくい。すなわち、図１のように、ガラス内のある層で、互いに垂直な電気変位ベクトルｄ１（ベクトル）、ｄ２（ベクトル）の偏光
【０１１５】
【数４７】
2004028843.tif 000046
【０１１６】
があり、次の層で、電気変位ベクトルおよび偏光がそれぞれｄ’１（ベクトル）、ｄ’２（ベクトル）、および
【０１１７】
【数４８】
2004028843.tif 000047
【０１１８】
に推移し、ｄ１（ベクトル）とｄ’１（ベクトル）のなす角度がΔθ≒０であるとすると、それらの光線の振幅は、
【０１１９】
【数４９】
2004028843.tif 000048
【０１２０】
となり、大部分がもとの電気変位ベクトルの方向とほぼ同じ方向の電気変位ベクトルの光線に推移すると考えられる。
【０１２１】
レンズ内部の光線追跡においては、屈折率が連続的に変化するので、ｓ（ベクトル）があらかじめ与えられたとして微小距離進行後のｓ（ベクトル）の変化を逐次計算すれば良いが、レンズ表面の屈折の場合は、ｓ（ベクトル）が屈折前後で急激に変化するために、ｓ（ベクトル）とｄ（ベクトル）の双方を同時に求める必要がある。
【０１２２】
そのため、前記第３の手段においては、一般的な屈折光線の方向余弦を求める式から出発する。すなわち、ｃ（ベクトル）およびｓ（ベクトル）を、それぞれ入射光および屈折光の波面法線方向の単位ベクトルとし、境界面の法線ベクトルをｑ（ベクトル）＝（ｑｘ，ｑｙ，ｑｚ）とし、屈折媒質に入射する前の媒質（空気等）の屈折率をｎｃ、入射後の媒質の屈折率をｎとすると、
【０１２３】
【数５０】
2004028843.tif 000049
【０１２４】
となり、
【０１２５】
【数５１】
2004028843.tif 000050
【０１２６】
となる。［三宅和夫：幾何光学（共立出版，１９８５）　ｐ．３０参照　］
【０１２７】
【数５２】
2004028843.tif 000051
【０１２８】
と置いて、
【０１２９】
【数５３】
2004028843.tif 000052
【０１３０】
が成り立つように、ｓ（ベクトル）を求める。屈折後の媒質が非等方性の場合、
【０１３１】
【数５４】
2004028843.tif 000053
【０１３２】
はｄ（ベクトル）の関数でもあるので、
【０１３３】
【数５５】
2004028843.tif 000054
【０１３４】
とし、（１４）式の
【０１３５】
【数５６】
2004028843.tif 000055
【０１３６】
と連立させて解く。この連立方程式を解くのには、Ｎｅｗｔｏｎ法を使用するのが好ましいが、他の公知の最適化手法を用いてもよい。Ｎｅｗｔｏｎ法を用いた用法について、以下に説明する。すなわち、ｄ（ベクトル）と同様にｓ（ベクトル）のさらに良い近似値を
【０１３７】
【数５７】
2004028843.tif 000056
【０１３８】
とすると、
【０１３９】
【数５８】
2004028843.tif 000057
【０１４０】
となるが、（６２）式、（６３）式の∂ｆｉ／∂ｓｉおよび∂ｇｉ／∂ｓｉ　（ｉ，ｊ＝ｙ，ｚ）　は、
【０１４１】
【数５９】
2004028843.tif 000058
【０１４２】
であることに注意すれば、未知数が５個の連立方程式に還元できる。このとき、（６２）式、（６３）式の微係数は（１６）式に加えて、
【０１４３】
【数６０】
2004028843.tif 000059
【０１４４】
が必要となる。なお、あまり精度が必要でない場合、波面法線楕円体を用いてｄ（ベクトル）をｓ（ベクトル）から求めることにより、単純に、
【０１４５】
【数６１】
2004028843.tif 000060
【０１４６】
から、Δｓ（ベクトル）のみを未知数にすればよいが、（６４）式の関係により、未知数を２個に還元できる。この場合、（６９）式の代わりに、
【０１４７】
【数６２】
2004028843.tif 000061
【０１４８】
とする必要がある。∂ｎ／∂ｓｊを含む項は、ｄ（ベクトル）をｓ（ベクトル）の関数とすることにより、新たに付加した項である。なお、（６２），（６３）式によりｓ（ベクトル）とｄ（ベクトル）を求めるときに、収束が遅い場合があるので、（７０）式でｓ（ベクトル）を求めて後、それらを初期値として、（６２），（６３）式でさらに高精度な計算をするとよい。このような計算手段は、前記第４の手段の１つの方法であることは明らかであり、前記第４の手段に含まれることは言うまでもない。
【０１４９】
【実施例】
以上のように述べた計算方法で、実施例として、外部応力によるガラス内の屈折率変動による収差変動を計算した。計算例として、図２のような球面収差のみ補正されたレンズ群に、波長２４８ｎｍ、半径５０ｍｍの光束を透過させた。このレンズ群の光学的諸元を表１に示す
【０１５０】
【表１】
2004028843.tif 000062
【０１５１】
外半径５２．５ｍｍのレンズＧ１（ガラス）の子午線方向（ｙ軸とする）の両端に、１０Ｎ／ｃｍ２の外部応力を加えた。計算に用いたガラスのヤング率は、０．７２×１０５Ｎ／ｍｍ２，ポアソン比は０．１６４，光弾性定数は、Ｃ１＝６．４９５×１０−７ｍｍ２／Ｎ，Ｃ２＝４．１７６×１０−６ｍｍ２／Ｎである。
【０１５２】
この場合のガラス内の応力分布は、図３のようになっている。図３（Ａ）はｙ軸方向の応力分布で、図３（Ｂ）はｚ軸方向（光軸およびｙ軸の両方に垂直な方向）の応力分布である。等高線の間隔は１Ｎ／ｃｍ２である。
【０１５３】
このような外力を加えた場合において、レンズＧ１に半径５０ｍｍの光束を透過させ、ｇｒａｄｉｅｎｔ　ｉｎｄｅｘ光線追跡をして波面収差を計算し、外力を加えない場合の波面収差との差をとり、波面収差変動を計算した。その結果、波面収差変動は図４のようになった。電気変位ベクトルがｙ軸方向の偏光を、前記のような応力分布のガラスに透過させて、像面上でｙ軸方向の偏光成分の波面収差変動を計算したのが図４（Ａ）であり、電気変位ベクトルがｚ軸方向の偏光を、同様な応力分布のガラスに透過させて、像面上でｚ軸方向の偏光成分の波面収差変動を計算したのが図４（Ｂ）である。図（Ａ）における等高線の間隔は０．００３２λ、図（Ｂ）における等高線の間隔は０．００４０λである。図４（Ａ）の差の最大値は、０．０３２λであり、図４（Ｂ）の差の最大値は、０．０４０λである。
【０１５４】
図３と図４とを比較すると、外力によるレンズの内部応力分布と偏光による波面収差変動の分布は、ほぼ対応するが、分布の形状がｙ軸とｚ軸とで互いに逆になっているのは、（３２）式において、Ｃ１よりもＣ２の絶対値が大きいことによるものと考えられる。
【０１５５】
なお、ｙ軸方向の偏光を同じ状態のガラスに透過させて、像面上でｚ軸方向の偏光成分を計算したものと、ｚ軸方向の偏光を同じ状態のガラスに透過させて、像面上でｙ軸方向の偏光成分を計算したものとは、図３（Ａ）および図３（Ｂ）の偏光よりも、かなり、強度が小さく無視できる程である。
【０１５６】
また、このような応力状態の場合、今回新たに発明した計算方法と、以前発明した波面法線楕円体による計算方法の精度を、（１１）式の残差ベクトルの絶対値の平均値を計算することにより、比較した。これを表２に示す。
【０１５７】
（表２）本発明による計算方法（Ａ）と先行論文による計算方法（Ｂ）による残差ベクトルの絶対値の平均値の比較
【０１５８】
【表２】
2004028843.tif 000063
【０１５９】
先行論文による計算方法（波面法線楕円体による計算方法）（Ｂ）の場合、外部応力が１０倍になると、残差ベクトルの絶対値がほぼ１００倍になるのに対し、本発明の計算方法（Ａ）によれば、外部応力が１０倍になっても、残差ベクトルの絶対値はほとんど変わらず、計算機の倍精度計算の誤差の範囲内に収まっていることが分かる。なお、波面法線楕円体による計算方法の場合、残差ベクトルがゼロにならない原因は、この方法がＭａｘｗｅｌｌ方程式を満足していないことにあると考えるよりも、平方根などの関数を使用する際の計算誤差にあると考える方がよい。
【０１６０】
【発明の効果】
以上説明したように、本発明の計算方法により、硝材内に不均一な応力分布がある場合の光線追跡において、よりＭａｘｗｅｌｌ方程式に近い式から計算精度を検証することにより、従来よりも高精度に計算することが可能になり、計算の信頼性を増した状態で光線追跡を行うことができる。
【図面の簡単な説明】
【図１】偏光の振幅の微小な推移を示す図である。
【図２】本発明の実施例に用いた球面収差のみ補正されたレンズの光路図である。
【図３】本発明の実施例におけるガラス中の応力分布の例を示す図である。
【図４】本発明の実施例においてガラスに応力を与えた場合の、波面収差変動の計算結果の例を示す図である。
【図１】
2004028843.tif 000067
【図２】
2004028843.tif 000068
【図３】
2004028843.tif 000069
【図４】
2004028843.tif 000070