【０００１】
【発明の属する技術分野】
本発明は、偏光を利用する光学系または、半導体素子製造用の投影光学系や、天体観測用の光学系等であって、特に光線通過部材の、加熱や外力による内部応力分布の変動による光線の偏光分離が問題となるような光学系における光線追跡方法に関するものである。
【０００２】
なお、本明細書においてはベクトルが多用されるが、電子出願の制約上、文章中においてベクトル記号を使用することができないので、数式中においてベクトル記号を用いて例えば
【０００３】
【数３】
【０００４】
と表されるものを、文章中においてはｘ（ベクトル）と表すことにする。
【０００５】
【従来の技術】
従来、複屈折がある硝材を光線が透過する場合の光線追跡については、例えば、１軸結晶を光線が入射する場合の屈折光線の屈折率と波面放線方向を求める方法について、Ｍ．Ｃ．Ｓｉｍｏｎによって、“Ｒａｙ ｔｒａｃｉｎｇ ｆｏｒｍｕｌａｓ ｆｏｒ ｍｏｎｏａｘｉａｌ ｏｐｔｉｃａｌ ｃｏｍｐｏｎｅｎｔｓ“という題名の文献で紹介された。（Ａｐｐｌ．Ｏｐｔ．（１９８３）ｐｐ３５４−３６０）に記載されている。
【０００６】
また、一般的に２軸結晶に光線が入射する場合の屈折光線の屈折率と波面法線ベクトルを求める方法の原理については、ボルン・ウオルフ著：「光学の原理ＩＩＩ」（邦訳：東海大学出版会）において、Ｆｒｅｓｎｅｌの公式（同文献、Ｐ９８０参照）を利用して求める方法が示されている。（同文献、Ｐ９９９〜Ｐ１０００参照）
【０００７】
【発明が解決しようとする課題】
しかし、ガラスの温度上昇分布に差がある場合や、ガラスの一部に外力を加えた場合、ガラス内で歪みが発生し、内部応力が発生するが、発生する内部応力は、通常方向異方性があるので、屈折率異方性があり、光線がガラスを透過すると、偏光になる。一方、内部応力は、ガラス内の位置により異なっていて、それに従って、屈折率も異なるので、ｇｒａｄｉｅｔ ｉｎｄｅｘ光線追跡をする必要がある。
【０００８】
内部応力のある媒質内でのｇｒａｄｉｅｎｔ ｉｎｄｅｘ光線追跡では、計算誤差が蓄積しないように、屈折率の計算は高精度に行う必要があるが、前記のＭ．Ｃ．Ｓｉｍｏｎの文献等で述べられているような屈折率を計算する方法では、ｇｒａｄｉｅｎｔ ｉｎｄｅｘ光線追跡に絶えられるような精度が保証できないことがわかった。そこで、本発明者は、波面法線楕円体によって、屈折率を計算する方法を開発した（橋本純夫：「内部応力分布があるガラス内の光線追跡」光学２６（１９９７） ６７７：以下この論文を「先行論文」という）。
【０００９】
この方法は、光学パラメータ及び光弾性パラメータを用いて屈折光の波面法線ベクトルを求め、これから「波面法線楕円体を考えて、波面法線楕円体の原点を通り、波面法線ベクトルと直交する平面と楕円体の光線は楕円になる。このとき、この楕円の主軸の長さは、位相速度の逆数に比例し、その主軸の方向は、電気変位ベクトルの振動方向に一致している。」という原理に基づいて、電気変位ベクトルを求め、さらに、この電気変位ベクトルと誘電テンソルの関係式により電気ベクトルを求めると共に、２つの電気変位ベクトルに直交する方向を磁気ベクトルの方向とし、このようにして求められた電気変位ベクトルと磁気ベクトルの外積からポインティングベクトル、すなわち光線の進行方法を求めるものである。
【００１０】
しかしながら、この方法での計算精度がどの程度であるかは、明確にされてはいなかった。また、ｇｒａｄｉｅｎｔ ｉｎｄｅｘ光線追跡によって、Ｍａｘｗｅｌｌ方程式を満足しているかどうかも明白でなかった。
【００１１】
本発明はこのような事情に鑑みてなされたもので、前記本発明者が開発した先行論文の方法よりも高精度で、またＭａｘｗｅｌｌ方程式を満足するような形で、光線追跡を行う方法を提供することを課題とする。
【００１２】
【課題を解決するための手段】
前記課題を解決するための第１の手段は、屈折率異方性のある媒質内での光線追跡であって、残差ベクトル
【００１３】
【数４】
【００１４】
（ただし、ｎは屈折率、ｅ（ベクトル）は電場ベクトル、ｄ（ベクトル）は電気変位ベクトル、ｓ（ベクトル）は屈折光の波面法線方向の単位ベクトル）を最小にすることにより、電気変位ベクトルｄ（ベクトル）を求める工程を有することを特徴とする光線追跡方法（請求項１）である。
【００１５】
前記課題を解決するための第２の手段は、前記第１の手段であって、残差ベクトルを最小にする方法が、電気変位ベクトルｄ（ベクトル）の成分を未知数とするニュートン法であることを特徴とするもの（請求項２）である。
【００１６】
前記課題を解決するための第３の手段は、屈折率異方性のある媒質に入射して屈折する光線の波面法線ベクトルと、屈折率を求める工程を有する光線追跡方法であって、請求項１に記載の残差ベクトルおよび、残差ベクトル
【００１７】
【数５】
【００１８】
（ただし、ｃ（ベクトル）およびｓ（ベクトル）は、それぞれ入射光および屈折光の波面法線方向の単位ベクトル、ｑ（ベクトル）は境界面の法線ベクトル、ｎｃは屈折媒質に入射する前の媒質（空気等）の屈折率、ｎは入射後の媒質の屈折率）
の両者を最小にすることにより、ｄ（ベクトル）とｓ（ベクトル）を求める工程を有することを特徴とする光線追跡方法（請求項３）である。
【００１９】
前記課題を解決するための第４の手段は、前記第３の手段であって、残差ベクトルを最小にする方法が、電気変位ベクトルｄ（ベクトル）の成分と、波面法線方向の単位ベクトルｓ（ベクトル）の成分を未知数とするニュートン法であることを特徴とするもの（請求項４）である。
【００２０】
本発明者は、先ず、屈折率異方性のある媒質内で、Ｍａｘｗｅｌｌ方程式を満足しているかどうかを検証するための、残差ベクトルを定義した。そして、その残差ベクトルを最小にするように、電気変位ベクトルｄ（ベクトル）の成分の値を求めることにより、光線追跡を行うことにした。ｄ（ベクトル）の成分を求めるには、これを未知数とするニュートン法により求めることが好ましい。
【００２１】
同様に、空気、または真空などの等方性媒質から屈折率異方性のある媒質内へ光線が入射する場合において、前記の残差ベクトルとともに、入射光の光線ベクトルと屈折光の波面法線ベクトルｓ（ベクトル）に関する残差ベクトルを定義して、両方の残差ベクトルを共に最小にするような、電気変位ベクトルｄ（ベクトル）の成分と、波面法線ベクトルｓ（ベクトル）の成分の値を求めることにより光線追跡を行うことにした。電気変位ベクトルｄ（ベクトル）の成分と、波面法線ベクトルｓ（ベクトル）の成分の値を求めるには、これらを未知数とするニュートン法により求めることが好ましい。
【００２２】
以下、これらの事項について詳細に説明する。
先ず、Ｍａｘｗｅｌｌ方程式から始めて、屈折率異方性のある媒質内での波面法線方向の単位ベクトルｓ（ベクトル）と屈折率ｎを計算する方法の精度を検証する方法について説明する。
【００２３】
光線の磁場、および電場、電気変位ベクトルをそれぞれ、
【００２４】
【数６】
【００２５】
とする。ただし、ｋ０は真空中の波数であり、φはｅｉｋｏｎａｌ
【００２６】
【数７】
【００２７】
である。ここで、ｒ（ベクトル）は光線の位置ベクトル、ｓ（ベクトル）は波面法線方向の単位ベクトル、ｎは光速ｃと位相速度ｖｐとによりｎ＝ｃ／ｖｐ と表される屈折率である。ｓ（ベクトル）は、例えば先行論文に詳述されているような方法により求めることができるが、前記第１の手段、第２の手段においては、これを求める方法については特に規定しておらず、他の公知の方法によって求めてもよい。これらをＭａｘｗｅｌｌ方程式に代入すると、
【００２８】
【数８】
となり、（４）式に（５）式を代入することにより、
【００２９】
【数９】
【００３０】
であるから、
【００３１】
【数１０】
【００３２】
となる。非等方媒質の場合、等方媒質の場合とは違って、ｄ（ベクトル）とｅ（ベクトル）は同じ方向のベクトルでないため、一般に、
【００３３】
【数１１】
【００３４】
である。よって、（６）式から、ｅｉｋｏｎａｌ方程式
【００３５】
【数１２】
【００３６】
を導出することはできない。しかし、ｎをｃ／ｖｐと表されるような屈折率とする限り、幾何光学的な波面の考察から、
【００３７】
【数１３】
【００３８】
が成立するので、この式の両辺を２乗すると、ｅｉｋｏｎａｌ方程式（１１）は成立する。（９）式に（１２）式を代入すると、
【００３９】
【数１４】
【００４０】
となる。そこで、
【００４１】
【数１５】
【００４２】
と置いて、残差ベクトル
【００４３】
【数１６】
【００４４】
を計算すると、先行論文において波面法線楕円体によりｄ（ベクトル）とｅ（ベクトル）を計算した方法の、計算精度がどの程度であるかがわかる。波面法線楕円体から計算された残差ベクトルは、内部応力が大きく成る程、大きくなる。
【００４５】
前記第１の手段においては、先行論文と異なり、（１４）式で定義される残差ベクトルが最小となるようにｄ（ベクトル）を定める。この方法は、先行文献と異なり、Ｍａｘｗｅｌｌの方程式を前提としているので、より正確にｄ（ベクトル）を求めることができる。
【００４６】
後述するように、ｎはｄ（ベクトル）の成分に関する２次式になるので、残差ベクトル
【００４７】
【数１７】
【００４８】
は、ｄ（ベクトル）の成分に関する非線形な式になる。そこで、
【００４９】
【数１８】
【００５０】
をゼロにするようなｄ（ベクトル）を求めるためには、Ｎｅｗｔｏｎ法を用いることが好ましい。［Ｔ．Ｒ．マッカーラ：計算機のための数値計算法概論（サイエンス社、１９８４） ｐ．７３参照］、以下、Ｎｅｗｔｏｎ法を用いた例について説明する。
すなわち、ｉ＝１〜３を、直交座標系のｘ、ｙ、ｚに対応させた場合、
【００５１】
【数１９】
【００５２】
の未知数が３個の連立方程式を満たす解Δｄｊ（ｊ＝１，３）を計算し、ｄ（ベクトル）のさらに良い近似値を
【００５３】
【数２０】
【００５４】
とし、この計算を繰り返して、高精度なｄ（ベクトル）の値を求める。
（１０）式において、ｓ（ベクトル）が与えられた場合の
【００５５】
【数２１】
【００５６】
の成分ｆｉを、ｄ（ベクトル）の成分ｄｊで微分すると、
【００５７】
【数２２】
【００５８】
となる。ここで、δｉｊはＫｒｏｎｅｃｋｅｒのδである。
【００５９】
以下、内部応力がある場合の屈折率の計算方法を示す。先ず、応力光学定数ｑ１１、ｑ１２と、光弾性定数Ｃ１、Ｃ２との関係は、
【００６０】
【数２３】
【００６１】
である。そして、楕円の主軸と座標軸とが一致しない一般の屈折率楕円体の式を、
【００６２】
【数２４】
【００６３】
とすると、βｉｊとｑ１１、ｑ１２の関係式は、応力が加わる前の屈折率ｎ０と、境界面と入射光線の交点上の応力分布σｉｊ（ｉ，ｊ＝ｘ，ｙ，ｚ）により、
【００６４】
【数２５】
【００６５】
となる。（２０）式が成り立つためには、応力光学定数ｑ１１，ｑ１２とｑ４４の関係において、
ｑ４４＝２（ｑ１１−ｑ１２） …（２１）
の関係式が成立する必要があるが、これは以下のようにして成立することが証明される。すなわち、応力によって変化した波面法線楕円体
【００６６】
【数２６】
【００６７】
において、立方晶系以上の対称性がある場合、
【００６８】
【数２７】
【００６９】
等となるが ［ボルン・ウオルフ：光学の原理ＩＩＩ（東海大学出版会、１９７６） ｐ．１０２６〜１０２９参照］、ガラスのように等方媒質の場合、応力のない場合の屈折率をｎ０とすると、
【００７０】
【数２８】
【００７１】
より、
【００７２】
【数２９】
【００７３】
となる。等方媒質は任意の座標変換に対して不変であるので、ｘｙ座標系をθ回転させて、ｘ’ｙ’座標系になったとすると、ｘ’ｙ’座標系での主応力σｘ’ｘ’，σｙ’ｙ’は、ｘｙ座標系での主応力σｘｘ，σｙｙおよびせん断応力σｘｙにより、
【００７４】
【数３０】
【００７５】
と表される。原点と、ｘ’軸と波面法線楕円体の交点の距離をｒ’とすると、（２２）式のｘ，ｙ，ｚをｘ’，ｙ’，ｚ’に置き換えた式に、
（ｘ’，ｙ’，ｚ’）＝（ｒ’，０，０）
を代入して、（２７）式、（２８）式より、
【００７６】
【数３１】
【００７７】
一方、（２２）式において、
（ｘ，ｙ，ｚ）＝（ｒｃｏｓθ，ｒｓｉｎθ，０）
を代入すると、（２４）式、（２５）式により、
【００７８】
【数３２】
【００７９】
ｘ軸をθ回転させたものがｘ’軸になることにより、（２９）式のｒ’と（３０）式のｒは同じなので、σｘｙの係数が等しいことにより、
ｑ４４＝２（ｑ１１−ｑ１２） …（３１）
となり、（２１）式が成立する。
【００８０】
また、波面法線楕円体を考慮することにより、
【００８１】
【数３３】
【００８２】
となる。また、
【００８３】
【数３４】
【００８４】
である。これを証明すると、以下のようになる。すなわち、一般の、主軸が座標軸と一致しない屈折率楕円体
【００８５】
【数３５】
【００８６】
が座標変換
【００８７】
【数３６】
【００８８】
によって楕円体の主軸と座標軸が一致するように直交変換した場合、すなわち、屈折率楕円体を
【００８９】
【数３７】
【００９０】
という形に直交変換した場合を考える。（３６）式に（３５）式を代入して、ｉとｊを入れ替えて、
【００９１】
【数３８】
【００９２】
（３４）式と比較して、
【００９３】
【数３９】
【００９４】
また、（３４）式の座標系に相当する電気ベクトル、および電気変位ベクトルの成分をｅｉ、ｄｉ、（３６）式の座標系に相当する電気ベクトルおよび電気変位ベクトルの成分を、ｅ’ｉ、ｄ’ｉとし、また、直交変換マトリックスはＣｉｊ −１＝Ｃｊｉであるので、
【００９５】
【数４０】
【００９６】
（３６）式のように、屈折率楕円体の主軸と座標軸が一致している場合、
β’ｉｉ＝１／ｎ’ｉ ２ …（４１）
であり、また、
ｎ’ｉ ２＝ε’ｉμ …（４２）
ｄ’ｉ＝ε’ｉｅ’ｉ …（４３）
ｅ’ｉ＝ｄ’ｉ／ε’ｉ＝μβ’ｉｉｄ’ｉ …（４４）
（４０）式に、（４４）式、（３９）式を代入して、
【００９７】
【数４１】
【００９８】
（３８）式より、
【００９９】
【数４２】
【０１００】
が求められ、これは（３３）式と同等となる。
【０１０１】
（１９）式、（２０）式、（３２）式、（３３）式等をｄｊで微分して後、（１６）式に代入して
∂ｆｉ／∂ｄｊを求め、（１５）式に代入することにより、この問題を解くことができる。
【０１０２】
そして、ｇｒａｄｉｅｎｔ ｉｎｄｅｘ光線追跡には、∂ｎ／∂ｘｊの値が必要であるが、それは以下のような方法で求められる。先ず、（１４）式に（３２）式を代入すると、その成分は、
【０１０３】
【数４３】
【０１０４】
となり、（１３）式の両辺をｘｊで微分すると、
∂ｆｉ／∂ｘｊ＝０…（４８）
となることから、
【０１０５】
【数４４】
【０１０６】
また、（３２）式より、
【０１０７】
【数４５】
【０１０８】
となり、これを（４９）式に代入することにより、∂ｄｋ／∂ｘｊに関する連立一次方程式が成立し、求めた∂ｄｋ／∂ｘｊの値を、（５０）式に代入することにより∂ｎ／∂ｘｊが求められる。
【０１０９】
以上のように、Ｎｅｗｔｏｎ法でｄ（ベクトル）を定めた後、（３３）式により、ｅ（ベクトル）を定める。光線の進行方向のベクトル（ポインティングベクトル）ｒ（ベクトル）は、ｐｏｉｎｔｉｎｇ ｖｅｃｔｏｒの方向なので、ｒ（ベクトル）は、ｄ（ベクトル）とｓ（ベクトル）の双方に垂直な磁気ベクトルｈ（ベクトル）により、
【０１１０】
【数４６】
【０１１１】
により求められる。以上の説明に置いては、ｄ（ベクトル）を定めるのにニュートン法を用いたが、（１４）式で定義される残差ベクトルを最小にするｄ（ベクトル）を求めるためには、他の公知の最適化方法を用いてもよい。
【０１１２】
以上を要約すると、前記第１の手段においては、公知の方法によりｓ（ベクトル）を求め、これから（１４）式で定義される残差ベクトルを最小にするｄ（ベクトル）を求め（（１４）式で必要とされるｅ（ベクトル）は（３３）式や、その他の公知の関係式を用いて求める）、求まったｄ（ベクトル）とｅ（ベクトル）に垂直なｈ（ベクトル）を求め、光線の進行方向を示すポインティングベクトルｒ（ベクトル）をｅ（ベクトル）とｈ（ベクトル）の外積として求める。
【０１１３】
微小距離進行後のｓ（ベクトル）の変化は、光線の微分方程式を差分近似するか、Ｒｕｎｇｅ−Ｋｕｔｔａ法 ［Ｔ．Ｒ．マッカーラ：計算機のための数値計算法概論（サイエンス社、１９８４） ｐ．２７７ 参照］ により求める。
【０１１４】
前記先行論文では、「ガラスを層状に分割した場合、各層に入射するごとに、電気変位ベクトルが互いに垂直な２つの偏光に分離するために、ガラス透過後の光線は光束となると考えられる」としたが、直角に近い角度へ電気変位ベクトルが推移するような偏光の成分はごく僅かであると考えられるので、無限に近い本数の光線になるとは考えにくい。すなわち、図１のように、ガラス内のある層で、互いに垂直な電気変位ベクトルｄ１（ベクトル）、ｄ２（ベクトル）の偏光
【０１１５】
【数４７】
【０１１６】
があり、次の層で、電気変位ベクトルおよび偏光がそれぞれｄ’１（ベクトル）、ｄ’２（ベクトル）、および
【０１１７】
【数４８】
【０１１８】
に推移し、ｄ１（ベクトル）とｄ’１（ベクトル）のなす角度がΔθ≒０であるとすると、それらの光線の振幅は、
【０１１９】
【数４９】
【０１２０】
となり、大部分がもとの電気変位ベクトルの方向とほぼ同じ方向の電気変位ベクトルの光線に推移すると考えられる。
【０１２１】
レンズ内部の光線追跡においては、屈折率が連続的に変化するので、ｓ（ベクトル）があらかじめ与えられたとして微小距離進行後のｓ（ベクトル）の変化を逐次計算すれば良いが、レンズ表面の屈折の場合は、ｓ（ベクトル）が屈折前後で急激に変化するために、ｓ（ベクトル）とｄ（ベクトル）の双方を同時に求める必要がある。
【０１２２】
そのため、前記第３の手段においては、一般的な屈折光線の方向余弦を求める式から出発する。すなわち、ｃ（ベクトル）およびｓ（ベクトル）を、それぞれ入射光および屈折光の波面法線方向の単位ベクトルとし、境界面の法線ベクトルをｑ（ベクトル）＝（ｑｘ，ｑｙ，ｑｚ）とし、屈折媒質に入射する前の媒質（空気等）の屈折率をｎｃ、入射後の媒質の屈折率をｎとすると、
【０１２３】
【数５０】
【０１２４】
となり、
【０１２５】
【数５１】
【０１２６】
となる。［三宅和夫：幾何光学（共立出版，１９８５） ｐ．３０参照 ］
【０１２７】
【数５２】
【０１２８】
と置いて、
【０１２９】
【数５３】
【０１３０】
が成り立つように、ｓ（ベクトル）を求める。屈折後の媒質が非等方性の場合、
【０１３１】
【数５４】
【０１３２】
はｄ（ベクトル）の関数でもあるので、
【０１３３】
【数５５】
【０１３４】
とし、（１４）式の
【０１３５】
【数５６】
【０１３６】
と連立させて解く。この連立方程式を解くのには、Ｎｅｗｔｏｎ法を使用するのが好ましいが、他の公知の最適化手法を用いてもよい。Ｎｅｗｔｏｎ法を用いた用法について、以下に説明する。すなわち、ｄ（ベクトル）と同様にｓ（ベクトル）のさらに良い近似値を
【０１３７】
【数５７】
【０１３８】
とすると、
【０１３９】
【数５８】
【０１４０】
となるが、（６２）式、（６３）式の∂ｆｉ／∂ｓｉおよび∂ｇｉ／∂ｓｉ （ｉ，ｊ＝ｙ，ｚ） は、
【０１４１】
【数５９】
【０１４２】
であることに注意すれば、未知数が５個の連立方程式に還元できる。このとき、（６２）式、（６３）式の微係数は（１６）式に加えて、
【０１４３】
【数６０】
【０１４４】
が必要となる。なお、あまり精度が必要でない場合、波面法線楕円体を用いてｄ（ベクトル）をｓ（ベクトル）から求めることにより、単純に、
【０１４５】
【数６１】
【０１４６】
から、Δｓ（ベクトル）のみを未知数にすればよいが、（６４）式の関係により、未知数を２個に還元できる。この場合、（６９）式の代わりに、
【０１４７】
【数６２】
【０１４８】
とする必要がある。∂ｎ／∂ｓｊを含む項は、ｄ（ベクトル）をｓ（ベクトル）の関数とすることにより、新たに付加した項である。なお、（６２），（６３）式によりｓ（ベクトル）とｄ（ベクトル）を求めるときに、収束が遅い場合があるので、（７０）式でｓ（ベクトル）を求めて後、それらを初期値として、（６２），（６３）式でさらに高精度な計算をするとよい。このような計算手段は、前記第４の手段の１つの方法であることは明らかであり、前記第４の手段に含まれることは言うまでもない。
【０１４９】
【実施例】
以上のように述べた計算方法で、実施例として、外部応力によるガラス内の屈折率変動による収差変動を計算した。計算例として、図２のような球面収差のみ補正されたレンズ群に、波長２４８ｎｍ、半径５０ｍｍの光束を透過させた。このレンズ群の光学的諸元を表１に示す
【０１５０】
【表１】
【０１５１】
外半径５２．５ｍｍのレンズＧ１（ガラス）の子午線方向（ｙ軸とする）の両端に、１０Ｎ／ｃｍ２の外部応力を加えた。計算に用いたガラスのヤング率は、０．７２×１０５Ｎ／ｍｍ２，ポアソン比は０．１６４，光弾性定数は、Ｃ１＝６．４９５×１０−７ｍｍ２／Ｎ，Ｃ２＝４．１７６×１０−６ｍｍ２／Ｎである。
【０１５２】
この場合のガラス内の応力分布は、図３のようになっている。図３（Ａ）はｙ軸方向の応力分布で、図３（Ｂ）はｚ軸方向（光軸およびｙ軸の両方に垂直な方向）の応力分布である。等高線の間隔は１Ｎ／ｃｍ２である。
【０１５３】
このような外力を加えた場合において、レンズＧ１に半径５０ｍｍの光束を透過させ、ｇｒａｄｉｅｎｔ ｉｎｄｅｘ光線追跡をして波面収差を計算し、外力を加えない場合の波面収差との差をとり、波面収差変動を計算した。その結果、波面収差変動は図４のようになった。電気変位ベクトルがｙ軸方向の偏光を、前記のような応力分布のガラスに透過させて、像面上でｙ軸方向の偏光成分の波面収差変動を計算したのが図４（Ａ）であり、電気変位ベクトルがｚ軸方向の偏光を、同様な応力分布のガラスに透過させて、像面上でｚ軸方向の偏光成分の波面収差変動を計算したのが図４（Ｂ）である。図（Ａ）における等高線の間隔は０．００３２λ、図（Ｂ）における等高線の間隔は０．００４０λである。図４（Ａ）の差の最大値は、０．０３２λであり、図４（Ｂ）の差の最大値は、０．０４０λである。
【０１５４】
図３と図４とを比較すると、外力によるレンズの内部応力分布と偏光による波面収差変動の分布は、ほぼ対応するが、分布の形状がｙ軸とｚ軸とで互いに逆になっているのは、（３２）式において、Ｃ１よりもＣ２の絶対値が大きいことによるものと考えられる。
【０１５５】
なお、ｙ軸方向の偏光を同じ状態のガラスに透過させて、像面上でｚ軸方向の偏光成分を計算したものと、ｚ軸方向の偏光を同じ状態のガラスに透過させて、像面上でｙ軸方向の偏光成分を計算したものとは、図３（Ａ）および図３（Ｂ）の偏光よりも、かなり、強度が小さく無視できる程である。
【０１５６】
また、このような応力状態の場合、今回新たに発明した計算方法と、以前発明した波面法線楕円体による計算方法の精度を、（１１）式の残差ベクトルの絶対値の平均値を計算することにより、比較した。これを表２に示す。
【０１５７】
（表２）本発明による計算方法（Ａ）と先行論文による計算方法（Ｂ）による残差ベクトルの絶対値の平均値の比較
【０１５８】
【表２】
【０１５９】
先行論文による計算方法（波面法線楕円体による計算方法）（Ｂ）の場合、外部応力が１０倍になると、残差ベクトルの絶対値がほぼ１００倍になるのに対し、本発明の計算方法（Ａ）によれば、外部応力が１０倍になっても、残差ベクトルの絶対値はほとんど変わらず、計算機の倍精度計算の誤差の範囲内に収まっていることが分かる。なお、波面法線楕円体による計算方法の場合、残差ベクトルがゼロにならない原因は、この方法がＭａｘｗｅｌｌ方程式を満足していないことにあると考えるよりも、平方根などの関数を使用する際の計算誤差にあると考える方がよい。
【０１６０】
【発明の効果】
以上説明したように、本発明の計算方法により、硝材内に不均一な応力分布がある場合の光線追跡において、よりＭａｘｗｅｌｌ方程式に近い式から計算精度を検証することにより、従来よりも高精度に計算することが可能になり、計算の信頼性を増した状態で光線追跡を行うことができる。
【図面の簡単な説明】
【図１】偏光の振幅の微小な推移を示す図である。
【図２】本発明の実施例に用いた球面収差のみ補正されたレンズの光路図である。
【図３】本発明の実施例におけるガラス中の応力分布の例を示す図である。
【図４】本発明の実施例においてガラスに応力を与えた場合の、波面収差変動の計算結果の例を示す図である。[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to an optical system using polarized light, a projection optical system for manufacturing a semiconductor element, an optical system for observing astronomical objects, and the like. The present invention relates to a ray tracing method in an optical system in which polarization separation becomes a problem.
[0002]
In this specification, vectors are frequently used, but vector symbols cannot be used in sentences due to restrictions of electronic applications.
[Equation 3]
Is expressed as x (vector) in the text.
[0005]
[Prior art]
Conventionally, as for ray tracing when a light beam passes through a glass material having birefringence, see, for example, a method of obtaining the refractive index and wavefront radiation direction of a refracted light beam when a light beam enters a uniaxial crystal. C. Introduced by Simon in a document entitled "Ray tracing formula for monoaxial optical components". (Appl. Opt. (1983) pp 354-360).
[0006]
In addition, the principle of a method for obtaining the refractive index and wavefront normal vector of a refracted ray when the ray is incident on a biaxial crystal is generally described by Born Wolff: "Principle of Optics III" (Japanese translation: Tokai University Press) ), A method of obtaining the value using Fresnel's formula (see the same document, p. 980) is shown. (Refer to the same document, P999 to P1000)
[0007]
[Problems to be solved by the invention]
However, when there is a difference in the temperature rise distribution of the glass, or when an external force is applied to a part of the glass, distortion occurs in the glass and internal stress is generated. Since the light has a refractive index anisotropy, the light becomes polarized when the light passes through the glass. On the other hand, since the internal stress varies depending on the position in the glass and the refractive index varies according to the position, it is necessary to perform ray index tracing.
[0008]
In the gradient index ray tracing in a medium having an internal stress, the calculation of the refractive index needs to be performed with high accuracy so that calculation errors do not accumulate. C. It has been found that the method of calculating the refractive index as described in the literature of Simon and the like cannot guarantee the accuracy that can be cut off by the gradient index ray tracing. Accordingly, the present inventor has developed a method of calculating the refractive index by using a wavefront normal ellipsoid (Sumio Hashimoto: "Racing a ray in glass with internal stress distribution" Optics 26 (1997) 677: "Preceding article").
[0009]
In this method, the wavefront normal vector of the refracted light is obtained using the optical parameters and the photoelastic parameters.From this, `` thinking the wavefront normal ellipsoid, passing through the origin of the wavefront normal ellipsoid, and orthogonal to the wavefront normal vector The length of the principal axis of the ellipse is proportional to the reciprocal of the phase velocity, and the direction of the principal axis coincides with the vibration direction of the electric displacement vector. Based on the principle, an electric displacement vector is obtained, an electric vector is obtained by a relational expression between the electric displacement vector and the dielectric tensor, and a direction orthogonal to the two electric displacement vectors is defined as a direction of the magnetic vector. The pointing vector, that is, the traveling method of the light beam, is obtained from the cross product of the electric displacement vector and the magnetic vector obtained as described above.
[0010]
However, the accuracy of the calculation in this method has not been clarified. Also, it was not clear from the gradient index ray tracing whether the Maxwell equation was satisfied.
[0011]
The present invention has been made in view of such circumstances, and provides a method of performing ray tracing with higher accuracy than the method of the previous paper developed by the inventor and satisfying the Maxwell equation. The task is to
[0012]
[Means for Solving the Problems]
A first means for solving the above-mentioned problem is ray tracing in a medium having a refractive index anisotropy, and includes a residual vector.
(Equation 4)
(Where n is the refractive index, e (vector) is the electric field vector, d (vector) is the electric displacement vector, and s (vector) is the unit vector of the refracted light in the normal direction of the wavefront). A ray tracing method (claim 1) comprising a step of obtaining a vector d (vector).
[0015]
A second means for solving the above-mentioned problem is the first means, wherein the method for minimizing a residual vector is a Newton method in which a component of an electric displacement vector d (vector) is unknown. (Claim 2).
[0016]
A third means for solving the above-mentioned problem is a ray tracing method including a step of obtaining a wavefront normal vector of a light ray incident on a medium having a refractive index anisotropy and refracted, and a refractive index. Item 1. The residual vector and the residual vector according to item 1.
(Equation 5)
(However, c (vector) and s (vector) is the unit vector of the wave normal direction of the respective incident and refracted light, q (vector) is the normal vector of the boundary surface, n c before entering the refractive medium (N is the refractive index of the medium after incidence)
A ray tracing method (claim 3) characterized by having a step of obtaining d (vector) and s (vector) by minimizing both.
[0019]
A fourth means for solving the above-mentioned problem is the third means, wherein the method for minimizing the residual vector includes a component of the electric displacement vector d (vector) and a unit vector in the wavefront normal direction. The present invention is characterized by the Newton method in which the component of s (vector) is an unknown number (claim 4).
[0020]
The inventor first defined a residual vector for verifying whether the Maxwell equation was satisfied in a medium having refractive index anisotropy. Then, ray tracing is performed by obtaining the value of the component of the electric displacement vector d (vector) so as to minimize the residual vector. In order to obtain the component of d (vector), it is preferable to obtain the component by the Newton method using this as an unknown number.
[0021]
Similarly, in the case where light is incident from an isotropic medium such as air or a vacuum into a medium having a refractive index anisotropy, together with the residual vector, the ray vector of the incident light and the wavefront normal of the refracted light. The value of the component of the electric displacement vector d (vector) and the value of the component of the wavefront normal vector s (vector) that define the residual vector with respect to the vector s (vector) and minimize both of the residual vectors Was determined to perform ray tracing. In order to obtain the value of the component of the electric displacement vector d (vector) and the value of the component of the wavefront normal vector s (vector), it is preferable to obtain the values by the Newton method using these as unknowns.
[0022]
Hereinafter, these matters will be described in detail.
First, starting from the Maxwell equation, a method for verifying the accuracy of the method of calculating the unit vector s (vector) and the refractive index n in the wavefront normal direction in a medium having refractive index anisotropy will be described.
[0023]
The magnetic field, electric field and electric displacement vector of the ray are
[0024]
(Equation 6)
[0025]
And Here, k 0 is the wave number in vacuum, and φ is eikonal
[0026]
(Equation 7)
[0027]
It is. Here, r (vector) is a position vector of a light ray, s (vector) is a unit vector in a normal direction of a wavefront, and n is a refractive index expressed as n = c / v p by a light speed c and a phase speed v p. is there. The s (vector) can be obtained by a method as detailed in a previous paper, for example, but the first means and the second means do not particularly define a method for obtaining this. May be determined by other known methods. Substituting these into the Maxwell equation gives
[0028]
(Equation 8)
And by substituting equation (5) into equation (4),
[0029]
(Equation 9)
[0030]
Because
[0031]
(Equation 10)
[0032]
It becomes. In the case of an anisotropic medium, unlike the case of an isotropic medium, d (vector) and e (vector) are not vectors in the same direction.
[0033]
[Equation 11]
[0034]
It is. Therefore, from equation (6), the eikonal equation
(Equation 12)
[0036]
Cannot be derived. However, as long as n is a refractive index expressed as c / v p , from consideration of the geometrical optical wavefront,
[0037]
(Equation 13)
[0038]
Holds, squaring both sides of this equation yields the eikonal equation (11). Substituting equation (12) into equation (9) gives
[0039]
[Equation 14]
[0040]
It becomes. Therefore,
[0041]
[Equation 15]
[0042]
And the residual vector
(Equation 16)
[0044]
, The degree of calculation accuracy of the method of calculating d (vector) and e (vector) using the wavefront normal ellipsoid in the previous paper can be understood. The residual vector calculated from the wavefront normal ellipsoid increases as the internal stress increases.
[0045]
In the first means, unlike the previous papers, d (vector) is determined so that the residual vector defined by equation (14) is minimized. Since this method is based on Maxwell's equation unlike the prior art, d (vector) can be obtained more accurately.
[0046]
As described later, since n is a quadratic expression relating to the component of d (vector), the residual vector
[Equation 17]
[0048]
Is a non-linear expression for the component of d (vector). Therefore,
[0049]
(Equation 18)
[0050]
It is preferable to use the Newton method in order to obtain d (vector) such that is zero. [T. R. McCarla: Introduction to Numerical Calculation Methods for Computers (Science, 1984) p. Hereinafter, an example using the Newton method will be described.
That is, when i = 1 to 3 correspond to x, y, and z in the rectangular coordinate system,
[0051]
[Equation 19]
[0052]
A solution Δd j (j = 1, 3) in which unknowns satisfy three simultaneous equations is calculated, and a better approximation of d (vector) is calculated as follows.
(Equation 20)
[0054]
This calculation is repeated to obtain a highly accurate value of d (vector).
In equation (10), when s (vector) is given,
(Equation 21)
[0056]
The component f i, is differentiated with component d j of d (vector),
[0057]
(Equation 22)
[0058]
It becomes. Here, δ ij is Kronecker's δ.
[0059]
Hereinafter, a calculation method of the refractive index when there is an internal stress will be described. First, the relationship between the stress optical constants q 11 and q 12 and the photoelastic constants C 1 and C 2 is as follows.
[0060]
(Equation 23)
[0061]
It is. Then, the formula of a general refractive index ellipsoid in which the main axis of the ellipse and the coordinate axis do not match,
[0062]
[Equation 24]
[0063]
Then, the relational expression between β ij and q 11 and q 12 is represented by the refractive index n 0 before the stress is applied and the stress distribution σ ij (i, j = x, y, z at the intersection of the boundary surface and the incident ray). )
[0064]
(Equation 25)
[0065]
It becomes. In order for the expression (20) to hold, in the relationship between the stress optical constants q 11 , q 12 and q 44 ,
q 44 = 2 (q 11 -q 12) ... (21)
Must be established, which is proved to be established as follows. That is, the wavefront normal ellipsoid changed by the stress
(Equation 26)
[0067]
In the case where there is more than cubic symmetry,
[0068]
[Equation 27]
[0069]
[Borne Wolf: Principles of Optics III (Tokai University Press, 1976) p. 1026 to 1029], in the case of an isotropic medium such as glass, assuming that the refractive index in the absence of stress is n 0 ,
[0070]
[Equation 28]
[0071]
Than,
[0072]
(Equation 29)
[0073]
It becomes. Since the isotropic medium is invariant with respect to arbitrary coordinate transformation, if the xy coordinate system is rotated by θ to become the x′y ′ coordinate system, the principal stress σ x′x in the x′y ′ coordinate system ' , Σ y'y' are expressed by principal stresses σ xx , σ yy and shear stress σ xy in the xy coordinate system.
[0074]
[Equation 30]
[0075]
It is expressed as Assuming that the distance between the origin and the intersection of the x 'axis and the wavefront normal ellipsoid is r', the expression (22) is obtained by replacing x, y, z with x ', y', z '.
(X ′, y ′, z ′) = (r ′, 0, 0)
And, from equations (27) and (28),
[0076]
(Equation 31)
[0077]
On the other hand, in equation (22),
(X, y, z) = (rcos θ, rsin θ, 0)
By substituting into Equations (24) and (25),
[0078]
(Equation 32)
[0079]
By rotating the x-axis by θ to become the x′-axis, r ′ in equation (29) and r in equation (30) are the same, so that the coefficients of σ xy are equal,
q 44 = 2 (q 11 -q 12) ... (31)
And the equation (21) is established.
[0080]
Also, by considering the wavefront normal ellipsoid,
[0081]
[Equation 33]
[0082]
It becomes. Also,
[0083]
[Equation 34]
[0084]
It is. Proof of this is as follows. That is, a general refractive index ellipsoid whose principal axis does not coincide with the coordinate axis.
(Equation 35)
[0086]
Is coordinate transformation.
[Equation 36]
[0088]
When the orthogonal transformation is performed so that the main axis of the ellipsoid coincides with the coordinate axis, that is, the refractive index ellipsoid is expressed as follows.
(37)
[0090]
Let us consider the case of orthogonal transformation to the form Substituting equation (35) into equation (36) and replacing i and j,
[0091]
[Equation 38]
[0092]
Compared with equation (34),
[0093]
[Equation 39]
[0094]
The electric vector corresponding to (34) of the coordinate system, and component e i of electric displacement vector, d i, the component of the electric vector and the electric displacement vectors corresponding to equation (36) of the coordinate system, e ' i , d ′ i, and the orthogonal transformation matrix is C ij −1 = C ji ,
[0095]
(Equation 40)
[0096]
When the main axis and the coordinate axis of the refractive index ellipsoid coincide with each other as in the equation (36),
β ′ ii = 1 / n ′ i 2 (41)
And also
n ′ i 2 = ε ′ i μ (42)
d 'i = ε' i e 'i ... (43)
e ′ i = d ′ i / ε ′ i = μβ ′ ii d ′ i (44)
Substituting equations (44) and (39) into equation (40),
[0097]
(Equation 41)
[0098]
From equation (38),
[0099]
(Equation 42)
[0100]
Is obtained, which is equivalent to the equation (33).
[0101]
(19), (20), (32), (33) after the type or the like is differentiated by d j, seeking ∂f i / ∂d j are substituted into (16), (15) This problem can be solved by substituting into the equation.
[0102]
Then, the gradient index ray tracing, it is necessary values of ∂n / ∂x j, it is obtained by the following method. First, when substituting equation (32) into equation (14), the component becomes
[0103]
[Equation 43]
[0104]
By differentiating both sides of equation (13) with x j ,
∂f i / ∂x j = 0 ... (48)
From
[0105]
[Equation 44]
[0106]
From equation (32),
[0107]
[Equation 45]
[0108]
By substituting this into equation (49), a simultaneous linear equation relating to ∂d k / ∂x j is established, and the obtained value of ∂d k / ∂x j is substituted into equation (50). ∂n / ∂x j is determined by.
[0109]
As described above, after d (vector) is determined by the Newton method, e (vector) is determined by Expression (33). Since the vector (pointing vector) r (vector) in the traveling direction of the ray is the direction of the pointing vector, r (vector) is represented by a magnetic vector h (vector) perpendicular to both d (vector) and s (vector).
[0110]
[Equation 46]
[0111]
Required by In the above description, Newton's method was used to determine d (vector). However, in order to find d (vector) that minimizes the residual vector defined by equation (14), other d (vector) must be determined. A known optimization method may be used.
[0112]
In summary, in the first means, s (vector) is obtained by a known method, and d (vector) that minimizes the residual vector defined by the equation (14) is obtained therefrom ((14) E (vector) required in the expression is obtained by using expression (33) or other known relational expressions), d (vector) and h (vector) perpendicular to e (vector) are obtained, A pointing vector r (vector) indicating the traveling direction of the light ray is obtained as an outer product of e (vector) and h (vector).
[0113]
The change in s (vector) after a short distance travel is obtained by differential approximation of the differential equation of the light beam or the Runge-Kutta method [T. R. McCarla: Introduction to Numerical Calculation Methods for Computers (Science, 1984) p. 277].
[0114]
According to the preceding paper, "when glass is divided into layers, every time it is incident on each layer, the electric displacement vector separates into two polarized lights that are perpendicular to each other, so the light beam transmitted through the glass is considered to be a light beam." However, since it is considered that the component of the polarized light whose electric displacement vector shifts to an angle close to a right angle is extremely small, it is unlikely that the number of rays becomes almost infinite. That is, as shown in FIG. 1, polarization of electric displacement vectors d 1 (vector) and d 2 (vector) perpendicular to each other in a certain layer in the glass.
[Equation 47]
[0116]
In the next layer, the electrical displacement vector and polarization are d ′ 1 (vector), d ′ 2 (vector), and
[Equation 48]
[0118]
If the angle between d 1 (vector) and d ′ 1 (vector) is Δθ ≒ 0, the amplitude of those rays is
[0119]
[Equation 49]
[0120]
It can be considered that most of the light changes to a light beam of the electric displacement vector in a direction substantially the same as the direction of the original electric displacement vector.
[0121]
In the ray tracing inside the lens, since the refractive index changes continuously, it is sufficient to calculate the change in s (vector) after traveling a minute distance assuming that s (vector) is given in advance. In the case of refraction, since s (vector) changes rapidly before and after refraction, it is necessary to obtain both s (vector) and d (vector) at the same time.
[0122]
Therefore, the third means starts from a formula for calculating the direction cosine of a general refracted ray. That, c (vector) and s (vector), respectively the unit vectors of the wavefront normal direction of the incident light and the refracted light, the normal vectors of the boundary surface q (vector) = (q x, q y , q z ), and the refractive index n c of the previous medium entering the refractive medium (such as air) and the refractive index of the medium after the incident is n,
[0123]
[Equation 50]
[0124]
Becomes
[0125]
(Equation 51)
[0126]
It becomes. [Kazuo Miyake: Geometrical Optics (Kyoritsu Shuppan, 1985) p. 30]
[0127]
(Equation 52)
[0128]
And put
[0129]
(Equation 53)
[0130]
S (vector) is obtained so that If the medium after refraction is anisotropic,
[0131]
(Equation 54)
[0132]
Is also a function of d (vector), so
[0133]
[Equation 55]
[0134]
And (14)
[Equation 56]
[0136]
And solve it. To solve this simultaneous equation, it is preferable to use the Newton method, but other known optimization methods may be used. The usage using the Newton method will be described below. That is, similar to d (vector), a better approximation of s (vector) is
[Equation 57]
[0138]
Then
[0139]
[Equation 58]
[0140]
Where ∂f i / ∂s i and ∂g i / ∂s i (i, j = y, z) in equations (62) and (63) are
[0141]
[Equation 59]
[0142]
Note that can be reduced to a system of five unknowns. At this time, the differential coefficients of the equations (62) and (63) are added to the equation (16),
[0143]
[Equation 60]
[0144]
Is required. In the case where very little precision is required, by obtaining d (vector) from s (vector) using a wavefront normal ellipsoid,
[0145]
[Equation 61]
[0146]
Therefore, only Δs (vector) may be set to an unknown value, but the unknown value can be reduced to two according to the relationship of equation (64). In this case, instead of equation (69),
[0147]
(Equation 62)
[0148]
It is necessary to The term including ∂n / ∂s j is a term newly added by making d (vector) a function of s (vector). In addition, when s (vector) and d (vector) are obtained by the equations (62) and (63), convergence may be slow. Therefore, after obtaining the s (vector) by the equation (70), they are initialized. As a value, it is preferable to perform a calculation with higher accuracy by the equations (62) and (63). Obviously, such a calculation means is one method of the fourth means, and it goes without saying that the calculation means is included in the fourth means.
[0149]
【Example】
With the above-described calculation method, as an example, aberration variation due to refractive index variation in glass due to external stress was calculated. As a calculation example, a light beam having a wavelength of 248 nm and a radius of 50 mm was transmitted through a lens group in which only spherical aberration was corrected as shown in FIG. Table 1 shows the optical data of this lens group.
[Table 1]
[0151]
An external stress of 10 N / cm 2 was applied to both ends of the lens G1 (glass) having an outer radius of 52.5 mm in the meridian direction (referred to as the y-axis). The Young's modulus of the glass used for the calculation is 0.72 × 10 5 N / mm 2 , the Poisson's ratio is 0.164, and the photoelastic constant is C 1 = 6.495 × 10 −7 mm 2 / N, C 2 = 4.176 × 10 −6 mm 2 / N.
[0152]
The stress distribution in the glass in this case is as shown in FIG. FIG. 3A shows a stress distribution in the y-axis direction, and FIG. 3B shows a stress distribution in the z-axis direction (a direction perpendicular to both the optical axis and the y-axis). The interval between the contour lines is 1 N / cm 2 .
[0153]
When such an external force is applied, a light beam having a radius of 50 mm is transmitted through the lens G1, and a gradient index ray tracing is performed to calculate a wavefront aberration. The variance was calculated. As a result, the fluctuation of the wavefront aberration was as shown in FIG. FIG. 4A shows the calculation of the wavefront aberration variation of the polarization component in the y-axis direction on the image plane by transmitting the polarization having the electric displacement vector in the y-axis direction through the glass having the stress distribution as described above. FIG. 4 (B) shows the calculation of the wavefront aberration variation of the polarization component in the z-axis direction on the image plane by transmitting the polarized light whose electric displacement vector is in the z-axis direction through glass having the same stress distribution. The interval between the contour lines in FIG. (A) is 0.0032λ, and the interval between the contour lines in FIG. (B) is 0.0040λ. The maximum value of the difference in FIG. 4A is 0.032λ, and the maximum value of the difference in FIG. 4B is 0.040λ.
[0154]
3 and 4, when the internal stress distribution of the lens due to the external force and the distribution of the wavefront aberration fluctuation due to the polarized light substantially correspond, the shapes of the distributions are opposite to each other in the y-axis and the z-axis. , in (32), it is believed to be due to the absolute value of C 2 than C 1 is large.
[0155]
Note that the polarization in the y-axis direction is transmitted through glass in the same state, the polarization component in the z-axis direction is calculated on the image plane, and the polarization in the z-axis direction is transmitted through glass in the same state. The above calculation of the polarization component in the y-axis direction is considerably smaller in intensity than the polarization in FIGS. 3A and 3B and can be ignored.
[0156]
In the case of such a stress state, the accuracy of the calculation method newly invented this time and the accuracy of the calculation method using the wavefront normal ellipsoid previously invented are calculated by calculating the average value of the absolute value of the residual vector of the equation (11). By doing so, they were compared. This is shown in Table 2.
[0157]
(Table 2) Comparison of the average value of the absolute value of the residual vector by the calculation method (A) according to the present invention and the calculation method (B) according to the preceding paper
[Table 2]
[0159]
In the case of the calculation method according to the preceding paper (calculation method using a wavefront normal ellipsoid) (B), when the external stress increases by 10 times, the absolute value of the residual vector increases by almost 100 times. According to (A), even when the external stress is increased by a factor of 10, the absolute value of the residual vector hardly changes and is within the error range of the double precision calculation by the computer. In the case of the calculation method using the wavefront normal ellipsoid, the reason why the residual vector does not become zero is that the method does not satisfy the Maxwell equation. It is better to think that there is a calculation error.
[0160]
【The invention's effect】
As described above, according to the calculation method of the present invention, in ray tracing in a case where a non-uniform stress distribution is present in the glass material, by verifying the calculation accuracy from an expression closer to the Maxwell equation, higher accuracy than before can be achieved. Calculations can be performed, and ray tracing can be performed with increased calculation reliability.
[Brief description of the drawings]
FIG. 1 is a diagram showing a small change in the amplitude of polarized light.
FIG. 2 is an optical path diagram of a lens used in an embodiment of the present invention, in which only spherical aberration is corrected.
FIG. 3 is a diagram showing an example of a stress distribution in glass in an example of the present invention.
FIG. 4 is a diagram showing an example of a calculation result of a wavefront aberration variation when a stress is applied to glass in an example of the present invention.
