１． はじめに
フィルムカメラの全盛期には、そこそこの写真マニアだったこともあり、レンズにもそれなりの興味はありました。しかし、持っている知識はせいぜい高校物理までのレベルです。そこで、自由時間をとり放題の今、もう少し深く学んでみたくなりました。手始めに、レンズの光路シミュレーションなどをやり始めてみましたが、世間では、まずは「近軸近似」なるものをマスターするのが先のようです。
しかし、「所詮 近似、大した役には立つまい」という先入観が強く全くヤル気が湧いてきません。高校物理の「薄いレンズ」で感じた安直さが今でも尾を引いています。ところが、記事のタイトルにもあるように、これは認識不足による大きな誤解でした。高級な複合レンズの焦点距離でも、ピッタリと正確に計算できる非常に役に立つ理論です。
参考になる文書やウェブページは沢山ありますが、説明内容に省略が多いうえに、さまざまな流儀が入り乱れて、認知機能が怪しくなってきた頭には優しさが足りません。そこで、自分にも理解できるように、これらを再構築して分かりやすい説明記事を書いてみました。
曖昧さを減らそうとしたため長文になってしまいました。要点だけが必要なときには下記の節さえ読んでいただければ十分と思います。
２.１ １つの境界面での屈折現象
２.３ １つの境界面での近軸基本式
３.１ 複数の境界面による屈折現象
３.２ u と h を使う計算方法
３.５ 光線行列解析 h,u 方式
４.２ 光路シミュレーションとの比較
２． １つの境界面での式の導出
２.１ １つの境界面での屈折現象
図１は、球面レンズの一つの境界面における光の屈折現象を示しています。この図の主要部分は、高校物理にも出てくる屈折の法則（スネルの法則）です。教科書では境界面が水平面として描かれ、その垂線を基準として入射角、出射角（屈折角）の関係で表現されています。しかし、この水平面と垂線を、レンズの接平面と法線とに置き換えて考えれば、スネルの法則はレンズの表面上でもそのまま成立します。なお、入射側と出射側で値が異なる変数については、出射側のシンボルに「$'$」を付けることで区別しています。
注意： 本記事では、とくに混同の恐れがない場合には、その現象の主体が「レンズ」なのか「境界面」なのかを明示せずに、「入射」「出射」を使うことがあります。どちらが主体かは、状況に応じて適切に解釈してください。
座標の原点は光軸と境界面との交点に置きますが、各座標軸の割り付け方法については文書によってマチマチです。ここでは、レンズの解析で最も一般的とされる方法に合わせて、光軸（横軸）を $\rm z$ 軸、縦軸を $\rm y$ 軸としています。$\rm x$ 軸は紙面の手前から奥に進む向きになりますが本記事では使用しません。（座標の採り方は便宜的なものであり、最終結果の計算式の形には関係しないため、あまり気にする必要はありません）
図１ 境界面での屈折の基本モード（流儀１）
図１が、レンズの入射側の面を示すものであれば左が空気、右がガラス、出射側の面であれば、その逆になります。左右の各領域は、それぞれの屈折率の材質で無限遠方まで満たされているものとします。図では、左に凸になった境界面の例を示していますが、凹/凸レンズの別や入射側/出射側の別により 右に凸になることもあります。これらの条件により多くの組み合わせが考えられ、その様子は図２に示すとおりです。
図２ 考えられる屈折の全モード
図２の上半分は左に凸、下半分は右に凸の境界面のモードです。また、左半分は 屈折率が「小 $\Rightarrow$ 大」、右半分は「大 $\Rightarrow$ 小」に変わる境界面のモードを示しています。全部で16種類のモードがありますが、〇で囲んだ同一番号のモードの a と b では数式が同じになるので、数式の違いで分類すれば全10モードです。
実は、図２の各図を、光軸を中心にして線対称に折り返すと、モード数はさらに倍増します。ところが、折り返しによる変形では、角度の全変数 $i$, $i'$, $u$, $u'$, $\varphi$ と、位置の変数 $h$ の符号が一斉に反転するため、数式の形に変化はありません。つまり、この変形では、$\cos$ を含む式には全く変化がなく、$\sin$, $\tan$ を含む式は両辺に「${-}1$」が掛かるだけなので、数式が表す意味は折り返し前と同じです。それゆえ、これらのモードは図２には含めていません。
さて、ここでは、全てのモードを統括する代表例として採用するものを「基本モード」、その他のモードを「異種モード」と呼ぶことにします。図１は、図２の ①a を基本モードとして選んだものです（流儀１）。しかし、場合によっては、図３のように、 ③a が基本モードに選ばれることもあります（流儀２）。また、稀には、④が基本モードに使われることもあります（流儀４）。さらに、説明の都合上、➄を基本モードにしたマイナーな流儀３も仲間に入れておきます。なお、流儀を区別するこれらの番号は、私が勝手に決めたもので一般的に通用するものではありません。
図３ 境界面での屈折の基本モード（流儀２）
結論から言えば、図２の中のどれを基本モードに採用しても、他のモードの動作まで含めて議論することができます。しかし、基本モードが変わると、結果として得られる数式の一部の形が異なってきます。そのため、利用者の混乱を招く原因にもなっていて少々困りものです。この記事では、まずは、図１を採用した 流儀１ の説明を優先し、最終的には他の流儀による結果にも言及することにします。
２.２ 基本モードだけで論じることの妥当性
世間では、直交座標の全ての象限で成り立つべき理論の証明を、第一象限の図だけを描いて済ませることが多いように感じます。例えば、高校数学の教科書では「$\cdots\cdots$ は、一般角に対しても成り立つことが知られている」と付記するだけで済ませ、具体的な確認方法にまでは言及していません。しかし、教える側、教えられる側の双方にとって、手数が省けることの方が歓迎されているためか、あまり苦情は聞いたことがありません。
本記事においても、多数のモード間の問題は、同様の論法で軽く乗り切りたいというのが本音ではあります。しかし、象限の場合よりもさらに複雑な様相を呈しているため、曖昧さが拭えずモヤモヤが残ってしまいそうです。そこで、基本モード図の中に書き込まれた数式だけで、異種モードの動作まで含めて議論できることを、もう少し掘り下げて確認しておこうと思います。
さて、図２の各モード図の中の数式の２行目以下の部分は、下記を前提として、三角形の幾何学的な性質を利用して導出したものです。
-
長さの符号
位置を表す変数 $\boldsymbol{r}$, $\boldsymbol{s}$, $\boldsymbol{s'}$, $\boldsymbol{\delta s}$, $\boldsymbol{h}$ の値は、横軸・縦軸の座標値をそのまま使用します。したがって、原点よりも右側/上側にある変数は正、左側/下側にある変数は負です。（余談： この結果、必然的に、境界面が左に凸の場合には $r$ は正、右に凸の場合には負となります。）
しかし、図中に含まれる三角形の辺のうち、立式に関係している辺の長さはすべて正の値となるように算段します。例えば、横軸を利用する場合には、軸の右側にある大きな数値から左側にある小さな数値を引くことで、辺の長さとして使う正の値を作ります。縦軸の場合もこれに準じます。 -
角度の符号
図中に表示されている角度 $\boldsymbol{i}$, $\boldsymbol{i'}$, $\boldsymbol{u}$, $\boldsymbol{u'}$, $\boldsymbol{\varphi}$ は、すべて正の値とみなしています。こうしておかないと、図と数式の直感的な対応づけができなくなるからです。
一方、図２の全モードについての包括的な議論をするためには、全図に共通する「数値の表現規則」が必要になります。上記の条件のうち、座標上の位置の表現方法については共通性があるので問題はありません。しかし、角度の符号の表現方法は各図にローカルなものであり、そのままでは使うことができません。そこで、角度の符号の共通的な扱いについては、別途、下記のように定めておくことにします。なお、ここで言う「角度」は鋭角側の角度を指しています。（光学関係では「時計回りが正」としている文書も多いようですが、ここでは数学的に順当な「反時計回りが正」のほうを採用しています。ここでどちらを採用したとしても、最終結果として得られる焦点距離などの計算式の形は変わらないので、それほど気にする必要はありません）
- $u$, $u'$, $\varphi$ については、「光軸」を基準にして、「光路」や「境界面の法線」が反時計回り側の角度にあるとき「global正」、時計回り側の角度にあるとき「global負」とします。
- $i$, $i'$ については、「法線」を基準にして、「光路」が反時計回り側の角度にあるとき「global正」、時計回り側の角度にあるとき「global負」とします。（余談： これにより、$i$ と $i'$ は必ず同符号になります。）
図１や図３を参照し、この対応関係をまとめると表１のようになります。表中に示す単なる「正/負」は、「各基本モードの都合（図示の角度がすべて正）に合わせてローカルに決めた正/負」を意味しています。
表１ 数式中の角度変数の符号の扱い
| 角度 変数 |
図１が基本モードのとき | 図３が基本モードのとき | ||
|---|---|---|---|---|
| global正 | global負 | global正 | global負 | |
| $i$ | 正 | 負 | 正 | 負 |
| $i'$ | 正 | 負 | 正 | 負 |
| $u$ | 負 | 正 | 正 | 負 |
| $u'$ | 負 | 正 | 正 | 負 |
| $\varphi$ | 負 | 正 | 負 | 正 |
結論から言うと、表１に示す符号の読み替えさえ行えば、基本モードの数式を使うだけで異種モードの動作をも代表できるようになります。これを具体的に確認してみます。
例えば、図１を基本モードに採用して、図２の ⑩ を異種モードとして議論するときには、⑩ 図から、$i$, $i'$, $u$, $u'$, $\varphi$ の角度の状態が、順に、global負,global負,global正,global負,global正 となっていることを読み取ります。これを、表１の「図１が基本モードのとき」の欄と照合させると、各変数に対応するローカルの符号は 負,負,負,正,負 となっていることが分かります。そこで、その符号に合わせて、⑩ 図中の数式の $i$, $i'$, $u$, $u'$, $\varphi$ を ${-}i$, ${-}i'$, ${-}u$, $u'$, ${-}\varphi$ と置き換えたうえで整理すれば、図１の中の数式と同じ形にすることができます（単調で冗長な変形なので詳細は省略します）。逆に、図１の数式の変数を上記と同様に置き換えても、⑩ の数式と同じ形にすることができます。
結局、global正/負を介することで、共通的な角度の整合性を保てていることが分かります。別の言い方をすれば、図１のモードでのローカル正/負の定義を、異種モードにもそのまま受け入れさせて式を書き替えれば、すべて、図１の中の式と同一になるということです。
以上では、一例だけを採り上げて説明しましたが、全 90 通りの組み合わせについて、どれも問題なく適合することが確認できています。これを手作業でやるのは面倒です。最後に確認プログラムを添付しておきます。（もっとスマートに証明したいのですが、シラミ潰しの方法しか思い浮かびません。）
それにしても、表１のような面倒なものが出てくると、読み進む気が失せてしまいそうです。しかし、心配は無用です。この後に続く式の変形の過程で $i$, $i'$, $\varphi$ は消去されてしまうため、これ以降で気に留めておくべきは、残りの $\boldsymbol{u}$, $\boldsymbol{u'}$ の符号だけになります。この結果、 「global正/負」＝「右上がり/右下がり」と 分かりやすい表現に言い替えることができるようにもなります。
２.３ １つの境界面での近軸基本式
さて、細部に拘ったため、議論が少々脇道に逸れてしまいました。ここから本題に戻ります。
この解析での最終課題は、複雑なレンズ群を通った光が、光軸上のどこに到達するかを求めることです。そのためには、それぞれの境界面において、光軸上の $z{=}s$ の点に向かうはずだった光が、$z{=}s'$ の点に向きを変える現象について定量的に把握できている必要があります。これを直接的に表せるのは、$s$ から $s'$ への関係式、あるいは、少し間接的にはなりますが、 $u$ から $u'$ への関係式です。まずは、これらを求めることから始めます。
図１から、次のような関係が成り立っていることが分かります。（ $\varphi{=}\cdots$ の式は使いません。$\varphi$ が小さいことさえ分かっていれば十分です。）
\begin{align}
&n\sin i=n'\sin i'\tag{1}\\
&r\sin i=(s-r)\sin u\tag{2}\\
&r\sin i'=(s'-r)\sin u'\tag{3}\\
&\delta s=r(1-\cos\varphi)\tag{4}\\
&h=(s-\delta s)\tan u=(s'-\delta s)\tan u'\tag{5}\\
\end{align}
(2), (3)式から $\sin i$ と $\sin i'$ の値を求めて、それらを(1)式に代入すると次式が得られます。
\begin{align}
n\frac{s-r}{r}\sin u=n'\frac{s'-r}{r}\sin u'\tag{6}
\end{align}
これによって、スネルの法則で重要な役割を担っていた変数 $i$ と $i'$ が都合よく消え、$s$, $s'$, $u$, $u'$ の関係が端的に表せるようになりました。しかし、$s$ と $s'$、$u$ と $u'$ の２つの関係式に分離するためには、式がもう一つ必要です。これには(5)式が使えそうなのですが、(5), (6)式にはともに三角関数が含まれており、そのままでは容易には解けそうもありません。しかし、屈折現象を光軸の近傍のごく狭い範囲（近軸）だけで考えれば、数式をもっと簡略化して扱うことができます。
(4), (5), (6)式の中で使われている三角関数は、角度の単位をラジアンで表したとき、マクローリン展開によって下記のように多項式として表現することができます。ここで、$\theta$ は $u$, $u'$, $\varphi$ と置き換えることができるシンポルです。
\begin{align}
&\sin\theta=\theta\color{#aaa}{-\frac{1}{6}\theta^3+\frac{1}{120}\theta^5+\cdots}\\
&\cos\theta=1\color{#aaa}{-\frac{1}{2}\theta^2+\frac{1}{24}\theta^4-\cdots}\\
&\tan\theta=\theta\color{#aaa}{+\frac{1}{3}\theta^3+\frac{2}{15}\theta^5+\cdots}\\
\end{align}
$\theta$ が非常に小さな値のときには、次数の高い $\theta^2$ 以上の項はさらに小さな値になるので無視することができます。$u$, $u'$, $\varphi$ のいずれも、近軸では非常に小さな値になるので、灰色の項を無視したうえ、(6), (4), (5)式の該当部分に適用すると、下記の(7)～(9)式に変形できます。厳密には、「$=$」 は 「$\simeq$」 とすべきですが、近軸だけに限定していることは既に宣言済みなので、そこまでは拘らないことにします。（ここからは、麻生さんもご満悦の脱三角関数の世界😁に入ります。これ以降も、$\boldsymbol{\tan\theta{=}\theta}$ の関係は随所で断りなく利用するので留意ください）。
\begin{align}
&n\frac{(s-r)u}{r}=n'\frac{(s'-r)u'}{r}\tag{7}\\[1ex]
&\delta s=r(1-1)=0\tag{8}\\[1ex]
&h=su=s'u'\tag{9}\\
\end{align}
(7)式の両辺に $r$ を掛けて、移項、整理すると、
\begin{align}
nsu-nru=n's'u'-n'ru'
\end{align}
\begin{align}
r(n'u'-nu)=n's'u'-nsu
\end{align}
これに(9)式を適用すると、さらに次のように変形していくことができます。
\begin{align}
r(n'u'-nu)=h(n'-n)
\end{align}
\begin{align}
n'u'=nu+\frac{h(n'-n)}{r}\tag{10}
\end{align}
\begin{align}
\bbox[white, 5pt, border: 1.5px solid black]{u'=\frac{n}{n'}u+\frac{h}{r}(1-\frac{n}{n'})}\tag{11}
\end{align}
これが、$u$ から $u'$ を求める関係式です。
さらに、(11)式に、(9)式から導かれる $u=h/s$ 、$u'=h/s'$ を代入すると、次のように変形できます。
\begin{align}
\frac{1}{s'}=\frac{n}{n'}{\cdot}\frac{1}{s}+\frac{1}{r}\left(1-\frac{n}{n'}\right)\tag{12}
\end{align}
これは、$s$ から $s'$ を求める関係式です。逆数としてではなく、$s'$ そのものを求めるのは次式です。
\begin{align}
\bbox[white, 5pt, border: 1.5px solid black]{s'=\frac{1}{\displaystyle\frac{n}{n'}{\cdot}\frac{1}{s}+\frac{1}{r}\left(1-\frac{n}{n'}\right)}}\tag{13}
\end{align}
３． 焦点距離と主点の計算
３.１ 複数の境界面による屈折現象
以上で、境界面が１面だけのときの計算式が決まりました。ここからは、複数の境界面を考えます。１枚のレンズであれば、入射面と出射面とで２面になります。複合レンズでは境界面はさらに増えますが、レンズの枚数の２倍になるとは限りません。例えば、２枚を張り合わせたレンズでは境界面は３面です。
図４ 複数の境界面での屈折
図４は、境界面が２面の場合の図です。ここからは、レンズ面は球面ではなく垂直な平面として図示しています。(8)式によって図１の $\delta s$ が $0$ ですから、これで何の問題もありません。レンズ面が球面であるという情報は、$r$ として既に(11), (13)式に含まれていますから、ただの平面として図示されていても、球面としての働きを無視した訳ではありません。
各変数には境界面番号を付記し、どの面についての情報かを区別できるようにします。図４では、第１,２面を描いているものとして、1, 2 の添字を付けています。この図からは、添字に 2 が付いた変数 $s$, $s'$ の原点は、第２面の位置 $z{=}z_2$ に移っていることも読み取れます。なお、これらの面を、任意の隣接した２つの境界面と捉えるのであれば、$1$ を $j$、$2$ を $j{+}1$ と読み替えることができます。
ここでは、第１面と第２面との間隔 $d$ の添字は、第２面が持つ属性と解釈して $d_2$ としています。屈折率が $n_2$ の区間の長さですから $d_2$ が自然と思います（しかし、文書によっては $d_1$ や $d'_1$ としていることもありますので混同しないようご注意ください）。また、どの添字の場合でも、$d$ の符号はすべて正として扱います。すなわち、$d_{j{+}1}{=}z_{j{+}1}{-}z_j$ あるいは、 $d_j{=}z_j{-}z_{j{-}1}$ です。
今までは、出射側を表す変数には「$'$」をつけて区別していました。しかし、図からも分かるように、これらのうち $n'$, $u'$ は、そのまま次の境界面の入射側の変数にもなっています。したがって、次のように、「$'$」なしでの表現も可能になります。
\begin{align}
n_{j+1}=n_j'\tag{14}\\
u_{j+1}=u_j'\tag{15}\\
\end{align}
ただし、横軸の原点は境界面ごとに異なるので、$s'$ については、そのままでは次の境界面の変数としては使えず、下記の変換が必要です。
\begin{align}
\bbox[white, 5pt, border: 1.5px solid black]{
s_{j{+}1}=s'_j-d_{j{+}1}\tag{16}
}
\end{align}
また、隣接する境界面間の $h$ の変換式は、図より次のように記述することができます。（$u'_j$ の部分は、本来は $\tan u'_j$ ですが、ここでも近軸での近似式を使っています。以降、この種の注記は省略します。）
\begin{align}
h_{j{+}1}=h_j-d_{j{+}1}u'_j
\end{align}
同じ値のものを異なった変数で表すと混同しやすいので、以降、$s$ 以外は「$'$」が付かない表現の方だけを採用することにします。そこで、上式に(15)式を代入すると、
\begin{align}
\bbox[white, 5pt, border: 1.5px solid black]{
h_{j{+}1}=h_j-d_{j{+}1}u_{j{+}1}\tag{17}
}\end{align}
最後に、(9)式に添字を追加したうえ(15)式を代入し、あとで使うことになる次の関係式を導いておきます。
\begin{align}
h_j=s_ju_j=s'_ju'_j=s'_ju_{j{+}1}\tag{18}
\end{align}
以上で準備は整いました。ここからは、１面だけで導出した(11)式 あるいは(13)式を使用し、多面を連続的に処理して、複合レンズの焦点距離と主点の位置を求める方法についてまとめます。なお、境界面の全面数は $k$ とします。
３.２ u と h を使う計算方法
焦点と主点の定義は図５のとおりです。すなわち、光軸と平行な「レンズ群への入射光」を想定したとき、「レンズ群からの出射光」（凹レンズのときはその延長線）が光軸と交わる点 $z{=}z_{\rm f}$ が焦点（Focal Point）、入射光の延長線と出射光の延長線が交わる点 $z{=}z_{\rm p}$ が主点（Principal Point）です。これら焦点と主点との間の長さが焦点距離（Focal Length）$f$ になります。なお、凸レンズの場合には $f$ の数値の符号は正、凹レンズの場合には負で表現することになっています。（主点を表す記号として、多くの文書でドイツ語 Hauputpunkt の「H」を使っています。しかし本記事では、少数派ながらWikipediaでも使っている英語流の「p」を使います。）
図５ 焦点距離と主点（後側、像側）
実は、焦点と主点は、各レンズにそれぞれ２つずつあります。レンズから見て、被写体が存在すると想定される側を「前側（物体側）」、フィルムや撮像素子が置かれる側を「後側（像側）」と呼びます。そこで、前方からの平行光線が後側に作る焦点や主点を「後側焦点」,「後側主点」、逆に、後方からの平行光線が前側に作るものを「前側焦点」,「前側主点」として区別しています。ここでは、まずは後側の焦点と主点について考えます。
さて、ここで(11)式にも境界面番号 $j$ を書き足し、(14), (15)式を代入したうえ、(17)式と組み合わせると、レンズの特性計算のための次の重要な基本式が得られます。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
&u_{j{+}1}=\frac{n_j}{n_{j{+}1}}u_j+\frac{h_j}{r_j}(1-\frac{n_j}{n_{j{+}1}})\\[1ex]
&h_{j{+}1}=h_j-d_{j{+}1}u_{j{+}1}
\end{align}\tag{19}
}
まずは、$\boldsymbol{j{=}1}$ として、(19)式の二つの式を上から順番に計算します。初期値の $u_1$ は、軸に平行に入射する光線を想定して $\boldsymbol{u_1{=}0}$ と置きます。また、$\boldsymbol{h_1}$ には任意の値を設定します。その他の値は、各境界面の属性として既知のはずですから、それらをそのまま代入します。すると、結果として $u_2$, $h_2$ が得られます。
次に、$j$ の値を一つ増やして $j{=}2$ とおき、先に得られた $u_2$, $h_2$ を入力として $u_3$, $h_3$ を計算します。あとは同様に $\boldsymbol{\boldsymbol{j}}$ の値を一つずつ増やしながら $\boldsymbol{j{=}k}$ になるまで繰り返します。ただし、最後の回の $\boldsymbol{j{=}k}$ の計算は１行目の式だけで終了します。この結果、$\boldsymbol{u_{k{+}1}}$ と $\boldsymbol{h_k}$ が得られます。
一方、図５より、次の関係があることが分かります。
\begin{align}
&h_k=(z_{\rm f}-z_k)u_{k{+}1}\tag{20}\\[1ex]
&h_1=h_{\rm p}=(z_{\rm f}-z_{\rm p})u_{k{+}1}\tag{21}
\end{align}
(20), (21)式を変形すると、
\begin{align}
&z_{\rm f}=\frac{h_k}{u_{k{+}1}}+z_k\tag{22}\\
&z_{\rm p}=-\frac{h_1}{u_{k{+}1}}+z_{\rm f}=\frac{h_k-h_1}{u_{k{+}1}}+z_k\tag{23}
\end{align}
ゆえに、焦点距離 $f$（ ${=}z_{\rm f}{-}z_{\rm p}$ ）は、
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f=\frac{h_1}{u_{k{+}1}}
\end{align}
}\tag{24}
また、最後面から後側主点までの距離 $s_{\rm p}$（ ${=}z_{\rm p}{-}z_k $ ）は次のようになります。主点が最後面よりも 前側/後側 にあるとき、$s_{\rm p}$ の符号はそれぞれ 負/正 になります。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm p}=\frac{h_k-h_1}{u_{k{+}1}}
\end{align}
}\tag{25}
なお、最前面を基準にして後側主点を表現するときには、$s_{\rm p}$ にレンズの全長 $\sum_{j{=}2}^{k}d_j$ を加えます（以下、全方式について同様）。
以上から分かるように、(19)式から得られる値のうち、焦点距離と主点の計算に必要となるのは、最終値の $u_{k{+}1}$ と $h_k$ の２つだけと少ないので便利です。しかし、最終的にはキャンセルされてしまうにも拘わらず、変数 $h_1$ に何らかの数値を仮定しないと計算が進められないのは厄介です。
レンズの前後とも空気など同一の屈折率の場合には、前側の焦点距離は後側のそれと同じなので再計算する必要はありません。しかし、前側の主点の位置は後側とは異なります。前側の主点位置や、レンズ前後の屈折率が異なる環境下での前側焦点距離を求めるためには、レンズ群全体を反転させたうえ、上記の方法に準じて再計算する必要があります。このとき、レンズデータの並びは下記の方法で変更します。
- $r$ は添字を置き換えて、符号を反転する。$r_j\ \Rightarrow\ {-}r_{k-j+1}$
- $n$, $d$ は添字を置き換える。$n_j\ \Rightarrow\ n_{k-j+2}$ 、$d_j\ \Rightarrow\ d_{k-j+2}$
しかし、この方法はミスを誘発しやすいため、後述の「光線行列解析」（３.５節、３.６節）を使うほうが安全で便利です。
３.３ s を使う計算方法
(13)式に境界面番号を適用して(14)式を代入したうえ、(16)式とまとめると次のようになります。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
&s'_j=\frac{1}{\displaystyle\frac{n_j}{n_{j{+}1}}{\cdot}\frac{1}{s_j}+\frac{1}{r_j}\left(1-\frac{n_j}{n_{j{+}1}}\right)}\\[1ex]
&s_{j{+}1}=s'_j-d_{j{+}1}
\end{align}\tag{26}
}
(26)式も、$j$ を $1$ から境界面の全数 $k$ まで一つずつ増やしながら、１行目、２行目の順に計算していきます。ただし、(19)式と同様、$k$ 回目については１行目の式だけで計算を終えます。初期値の $s_1$ は、軸に平行な光線を想定して$s_1{=}\infty$ と置きます。
焦点距離 $f$ は、各境界面で得られた $s$, $s'$ の値から次の式で求められます（導出過程については後述します）。ただし、レンズ系の途中に光軸と平行な光路ができる場合には、この部分で $s$ や $s'$ の値が $\infty$ になることがありそうなので、計算エラーにならないような注意が必要です。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f=\frac{s'_k{\cdot} s'_{k{-}1}{\cdot} s'_{k{-}2}\cdots\cdots s'_2 {\cdot} s'_1}{s_k{\cdot} s_{k{-}1}{\cdot} s_{k{-}2}\cdots s_3 {\cdot} s_2}=\frac{\prod_{j=1}^{k}s'_j}{\prod_{j=2}^{k}s_j}
\end{align}
}\tag{27}
図５より、最後面から後側主点までの距離 $s_{\rm p}$（ $ {=}z_{\rm p}{-}z_k$ ）は、
s_{\rm p}=(z_{\rm p}-z_{\rm f})+(z_{\rm f}-z_k)=-f+s'_k
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm p}=s'_k-f
\end{align}
}\tag{28}
「$h$ と $u$ を使う計算方法」の場合とは異なり、計算開始時に $h_1$ のような暫定的な数値を仮定する必要はありません。しかし、焦点距離と主点の計算には、(26)式の各計算ステップで得られる全ての $s$, $s'$ の値が必要になります。
さて、先に結果だけを示した「焦点距離を計算する(27)式」は、次のようにして求められます。
(18)式より、隣接する境界面の各 $u$ の値の間には次の関係があることが分かります。
\begin{align}
u_{j{+}1}=\frac{s_j}{s'_j}u_j
\end{align}
$j$ は、$k$ 以下の任意の正整数をとれる添字なので、隣接面どうしの関係であれば他の多くの表現も可能です。
\begin{align}
&u_j=\frac{s_{j{-}1}}{s'_{j{-}1}}u_{j{-}1}\hspace{2em}
u_{j{-}1}=\frac{s_{j{-}2}}{s'_{j{-}2}}u_{j{-}2}\hspace{2em}
u_{j{-}2}=\frac{s_{j{-}3}}{s'_{j{-}3}}u_{j{-}3}\\[1ex]
&\cdots\cdots\cdots\hspace{2em}
u_4=\frac{s_3}{s'_3}u_3\hspace{2em}
u_3=\frac{s_2}{s'_2}u_2\hspace{2em}
u_2=\frac{s_1}{s'_1}u_1\hspace{2em}
\end{align}
これらの式の右辺にある $u$ に、後続の式の左辺の $u$ を順次代入していき、(18)式より $h_1{=}s_1u_1$ であることを考慮すると次の式が得られます。
\begin{align}
u_{j{+}1}&=\frac{s_j}{s'_j}{\cdot}\frac{s_{j{-}1}}{s'_{j{-}1}}{\cdot}\frac{s_{j{-}2}}{s'_{j{-}2}}\cdots\cdots\frac{s_3}{s'_3}{\cdot}\frac{s_2}{s'_2}{\cdot}\frac{s_1}{s'_1}u_1\\[1.5ex]
&=\frac{s_j{\cdot} s_{j{-}1}{\cdot} s_{j{-}2}\cdots s_3 {\cdot} s_2}{s'_j{\cdot} s'_{j{-}1}{\cdot} s'_{j{-}2}\cdots\cdots s'_3 {\cdot} s'_2 {\cdot} s'_1}h_1
\end{align}
(24)式の $f=h_1/u_{k{+}1}$ の分母に、上式の $j$ を $k$ と置き換えたものを代入すれば、(27)式が得られます。
３.４ nu と h を使う計算方法
前述の「$u$ と $h$ を使う方法」で提示した(19)式は、(11)式をもとにして作られています。(11)式は、求めるべき値 $u'$ が左辺に単独で明示されており、素人目には自然で好ましい形です。しかし、世間ではこの形は主流ではなく、$n$ と $u$ の積 $nu$ を一纏まりの変数として扱う方が一般的なようです。この場合、(11)式の１ステップ前の(10)式を採用することになり、(19)式は次のように書き換えられます。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
&(nu)_{j{+}1}=(nu)_j+\frac{h_j}{r_j}(n_{j{+}1}-n_j)\\[1ex]
&h_{j{+}1}=h_j-d_{j{+}1}\frac{(nu)_{j{+}1}}{n_{j{+}1}}
\end{align}\tag{29}
}
計算の手順は(19)式に準じます。すなわち、まずは $j{=}1$ とし、$u_1{=}0$（ゆえに $(nu)_1{=}0$）、任意の $h_1$ を初期値として計算を開始し、$j{=}k$ 回の１行目まで計算して $u_{k{+}1}$（ ${=}(nu)_{k{+}1}/n_{k{+}1}$ ）と $h_k$ を求めます。焦点距離や主点の計算には(24), (25)式がそのまま使えます。
(19)式をこの形式に置き換えることにより、１行目の計算での乗除算は各１回ずつ減りますが、２行目の計算では $u_{j{+}1}$ が必要になるため、$u_{j{+}1}{=}(nu)_{j{+}1}/n_{j{+}1}$ の除算が１回だけ増えます。結局、１つの $j$ について１回だけ乗算が節約できることになります。しかし、手回しのタイガー計算機の時代ならともかく、このコンピュータの時代に、これがどれだけのメリットになるのでしょうか。ただし、素人が気づけないだけで、これ以外にも何か良いことがあるのかもしれません。
３.５ 光線行列解析 h,u 方式
「$s$ を使う計算方法」の(26)式には、$s$ が逆数の形で入っていて複雑ですが、「$u$ と $h$ を使う計算方法」の(19)式は、$u$, $h$ を変数とする綺麗な線形の多項式になっています。したがって、$j$ の次ステップへの変換は行列で表現することができるはずです。しかし、(19)式は連立方程式ではないため、変数の係数をそのまま並べただけの行列にはなりません。
(19)式は、初期値の $h_1$ と $u_1$ から、焦点距離や主点の計算に使う $h_k$ と $u_{k{+}1}$ を求めるための計算式でした。これを変換行列で表せば、最終的に次のような形になるはずです。この式を用いて計算する方法を光線行列解析あるいは ABCD行列解析と呼びます。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
\left[
\begin{array}{c}
h_k\\
u_{k{+}1}
\end{array}
\right]
=
\left[
\begin{array}{cc}
A & B\\
C & D
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
=
\boldsymbol{M}
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
\end{align}\tag{30}
}
$h$ と $u$ の配置の順序が(19)式とは逆になっていますが、殆ど全ての文書でこの並びが採用されています。これを変えると $A$, $B$, $C$, $D$ の各要素が持つ意味合いが変わってしまうので、本記事でも、この慣例に従って書くことにします。
それでは、$A$, $B$, $C$, $D$ の値を求めていきましょう。$[h_1,\ u_1]^{\rm T}$ を入力ベクトルとして、(19)式の $j{=}1$ での１行目の処理を行列で表すと、次のようになります。(19)式の１行目では $u_2$ を計算しているだけなので、 $h_1$ は何も変化せず、そのままの値が出力されます。この行列を $\boldsymbol{R}_1$ と置きます。
\begin{align}
\left[
\begin{array}{c}
h_1\\
u_2
\end{array}
\right]
=
\left[
\begin{array}{cc}
1 & 0\\
(1{-}n_1/n_2)/r_1 & n_1/n_2
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
=
\boldsymbol{R}_1
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
\end{align}\tag{31}
次に、上式の出力ベクトルを入力として、(19)式の２行目の式を行列で表すと次のようになります。ここでは、$u_2$ は変化せず $h_2$ の計算だけが行われます。この行列を $\boldsymbol{T}_1$ と置きます。
\begin{align}
\left[
\begin{array}{c}
h_2\\
u_2
\end{array}
\right]
=
\left[
\begin{array}{cc}
1 & -d_2\\
0 & 1
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
u_2
\end{array}
\right]
=
\boldsymbol{T}_1
\left[
\begin{array}{c}
h_1\\
u_2
\end{array}
\right]
\end{align}\tag{32}
これを $j{=}k{-}1$ まで繰り返し、最後の $j{=}k$ で(19)式の１行目だけを計算すると次のようになります。
\begin{align}
\left[
\begin{array}{c}
h_k\\
u_{k{+}1}
\end{array}
\right]
=
\boldsymbol{R}_k
\boldsymbol{T}_{k{-}1}
\boldsymbol{R}_{k{-}1}
\cdots
\boldsymbol{T}_2
\boldsymbol{R}_2
\boldsymbol{T}_1
\boldsymbol{R}_1
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
\end{align}\tag{33}
(30)式と(33)式を対比させ、$A$, $B$, $C$, $D$ の値は次の式から求められます。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
&\boldsymbol{M}=
\left[
\begin{array}{cc}
A & B\\
C & D
\end{array}
\right]
=
\boldsymbol{R}_k
\boldsymbol{T}_{k{-}1}
\boldsymbol{R}_{k{-}1}
\cdots
\boldsymbol{T}_2
\boldsymbol{R}_2
\boldsymbol{T}_1
\boldsymbol{R}_1\\[1ex]
&\boldsymbol{R}_j
=
\left[
\begin{array}{cc}
1 & 0\\
(1{-}n_j/n_{j{+}1})/r_j & n_j/n_{j{+}1}
\end{array}
\right]\\[1.5ex]
&\boldsymbol{T}_j
=
\left[
\begin{array}{cc}
1 & -d_{j{+}1}\\
0 & 1
\end{array}
\right]
\end{align}\tag{34}
}
したがって、焦点距離 $f$ は、(24)式と(30)式から、
\begin{align}
f=\frac{h_1}{u_{k{+}1}}
=
\frac{h_1}{Ch_1+Du_1}
\end{align}
ここで、レンズへの入射光は光軸と平行なので、(19)式の計算のときと同様に $u_1{=}0$ とおくと、
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f=\frac{1}{C}
\end{align}\tag{35}
}
最後面から後側主点までの距離 $s_{\rm p}$ は、(25)式を(30)式で変形して、
s_{\rm p}=\frac{h_k-h_1}{u_{k{+}1}}=\frac{Ah_1+Bu_1-h_1}{Ch_1+Du_1}
ここでも $u_1{=}0$ とおくと、
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm p}=\frac{A-1}{C}
\end{align}\tag{36}
}
いままで求めてきた焦点距離と主点は、すべてレンズの後側（posterior、像側）のものです。ここからは、$\color{#aaa}{\boldsymbol{A}}$,$\color{#aaa}{\boldsymbol{B}}$,$\boldsymbol{C}$,$\boldsymbol{D}$ を使って前側（anterior、物体側）の焦点距離と主点を求めていきます。最終段の第 $k$ 面からの出射光が光軸に平行になるように、第１面への入射光を調整し、このときの入出射光の高さと入射角を求めれば、これらの値から計算することができます（図６）。図６は、光線が前方から後方に進むように描かれています、しかし、光は逆方向に進んでも同一の光路を通る性質があるので、後方からの平行光線が前側の焦点に収束している図であるとみなすこともできます。
図６ 焦点距離と主点（前側、物体側）
前側の特性を計算するために、逆行列を利用して(30)式の入出力を入れ替えます。
\begin{align}
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
&=
\boldsymbol{M}^{-1}
\left[
\begin{array}{c}
h_k\\
u_{k{+}1}
\end{array}
\right]
=
\frac{{\rm adj}(\boldsymbol{M})}{|\boldsymbol{M}|}
\left[
\begin{array}{c}
h_k\\
u_{k{+}1}
\end{array}
\right]\\[2ex]
&=
\frac{
\left[
\lower0.3ex{
\begin{array}{cc}
D & {-}B\\[-1ex]
{-}C & A
\end{array}
}
\right]
}{\lower0.8ex{|\boldsymbol{M}|}}
\left[
\begin{array}{c}
h_k\\
u_{k{+}1}
\end{array}
\right]
\end{align}\tag{37}
ここで、(34)式より、$|\boldsymbol{R}_j|=n_j/n_{j{+}1}$、$|\boldsymbol{T}_j|=1$ であることを考慮すると、$|\boldsymbol{M}|$ （${=}AD{-}BC$）の値は、
\begin{align}
|\boldsymbol{M}|
&=
|\boldsymbol{R}_k|
|\boldsymbol{T}_{k{-}1}|
|\boldsymbol{R}_{k{-}1}|
\cdots
|\boldsymbol{T}_2|
|\boldsymbol{R}_2|
|\boldsymbol{T}_1|
|\boldsymbol{R}_1|\\
&=
\frac{\cancel{n_k}}{n_{k{+}1}}\frac{\cancel{n_{k{-}1}}}{\cancel{n_k}}
\cdots
\frac{\cancel{n_2}}{\cancel{n_3}}\frac{n_1}{\cancel{n_2}}
=\frac{n_1}{n_{k{+}1}}
\end{align}\tag{38}
したがって、
\begin{align}
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
=
\frac{n_{k{+}1}}{n_1}
\left[
\begin{array}{cc}
D & {-}B\\
{-}C & A
\end{array}
\right]
\left[
\begin{array}{c}
h_k\\
u_{k{+}1}
\end{array}
\right]
\end{align}\tag{39}
一般には、レンズを空気中で使うことだけを想定して、$\boldsymbol{n_1{=}n_{k{+}1}}$ とし、早々と $\boldsymbol{n_{k{+}1}/n_1{=}1}$ と置いてしまうことも多いようです。しかし、この記事では、特殊環境にも対応できるように、(39)式のまま説明を続けることにします。
その前に注意しておくことがあります。図６では入射角が「右上がり」となっており、基本モード（図１）の条件「右下がりがローカル正」とは逆です。そこで、この図では、注意喚起のために、角度 $u_1$ の文字を赤色に変えています。この場合、$u_1$ の中身の数値は負の値になっているはずです。しかし、各モードでの立式のときの既述の規則（２.２節）では、「角度はすべて正の値」と規定していました。ここでもそれに準じた扱いが必要になってきます。そこで、赤色の部分の $u_1$ は ${-}u_1$ と置き換えることで正の値に変換します。
すると、
\begin{align}
&h_1=(0-z_{\rm fa})(-u_1)=z_{\rm fa}u_1\\[1ex]
&h_k=h_{\rm pa}=(z_{\rm pa}-z_{\rm fa})(-u_1)=(z_{\rm fa}-z_{\rm pa})u_1
\end{align}
両式から $z_{\rm fa}$ と $z_{\rm pa}$ を求めて(39)式を代入し、出射光が光軸と平行であることから $u_{k{+}1}{=}0$ とおくと、
\begin{align}
z_{\rm fa}&=\frac{h_1}{u_1}=\frac{\frac{n_{k{+}1}}{n_1}(Dh_k-Bu_{k{+}1})}{\frac{n_{k{+}1}}{n_1}(-Ch_k+Au_{k{+}1})}=-\frac{D}{C}\\[2ex]
z_{\rm pa}&=-\frac{h_k}{u_1}+z_{\rm fa}\\
&=-\frac{h_k}{\frac{n_{k{+}1}}{n_1}(-Ch_k+Au_{k{+}1})}+z_{\rm fa}\\
&=\frac{n_1}{n_{k{+}1}C}-\frac{D}{C}=\frac{{n_1}/{n_{k{+}1}}-D}{C}
\end{align}
以上より、前側焦点距離 $f_{\rm a}$（ $ {=}z_{\rm pa}{-}z_{\rm fa}$ ）および 最前面から前側主点までの距離 $s_{\rm pa}$（ ${=}z_{\rm pa}{-}0$ ）は次のようになります。主点が最前面よりも 前側/後側 にあるとき、$s_{\rm pa}$ の符号はそれぞれ 負/正 になります。なお、上記 $f_{\rm a}$ の引き算の順序が後側焦点距離の計算時とは逆になっているのは、凸/凹レンズの $f$ の符号を 正/負 に合わせるための意図的な操作です。（一部には、前側焦点距離の符号を後側とは逆にしている文書も見られますが、どうも主流ではなさそうなので、ここでは見習わないことにしました。）
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f_{\rm a}=\frac{n_1}{n_{k{+}1}C}
\end{align}\tag{40}
}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm pa}=\frac{{n_1}/{n_{k{+}1}}-D}{C}
\end{align}\tag{41}
}
光線行列解析は、計算開始時の暫定的な変数の仮定も必要なく、特性の計算に使う数値は最終値の $A$, $C$, $D$ だけで済むうえ、前側の特性の計算も簡単にできて便利です。しかし、計算量はかなり増えそうです。そうは言っても、行列計算が簡単にできるツールがあればそれほど気にはならないでしょう。
３.６ 光線行列解析 h,nu 方式
前３.５節の(30)式は、３.２節の「$u$ と $h$ を使う計算方法」の応用でした。これと同様に、３.４節の「$nu$ と $h$ を使う計算方法」を利用した光線行列解析も可能です。この場合、(30)式は次のように変わります。一般には、どちらの方法を使っている場合でも $A$, $B$, $C$, $D$ の文字をそのまま使いますが、本記事では、後者の場合には「$'$」をつけて区別することにします。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
\left[
\begin{array}{c}
h_k\\
(nu)_{k{+}1}
\end{array}
\right]
=
\left[
\begin{array}{cc}
A' & B'\\
C' & D'
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
(nu)_1
\end{array}
\right]
=
\boldsymbol{M}'
\left[
\begin{array}{c}
h_1\\
(nu)_1
\end{array}
\right]
\end{align}\tag{42}
}
(29)式を参照し、３.５節と同様にして $\boldsymbol{R}'_j$ , $\boldsymbol{T}'_j$ 行列を作れば、 $A'$, $B'$, $C'$, $D'$ の値は次の計算で求められます。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
&\boldsymbol{M}'=
\left[
\begin{array}{cc}
A' & B'\\
C' & D'
\end{array}
\right]
=
\boldsymbol{R}'_k
\boldsymbol{T}'_{k{-}1}
\boldsymbol{R}'_{k{-}1}
\cdots
\boldsymbol{T}'_2
\boldsymbol{R}'_2
\boldsymbol{T}'_1
\boldsymbol{R}'_1\\[1ex]
&\boldsymbol{R}'_j
=
\left[
\begin{array}{cc}
1 & 0\\
(n_{j{+}1}{-}n_j)/r_j & 1
\end{array}
\right]\\[1.5ex]
&\boldsymbol{T}'_j
=
\left[
\begin{array}{cc}
1 & -d_{j{+}1}/n_{j{+}1}\\
0 & 1
\end{array}
\right]
\end{align}\tag{43}
}
焦点距離や主点の位置を表す数式も、前節の手順と同様にして導くことはできます。しかし、似たような変形ばかりを繰り返すのも面白くないので、別の方法でやってみましょう。
$A'$, $B'$, $C'$, $D'$ と前節の $A$, $B$, $C$, $D$ との関係を調べるために式を変形してみます。(30)式の入出力と(42)式の入出力の間には次の関係があります。
\begin{align}
\left[
\begin{array}{c}
h_1\\
(nu)_1
\end{array}
\right]
=
\left[
\begin{array}{cc}
1 & 0\\
0 & n_1
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
\end{align}
\begin{align}
\left[
\begin{array}{c}
h_k\\
(nu)_{k{+}1}
\end{array}
\right]
=
\left[
\begin{array}{cc}
1 & 0\\
0 & n_{k{+}1}
\end{array}
\right]
\left[
\begin{array}{c}
h_k\\
u_{k{+}1}
\end{array}
\right]
\end{align}
これらを(42)式に代入して変形していき、最後に(30)式と対比させます。
\begin{align}
\left[
\begin{array}{cc}
1 & 0\\
0 & n_{k{+}1}
\end{array}
\right]
\left[
\begin{array}{c}
h_k\\
u_{k{+}1}
\end{array}
\right]
=
\left[
\begin{array}{cc}
A' & B'\\
C' & D'
\end{array}
\right]
\left[
\begin{array}{cc}
1 & 0\\
0 & n_1
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
\end{align}
\begin{align}
\left[
\begin{array}{c}
h_k\\
u_{k{+}1}
\end{array}
\right]
=
\left[
\begin{array}{cc}
1 & 0\\
0 & n_{k{+}1}
\end{array}
\right]^{{-}1}
\left[
\begin{array}{cc}
A' & B'\\
C' & D'
\end{array}
\right]
\left[
\begin{array}{cc}
1 & 0\\
0 & n_1
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
\end{align}
\begin{align}
\left[
\begin{array}{c}
h_k\\
u_{k{+}1}
\end{array}
\right]
&=
\left[
\begin{array}{cc}
A' & n_1B'\\
(1/n_{k{+}1})C' & (n_1/n_{k{+}1})D'
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]\\[1ex]
&=
\left[
\begin{array}{cc}
A & B\\
C & D
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
\end{align}
この結果から、$A=A'$ 、$B=n_1B'$ 、$C=(1/n_{k{+}1})C'$ 、$D=(n_1/n_{k{+}1})D'$ なる対応関係があることが分かります。これを前節の(35), (36), (40), (41)式に代入すれば、後側の焦点距離 $f$、主点位置 $s_{\rm p}$、前側の焦点距離 $f_{\rm a}$、主点位置 $s_{\rm pa}$ が一気に求まります。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f=\frac{n_{k{+}1}}{C'}
\end{align}\tag{44}
}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm p}=\frac{n_{k{+}1}(A'-1)}{C'}
\end{align}\tag{45}
}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f_{\rm a}=\frac{n_1}{C'}
\end{align}\tag{46}
}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm pa}=\frac{n_1(1-D')}{C'}
\end{align}\tag{47}
}
３.７ 流儀２に対応した読み替え
図３から、流儀２での基本式は次のようになり、流儀１（図１）のものとは若干異なります（赤色部分で項の順序が入れ替わっています）。なお、本節での数式番号は、対応する流儀１の式と同番号とし「$'$」を付けて区別しています。
\begin{align}
&n\sin i=n'\sin i'\tag{1'}\\
&r\sin i=(\color{red}{r-s})\sin u\tag{2'}\\
&r\sin i'=(\color{red}{r-s'})\sin u'\tag{3'}\\
&\delta s=r(1-\cos\varphi)\tag{4'}\\
&h=(\color{red}{\delta s-s})\tan u=(\color{red}{\delta s-s'})\tan u'\tag{5'}\\
\end{align}
上記の違いによって、後続の式も影響を受けていきます。しかし、大半の式は、両辺に ${-}1$ が掛かる程度の変化で済み、結果的に、対応する流儀１の式と等価のままで維持されます。
以降は、流儀１から実質的に変化する部分だけを重点的に拾い上げていきます。
(5)式が(5')式に変わったことにより、近軸として近似した後の(9), (10), (11)式も次のように変化します。
h=\color{red}{\boldsymbol{-}}su=\color{red}{\boldsymbol{-}}s'u'\tag{9'}
\begin{align}
n'u'=nu\color{red}{\boldsymbol{-}}\frac{h(n'-n)}{r}\tag{10'}
\end{align}
\begin{align}
\bbox[white, 5pt, border: 1.5px solid black]{
u'=\frac{n}{n'}u\color{red}{\boldsymbol{-}}\frac{h}{r}(1-\frac{n}{n'})}\tag{11'}
\end{align}
しかし、(11')式に、(9')式からの $u{=}{-}h/s$ 、$u'{=}{-}h/s'$ を代入すると、これらの変化は相殺されて、$s'$ を計算する(12), (13)式はそのまま維持されます。角度の変数を含まない数式が持つ強みです。
次に変わるのは、境界面間での $h$ の変化を表す(17)式です。流儀１では「右下がりが正」の角度だったため、光路が先に進むほど $h$ が小さくなりました。しかし、流儀２では「右上がりが正」なので、変化分が逆になり符号が変わります。
\begin{align}
\bbox[white, 5pt, border: 1.5px solid black]{
h_{j{+}1}=h_j\color{red}{\boldsymbol{+}}d_{j{+}1}u_{j{+}1}\tag{17'}
}\end{align}
また、(9)式 が (9')式 に変わったことにより、(18)式 は次のようになります。
\begin{align}
h_j=\color{red}{\boldsymbol{-}}s_ju_j=\color{red}{\boldsymbol{-}}s'_ju'_j=\color{red}{\boldsymbol{-}}s'_ju_{j{+}1}\tag{18'}
\end{align}
=====「$\boldsymbol{u}$ と $\boldsymbol{h}$ を使う計算方法」の部分 =====
(11')式と(17')式より、
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
&u_{j{+}1}=\frac{n_j}{n_{j{+}1}}u_j\color{red}{\boldsymbol{-}}\frac{h_j}{r_j}(1-\frac{n_j}{n_{j{+}1}})\\[1ex]
&h_{j{+}1}=h_j\color{red}{\boldsymbol{+}}d_{j{+}1}u_{j{+}1}
\end{align}\tag{19'}
}
図５の $u_{k{+}1}$ は、流儀１ではローカル正でしたが、流儀２ではローカル負になります。そこで、流儀１の図６において $u_1$ を ${-}u_1$ と置き換えたのと同様に、(20), (21)式の $u_{k{+}1}$ を ${-}u_{k{+}1}$ に置き換えると、
\begin{align}
&h_k=(z_{\rm f}-z_k)(\color{red}{\boldsymbol{-}}u_{k{+}1})\tag{20'}\\[1ex]
&h_1=h_{\rm p}=(z_{\rm f}-z_{\rm p})(\color{red}{\boldsymbol{-}}u_{k{+}1})\tag{21'}
\end{align}
\begin{align}
&z_{\rm f}=\color{red}{\boldsymbol{-}}\frac{h_k}{u_{k{+}1}}+z_k\tag{22'}\\
&z_{\rm p}=\color{red}{\boldsymbol{+}}\frac{h_1}{u_{k{+}1}}+z_{\rm f}=\frac{\color{red}{\boldsymbol{-}}h_k\color{red}{\boldsymbol{+}}h_1}{u_{k{+}1}}+z_k\tag{23'}
\end{align}
これより、 $f$（ ${=}z_{\rm f}{-}z_{\rm p}$ ）と $s_{\rm p}$（ ${=}z_{\rm p}{-}z_k$ ）は、
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f=\color{red}{\boldsymbol{-}}\frac{h_1}{u_{k{+}1}}
\end{align}
}\tag{24'}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm p}=\color{red}{\boldsymbol{-}}\frac{h_k-h_1}{u_{k{+}1}}
\end{align}
}\tag{25'}
=====「$\boldsymbol{s}$ を使う計算方法」の 部分 =====
先に示したように $s'$ を計算する(12), (13)式は変化しないので、(26), (27), (28)式はそのまま維持されます。
=====「$\boldsymbol{nu}$ と $\boldsymbol{h}$ を使う計算方法」の部分 =====
(19')式と同様の符号の変化が生じます。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
&(nu)_{j{+}1}=(nu)_j\color{red}{\boldsymbol{-}}\frac{h_j}{r_j}(n_{j{+}1}-n_j)\\[1ex]
&h_{j{+}1}=h_j\color{red}{\boldsymbol{+}}d_{j{+}1}\frac{(nu)_{j{+}1}}{n_{j{+}1}}
\end{align}\tag{29'}
}
=====「光線行列解析 $\boldsymbol{h},\,\boldsymbol{u}$ 方式」の部分 =====
後側の特性計算
(19)式が(19')式に変わったことにより、次のように変化します。
\begin{align}
\left[
\begin{array}{c}
h_1\\
u_2
\end{array}
\right]
=
\left[
\begin{array}{cc}
1 & 0\\
\color{red}{\boldsymbol{-}}(1{-}n_1/n_2)/r_1 & n_1/n_2
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
=
\boldsymbol{R}_1
\left[
\begin{array}{c}
h_1\\
u_1
\end{array}
\right]
\end{align}\tag{31'}
\begin{align}
\left[
\begin{array}{c}
h_2\\
u_2
\end{array}
\right]
=
\left[
\begin{array}{cc}
1 & \color{red}{\boldsymbol{+}}d_2\\
0 & 1
\end{array}
\right]
\left[
\begin{array}{c}
h_1\\
u_2
\end{array}
\right]
=
\boldsymbol{T}_1
\left[
\begin{array}{c}
h_1\\
u_2
\end{array}
\right]
\end{align}\tag{32'}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
&\boldsymbol{M}=
\left[
\begin{array}{cc}
A & B\\
C & D
\end{array}
\right]
=
\boldsymbol{R}_k
\boldsymbol{T}_{k{-}1}
\boldsymbol{R}_{k{-}1}
\cdots
\boldsymbol{T}_2
\boldsymbol{R}_2
\boldsymbol{T}_1
\boldsymbol{R}_1\\[1ex]
&\boldsymbol{R}_j
=
\left[
\begin{array}{cc}
1 & 0\\
\color{red}{\boldsymbol{-}}(1{-}n_j/n_{j{+}1})/r_j & n_j/n_{j{+}1}
\end{array}
\right]\\[1.5ex]
&\boldsymbol{T}_j
=
\left[
\begin{array}{cc}
1 & \color{red}{\boldsymbol{+}}d_{j{+}1}\\
0 & 1
\end{array}
\right]
\end{align}\tag{34'}
}
さらに、(24), (25)式が (24'), (25')式に変わったことにより、
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f=\color{red}{\boldsymbol{-}}\frac{1}{C}
\end{align}\tag{35'}
}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm p}=\color{red}{\boldsymbol{-}}\frac{A-1}{C}
\end{align}\tag{36'}
}
前側の特性計算
流儀１では、$u_1$ の値がローカル負だったため、$u_1$ を ${-}u_1$ に置き換えました。しかし、流儀２においては、$u_1$ はローカル正の状態なので置き換えの必要はなく、無採番の数式部分は次のように変わります。
\begin{align}
&h_1=\color{red}{\boldsymbol{-}}z_{\rm fa}u_1\\[1ex]
&h_k=h_{\rm pa}=\color{red}{\boldsymbol{-}}(z_{\rm fa}-z_{\rm pa})u_1
\end{align}
\begin{align}
z_{\rm fa}&=\color{red}{\boldsymbol{-}}\frac{h_1}{u_1}=\color{red}{\boldsymbol{-}}\frac{\frac{n_{k{+}1}}{n_1}(Dh_k-Bu_{k{+}1})}{\frac{n_{k{+}1}}{n_1}(-Ch_k+Au_{k{+}1})}=\color{red}{\boldsymbol{+}}\frac{D}{C}\\[2ex]
z_{\rm pa}&=\color{red}{\boldsymbol{+}}\frac{h_k}{u_1}+z_{\rm fa}\\
&=\color{red}{\boldsymbol{+}}\frac{h_k}{\frac{n_{k{+}1}}{n_1}(-Ch_k+Au_{k{+}1})}+z_{\rm fa}\\
&=\color{red}{\boldsymbol{-}}\frac{n_1}{n_{k{+}1}C}\color{red}{\boldsymbol{+}}\frac{D}{C}=\color{red}{\boldsymbol{-}}\frac{{n_1}/{n_{k{+}1}}-D}{C}
\end{align}
これにより、$f_{\rm a}$（ ${=}z_{\rm pa}{-}z_{\rm fa}$ ）と $s_{\rm pa}$（ ${=}z_{\rm pa}{-}0$ ）は、
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f_{\rm a}=\color{red}{\boldsymbol{-}}\frac{n_1}{n_{k{+}1}C}
\end{align}\tag{40'}
}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm pa}=\color{red}{\boldsymbol{-}}\frac{{n_1}/{n_{k{+}1}}-D}{C}
\end{align}\tag{41'}
}
=====「光線行列解析 $\boldsymbol{h},\,\boldsymbol{nu}$ 方式」の部分 =====
(29)式が(29')式に変わったことにより、(43)式は次のように変わります。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
&\boldsymbol{M}'=
\left[
\begin{array}{cc}
A' & B'\\
C' & D'
\end{array}
\right]
=
\boldsymbol{R}'_k
\boldsymbol{T}'_{k{-}1}
\boldsymbol{R}'_{k{-}1}
\cdots
\boldsymbol{T}'_2
\boldsymbol{R}'_2
\boldsymbol{T}'_1
\boldsymbol{R}'_1\\[1ex]
&\boldsymbol{R}'_j
=
\left[
\begin{array}{cc}
1 & 0\\
\color{red}{\boldsymbol{-}}(n_{j{+}1}{-}n_j)/r_j & 1
\end{array}
\right]\\[1.5ex]
&\boldsymbol{T}'_j
=
\left[
\begin{array}{cc}
1 & \color{red}{\boldsymbol{+}}d_{j{+}1}/n_{j{+}1}\\
0 & 1
\end{array}
\right]
\end{align}\tag{43'}
}
(35'), (36'), (40'), (41')式に、流儀１と同様に $A=A'$ 、$B=n_1B'$ 、$C=(1/n_{k{+}1})C'$ 、$D=(n_1/n_{k{+}1})D'$ なる変換を行えば、流儀２における焦点距離や主点が次のように求まります。
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f=\color{red}{\boldsymbol{-}}\frac{n_{k{+}1}}{C'}
\end{align}\tag{44'}
}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm p}=\color{red}{\boldsymbol{-}}\frac{n_{k{+}1}(A'-1)}{C'}
\end{align}\tag{45'}
}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
f_{\rm a}=\color{red}{\boldsymbol{-}}\frac{n_1}{C'}
\end{align}\tag{46'}
}
\bbox[white, 5pt, border: 1.5px solid black]{
\begin{align}
s_{\rm pa}=\color{red}{\boldsymbol{-}}\frac{n_1(1-D')}{C'}
\end{align}\tag{47'}
}
３.８ その他の流儀への対応
２.２節の最後で「最終的に気に留めておくべきは、$u$, $u'$ の符号だけ」であることを述べました。図２には 16 種類のモードがありますが、$u$ / $u'$ それぞれの「右上がり」や「右下がり」という傾斜だけに着目すれば、下がり/下がり（流儀１）、上がり/上がり（流儀２）、下がり/上がり（流儀３）、上がり/下がり（流儀４） の４種類に集約できます。このうち、流儀１,２については既に最終式を導出済みなので、ここでは残りの流儀３,４について考えることにします。
結論から言えば、流儀３,４の検討結果は、それぞれが流儀１,２と同一の式に帰着し、最終的に残る計算式の形は、流儀１,２から導かれる２種類だけになることが分かります。
各流儀における数式を比較して表２に示しています。最上段には、入射光角 $u$ と出射光角 $u'$ について、ローカル正とみなす傾きの方向も併記しています。行 [01], [02], [03], [04] は、既述の (1), (2), (3), (5)式に相当する、各流儀における基本式を近軸近似して導かれた式です。[01] には、後続式を代入しやすいように、変形した式も併記しています。[01] に [02], [03] を代入して [05] が得られます。[04] を利用してこれを変形すると [06] になります。これを $u'$ について解くと [07] になります。これは、既述の(11)式に相当する式です。
既述の(19)式（[11], [12] に相当）は、流儀１,２だけに通用する $u'_j=u_{j{+}1}$ の関係を利用して作られた式です。条件が異なる流儀３,４まで含めて検討する場合には、一旦、[07] から [08], [09] を 経て、[10] の条件を加味して [11], [12] に至る必要があります。この結果、流儀３は流儀１と、流儀４は流儀２と同一の式になります（式が同一の枠内は、背景色も同じにしています）。
表２ 各流儀での数式の違い
以下、詳細な説明は省きますが、さらに検討を進めていくと、[15], [16] においても、流儀３,４はそれぞれ流儀１,２と同一の式になります。これらの [11], [12], [15], [16] は「$u$ と $h$ を使う計算方法」の基本式です。他の「$nu$ と $h$ を使う計算方法」や「光線行列解析」は、これらの式を元にした計算方法ですから、どの方法もすべて流儀１,２の式に集約されることになります。
残りの「$s$ を使う計算方法」についても、基本を押さえながら変形していけば、[20] に示すように全ての流儀において完全に同一の式になります。
４． 実際のレンズによる確認
４.１ 近軸近似による計算結果
これまでの検討で、計算の流儀は２種類、さらに、各流儀には計算方法が５種類もあることが分かりました。理論式ばかりをこんなに多量に並べられても目移りはするし、一部には誤記もありそうで不信感も湧いてきます。そこで、これらの不満や不安を解消するために、実際のレンズデータをモデルにして、焦点距離や主点位置を計算し、各計算方法の結果を評価してみようと思います。
これらのレンズデータは、書籍「レンズ設計のすべて - 光学設計の神髄を探る」（辻貞彦 著、電波新聞社）から借用した５種類 と 自作の２種類 です。詳細なデータは第６章に添付したプログラムの中に記述しています。
レンズ１: No.58 Tessar型 公称 f50mm F2.8
レンズ２: No.72 Ernostar型 公称 f50mm F2
レンズ３: No.84 Sonnar型 公称 f50mm F2
レンズ４: No.93 Gauss型 公称 f50mm F2
レンズ５: No.113 Gauss型 公称 f50mm F2
レンズ６: 自作単凸レンズ f51.6mm
レンズ７: 自作単凹レンズ f51.6mm
このプログラムについては、本記事中の数式を参照しながら、そのままの形で間違いのないようにコード化し、チェックにも念を入れています。これでおかしな計算結果さえ出なければ、記事中の数式も信用できることになります。この計算出力の最初の部分だけを次に示しています。どの計算方法においても同一の数値が得られていることが分かります。
No.58 Tessar型 公称50mm, F2.8
【流儀１】uとhを使う計算方法
前/後 焦点距離 51.5995/51.5995 前/後 主点位置 6.8399/-11.9054 (7.5946)
【流儀１】sを使う計算方法
前/後 焦点距離 51.5995/51.5995 前/後 主点位置 6.8399/-11.9054 (7.5946)
【流儀１】nuとhを使う計算方法
前/後 焦点距離 51.5995/51.5995 前/後 主点位置 6.8399/-11.9054 (7.5946)
【流儀１】光線行列解析 h,u 方式
前/後 焦点距離 51.5995/51.5995 前/後 主点位置 6.8399/-11.9054 (7.5946)
A B = 0.769273 -17.1672
C D 0.01938 0.867442
【流儀１】光線行列解析 h,nu 方式
前/後 焦点距離 51.5995/51.5995 前/後 主点位置 6.8399/-11.9054 (7.5946)
A B = 0.769273 -17.1672
C D 0.01938 0.867442
【流儀２】uとhを使う計算方法
前/後 焦点距離 51.5995/51.5995 前/後 主点位置 6.8399/-11.9054 (7.5946)
【流儀２】sを使う計算方法
前/後 焦点距離 51.5995/51.5995 前/後 主点位置 6.8399/-11.9054 (7.5946)
【流儀２】nuとhを使う計算方法
前/後 焦点距離 51.5995/51.5995 前/後 主点位置 6.8399/-11.9054 (7.5946)
【流儀２】光線行列解析 h,u 方式
前/後 焦点距離 51.5995/51.5995 前/後 主点位置 6.8399/-11.9054 (7.5946)
A B = 0.769273 17.1672
C D -0.01938 0.867442
【流儀２】光線行列解析 h,nu 方式
前/後 焦点距離 51.5995/51.5995 前/後 主点位置 6.8399/-11.9054 (7.5946)
A B = 0.769273 17.1672
C D -0.01938 0.867442
他のレンズについても、すべての計算方法で同一の結果が得られて表３のようになります。
表３ 各レンズの特性計算結果
| レンズ の種類 |
焦点距離 [mm] | 主点位置 [mm] | |||
|---|---|---|---|---|---|
| 前側 |
後側 |
前側 (前面基準) |
後側 | ||
| (後面基準) | (前面基準) | ||||
| No.58 Tessar型 | 51.5995 | 51.5995 | 6.8399 | -11.9054 | 7.5946 |
| No.72 Ernostar型 | 51.5999 | 51.5999 | 6.2542 | -20.6015 | 4.6596 |
| No.84 Sonnar型 | 51.5998 | 51.5998 | 4.5749 | -25.5508 | 8.9392 |
| No.93 Gauss型 | 51.6000 | 51.6000 | 23.9329 | -22.1485 | 13.3515 |
| No.113 Gauss型 | 51.6001 | 51.6001 | 34.2580 | -25.2506 | 19.4894 |
| 単凸レンズ | 51.6000 | 51.6000 | 0.8759 | -1.3860 | 2.1954 |
| 単凹レンズ | -51.6000 | -51.6000 | 0.2362 | -0.3870 | 0.6130 |
| $\hspace{7em}$ | $\hspace{4em}$ | $\hspace{4em}$ | $\hspace{4em}$ | $\hspace{4em}$ | $\hspace{4em}$ |
上記での計算条件は、レンズは空気中に存在するものとして前/後とも同一の屈折率としています。念のため、前/後が別々の屈折率の条件でも計算してみました。この場合、主点の位置は変わるし、焦点距離も前/後で異なる値になります。しかし、どの計算方法においても同一の結果が得られることには変わりがありません（結果の詳細は割愛）。
ところで、公称 50mm のレンズなのに、近軸近似での焦点距離が 51.6mm という結果になってしまいました。 「1.6mmもの誤差があるのに "精確そのもの" などと強弁するのは詐欺ではないか！」 とお叱りを受けそうです。しかし、ご安心ください。前記の参考書籍の 18 ページには「なお50mmのレンズはLeitz社のレンズに合わせて実焦点距離51.6mmとすることが多いので、以下の設計例でもそのように扱う」とあります。この 51.6mm は、設計者の意図どおりのレンズの性能です。これで疑いは晴れ、あらためて "精確そのもの" を実感することができました。
以上により、本記事中に記載された計算式はすべて正確であり、どの計算方法を使っても同一の結果が精度よく得られることが分かりました。しかし、どの流儀のどの計算方法を使うべきかは未だ決め兼ねます。これについては、第５章の「おわりに」で考えてみることにします。
４.２ 光路シミュレーションとの比較
前節までで、精度の良い計算結果が得られそうなことは分かってきました。しかし、表面的な数値だけを多量に示されても十分には納得できないでしょう。そこで、収差を抑えたレンズであれば、「光軸から離れた光束」も全て「近軸近似の焦点」に収束することをシミュレーションによる画像で最終確認してみようと思います。
図７～１１に複合レンズのシミュレーション結果を示しています。各図の上側の図は、レンズを本来の向きで使い、前側からの平行光線が後側へ収束する様子を示しています。ただし、小さな枠内の図は、焦点付近を100倍にも拡大した苛酷なものですから、まだ見なかったことにしておいてください。そうすれば、収差を補正した複合レンズであれば、近軸近似で計算した後側焦点（赤色の点）に全ての光線が正確に集束していることが分かります。これで、計算結果が正しいことを視覚的にも確認することができました。なお、緑色の点は近軸近似で計算した後側の主点です。この主点から上下に伸びる黒色の線は主面と呼ばれるもので、入射光と出射光の延長線の交点群をプロットすることで得られる面です。今回の記事では脇役ですが、灰色の縦線は解放状態の「絞り」です。
一方、下側の図は、平行光線をレンズの後側から照射した場合の様子です。普通の使い方とは逆ですが、近軸近似による前側の焦点と主点の計算結果を確認するためのシミュレーションです。残念ながら、前側焦点への収束状態は後側ほどには正確ではなく少しずれているようです。拡大図によれば、後側のズレは 0.15[mm] 前後でしたが、前側のズレはレベルが違います。しかし、前側の拡大図を細かく観察すれば、光軸に近い光線ほど焦点近くに収束している様子が確認できます。近軸近似の精度が悪いというよりも、設計者にとっては後側の収差の補正が最優先課題で、前側の収差にまでは手が回り兼ねているというのが実状のような気もします。
図１２,１３は単レンズのシミュレーション結果です。収差だらけの単レンズでは評価の余地はありません。あらためて複合レンズの有難さが分かります。
図７
図８
図９
図１０
図１１
図１２
図１３
５． おわりに
以上の検討結果から我流で判断すれば、「収差を補正する」という言葉は「光軸から離れた光線でも近軸近似の焦点に収束させる」と言い替えることができそうです。そうであれば、収差が補正された上等なレンズほど、その焦点距離は近軸近似での計算値に近づくことになります。したがって、「精確そのもの」という表現は、「近軸近似の精度が凄い」というよりも、「レンズ屋の収差の補正のレベルが凄い」と解釈する方が当たっているような気がします。
これだけ精度が良いことが分かると、近軸近似という言葉には違和感が湧いてきます。そのような視点から眺め直してみると、「近軸近似」という用語は、アマチュアや素人マニア向けの文書で多く用いられているように思えます。専門家向けの文書では「近軸追跡」,「近軸光線追跡」,「近軸解析」などの用語が使われることが多いようです。
それにしても、レンズの計算に限らず、多くの分野で流儀が乱立してるのは本当に迷惑です。なんとかならないものでしょうか？ その一助となるべく、私の好みと独断により、近軸近似の多くの計算方法の中からお薦めの方法を挙げておきます。迷うことなく、流儀１の３.５節「光線行列解析 $\boldsymbol{h}$,$\boldsymbol{u}$ 方式」あるいは、行列の計算を避けたければ３.２節「$\boldsymbol{u}$ と $\boldsymbol{h}$ を使う計算方法」を推します。
そもそも流儀２では、光の収束具合を議論しようとしているときに、発散を印象づける図３を使っているため素直な理解を妨げています。反時計回りの角度を正方向として扱うことだけは数学的に納得できますが、結果として得られる焦点距離などの式に、不自然な「$-$」符号が付くことも美しくありません。
「$s$ を使う計算方法」では、計算の途中に得られた多くの変数を、計算が終わるまで保持しておく必要があり面倒です。また、光路の途中に水平部分が生じたときの例外処理が必要です。添付のプログラムでもこの処理は考慮に入れておりません。「$nu$ と $h$ を使う計算方法」は、計算量を節減できるとは言っても、コンピュータが相手では効果は軽微です。解釈が困難な $nu$ のような複合物理量を扱っているのも馴染めません。また、この方法から発展した「光線行列解析 $h$,$nu$ 方式」も同様です。
６． プログラム
すべて MATLAB で書いたプログラムです。使用しているバージョンは、かなり古めの R2019a です。６.１節だけは、Symbolic Math Toolbox も利用しています。
６.１ モード間での式の等価性の確認
２.２節で紹介したプログラムです。実行結果の一部は下記のようになります。Formula の欄に 0, Inf, NaN 以外の数値が並べば、等価性があることが証明できます。一見、有難みが分からない結果ですが、プログラム中に登録された数式内の符号を一つでも書き換えてみれば、NaN が頻発して等価性判定の実力が実感できます。
ans =
90×4 table
Base Else Formula_1234 Formula_5678
____ ____ ____________________ ____________________
1 2 -1 -1 -1 1 1 1 1 1
1 3 1 1 1 1 1 1 1 1
1 4 1 1 1 1 1 1 1 1
1 5 1 1 1 1 1 1 1 1
1 6 -1 1 1 -1 -1 1 1 1
1 7 1 -1 -1 -1 -1 1 1 1
1 8 -1 1 1 -1 -1 1 1 1
1 9 -1 1 1 -1 -1 1 1 1
1 10 -1 1 1 -1 -1 1 1 1
2 1 -1 -1 -1 1 1 1 1 1
2 3 -1 -1 -1 1 1 1 1 1
プログラム（折り畳んであります）
% Lens50.m
% 角度変数の符号をうまく読み替えることで、
% 近軸近似の全10モードの動作が、
% 一つのモードの図に基づく一組の式だけで証明できることの裏付け。
% （全10モード： Qiitaの記事中の図２の全16モードのうち、数式の形が
% 同じものについては１つとして数えると10種類になる）
%
% ２つの等式 f(x)=g(x) と u(x)=v(x) が等価であることの確認方法。
% まずは下記のように変形する。
% f(x)-g(x)=0 と u(x)-v(x)=0
% 両式に次の関係があれば、文句なしの等価であるはず。
% {f(x)-g(x)}/{u(x)-v(x)} = [0や∞ではない数値]
% 他にも、等価な条件はあるかも知れないが、上記だけで等価であることが
% 分かれば、それ以上の余計なことを考える必要はない。ダメなときは考え
% 直せば良い。
% 結論としては、今回の課題については、上記だけで十分だった。
clear
close all
% シンボル変数の宣言
% 座標の原点は、球面と光軸との交点にあるものとする。
syms i id u ud ph n nd r s sd ds h
% i: 入射角
% id: 出射角
% u: 入射光と光軸との交角
% ud: 出射光と 〃
% ph: 入射点での球面の法線と光軸との交角
% n: 入射側の屈折率
% nd: 出射側 〃
% r: 球面の中心の座標
% s: 入射光と光軸との交点の座標
% sd: 出射光と 〃
% ds: 球面上の入射点の横軸座標
% h: 〃 縦軸座標
% 角度 i,id,u,ud,ph のグローバル符号の定義。
% 各モードの図の中に示す式は、図中の全ての角度が正である
% と仮定（ローカル符号）して導出されたもの。グローバル符号は
% これとは別物で、i,id,u,ud,ph のそれぞれに、全モードに共通の次の
% 判断基準によって +1/-1 を定義したもの。
% ただし、法線や水平線を回転させる方向は、
% 90度未満の回転角で、目的の線に重ねられる方向とする。
% （注： Qiita記事中では、 +1/-1 = global正/global負）
%
% i,id : 法線を反時計回りに回転した傾きのとき +1。
% 〃 時計回り 〃 -1。
% u,ud,ph : 水平線（z 軸）を反時計回りに回転した傾きのとき +1。
% 〃 時計回り 〃 -1。
% i,id,u,ud,ph のグローバル符号（行番号 = モード番号）
% 各モードの図と上記の基準とを対比させて、人為的に仕分けた結果。
global_sign=[
+1 +1 -1 -1 -1; % 図２の①
-1 -1 -1 -1 -1; % 〃 ②
+1 +1 +1 +1 -1; % 〃 ③
+1 +1 +1 -1 -1; % 〃 ④
+1 +1 -1 +1 -1; % 〃 ⑤
-1 -1 +1 +1 +1; % 〃 ⑥
+1 +1 +1 +1 +1; % 〃 ⑦
-1 -1 -1 -1 +1; % 〃 ⑧
-1 -1 -1 +1 +1; % 〃 ⑨
-1 -1 +1 -1 +1; % 〃 ⑩
];
% 各モードの図中に書かれている数式。
% formulaの第3引数がモード番号。
% 第1列が式の左辺、第2列が式の右辺。
formula(:,:,1)=[
n*sin(i) nd*sin(id)
(r-0)*sin(i) (s-r)*sin(u)
(r-0)*sin(id) (sd-r)*sin(ud)
ph i+u
ph id+ud
ds (r-0)*(1-cos(ph))
h (s-ds)*tan(u)
h (sd-ds)*tan(ud)
];
formula(:,:,2)=[
n*sin(i) nd*sin(id)
(r-0)*sin(i) (r-s)*sin(u)
(r-0)*sin(id) (r-sd)*sin(ud)
ph -i+u
ph -id+ud
ds (r-0)*(1-cos(ph))
h (s-ds)*tan(u)
h (sd-ds)*tan(ud)
];
formula(:,:,3)=[
n*sin(i) nd*sin(id)
(r-0)*sin(i) (r-s)*sin(u)
(r-0)*sin(id) (r-sd)*sin(ud)
ph i-u
ph id-ud
ds (r-0)*(1-cos(ph))
h (ds-s)*tan(u)
h (ds-sd)*tan(ud)
];
formula(:,:,4)=[
n*sin(i) nd*sin(id)
(r-0)*sin(i) (r-s)*sin(u)
(r-0)*sin(id) (sd-r)*sin(ud)
ph i-u
ph id+ud
ds (r-0)*(1-cos(ph))
h (ds-s)*tan(u)
h (sd-ds)*tan(ud)
];
formula(:,:,5)=[
n*sin(i) nd*sin(id)
(r-0)*sin(i) (s-r)*sin(u)
(r-0)*sin(id) (r-sd)*sin(ud)
ph i+u
ph id-ud
ds (r-0)*(1-cos(ph))
h (s-ds)*tan(u)
h (ds-sd)*tan(ud)
];
formula(:,:,6)=[
n*sin(i) nd*sin(id)
(0-r)*sin(i) (r-s)*sin(u)
(0-r)*sin(id) (r-sd)*sin(ud)
ph i+u
ph id+ud
ds (r+(0-r)*cos(ph))
h (ds-s)*tan(u)
h (ds-sd)*tan(ud)
];
formula(:,:,7)=[
n*sin(i) nd*sin(id)
(0-r)*sin(i) (s-r)*sin(u)
(0-r)*sin(id) (sd-r)*sin(ud)
ph -i+u
ph -id+ud
ds (r+(0-r)*cos(ph))
h (ds-s)*tan(u)
h (ds-sd)*tan(ud)
];
formula(:,:,8)=[
n*sin(i) nd*sin(id)
(0-r)*sin(i) (s-r)*sin(u)
(0-r)*sin(id) (sd-r)*sin(ud)
ph i-u
ph id-ud
ds (r+(0-r)*cos(ph))
h (s-ds)*tan(u)
h (sd-ds)*tan(ud)
];
formula(:,:,9)=[
n*sin(i) nd*sin(id)
(0-r)*sin(i) (s-r)*sin(u)
(0-r)*sin(id) (r-sd)*sin(ud)
ph i-u
ph id+ud
ds (r+(0-r)*cos(ph))
h (s-ds)*tan(u)
h (ds-sd)*tan(ud)
];
formula(:,:,10)=[
n*sin(i) nd*sin(id)
(0-r)*sin(i) (r-s)*sin(u)
(0-r)*sin(id) (sd-r)*sin(ud)
ph i+u
ph id-ud
ds (r+(0-r)*cos(ph))
h (ds-s)*tan(u)
h (sd-ds)*tan(ud)
];
BBB=[]; % 結果を収納するための行列
old=[i id u ud ph]; % 符号を読み替える前のそのままの角度変数
for nmb=1:10 % 基本モードとして採用するモード番号
fomb=formula(:,1,nmb)-formula(:,2,nmb);
% 基本モードの8行の式の全項を左辺へ移項した式。
for nme=1:10 % 異種モードとして扱うモード番号
if nmb==nme % 同一モードは元々等価なのでスキップ
continue;
end
fomo=formula(:,1,nme)-formula(:,2,nme);
% 異種モードの8行の式の全項を左辺へ移項した式。
Sr=global_sign(nmb,:).*global_sign(nme,:);
% モード間の、グローバル符号の相対的な差異。
% 基本モードと異種モードで同符号の部分は +1、
% 異符号の部分は -1。
% この結果が -1 となった角度変数は、次ステップ
% で負の値に置き換え、+1 となった角度変数は
% 正の値のままとする。
new=old.*Sr; % 読み替えたあとの符号に変換した角度変数
fombx=subs(fomb,old,new);
% 移項後の基本モードの式に含まれる全角度変数を、
% 読み替え後の角度変数と置き換える。
quo=simplify(fomo./fombx)';
% fomo./fombx が0,∞以外の数値だけの配列に
% なれば、異種モードの式と、置き換え後の基本
% モードの式とは等価である。
AAA=zeros([1 8]);
for n=1:8 %
try
dummy=double(quo(n)); % quoが数値ではなく変数のときは、
% 本来ならばエラーになるはず。
catch % エラーになることが分かったとき。
AAA(n)=NaN;
disp(['基本モード ' num2str(nmb) ...
' と 異種モード ' num2str(nme) ' 間で、第 ' ...
num2str(n) ' 式が不整合（変数あり）']);
end
if ~isnan(AAA(n))
AAA(n)=double(quo(n)); % 数値化
if AAA(n)==0 | isinf(AAA(n))
disp(['基本モード ' num2str(nmb) ...
' と 異種モード ' num2str(nme) ' 間で、第 ' ...
num2str(n) ' 式が不整合（0,Inf あり）']);
end
end
end
% 一括して一つの行列にまとめる。
BBB=[BBB;nmb nme AAA];
% nmb,nme は基本モード番号と異種モード番号。
end
end
% 裏付けがとれたことのテーブル表示。
% Formula_1234, Formula_5678 欄に NaN や 0 や Inf が皆無であれば、
% 等価であることが証明できたことになる。
Base=BBB(:,1);
Else=BBB(:,2);
Formula_1234 =BBB(:,[3:6]); % 一つの欄内の列数が多すぎると表示できな
Formula_5678=BBB(:,[7:10]); % い仕様のようなので、２つに分割した。
table(Base,Else,Formula_1234,Formula_5678)
６.２ 近軸近似計算方法、10 種類での結果の比較
４.１節で使ったプログラムです。180 行目あたりの % （レンズの前後の媒質が空気以外のときを想定） の下にある２行分の「%」を外せば、レンズの前後が空気以外の場合の計算例も表示できます。
プログラム（折り畳んであります）
% Lens51.m
% レンズの「各種の近軸近似式」による焦点距離と主点の計算。
%
% ■モデルとして、プロが設計した下記の５種類の複合レンズと、
% 自作の凸/凹単レンズ各１種を利用する。
%
% 「レンズ設計のすべて - 光学設計の神髄を探る」（辻貞彦 著、電波新聞社）
% には、公称50mm,F2のレンズの34種の設計例が挙げられているが、その中から
% No.(58),72,84,93,113のレンズを選んだ。（No.58はF2.8）
%
% No.93を例として、文献中の生のデータの内容を示せば次のとおり。
%
% 面 有効径 r(曲率半径) d(肉厚, nd(硝材 vd(ア
% 空気空隙) 屈折率) ッベ数)
% 1 32.4 28.45336 5.00000 1.65844 50.88
% 2 31.0 81.93668 0.10000 1.00000
% 3 24.9 17.86823 6.50000 1.64850 53.02
% 4 22.0 244.84648 1.30000 1.64769 33.79
% 5 16.9 11.87418 6.50000 1.00000
% 6 P(絞り面) 15.9 0.00000 6.50000 1.00000
% 7 14.6 -13.85226 1.00000 1.64769 33.79
% 8 17.1 -49.52634 4.00000 1.64850 53.02
% 9 19.8 -19.13128 0.10000 1.00000
% 10 23.9 667.74212 4.50000 1.65844 50.88
% 11 25.0 -29.07016 0.00000 1.00000
%
% 【読み方】
% ●「面」はレンズの面番号。独立した1枚のレンズは、表と裏に1面ずつ、
% 合わせて2つの面を持つ。ただし、レンズ同士が張り合わされてできた
% 接着面は、2面ではなく、1つの面として数える。
% なお、これ以外にも、絞りが配置される面も1つの面として数える。
% ●「有効径」はレンズ面の直径[mm]。あるいは、開放絞りの直径を表す。
% ●「r(曲率半径)」はレンズ面の曲率半径[mm]。
% 正の数のときは被写体側に凸、負の数のときはフィルム側に凸を示す。
% 上例のレンズには無いが、レンズ面が平面のときにはr=0.00000と表す。
% また、絞り面もr=0.00000と表す。
% ●「d(肉厚,空気空隙)」は、その面から次の面までの、光軸上の距離[mm]
% を表す。
% ●「nd(硝材屈折率)」は、その面から次の面の間の素材の屈折率を表す。
% ●「vd(アッベ数)」は、色収差の起こり難さを示すガラス材の定数。
% νd=(nd-1)/(nF-nC)。ただし、nF,nd,nCは、波長486.1, 587.6, 656.3[nm]
% （青緑,黄,赤）における屈折率。
% 本プログラムにおいては、上記の表を次のように改変して利用する。
% ・レンズの平面や絞り面は、r=0.00000ではなく、r=Infとする。
% 平面の曲率半径は 0 ではなく∞だから、それが当然。
% ・絞り面も無視せず、その前後が空気層からなる1つのレンズ平面とみなす。
% ・行列中に空欄はまずいので、空気のアッベ数はInfとする。
% しかし、ここでは色収差は無視するので、アッベ数は利用しない。
% ・変数の添字の決め方
% 面 : 物体側の面から j=1,2,3, ..., k。 全 k 面。
% 面の曲率半径 r : その面の添字と同じ
% 屈折率 n : その媒質の後側にある面の添字と同じ
% 面の間隔 d : その区間の後側にある面の添字と同じ
clear
close all
% ==================================================================
% レンズの仕様を設定（すべて、公称の焦点距離50mm、実焦点距離51.6mm）
% ==================================================================
% No.58 Tessar型 F2.8
% No.72 Ernostar型 F2
% No.84 Sonnar型 F2
% No.93 Gauss型 F2
% No.113 Gauss型 F2
% 自作の単凸レンズ
% 自作の単凹レンズ
% 有効径 曲率半径 肉厚d 屈折率 アッベ数
% No.058
FFo(1)=2.8; % 開放F値
Lnamex{1}=['No.58 Tessar型 公称50mm, F2.8'];
Surx{1}=[
19.1 19.48717 4.50000 1.77250 49.60
17.3 91.06909 4.00000 1.00000 Inf
15.2 -62.93988 2.00000 1.72825 28.46
14.1 17.94764 2.00000 1.00000 Inf
14.2 Inf 2.00000 1.00000 Inf
14.6 141.35119 1.00000 1.59270 35.31
14.7 24.49046 4.00000 1.80440 39.59
14.7 -37.36912 0.00000 1.00000 Inf
];
% No.072
FFo(2)=2; % 開放F値
Lnamex{2}=['No.72 Ernostar型 公称50mm, F2'];
Surx{2}=[
29.0 22.09251 6.00000 1.65844 50.88
27.5 123.19285 0.10000 1.00000 Inf
22.4 21.27789 4.00000 1.63854 55.38
20.1 32.11702 2.50000 1.00000 Inf
19.2 -424.31820 1.20000 1.76182 26.52
17.0 15.47789 4.00000 1.00000 Inf
16.9 Inf 4.96104 1.00000 Inf
16.6 58.13366 2.50000 1.72342 37.93
16.4 -52.72323 0.00000 1.00000 Inf
];
% No.084
FFo(3)=2; % 開放F値
Lnamex{3}=['No.84 Sonnar型 公称50mm, F2'];
Surx{3}=[
28.7 32.73326 3.80000 1.69680 55.53
27.4 97.63145 0.19000 1.00000 Inf
23.5 17.57712 4.50000 1.69680 55.53
22.0 36.68896 3.30000 1.48749 70.21
20.5 2255.56502 1.00000 1.64769 33.79
16.8 12.10293 4.50000 1.00000 Inf
16.5 Inf 1.20000 1.00000 Inf
16.2 -108.10542 6.00000 1.53172 48.84
15.3 27.87691 10.00000 1.77250 49.60
19.8 -64.34508 0.00000 1.00000 Inf
];
% No.093
FFo(4)=2; % 開放F値
Lnamex{4}=['No.93 Gauss型 公称50mm, F2'];
Surx{4}=[
32.4 28.45336 5.00000 1.65844 50.88
31.0 81.93668 0.10000 1.00000 Inf
24.9 17.86823 6.50000 1.64850 53.02
22.0 244.84648 1.30000 1.64769 33.79
16.9 11.87418 6.50000 1.00000 Inf
15.9 Inf 6.50000 1.00000 Inf
14.6 -13.85226 1.00000 1.64769 33.79
17.1 -49.52634 4.00000 1.64850 53.02
19.8 -19.13128 0.10000 1.00000 Inf
23.9 667.74212 4.50000 1.65844 50.88
25.0 -29.07016 0.00000 1.00000 Inf
];
% No.113
FFo(5)=2; % 開放F値
Lnamex{5}=['No.113 Gauss型 公称50mm, F2'];
Surx{5}=[
38.5 41.22676 5.50000 1.80610 40.92
36.7 141.20948 0.10000 1.00000 Inf
28.7 20.82621 9.00000 1.71700 47.92
22.9 -1427.55747 1.32000 1.71736 29.51
16.2 13.55312 6.45000 1.00000 Inf
15.1 Inf 6.45000 1.00000 Inf
15.1 -18.95597 1.32000 1.71736 29.51
19.4 86.65094 9.00000 1.71700 47.92
25.6 -24.97014 0.10000 1.00000 Inf
33.0 70.76937 5.50000 1.80610 40.9
33.7 -86.13304 0.00000 1.00000 Inf
];
% 単凸レンズ
FFo(6)=2;
Lnamex{6}=['単凸レンズ 51.6mm'];
Surx{6}=[
25.0 50.0000 3.58140 1.60000 Inf
25.0 -79.1187 0.00000 1.00000 Inf
];
% 単凹レンズ
FFo(7)=2;
Lnamex{7}=['単凹レンズ 51.6mm'];
Surx{7}=[
25.0 -50.0000 1.00000 1.60000 Inf
25.0 81.9123 0.00000 1.00000 Inf
];
for nc=1:7 % ７種類のレンズについて繰り返す
disp(Lnamex{nc})
% レンズデータをQiita記事の本文の添字に合わせて変換
Sur=Surx{nc}(:,[2 3 4]);
k=size(Sur,1); % 面数
r=Sur(:,1); % 曲率半径 r_j (j=1～k)
d=[NaN ; Sur([1:end-1],2)]; % 境界面の間隔 d_j (j=1～k, d_1はNaN)
n=[1.00 ; Sur([1:end],3)]; % 屈折率 n_j (j=1～k+1)
% (j=1, k+1 はレンズの前後の空気層)
% ■■■■■ テスト用（常時はコメント化）■■■■■
% （レンズの前後の媒質が空気以外のときを想定）
%n=[1.20 ; Sur([1:end],3)]; % 前方の屈折率を 1.20 に置き換え
%n(end)=1.10; % 後方の屈折率を 1.10 に置き換え
% ■■ 流儀１ ■■
% ======================
% u と h を使う計算方法（流儀１）
% ======================
disp(' 【流儀１】uとhを使う計算方法')
% == 像側 posterior ==
u(1)=0;
h(1)=1.0;
for j=1:k
u(j+1)=n(j)*u(j)/n(j+1)+h(j)*(1-n(j)/n(j+1))/r(j);
if j~=k
h(j+1)=h(j)-d(j+1)*u(j+1);
end
end
fp=h(1)/u(k+1);
pp=(h(k)-h(1))/u(k+1); % 後面を基準に測った後側主点位置
% （基準より右側=正,左側=負）
ppx=sum(d(2:end))+pp; % 前面を 〃
% == 物体側 anterior ==
% レンズの左右反転に合わせて、データを並べ替え編集
rx=-flipud(r);
dx=[NaN ; flipud(d)];
dx=dx(1:end-1);
nx=flipud(n);
% ここで、j=1 は、反転前には最右だった面
u(1)=0;
h(1)=1.0;
for j=1:k
u(j+1)=nx(j)*u(j)/nx(j+1)+h(j)*(1-nx(j)/nx(j+1))/rx(j);
if j~=k
h(j+1)=h(j)-dx(j+1)*u(j+1);
end
end
fa=h(1)/u(k+1);
pa=-(h(k)-h(1))/u(k+1); % 前面を基準に測った前側主点位置
% （基準より右側=正,左側=負）
% 反転後に最右になった面から見た主点位置は、
% pax = (h(k)-h(1))/u(k+1) （これは pp と同じ式）
% 再反転させて元に戻し、本来の最左だった面から見た主点位置は、
% pa = -pax = -(h(k)-h(1))/u(k+1)
disp([' 前/後 焦点距離 ' ...
num2str(fa,'%7.4f') '/' num2str(fp,'%7.4f') ...
' 前/後 主点位置 ' ...
num2str(pa,'%7.4f') '/' num2str(pp,'%7.4f') ...
' (' num2str(ppx,'%7.4f') ')']);
% ======================
% s を使う計算方法（流儀１）
% ======================
disp(' 【流儀１】sを使う計算方法')
% == 像側 posterior ==
s(1)=Inf;
for j=1:k
sd(j)=n(j)/(n(j+1)*s(j))+(1-n(j)/n(j+1))/r(j);
sd(j)=1/sd(j);
if j~=k
s(j+1)=sd(j)-d(j+1);
end
end
fp=prod(sd(1:k))/prod(s(2:k));
pp=sd(k)-fp;
ppx=sum(d(2:end))+pp;
s=[];
sd=[];
% == 物体側 anterior ==
% レンズの左右反転に合わせて、データを並べ替え編集
rx=-flipud(r);
dx=[NaN ; flipud(d)];
dx=dx(1:end-1);
nx=flipud(n);
% ここで、j=1 は、反転前には最右だった面
s(1)=Inf;
for j=1:k
sd(j)=nx(j)/(nx(j+1)*s(j))+(1-nx(j)/nx(j+1))/rx(j);
sd(j)=1/sd(j);
if j~=k
s(j+1)=sd(j)-dx(j+1);
end
end
fa=prod(sd(1:k))/prod(s(2:k));
pa=-sd(k)+fa;
disp([' 前/後 焦点距離 ' ...
num2str(fa,'%7.4f') '/' num2str(fp,'%7.4f') ...
' 前/後 主点位置 ' ...
num2str(pa,'%7.4f') '/' num2str(pp,'%7.4f') ...
' (' num2str(ppx,'%7.4f') ')']);
% ======================
% nu と h を使う計算方法（流儀１）
% ======================
disp(' 【流儀１】nuとhを使う計算方法')
% == 像側 posterior ==
nu(1)=0;
h(1)=1.0;
for j=1:k
nu(j+1)=nu(j)+h(j)*(n(j+1)-n(j))/r(j);
if j~=k
h(j+1)=h(j)-d(j+1)*nu(j+1)/n(j+1);
end
end
fp=h(1)/(nu(k+1)/n(k+1));
pp=(h(k)-h(1))/(nu(k+1)/n(k+1));
ppx=sum(d(2:end))+pp;
% == 物体側 anterior ==
% レンズの左右反転に合わせて、データを並べ替え編集
rx=-flipud(r);
dx=[NaN ; flipud(d)];
dx=dx(1:end-1);
nx=flipud(n);
% ここで、j=1 は、反転前には最右だった面
nu(1)=0;
h(1)=1.0;
for j=1:k
nu(j+1)=nu(j)+h(j)*(nx(j+1)-nx(j))/rx(j);
if j~=k
h(j+1)=h(j)-dx(j+1)*nu(j+1)/nx(j+1);
end
end
fa=h(1)/(nu(k+1)/nx(k+1));
pa=-(h(k)-h(1))/(nu(k+1)/nx(k+1));
disp([' 前/後 焦点距離 ' ...
num2str(fa,'%7.4f') '/' num2str(fp,'%7.4f') ...
' 前/後 主点位置 ' ...
num2str(pa,'%7.4f') '/' num2str(pp,'%7.4f') ...
' (' num2str(ppx,'%7.4f') ')']);
% ======================
% 光線行列解析 h,u 方式（流儀１）
% ======================
disp(' 【流儀１】光線行列解析 h,u 方式')
M=1;
for j=1:k
M=[1 0 ; (1-n(j)/n(j+1))/r(j) n(j)/n(j+1)]*M;
if j~=k
M=[1 -d(j+1) ; 0 1]*M;
end
end
A=M(1,1);
B=M(1,2);
C=M(2,1);
D=M(2,2);
fp=1/C;
pp=(A-1)/C;
ppx=sum(d(2:end))+pp;
fa=(n(1)/n(k+1))/C;
pa=(n(1)/n(k+1)-D)/C;
disp([' 前/後 焦点距離 ' ...
num2str(fa,'%7.4f') '/' num2str(fp,'%7.4f') ...
' 前/後 主点位置 ' ...
num2str(pa,'%7.4f') '/' num2str(pp,'%7.4f') ...
' (' num2str(ppx,'%7.4f') ')']);
% ABCD行列の表示
disp([[' A B = ' ; ' C D '] num2str([A B;C D])])
% ======================
% 光線行列解析 h,nu 方式（流儀１）
% ======================
disp(' 【流儀１】光線行列解析 h,nu 方式')
M=1;
for j=1:k
M=[1 0 ; (n(j+1)-n(j))/r(j) 1]*M;
if j~=k
M=[1 -d(j+1)/n(j+1) ; 0 1]*M;
end
end
A=M(1,1);
B=M(1,2);
C=M(2,1);
D=M(2,2);
fp=n(k+1)/C;
pp=n(k+1)*(A-1)/C;
ppx=sum(d(2:end))+pp;
fa=n(1)/C;
pa=n(1)*(1-D)/C;
disp([' 前/後 焦点距離 ' ...
num2str(fa,'%7.4f') '/' num2str(fp,'%7.4f') ...
' 前/後 主点位置 ' ...
num2str(pa,'%7.4f') '/' num2str(pp,'%7.4f') ...
' (' num2str(ppx,'%7.4f') ')']);
% ABCD行列の表示
disp([[' A B = ' ; ' C D '] num2str([A B;C D])])
% ■■ 流儀２ ■■
% ======================
% u と h を使う計算方法（流儀２）
% ======================
disp(' 【流儀２】uとhを使う計算方法')
% == 像側 posterior ==
u(1)=0;
h(1)=1.0;
for j=1:k
u(j+1)=n(j)*u(j)/n(j+1)-h(j)*(1-n(j)/n(j+1))/r(j);
if j~=k
h(j+1)=h(j)+d(j+1)*u(j+1);
end
end
fp=-h(1)/u(k+1);
pp=-(h(k)-h(1))/u(k+1);
ppx=sum(d(2:end))+pp;
% == 物体側 anterior ==
% レンズの左右反転に合わせて、データを並べ替え編集
rx=-flipud(r);
dx=[NaN ; flipud(d)];
dx=dx(1:end-1);
nx=flipud(n);
% ここで、j=1 は、反転前には最右だった面
u(1)=0;
h(1)=1.0;
for j=1:k
u(j+1)=nx(j)*u(j)/nx(j+1)-h(j)*(1-nx(j)/nx(j+1))/rx(j);
if j~=k
h(j+1)=h(j)+dx(j+1)*u(j+1);
end
end
fa=-h(1)/u(k+1);
pa=(h(k)-h(1))/u(k+1);
% 反転後に最右になった面から見た主点位置は、
% pax = -(h(k)-h(1))/u(k+1) （これは pp と同じ式）
% 再反転させて元に戻し、本来の最左だった面から見た主点位置は、
% pa = -pax = (h(k)-h(1))/u(k+1)
disp([' 前/後 焦点距離 ' ...
num2str(fa,'%7.4f') '/' num2str(fp,'%7.4f') ...
' 前/後 主点位置 ' ...
num2str(pa,'%7.4f') '/' num2str(pp,'%7.4f') ...
' (' num2str(ppx,'%7.4f') ')']);
% ======================
% s を使う計算方法（流儀２）
% ======================
disp(' 【流儀２】sを使う計算方法')
% == 像側 posterior ==
s(1)=Inf;
for j=1:k
sd(j)=n(j)/(n(j+1)*s(j))+(1-n(j)/n(j+1))/r(j);
sd(j)=1/sd(j);
if j~=k
s(j+1)=sd(j)-d(j+1);
end
end
fp=prod(sd(1:k))/prod(s(2:k));
pp=sd(k)-fp;
ppx=sum(d(2:end))+pp;
s=[];
sd=[];
% == 物体側 anterior ==
% レンズの左右反転に合わせて、データを並べ替え編集
rx=-flipud(r);
dx=[NaN ; flipud(d)];
dx=dx(1:end-1);
nx=flipud(n);
% ここで、j=1 は、反転前には最右だった面
s(1)=Inf;
for j=1:k
sd(j)=nx(j)/(nx(j+1)*s(j))+(1-nx(j)/nx(j+1))/rx(j);
sd(j)=1/sd(j);
if j~=k
s(j+1)=sd(j)-dx(j+1);
end
end
fa=prod(sd(1:k))/prod(s(2:k));
pa=-sd(k)+fa;
disp([' 前/後 焦点距離 ' ...
num2str(fa,'%7.4f') '/' num2str(fp,'%7.4f') ...
' 前/後 主点位置 ' ...
num2str(pa,'%7.4f') '/' num2str(pp,'%7.4f') ...
' (' num2str(ppx,'%7.4f') ')']);
% ======================
% nu と h を使う計算方法（流儀２）
% ======================
disp(' 【流儀２】nuとhを使う計算方法')
% == 像側 posterior ==
nu(1)=0;
h(1)=1.0;
for j=1:k
nu(j+1)=nu(j)-h(j)*(n(j+1)-n(j))/r(j);
if j~=k
h(j+1)=h(j)+d(j+1)*nu(j+1)/n(j+1);
end
end
fp=-h(1)/(nu(k+1)/n(k+1));
pp=-(h(k)-h(1))/(nu(k+1)/n(k+1));
ppx=sum(d(2:end))+pp;
% == 物体側 anterior ==
% レンズの左右反転に合わせて、データを並べ替え編集
rx=-flipud(r);
dx=[NaN ; flipud(d)];
dx=dx(1:end-1);
nx=flipud(n);
% ここで、j=1 は、反転前には最右だった面
nu(1)=0;
h(1)=1.0;
for j=1:k
nu(j+1)=nu(j)-h(j)*(nx(j+1)-nx(j))/rx(j);
if j~=k
h(j+1)=h(j)+dx(j+1)*nu(j+1)/nx(j+1);
end
end
fa=-h(1)/(nu(k+1)/nx(k+1));
pa=(h(k)-h(1))/(nu(k+1)/nx(k+1));
disp([' 前/後 焦点距離 ' ...
num2str(fa,'%7.4f') '/' num2str(fp,'%7.4f') ...
' 前/後 主点位置 ' ...
num2str(pa,'%7.4f') '/' num2str(pp,'%7.4f') ...
' (' num2str(ppx,'%7.4f') ')']);
% ======================
% 光線行列解析 h,u 方式（流儀２）
% ======================
disp(' 【流儀２】光線行列解析 h,u 方式')
M=1;
for j=1:k
M=[1 0 ; -(1-n(j)/n(j+1))/r(j) n(j)/n(j+1)]*M;
if j~=k
M=[1 d(j+1) ; 0 1]*M;
end
end
A=M(1,1);
B=M(1,2);
C=M(2,1);
D=M(2,2);
fp=-1/C;
pp=-(A-1)/C;
ppx=sum(d(2:end))-(A-1)/C;
fa=-(n(1)/n(k+1))/C;
pa=-(n(1)/n(k+1)-D)/C;
disp([' 前/後 焦点距離 ' ...
num2str(fa,'%7.4f') '/' num2str(fp,'%7.4f') ...
' 前/後 主点位置 ' ...
num2str(pa,'%7.4f') '/' num2str(pp,'%7.4f') ...
' (' num2str(ppx,'%7.4f') ')']);
% ABCD行列の表示
disp([[' A B = ' ; ' C D '] num2str([A B;C D])])
% ======================
% 光線行列解析 h,nu 方式（流儀２）
% ======================
disp(' 【流儀２】光線行列解析 h,nu 方式')
M=1;
for j=1:k
M=[1 0 ; -(n(j+1)-n(j))/r(j) 1]*M;
if j~=k
M=[1 d(j+1)/n(j+1) ; 0 1]*M;
end
end
A=M(1,1);
B=M(1,2);
C=M(2,1);
D=M(2,2);
fp=-n(k+1)/C;
pp=-n(k+1)*(A-1)/C;
ppx=sum(d(2:end))+pp;
fa=-n(1)/C;
pa=-n(1)*(1-D)/C;
disp([' 前/後 焦点距離 ' ...
num2str(fa,'%7.4f') '/' num2str(fp,'%7.4f') ...
' 前/後 主点位置 ' ...
num2str(pa,'%7.4f') '/' num2str(pp,'%7.4f') ...
' (' num2str(ppx,'%7.4f') ')']);
% ABCD行列の表示
disp([[' A B = ' ; ' C D '] num2str([A B;C D])])
end
６.３ 近軸近似計算と光路シミュレーションの比較
この記事に掲載している画像は、全て MATLAB で作成しました。この中から、図７～図１３の作図用のプログラムのみを添付します。なお、画像ファイル出力の print 部分はコメント化しています。必要に応じてコメントを外してください。また、光路シミュレーションに汎用的に利用できそうなローカル関数もいくつか含んでいます。
プログラム（折り畳んであります）
% Lens57.m
% 近軸近似の結果 と 光路シミュレーション の合成
% 図７～図１３の作成用
clear
close all
% 標準線色
ao=[0 0.4470 0.7410]; % '#0072BD'
aka=[0.8500 0.3250 0.0980]; % '#D95319'
daidai=[0.9290 0.6940 0.1250]; % '#EDB120'
murasaki=[0.4940 0.1840 0.5560]; % '#7E2F8E'
midori=[0.4660 0.6740 0.1880]; % '#77AC30'
sorairo=[0.3010 0.7450 0.9330]; % '#4DBEEE'
chairo=[0.6350 0.0780 0.1840]; % '#A2142F'
haiiro=[0.5 0.5 0.5];
% レンズデータ
% 有効径 曲率半径 肉厚 屈折率 アッベ数
% No.058
FFo(1)=2.8; % 開放F値
Lnamex{1}=['No.58 Tessar型 公称50mm, F2.8'];
Surx{1}=[
19.1 19.48717 4.50000 1.77250 49.60
17.3 91.06909 4.00000 1.00000 Inf
15.2 -62.93988 2.00000 1.72825 28.46
14.1 17.94764 2.00000 1.00000 Inf
14.2 Inf 2.00000 1.00000 Inf
14.6 141.35119 1.00000 1.59270 35.31
14.7 24.49046 4.00000 1.80440 39.59
14.7 -37.36912 0.00000 1.00000 Inf
];
% No.072
FFo(2)=2; % 開放F値
Lnamex{2}=['No.72 Ernostar型 公称50mm, F2'];
Surx{2}=[
29.0 22.09251 6.00000 1.65844 50.88
27.5 123.19285 0.10000 1.00000 Inf
22.4 21.27789 4.00000 1.63854 55.38
20.1 32.11702 2.50000 1.00000 Inf
19.2 -424.31820 1.20000 1.76182 26.52
17.0 15.47789 4.00000 1.00000 Inf
16.9 Inf 4.96104 1.00000 Inf
16.6 58.13366 2.50000 1.72342 37.93
16.4 -52.72323 0.00000 1.00000 Inf
];
% No.084
FFo(3)=2; % 開放F値
Lnamex{3}=['No.84 Sonnar型 公称50mm, F2'];
Surx{3}=[
28.7 32.73326 3.80000 1.69680 55.53
27.4 97.63145 0.19000 1.00000 Inf
23.5 17.57712 4.50000 1.69680 55.53
22.0 36.68896 3.30000 1.48749 70.21
20.5 2255.56502 1.00000 1.64769 33.79
16.8 12.10293 4.50000 1.00000 Inf
16.5 Inf 1.20000 1.00000 Inf
16.2 -108.10542 6.00000 1.53172 48.84
15.3 27.87691 10.00000 1.77250 49.60
19.8 -64.34508 0.00000 1.00000 Inf
];
% No.093
FFo(4)=2; % 開放F値
Lnamex{4}=['No.93 Gauss型 公称50mm, F2'];
Surx{4}=[
32.4 28.45336 5.00000 1.65844 50.88
31.0 81.93668 0.10000 1.00000 Inf
24.9 17.86823 6.50000 1.64850 53.02
22.0 244.84648 1.30000 1.64769 33.79
16.9 11.87418 6.50000 1.00000 Inf
15.9 Inf 6.50000 1.00000 Inf
14.6 -13.85226 1.00000 1.64769 33.79
17.1 -49.52634 4.00000 1.64850 53.02
19.8 -19.13128 0.10000 1.00000 Inf
23.9 667.74212 4.50000 1.65844 50.88
25.0 -29.07016 0.00000 1.00000 Inf
];
% No.113
FFo(5)=2; % 開放F値
Lnamex{5}=['No.113 Gauss型 公称50mm, F2'];
Surx{5}=[
38.5 41.22676 5.50000 1.80610 40.92
36.7 141.20948 0.10000 1.00000 Inf
28.7 20.82621 9.00000 1.71700 47.92
22.9 -1427.55747 1.32000 1.71736 29.51
16.2 13.55312 6.45000 1.00000 Inf
15.1 Inf 6.45000 1.00000 Inf
15.1 -18.95597 1.32000 1.71736 29.51
19.4 86.65094 9.00000 1.71700 47.92
25.6 -24.97014 0.10000 1.00000 Inf
33.0 70.76937 5.50000 1.80610 40.9
33.7 -86.13304 0.00000 1.00000 Inf
];
% 単凸レンズ
FFo(6)=2;
Lnamex{6}=['単凸レンズ 51.6mm'];
Surx{6}=[
25.0 50.0000 3.58140 1.60000 Inf
25.0 -79.1187 0.00000 1.00000 Inf
];
% 単凹レンズ
FFo(7)=2;
Lnamex{7}=['単凹レンズ 51.6mm'];
Surx{7}=[
25.0 -50.0000 1.00000 1.60000 Inf
25.0 81.9123 0.00000 1.00000 Inf
];
% 光源（z,y平面内で光軸(z軸)に平行な光線群）
Piy=linspace(-18,18,53); % 光線の本数分の要素数の行ベクトル
Pi0=[0*Piy' Piy' 0*Piy'-10]; % 個々の光線の始点の位置座標
Vi0=[0*Piy' 0*Piy' 0*Piy'+1]; % 個々の光線の進行方向ベクトル
for nc=1:7 % 全てのレンズについて繰り返す
Sur=Surx{nc}; % レンズデータの取り込み
% ローカル関数の呼び出し（近軸近似での諸量の計算）
% fa: 前側焦点距離、 pa: 前側主点位置
% fp: 後側焦点距離、 pp: 後側主点位置（後面基準）
% ppx: 後側主点位置（前面基準） Lt: レンズ長
[fa pa fp pp ppx Lt] = Lens_spec_calc(Sur);
% 作図位置の調整
Zbias=-Lt/2; % 横方向の 0 位置はレンズ長の中心とする
Ybias=[20 -20]; % [前側光路図,後側光路図]の縦方向の基準位置
% レンズデータをQiita記事の本文の添字に合わせて変換
kk=size(Sur,1); % 面数
DD=Sur(:,1); % 口径
rr=Sur(:,2); % 曲率半径 r_j (j=1～k)
dd=[NaN ; Sur([1:end-1],3)]; % 境界面の間隔 d_j (j=1～k, d_1はNaN)
nn=[1.00 ; Sur([1:end],4)]; % 屈折率 n_j (j=1～k+1)
% (j=1, k+1 はレンズの前後の空気層)
for dir=1:2 % dir=1: 後側特性、2: 前側特性
if dir==2 % レンズの前後を反転させたとき
% レンズの左右反転に合わせて、データを並べ替え編集
DD=flipud(DD);
rr=-flipud(rr);
dd=[NaN ; flipud(dd)];
dd=dd(1:end-1);
nn=flipud(nn);
% ここで、j=1 は、反転前には最右だった面
end
Zc=[0;cumsum(dd([2:end]))]; % 各境界面と光軸との交点のz座標
% 要素数が面数の列ベクトル
Pi=Pi0; % 光路の始点の初期値
Vi=Vi0; % 光の進行方向の初期値
% 全ての交点位置を記録する三次元行列の準備
% 次元 ＝ 光線番号（行番号）
% x,y,z成分（列番号）
% 境界面番号（層番号）
% 層数 kk+3 ＝
% 光源始点1 + 境界面数kk
% + 最終投射面1 + データ区切り用NaN面1
Points=ones(size(Pi0,1),3,kk+3);
Points(:,:,1)=Pi0; % 光源の始点の位置情報を記録
for n=1:kk
% ローカル関数の呼び出し（入射点 Po と出射方向 Vo の計算）
[Po,Vo,~] = Lens_path02 (Zc(n),DD(n),rr(n),nn(n),nn(n+1),Pi,Vi);
Points(:,:,n+1)=Po; % 境界面との交点の位置情報を記録
Pi=Po;
Vi=Vo;
end
% 最終投影面との交点
if dir==1 % レンズの前後の配置は通常状態
if fp>0 % 凸レンズ
Ze=ppx+fp+20; % 最終投影面の位置
else % 凹レンズ
Ze=Lt+20;
end
else % レンズの前後を反転させた状態
if fp>0
Ze=fa-pa+20+Lt;
else
Ze=Lt+20;
end
end
% ローカル関数の呼び出し（入射点 Po と出射方向 Vo の計算）
[Po,Vo,~] = Lens_path02 (Ze,50,Inf,1,1,Pi,Vi);
Points(:,:,kk+2)=Po;
% 光線ごとのデータ区切り用の NaN 面
Points(:,:,kk+3)=Pi0*NaN;
% 多次元行列の次元の入れ替え。
% 入れ替え後の次元 = x,y,z成分（行番号）
% 境界面番号（列番号）
% 光線番号（層番号）
Q=permute(Points,[2 3 1]);
% 全光路を NaN を介して１本に繋げる
Dx=[];
Dy=[];
Dz=[];
for n=1:size(Q,3)
Dx=[Dx Q(1,:,n)];
Dy=[Dy Q(2,:,n)];
Dz=[Dz Q(3,:,n)];
end
% figure の大き目の画面を準備
hf=figure(nc);
hf.Position=[hf.Position(1) hf.Position(2)-70 ...
hf.Position(3)*1.2 hf.Position(4)*1.2];
hf.Color='w';
plot([-80 60],[1 1]*Ybias(dir),'Color','#aaa'); % 光軸の描画
hold on
% 光束の描画
if fp>0 % 凸レンズのとき
% 実在の光束の描画（濃色）
if dir==1 % レンズの前後の配置は通常状態
hL1=plot(Dz+Zbias,Dy+Ybias(dir),'Color',ao);
else % レンズの前後を反転させた状態
hL2=plot(-Dz+Lt+Zbias,Dy+Ybias(dir),'Color',ao);
end
else % 凹レンズのとき
% 実在の光束の描画（淡色）
if dir==1 % レンズの前後の配置は通常状態
plot(Dz+Zbias,Dy+Ybias(dir),'Color','#aaa');
else % レンズの前後を反転させた状態
plot(-Dz+Lt+Zbias,Dy+Ybias(dir),'Color','#aaa');
end
% === 虚像を作る仮想光束の描画 ===
% 虚像の投影面のz座標
if dir==1 % レンズの前後の配置は通常状態
Ze=ppx+fp-20;
else % レンズの前後を反転させた状態
Ze=fa-pa-20+Lt;
end
% 虚像を作る仮想光束の計算
[x,y]=intersection_of_lines02 ...
([Points(:,1,1)*0+Ze Points(:,1,1)*0-50], ...
[Points(:,1,1)*0+Ze Points(:,1,1)*0+50], ...
Points(:,[3 2],end-2), Points(:,[3 2],end-1));
% 全光線を NaN を介して１本に繋げる
Ex=[];
Ey=[];
Ez=[];
for n=1:size(Points,1)
Ex=[Ex 0 0 NaN];
Ey=[Ey Points(n,2,end-2) y(n) NaN];
Ez=[Ez Points(n,3,end-2) x(n) NaN];
end
% 仮想光束の描画（濃色）
if dir==1 % レンズの前後の配置は通常状態
hL3=plot(Ez+Zbias,Ey+Ybias(dir),'Color',ao);
else % レンズの前後を反転させた状態
hL4=plot(-Ez+Lt+Zbias,Ey+Ybias(dir),'Color',ao);
end
end % end of "if fp>0"
axis equal
axis([-80 60 -42 48])
grid on
% ローカル関数の呼び出し（レンズの断面図の描画）
draw_cross_section(Sur,1,[Zbias Ybias(dir)])
% 主面の計算と描画
% ローカル関数の呼び出し（入射光・出射光の延長線の交点の検出）
[x,y]=intersection_of_lines02 ...
(Points(:,[3 2],1), Points(:,[3 2],2), ...
Points(:,[3 2],end-2), Points(:,[3 2],end-1));
if dir==1 % レンズの前後の配置は通常状態
hp0=plot(x+Zbias,y+Ybias(dir),'k');
else % レンズの前後を反転させた状態
plot(-x+Lt+Zbias,y+Ybias(dir),'k');
end
% 近軸近似での焦点と主点の位置を描画
hp1=plot(pa+Zbias,0+Ybias(dir),'.','MarkerSize',15, ...
'Color','#bb0'); % 前側主点
hp2=plot(ppx+Zbias,0+Ybias(dir),'.','MarkerSize',15, ...
'Color','#0aa'); % 後側主点
hp3=plot(pa-fa+Zbias,0+Ybias(dir),'.','MarkerSize',15, ...
'Color','#f0f'); % 前側焦点
hp4=plot(ppx+fp+Zbias,0+Ybias(dir),'.','MarkerSize',15, ...
'Color','#f00'); % 後側焦点
% レンズ名の表示
text(-75,43,Lnamex{nc},'FontSize',12,'FontName','メイリオ');
haxmain = gca; % 主座標のハンドルを取得
% 補助座標に拡大図を描画
Wz=30; % 補助座標の幅
Hz=15; % 補助座標の高さ
ka=100; % 拡大倍率（拡大座標の横幅は Wz/ka [mm] 相当となる）
if dir==1
Xf=hp4.XData;
Yf=hp4.YData;
if nc~=7
[ax1] = add_sub_axes(haxmain,[28 -3 Wz Hz]);
copyobj([hp4 hL1],ax1);
else
[ax1] = add_sub_axes(haxmain,[-78 -2 Wz Hz]);
copyobj([hp4 hL3],ax1);
end
else
Xf=hp3.XData;
Yf=hp3.YData;
if nc~=7
[ax2] = add_sub_axes(haxmain,[-78 -12 Wz Hz]);
copyobj([hp3 hL2],ax2);
else
[ax2] = add_sub_axes(haxmain,[28 -12 Wz Hz]);
copyobj([hp3 hL4],ax2);
end
end
% 焦点の位置を補助座標の中心に置いて拡大する
axis([Xf-Wz/(2*ka) Xf+Wz/(2*ka) Yf-Hz/(2*ka) Yf+Hz/(2*ka)]);
% 補助座標の白色背景
hh=axis
hpatch = patch([hh(1) hh(2) hh(2) hh(1)], ...
[hh(3) hh(3) hh(4) hh(4)],'w','EdgeColor','k');
uistack(hpatch,'bottom');
axis off
axes(haxmain); % アクティブ座標をサブからメインに切り替え
end % end of "for dir=1:2"
axes(ax1); % 隠れている拡大図を前面に出す
axes(ax2);
% ローカル関数の呼び出し（figure のトリミング）
% （legend 付きの figure ではエラーになるので、legend 前に実行）
trim_and_expand_figure_margin(hf,[-85 -85 -70 -70]);
% 凡例の表示
% Position属性を使うときは、legendのハンドルを作り、しかも
% その内容を表示させないと位置や枠の大きさが思いどおりにならない。
% ■■■ MATLAB のバグか？ ■■■
if nc~=7
hleg=legend([hp3 hp4 hp1 hp2 hp0], ...
{'前側焦点','後側焦点','前側主点','後側主点','主面'}, ...
'FontSize',12,'Position',[0.09 0.65 0.2 0.2])
else
hleg=legend([hp3 hp4 hp1 hp2 hp0], ...
{'前側焦点','後側焦点','前側主点','後側主点','主面'}, ...
'FontSize',12,'Position',[0.09 0.15 0.2 0.2])
end
axes(ax1); % 隠れている拡大図を前面に出す
axes(ax2);
% figureの画像ファイル出力
basename=mfilename('fullpath'); % このプログラムのファイル名
% print(gcf,'-dpng',[basename '_0' num2str(nc) '.png'],'-r600')
end % end of "for nc=1:7"
% ==============================================================
% ==============================================================
% 以下、ローカル関数群
% ==============================================================
% ==============================================================
% ====================================
% レンズの主点と焦点の計算（流儀１の光線行列解析 h,u 方式による）
% ====================================
function [fa pa fp pp ppx Lt] = Lens_spec_calc(Sur)
% 【入力】
% Sur: レンズデータの行列
% 境界面の数だけの行数 × ５列
% ５列=[口径,曲率半径,肉厚,屈折率,アッベ数]
%
% 【出力】
% fa: 前側焦点距離
% pa: 前側主点（最前面を基準にして、前側は-、後側は+）
% fp: 後側焦点距離
% pp: 後側主点（最後面を基準にして、前側は-、後側は+）
% ppx: 〃 （最前面 〃 ）
% Lt: 光軸上のレンズ長
k=size(Sur,1); % 面数
r=Sur(:,2); % 曲率半径 r_j (j=1～k)
d=[NaN ; Sur([1:end-1],3)]; % 境界面の間隔 d_j (j=1～k, d_1はNaN)
n=[1.00 ; Sur([1:end],4)]; % 屈折率 n_j (j=1～k+1)
% (j=1, k+1 はレンズの前後の空気層)
M=1;
for j=1:k
M=[1 0 ; (1-n(j)/n(j+1))/r(j) n(j)/n(j+1)]*M;
if j~=k
M=[1 -d(j+1) ; 0 1]*M;
end
end
A=M(1,1);
B=M(1,2);
C=M(2,1);
D=M(2,2);
fp=1/C;
pp=(A-1)/C;
Lt=sum(d(2:end));
ppx=Lt+pp;
fa=(n(1)/n(k+1))/C;
pa=(n(1)/n(k+1)-D)/C;
end
% ====================================
% 平面上の2直線（含む:延長線）の交点の座標の取得
% ====================================
function [x,y]=intersection_of_lines02(xy1a,xy1b,xy2a,xy2b)
% 平面上の2直線（含む:延長線）の交点の座標の取得（多重一括処理）。
% 2本の線分を1組として、それが複数組あるとき、
% 全組の交点を一度の処理で取得する。
%
% 【入力】
% xy1a: 1本目の直線のa端のx,y座標を [x1a y1a] とし、
% これを縦方向に組数分だけ並べた、多行2列の行列
% xy1b: 同じく、b端の [x1b y1b] を上記と同様に並べた行列
% xy2a: 2本目の a端の [x2a y2a] 〃
% xy2b: 同じく、b端の [x2b y2b] 〃
% 【出力】
% x,y: 交点の座標。それぞれが組数分の長さの列ベクトル
Nc=size(xy1a,1);
if (size(xy1b,1)~=Nc | size(xy2a,1)~=Nc | size(xy2b,1)~=Nc | ...
size(xy1a,2)~=2 | size(xy1b,2)~=2 | size(xy2a,2)~=2 | ...
size(xy2b,2)~=2)
error('引数の形式や値が不正です。');
end
% 入力された x,y の２次元表現を x,y,z の３次元表現に変える。
% ３要素にしないとcrossコマンドが使えないため。z 成分は 0 に固定。
P1=[xy1a 0*ones(Nc,1)]; % 1本目の線の両端座標
P2=[xy1b 0*ones(Nc,1)];
P3=[xy2a 0*ones(Nc,1)]; % 2本目の線の両端座標
P4=[xy2b 0*ones(Nc,1)];
Q=cross(P2-P1,P4-P3,2);
D = Q(:,1).^2 + Q(:,2).^2 + Q(:,3).^2;
D(D<eps(1))=NaN; % 平行線同士で交点が無い組はNaN
cA=cross(P4-P3,P3-P1,2);
cB=cross(P4-P3,P2-P1,2);
k=cA(:,3).*cB(:,3)./D;
P21=P2-P1;
XY=P1+[k.*P21(:,1) k.*P21(:,2) k.*P21(:,3)];
x=XY(:,1);
y=XY(:,2);
end
% ==================================================
% レンズの断面図の描画
% ==================================================
function draw_cross_section(Sur,Fc,Bxy)
% 注： 計算部ではレンズ光学の慣例に合わせて、横軸/縦軸を z/y と
% しているが、描画部では一般常識の x/y としているので、
% 混同なきよう要注意。
% 【入力】
% Sur: レンズデータの行列
% 境界面の数だけの行数 × ５列
% ５列=[口径,曲率半径,肉厚,屈折率,アッベ数]
% Fc: 設定絞り値/解放絞り値（絞りが解放のときは1.0）
% Bxy: 描画位置調整用の行ベクトル（[横 縦]方向への移動量）
% バイアスが[0 0]のとき、描画のための
% 縦軸原点は水平光軸上、
% 横軸原点はレンズ最前面と光軸の交点
%
% 【出力】
% なし
if (size(Sur,2)~=5 | Fc>1.0 | size(Bxy,1)~=1 | size(Bxy,2)~=2)
error('引数の形式や値が不正です。');
end
maxD=max(Sur(:,1)); % 最大口径
Zt=0; % 最前面の境界と光軸の交点のz座標
for nn=1:size(Sur,1) % 合成レンズの1面ずつを処理。
% 現面の有効径Dc
if nn>=2 & isinf(Sur(nn,2)) & ...
Sur(nn-1,4)==Sur(nn,4) % 現面が絞り面のとき
Dc=Sur(nn,1)*Fc;
else % 現面がレンズ面のとき
Dc=Sur(nn,1);
end
Rs=Sur(nn,2); % 現面の曲率半径
Zo=Rs+Zt; % 現面の曲率半径の中心点のz座標
% 描画準備
B=[-Dc/2,Dc/2]; % レンズの両端のx座標
if ~isinf(Rs) % 面が球面のとき
an=asin((Dc/2)/abs(Rs));
if Rs>0
A=linspace(pi+an,pi-an,100);
elseif Rs<0
A=linspace(-an,an,100);
else
error('Rsの誤指定です。')
end
yy1=abs(Rs)*sin(A); % 現面の断面の円弧のデータ
zz1=abs(Rs)*cos(A)+Zo;
else % 現面が平面のとき
yy1=B; % 平面の断面のデータ
zz1=B*0+Zt;
end
% 描画
if nn>=2
if Sur(nn-1,4)>=1.0003
% 現面より被写体側にある媒質がガラスのとき
hp=patch([zz1 fliplr(zz1_old)]+Bxy(1), ...
[yy1 fliplr(yy1_old)]+Bxy(2),'m','EdgeColor','r');
hp.FaceAlpha=0.1;
end
if isinf(Sur(nn,2)) & Sur(nn-1,4)==Sur(nn,4)
% 現面が絞り面のとき
plot([1 1 NaN 1 1]*Zt+Bxy(1), ...
[-1.2*maxD/2 yy1(1) NaN yy1(end) 1.2*maxD/2]+Bxy(2), ...
'Color',[0.5 0.5 0.5])
end
end
yy1_old=yy1;
zz1_old=zz1;
Zt=Zt+Sur(nn,3); % 次に処理する境界面と光軸の交点のz座標
end
end
% ==================================================
% 球面レンズの１つの境界面での屈折の計算
% ==================================================
function [Po,Vo,Err] = Lens_path02 (Zc,Dc,Rs,ni,no,Pi,Vi)
% 両側の媒体の屈折率が異なる球面状の境界面での屈折の計算。
% 具体的には・・・・
% 「入射光（直線）と境界面（球面）との交点の座標」の計算と
% 「屈折後の光の進路方向の単位ベクトル」の計算。
%
% レンズは２つの面を持つので、例えば、１枚のレンズの計算では、
% 裏と表でこの関数をそれぞれ１回ずつ呼び出す必要がある。
%
% フィルム側から被写体側を見たとき、
% x軸は右向き、y軸は上向き、z軸は被写体からフィルムへの向き。
% 光の進路は、z軸の負方向から正方向に進むものとする。
% 長さの単位は[mm]、位置や方向のベクトルは三次元。
%
% 【入力】
% Zc: 境界面と光軸との交点のz座標値。
% Dc: 境界面の口径
% Rs: 境界面の曲率半径。入射方向に凸のときは「+」、
% 出射方向に凸のときは「-」、平面のときは「Inf」。
% ni: 入射側の媒質の屈折率
% no: 出射側 〃
% Pi: 入射光（直線）の「始点の座標」の行ベクトルを、
% 光線の本数分だけ縦に並べた行列。
% Vi: 入射光の進行方向を示す行ベクトルを、光線の本数分だけ縦に
% 並べた行列。（内部で直ちに単位化するので、長さは任意）
%
% 【出力】
% Po: 入射光と境界面との「交点の座標」の行ベクトルが、
% 光線の本数分だけ縦に並んた行列。
% Vo: 屈折後の光の進行方向を示す単位行ベクトルが、
% 光線の本数分だけ縦に並んた行列。
% Err: エラーコード。光線の本数分の要素の列ベクトル。
% 0: エラーなし
% 1: 光線がレンズ球面の全球とさえも交差しない。
% 2: 交差はするが、交点がレンズの口径外に外れる。
% 3: ガラスから空気への境界面で全反射が起こる。
% ■■計算根拠
%
% ■球面（レンズ面）と直線（光路）の交点
% 球面の式 ||x-Cs||^2=Rs^2
% 直線の式 x=Pi+t*Vi
% x:三次元位置ベクトル、Cs:球の中心のベクトル、Rs:球の半径
% Pi:直線上の任意の一点のベクトル、Vi:光路方向の単位ベクトル
% t:媒介変数
% 以上より、||Pi+t*Vi-Cs||^2=Rs^2 となるtの値txを求めると、
% t^2 + 2Vi・(Pi-Cs)t + (Pi-Cs)・(Pi-Cs) - Rs^2 = 0 より、
% tx = -Vi・(Pi-Cs) ± sqrt( {Vi・(Pi-Cs)}^2 - ||Pi-Cs||^2 + Rs^2 )
% 交点PoはPi+tx*Viとなる。txの式の中の±は状況に応じて+/-を選択。
%
% ■屈折の計算
% 固定座標上での任意の向きの計算は複雑になるので、一旦、都合の良い
% 暫定座標上へ変換し、屈折の計算をした後に元の座標に戻す。
% ① 曲面上の交点を暫定原点とし、固定座標上で被写体側を向いた
% 単位法線ベクトルNorを暫定y軸に選び、
% Vy=[0 1 0]とおく。
% ② 上記法線と入射光線が作る平面を暫定xy平面とし、この平面を、
% 入射光が第１象限にくるような方向から観察して、暫定原点から
% 右向きに暫定x軸を設定する。
% 入射角をθiとすると、固定座標上で光路と逆向き
% のベクトル-Viは、暫定座標上では
% Vv=[cos(pi/2-θi) sin(pi/2-θi) 0]に置き換わる。
% ただし、θi=acos((-Vi)・Nor)
% ③ 固定座標上で表された光路面に垂直な単位ベクトル
% Vp={(-Vi)×Nor}/norm{(-Vi)×Nor}は
% 暫定z軸のVz=[0 0 1]に置き換わる。
% 上記の変換に必要な行列をMとすると、
% [Vv' Vy' Vz']=M*[-Vi' Nor' Vp']より、
% M=[Vv' Vy' Vz']*inv([-Vi' Nor' Vp'])
% 屈折後の暫定xy平面上での光の逆進路Voxは、
% Vox=[cos(pi/2-θo) sin(pi/2-θo) 0]
% ただし、θo=asin(ni*sinθi/no) （スネルの法則）
% これを元の座標に戻したうえ向きを反転すると、屈折後の光の進路Voは
% Vo=(-inv(M)*Vox')'
%
% 注： 上記の計算手順は１本の光線について記したもの。
% 実際には複数の光線を同時に処理するため、プログラム上の
% 計算式は少し異なる部分がある。
if (Dc<=0 | Rs==0 | ni<1 | no<1 ...
| size(Pi,2)~=3 | size(Vi,2)~=3 |size(Pi,1)~=size(Vi,1))
error('引数の形式や値が不正です。');
end
% 光の進行方向ベクトルを単位化
Lvi = sqrt(Vi(:,1).^2 + Vi(:,2).^2 + Vi(:,3).^2);
Vi = [Vi(:,1)./Lvi Vi(:,2)./Lvi Vi(:,3)./Lvi];
Err = zeros(size(Pi,1),1); % 光線ごとのエラー内容の記録用
if ~isinf(Rs) % 境界面が球面のとき
Cs = Rs + Zc; % 球の中心のz座標
Pc = [0 0 Cs]; % 球の中心点の座標
% 球面と直線の交点の有無を確認する判別式
Pic = Pi - Pc;
q = dot(Vi,Pic,2).^2 ...
- (Pic(:,1).^2 + Pic(:,2).^2 + Pic(:,3).^2) + Rs^2;
Err(q<0) = 1; % エラー1: 光線が球面と交差しない。
q(q<0) = NaN;
if Rs>0
tx = -dot(Vi,Pic,2) - sqrt(q); % 交点での媒介変数の値
else
tx = -dot(Vi,Pic,2) + sqrt(q);
end
% 境界面と光線の交点の座標
Po = Pi + [tx.*Vi(:,1) tx.*Vi(:,2) tx.*Vi(:,3)];
% その交点で、入射側を向く「球面の法線」のベクトル
Nor = sign(Rs)*(Po - Pc);
% これを単位ベクトル化
Lnor = sqrt(Nor(:,1).^2 + Nor(:,2).^2 + Nor(:,3).^2);
Nor = [Nor(:,1)./Lnor Nor(:,2)./Lnor Nor(:,3)./Lnor ];
else % 境界面が平面のとき
tx = (Zc - Pi(:,3))./Vi(:,3); % 交点での媒介変数の値
% 境界面と光線の交点の座標
Po = [tx.*Vi(:,1)+Pi(:,1) ...
tx.*Vi(:,2)+Pi(:,2) ...
Zc*ones(size(Pi,1),1)];
Nor = repmat([0 0 -1],size(Pi,1),1); % 入射側を向く単位法線
end
Err(Po(:,1).^2+Po(:,2).^2>(Dc/2)^2) = 2;
% エラー2: 交点がレンズの口径外に外れた
Cnv = cross(Nor,-Vi,2);
Ctir = sqrt(Cnv(:,1).^2 + Cnv(:,2).^2 + Cnv(:,3).^2); % 全反射条件
Err(Ctir>no/ni) = 3; % エラー3: ガラスから空気の面で全反射
% 屈折の計算（光線の本数だけ繰り返し処理）
Vo=ones(size(Pi,1),3)*NaN; % 出射光の進路方向の単位ベクトルの初期値
for n=1:size(Pi,1)
if Err(n)~=0
continue;
end
Norx=Nor(n,:);
Vix=Vi(n,:);
if norm(cross(Norx,-Vix,2))>1e-8
% 入射角が0ではないと見なせるとき
thi=acos(dot(Norx,-Vix)); % 入射角θi
Vy=[0 1 0]; % 暫定y軸方向の単位ベクトル
Vv=[cos(pi/2-thi) sin(pi/2-thi) 0]; % 暫定xy面での逆光路ベクトル
Vz=[0 0 1]; % 暫定z軸方向の単位ベクトル
Vp=cross(-Vix,Norx,2);
Vp=Vp/norm(Vp); % 元座標表現で光路面に垂直な単位ベクトル
M=[Vv' Vy' Vz']*inv([-Vix' Norx' Vp']);
% 座標変換用の回転行列
tho=asin(ni*sin(thi)/no); % 出射角θo
Vox = [cos(pi/2-tho) sin(pi/2-tho) 0];
Vox = (-inv(M)*Vox')'; % 出射光の進路方向の単位ベクトル
else
% 入射角がほとんど0のとき
Vox = Vix;
end
Vo(n,:)=Vox;
end % end of "for n=1:size(Pi,1)"
end % end of "function Lens_path02"
% ====================================
% figureの上下左右の余白調整
% ====================================
function trim_and_expand_figure_margin(hf,KK)
% 注： figure に legend が含まれているとエラーになる。
% legend は、この関数を実行後に作成すること。
% 【入力】
% hf: 余白調整したいfigureのハンドル
% KK: 余白の変更量から成る配列。 単位は[pixel]。
% （[上余白 下余白 左余白 右余白]。増量は+、減量は-）
% 【出力】
% 返り値は無し
posf=hf.Position; % 余白変更前のfigureのposition
% [左端 下端 幅 高さ]、単位は[pixel]
naxe=size(hf.Children,1); % figure内のChildren（axesやcolorbar）
% などの総数。
% 全Childrenのposition行列を作成
% 各行がそれぞれのChildrenに対応。
% 各列は[左端 下端 幅 高さ]、単位は[normalized]。
Posa=[]; % position行列の初期値
for n=1:naxe
axes(hf.Children(n))
hn(n)=gca;
Posa=[Posa ; hn(n).Position];
end
% 各Childrenの、余白変更前のfigure内でのpixel単位でのPositionを計算
Posp(:,[1 3])=Posa(:,[1 3])*posf(3);
Posp(:,[2 4])=Posa(:,[2 4])*posf(4);
% 余白変更後のfigureのpixel単位でのPositionを計算。
posfx=posf+[-KK(3) -KK(2) KK(3)+KK(4) KK(1)+KK(2)];
% 余白変更後のfigure内での各Childrenのpixel単位でのPositionを計算
Pospx(:,1)=Posp(:,1)+KK(3);
Pospx(:,2)=Posp(:,2)+KK(2);
Pospx(:,3)=Posp(:,3);
Pospx(:,4)=Posp(:,4);
% 余白変更後の各Childrenのnormalized単位でのPositionを計算する。
Posax(:,[1 3])=Pospx(:,[1 3])/posfx(3);
Posax(:,[2 4])=Pospx(:,[2 4])/posfx(4);
% figureとChildrenのpositionを更新
hf.Position=posfx;
for nn=1:naxe
hn(nn).Position=Posax(nn,:);
end
end
% ====================================
% 現在の座標上に補助的な座標軸を追加する
% ====================================
function [hax2] = add_sub_axes(hax1,Pos)
% 備忘： _px ...ピクセル単位、_nr ...normalized単位
% _ax ...hax1_axesの座標スケール単位
% 【入力】
% hax1: 現在の座標のハンドル
% Pos: ベクトル [X_ax Y_ax W_ax H_ax]
% （単位: hax1の座標スケール）
% ただし
% X_ax: 追加する補助座標面の左下端のx座標
% Y_ax: 〃 のy座標
% W_ax: 追加する補助座標面の幅
% H_ax: 〃 の高さ
% 【出力】
% hax2: 追加した補助座標面のハンドル
% 補助座標の位置と大きさ（単位: hax1の座標スケール）
X_ax=Pos(1); Y_ax=Pos(2); W_ax=Pos(3); H_ax=Pos(4);
Pf=hax1.Parent.Position; % モニタ画面内の figure の位置と大きさ
% （単位: pixel）
% [左下端x 左下端y 幅 高さ]
% Figure の大きさ（単位: pixel）
Wf_px=Pf(3); Hf_px=Pf(4);
Pa=hax1.Position; % figure内のhax1の位置
% （単位: normalized）
% [左下端x 左下端y 幅 高さ]
% 注： axis equal コマンドによって axes の実質的な高さや
% 幅が圧縮されてしまった後でも、上記コマンドで
% 得られる Pa の値は axis equal を指定する前の値に
% なる。（後の処理が面倒になって困る！！）
% figure 内の hax1 の位置と大きさ（単位: normalized）
Xa_nr=Pa(1); Ya_nr=Pa(2); Wa_nr=Pa(3); Ha_nr=Pa(4);
% hax1 の大きさ（単位: pixel）
Wa_px=Wf_px*Wa_nr; % 幅
Ha_px=Hf_px*Ha_nr; % 高さ
axes(hax1);
Pb=axis; % hax1_axes内のx,y座標軸の範囲
% [Xmin_ax Xmax_ax Ymin_ax Ymax_ax]
% hax1 の座標軸の範囲（単位: hax1の座標スケール）
Xmin_ax=Pb(1); Xmax_ax=Pb(2); Ymin_ax=Pb(3); Ymax_ax=Pb(4);
% アスペクト比
K0=Wa_px/Ha_px; % axis equal 前の hax1 のアスペクト比
Wa_ax = Xmax_ax - Xmin_ax;
Ha_ax = Ymax_ax - Ymin_ax;
Keq = Wa_ax/Ha_ax; % axis equal 後の hax1 の実質アスペクト比
K = Keq/K0; % axis equal 前後のアスペクト比の変化倍率
if K >= 1.0
X_nr = Wa_nr*(X_ax - Xmin_ax)/Wa_ax + Xa_nr;
Y_nr = (Ha_nr/K)*(Y_ax - Ymin_ax)/Ha_ax + Ya_nr + (1-1/K)*Ha_nr/2;
W_nr = Wa_nr*W_ax/Wa_ax;
H_nr = (Ha_nr/K)*H_ax/Ha_ax;
else
X_nr = Wa_nr*K*(X_ax - Xmin_ax)/Wa_ax + Xa_nr + (1-K)*Wa_nr/2;
Y_nr = Ha_nr*(Y_ax - Ymin_ax)/Ha_ax;
W_nr = Wa_nr*K*W_ax/Wa_ax;
H_nr = Ha_nr*H_ax/Ha_ax;
end
hax2=axes('Position',[X_nr Y_nr W_nr H_nr]);
end