統合趣味環境

その都度、気になったことを綴ります。

ゼロから始める八丈島旅行計画

なぜ八丈島なのか

非日常を，味わいたくはないか？

都会にいると，ビル群に囲まれて鬱になる．

なんか，こう，ヒト以前に生命体として大事な何かを失った気がする．

あーあ，何も考えなくていい場所に行きたい．

すると，おぼろげながら浮かんできたんです，島嶼部という語が...

それに，なんか学生っぽくてイイ．

近すぎるのも遠すぎるのもイヤだし，ここは八丈島で手を打とうじゃないか．

というわけで，

海！空！島！

みたいな理想の風景を思い描いている．

とはいえ曇りの多い島であるらしい．過度な期待はしないでおこう．

八丈島へ行こう

船の予約

飛行機なんて乗りません．船です．

絶海の孤島を味わうには，ぜひとも船で行くべきです．

東海汽船の客船に乗って，竹芝から10時間半くらいで行けるみたいです．

学割あります．私は電話で予約．

最初は自動音声ですが，途中からオペレータさんに切り替わります．

極貧学生の気分なので，二等の雑魚寝席で行こうと思います．

船旅の概要

竹芝 22:30 → 三宅島 5:00 → 御蔵島 6:00 → 八丈島 8:55

youtubeの八丈旅行記を見ていると，最初の寄港地である三宅島に到着する前の4:00くらいに，室内灯が全部つくらしい．まあ，安い席を買うからなのでしょうけど...

寝れるのだろうか．

それと，普通の人間は寝る前にご飯を食べて風呂に入ったりすると思うが，船内にもレストランがあるし，共用シャワーも完備だ．お金を積むと部屋にシャワーが付いたりする．

しかし出発も遅いので，夕飯は陸で食べることにし，風呂は家で浸かって来ようと思う．

実は前にサークル関係のイベントで，伊豆大島に行ったことがある．その際には竹芝ターミナル近くの「ゆで太郎」で済ませた記憶がある．

八丈島に泊まろう

民宿かホテルか，それが問題だ．

ここで私は地元の店を開拓したかったので，ホテル素泊まりを選択．

八丈島を車で回ろう

極貧学生ゆえ，最初は徒歩強行軍を模索...

うーん，思ったより島が広いぞ．

広いというか，2つある山が曲者．

道がグネグネしてるし起伏はありそうだし，車しかないな．

レンタカーを予約．

八丈島でご飯を食べよう

島寿司，というものがあるらしい．

これは，ヅケにされたタネ，そしてからしの和えられた特殊なシャリに特徴づけられる島の伝統料理らしい．

いろいろ店はあるけど，どこも電話予約が必須．

昼に寿司を食べ，夜は居酒屋に乗り込む，かな．

八丈島を観光しよう

ダイビング目当ての人も多いみたい．

でも私は泳ぎませんよ．それにまだ2月だし．

今のところ行った方が良さそうなポイントは，ここだ．

そして航空ファンとしては八丈島ORSRも見逃せない．本邦航空局はこれによって上空15,000ft以上の監視能力を獲得しているらしく，航空機同士の間隔を正確に知るためなどに利用される重要設備なのだ．

特に八丈島はSSRモードSの高機能レーダーに改修済みらしく，見ても分からないだろうが一見の価値があると思う．

ちなみにGoogleマップにすら場所が載っていない超極秘施設（笑）だが，航空写真で見る限り，えこ・あぐりまーとの近所だと判明している．

近所にRCAG施設もあると睨んでいる．

八丈島考

事前学習は欠かせない．

ここでは江戸時代の八丈島をみてみよう．

早速だが八丈島は歴史的にかなり明暗のある島に思える．

以下はネットの情報と，大学で借りた『江戸時代の八丈島――孤島苦の究明――』によるところが多い．

注意：以下，とがった内容が多いので折り畳みとした．

詳細はこちら

名称由来

ちまたで有力なのは八丈絹に因むという説だろう．

それでは絹の「八丈」というのに別の由来があるのではないかと疑うかもしれないが， 八丈はサイズのことであって，織物固有の名詞ではないから，「八丈」の由来は「八丈絹」ただこれに帰着される．

後に記すように，かつての島における黄八丈の重要性は筆舌に尽くしがたいほどであり，一見して安直なネーミングではあるものの，この島の名が意味する真の重みは計り知れないと感じた．

政治支配

室町のころには奥山宗麟が，のち小田原北条氏が支配したようだ．江戸にはもちろん徳川の支配下に置かれた．

しかし幕府が送ってくる島役人はカタチだけであり，実質的に支配者的存在だったのは，地役人こと島の有力者だった．

支配されれば年貢を納めねばならないが...八丈島は絶海の孤島で，作物もよく育たない．そこで，いずれの時代にも，八丈絹を米の代わりに納めていたようだ．

もっと言えば，逆に織物を納めて米をねだるという慢性的な食糧不足ぶりだった．

流人

「鳥も通わぬ」絶海の孤島といった文言に加え，「八丈島と言えば流罪人の島」とも文献に書いてあった．

前者はともかく後者に関してはその通りで，江戸時代に相当数の流罪人が送られたのは事実だ．

最初に送られた流罪人は宇喜多秀家といわれる．

五代将軍：綱吉のころから頻繁に罪人が送られるようになり，食糧問題はさらに深刻な問題となった．

どんな罪で送られたのだと問う人がいるかもしれないが，これには少々微妙な点がある．内訳をみるに女犯の僧侶が多いようだが，これを素直に受け取るのは危険だ．

どうやら当時，日蓮宗の一派が弾圧を受けており，必ずしも罪状通りの罪によって遠島されたとは言えないらしい．

流人が島に行くまでの流れ

八丈島版：流人道中記，ではないが文献に詳しく書いてあった．

  • 前日に牢屋見廻りが囚人一同を呼び出す
  • 手鎖・腰縄で連れ出される
  • 市内の髪結に連れていかれる
  • 入用のものを渡される
  • 身分に応じて駕籠なり羽交い絞めなりで与力に引き渡され，罪状改め
  • 伊豆代官手付の者に渡される
  • 永代橋近くの三宅島行きの船で出帆
  • 途中，浦賀の番所の流人改めを受ける
  • 三宅島へ
  • 三宅島役所の責任で船を雇われ，数か月後の順風を待って八丈へ送られる

八丈に着くと，くじ引きなどでどの村に配置されるかが決められた．文献には明記されていなかったが，おそらくこれは流人が多くなって少し厄介者のようになってきた江戸後期の話で，まだ流人がまれにしか送られてこなかった時代には割と歓待されていたようだ．

流人の生活

武士や僧侶などの流人は島人の教育係となったり，大工は神社再興だとか江戸風の民家を建てるだとかで役に立っていたそうだ．

また八丈は無医の島であったから，平時には素人療法・あやしげな加持祈祷が横行した．従って医者，ないし薬の調合を多少知っているだけの者でも流されてくると重宝された．もっとも，罪人の中には遊び人だっているわけで，なかなか島の風紀を乱してくれたらしい．

また罪人の結婚に関してだが，独身で過ごすのもしんどいので，大抵は水汲み女という名称で，島の女性を妻として生きたらしい．ただ女性問題で島人と紛争を起こした者は相当あったようだ．

何にせよ，流人もまた食糧難の被害者であった．

宇喜多秀家は83歳で亡くなる際，「せめて米の粥を食って死にたい」と言った，という話が残されている．

また永見大蔵は米を求めても売ってくれる人がおらず，千両箱を枕に餓死したらしいという話が残されている．

彼らは比較的恵まれた流人だと思われるが，生活に根の張っていない並みの流人は悲惨というほかないような状態であり，流人が凶作における一番の被害者であったといってもよい．

ちなみに千両箱のくだりだが，島では物々交換が長く行われており，そもそも金貨銀貨などが意味をなさなかったのではないか，とも言われているようだ．

刑罰

そんな状態に耐え切れず，飢饉の際に流人はまれに島の脱出を企てた．これは事前にばれたり，追いかけて行って見つかるなどすることが多かったようだ．

問題はその後の刑罰だ．

脱走犯が見つかった折には死刑となるわけだが，島にはもともと刑場などなかった．

ここが少々興味深いところだが，島では古くより血のケガレを嫌い，幕府に上納する織物を織るという点でもケガレを忌避する動きがあったそうだ．従って死刑の場合には海中へ突き落として刑を執り行っていたらしい．

とはいえ，その後，底土が浜にて刑が行われるようになった．

この血のケガレについては，後に記す「他屋」の節にも関連するので，興味があればそちらも見てほしい．

一般の島人たちも犯罪を犯すことがあるわけだが，この場合には隣の八丈小島に流されるなどされたらしい．

産業

田畑，そして八丈絹．これに尽きる．

今でこそ，葉物の市場とか，観光とかが栄えているが，江戸当時は漁業とか牧畜を盛んにしようという発想には至らなかったようだ．

背景には，流刑地でもあることから船の製造が規制されたこと，また「四つ足」と言われる動物を食らうのを忌避したことが挙げられるそうだ．

すなわち島の食料は，もはや八丈絹による輸入頼みでしかなく，黄八丈の生産こそが八丈島の生命線だった．

これを織る女性の社会的地位は高かったといわれ，当時の人口比率でみても女性の方が若干多い．だからといって男子の間引きだとかは文献として残されていないようだが，人口の増えすぎは問題となるので，出生管理は少々行われていたらしい．これに関連した話では，内地（本州）への移民政策などもあったようで，特に未開墾の土地が多かった茨城に移り住むことがあったようだ．しかし現地の受け入れ態勢がままならず，成功例は多くなかったという．

産業の話に戻って，今でも島での重要な産品として知られる芋焼酎について書いておこう．

不作の島ゆえ青木昆陽先生に有名なサツマイモの生産が，島でも行われていたらしい．本格的な生産に至るまでは相当な年月があったようだが，飢饉の際にはこれに救われた島民も多かったという．

さて，島には娯楽も少ないゆえ，酒は格好の嗜好品だった．古くから伝わっていたのは米による濁り酒らしいが，食糧難の島で米がどれだけ貴重なものかは想像するに余りある．そこで芋による酒造りを思い立ったのは当然の帰結だろう．

伝説

ざっくりしか書かないが，女護ヶ島伝説というのがある．

出典は覚えていないが，その昔，八丈島には女が住み，隣の青ヶ島には男が住んでいたという．彼ら青ヶ島の男たちは年に一度，八丈島へ向かう．八丈島の女たちは浜にわらじを並べてこれを待った，などなど．

真偽不明ではあるが，前節に書いた通り女性の地位は高く，女護ヶ島伝説の背景に人口比率があるのは間違いないのではなかろうか．文献では，島に漂着した難破船の乗組員が島の女性たちに手厚くもてなされ，内地に帰った際に尾ひれがついて話が広まったのではないかと考察している．

なお源実朝が八丈島に来た折，男女で島を住み分けるこの風習を止めさせたなども語られているところではあるらしい．

他屋

八丈島の昔の慣習のなかでもひときわ珍しく，そして悪しきものは他屋（たや）だ．

これは，女性の出産ないし月経時に血を避けて，別の家屋に移ることなのだが，出産時はともかく後者の場合にはその小屋に酒や肴を持った男が押し入り，夜通しランチキ騒ぎとなったという．

これは明治の世になって徹底的に取り締まられた．

まとめ

江戸時代の八丈島に関する文献のみを見ていたため，かなり尖った内容になってしまった．最近の島では観光業を中心に，農業では葉っぱや切り花の生産が，漁業ではくさやなどの生産が盛んだと聞く．こうしたデータに関しては東京都の八丈支庁が出している事業計画などが詳しい．

中でも気になったのは高速マイムマイムだ．普通に面白そう．

現在では羽田空港から八丈島空港まで55分で行けるようになり，アクセスに関してはもはや絶海の孤島ではない．日常生活に関しても不自由はほぼ無いといってよさそうで，携帯電話も大半のエリアで使用可能だ．現在の島の様子，ぜひ実際にこの目で確かめたいと思う．

技術資料

レンズ解析系の技術まとめ

まだ仮開設

非球面の式

いろいろある．その点は設計の自由度なんじゃないだろうか． たぶん一般的な非球面式は

\begin{aligned} x = \frac{cy^2}{1 + \sqrt{1 - (1 + k) c^2 y^2}} + a_4y^4 + a_6y^6 + \cdots \end{aligned}

で，第二項以降の次数は2次からの人もいれば奇数次を入れる人もいたりする*1．

中川『レンズ設計工学』では

\begin{aligned} x = \frac{cy^2}{1 + \sqrt{1 - c^2 y^2k}} + a_4y^4 + a_6y^6 + \cdots \end{aligned}

のようにコーニック係数kが違うという奇天烈な状況である．いずれにせよレンズ頂点を0としたときのサグ量を求める式で非球面を定義するのが普通と思われる．

非球面式の1階微分

非球面式を

\begin{aligned} x = \frac{cy^2}{1 + \sqrt{1 - (1 + k) c^2 y^2}} + a_4y^4 + a_6y^6 + \cdots \end{aligned}

とした場合，1階微分は

\begin{aligned} \frac{dx}{dy} = cy\left\{1 + \frac{c^2y^2 (k+1 )}{\sqrt{1 - c^2 y^2(1 + k) } + (1 - c^2 y^2(1 + k) )}\right\} + 4a_4y^3 + 6a_6y^5 + \cdots \end{aligned}

でまとめられる．

法線ベクトル

レンズ曲面は回転体だとすれば

\begin{aligned} x = u\cos v \qquad y = u\sin v \qquad z = f(u) \end{aligned}

で表現でき，法線ベクトルは

\begin{aligned} \frac{1}{\sqrt{1 + \{f'(u)\}^2}} \left( \begin{array}{c} -f'(u)\cos v\\ -f'(u)\sin v\\ 1 \end{array} \right) \end{aligned}

となる．ただし単位ベクトル化しておいた．

非球面式の2階微分

今のところ1度も使っていないが，せっかく計算してあるので供養する． 非球面式を

\begin{aligned} x = \frac{cy^2}{1 + \sqrt{1 - (1 + k) c^2 y^2}} + a_4y^4 + a_6y^6 + \cdots \end{aligned}

とした場合，2階微分は

\begin{aligned} \frac{d^2x}{dy^2} = \left(c + \frac{c^3y^2(k+1)}{1-c^2y^2(k+1)}\right) \left\{1 + \frac{c^2y^2 (k+1 )}{\sqrt{1 - c^2 y^2(1 + k) } + (1 - c^2 y^2(1 + k) )}\right\} + \sum_{i=1} 2(i+1)(2i+1)a_{2(i+1)}y^{2i} \end{aligned}

でまとめられる．和記号の部分はグチャっているが，i=1,2,3,\cdotsと入れてもらえれば係数a_4, a_6, a_8,\cdotsなど偶数次の係数が得られるはず．

*1:zemaxのサイトでは拡張多項式と呼んでいた

https://support.zemax.com/hc/ja/articles/1500005489981-ZOS-API-%E3%83%A6%E3%83%BC%E3%82%B6%E3%83%BC%E6%8B%A1%E5%BC%B5%E3%82%92%E4%BD%9C%E6%88%90%E3%81%97%E3%81%A6%E3%83%81%E3%82%A7%E3%83%93%E3%82%B7%E3%82%A7%E3%83%95%E3%81%8B%E3%82%89%E6%8B%A1%E5%BC%B5%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%81%AB%E5%A4%89%E6%8F%9B%E3%81%99%E3%82%8B%E6%96%B9%E6%B3%95

tkinterで画面サイズがうまく得られない案件

内容：画面サイズが正確に得られない

題名の通り．tkinterのwinfo_screenwidth()やwinfo_screenheight()で，カタログスペックどおりに画面サイズが得られない！

Python標準ライブラリのctypesを使って正確なカタログスペックの値が得られたよ．

Pythonを使って画像の表示位置を調整したい人，例えば画面のど真ん中に画像を出してオープニング画面（スプラッシュ画面）を出したい人などに役立つといいな．

きっかけ　Python 3.8 -> Python 3.10にしたとき

同じ症状で悩んでいる人もいるのではないかと思うので，きっかけとなった問題点から順に書いていくことにする（この問題点の方は未だ改善されず）．

引き続き，Pythonでレンズ解析ソフトを組んでいます．

やさしいコードを目指し「なるべくPythonにデフォで付いてくるライブラリを使う」という方針で組んでいますが*1，そうすると使えるGUIパッケージはtkinterの一択．

ここでPython3.8のtkinterで不具合があり，表を書くときに使うtreeviewの行の色を変えるときに

(treeview名).tag_configure((tagの名前), foreground=(色), background=(色))

だけだとなぜか動作せず，その前に

def fixed_map(option):
        return [elm for elm in style.map('Treeview', query_opt=option) if elm[:2] != ('!disabled', '!selected')]

style.map('Treeview', foreground=fixed_map('foreground'), background=fixed_map('background'))

というオマジナイを書かねばならない*2．ただしstyleはttk.Style()で得るものとする．例えばstyle=ttk.Style()．

これだけだと行の色が変わる．

が，表作成の時に行数を過剰にしてしまった余分な領域は白塗りのままとなるから，不格好なのだ．

これを避けるには事前にstyle.theme_use("alt")とかを用いて（altの他にも"clam"などが動作する），

style.configure("Treeview", background=(色), fieldbackground=(色), foreground=(色))

を入れる．

これだと前のオマジナイとかは必要なく，一括して「行の背景色」，「余白の色」，「行の文字色」が設定できる．はずだった．

実際はそうならず（泣），余白の色だけ変わって残りは元のまま白地に黒字である．前のオマジナイを入れてもダメ．

このとき使っていたPythonのバージョンが3.8だった．「これを最新の3.10にすれば既に改善されているかしら（全然バグの報告が無いから皆さん動いてるみたいだし）」？と思ったのがきっかけである．

アップデートして，パッケージ全部入れ直して，いざ動作させると...

あれ？オープニング画面が真ん中に表示されない！

tkinterで自称：画面サイズが得られる方法と実際，その解決法

全体的にこれが良い解答だと思う（→）*3．

冒頭に書いた通り，tkinterのwinfo_screenwidth()やwinfo_screenheight()で画面サイズが得られるはずである．

ただ，高解像度スクリーンではシステム側で解像度を調整しているらしく*4，この関数では調整後のニセの値が帰ってくる．

つまり，システム上どんな解像度で動いているかを知るには十分なのだが，そもそもの画面サイズが必要になる1280×700pxの画像を画面のど真ん中に出したいなどの要求では無理が生じる．

そこでPython標準ライブラリのctypesを用い，

import ctypes
user32 = ctypes.windll.user32
user32.SetProcessDPIAware()
[scw, sch] = [user32.GetSystemMetrics(0), user32.GetSystemMetrics(1)]

とすれば変数scwとschに正しい値が得られる．

他の疑問

なぜ3.8できちんと動いていたんだ？

また，画像サイズがなぜか前より小さくなった．

どうも調べてみると前のサイズがデカかったらしく，今までは上手に画像の解像度を読み替えて何となく動いていたらしい．

*1:さすがにグラフ描画はmatplotlib，計算処理はnumpyを使っています．

*2:https://ja.stackoverflow.com/questions/64095/python-ttk-treeview-python3-7%E3%81%A7%E3%83%AA%E3%82%B9%E3%83%88%E3%81%AB%E5%89%B2%E3%82%8A%E5%BD%93%E3%81%A6%E3%81%9Ftag%E3%81%AB%E5%AF%BE%E3%81%97%E3%81%A6%E8%89%B2%E3%82%92%E8%A8%AD%E5%AE%9A%E3%81%99%E3%82%8B%E6%96%B9%E6%B3%95

*3:https://qastack.jp/programming/3129322/how-do-i-get-monitor-resolution-in-python

*4:https://learn.microsoft.com/en-us/previous-versions/windows/internet-explorer/ie-developer/platform-apis/aa770067(v=vs.85)?redirectedfrom=MSDN

リモートスイッチ：Canon RS-60E3レビュー

購入のきっかけ

比較的長時間（たぶん25 s 程度だが）の露光を行って夜景が撮りたかった．

使用条件

Canon EOS 80D + RS-60E3

ライブビューで撮影

最初に手こずった点

半押ししてもピントが動かず，全押しすると突然シャッターが切れる現象に遭遇．

課題解決方法

原因は端子の接触不良による．奥まで，けっこうシッカリと差し込むことが大事．

どこかのサイトで「手振れ補正がONだとまずい」みたいな記載があった気がするが，それは正しくない．

あくまで本体のシャッターと同様に使用できる．できていない場合は本体設定か接触を見直そう．

本体のボタン設定を変えると，たぶんリモートスイッチの方も追従して変わるだろうと思われる．

今後の使用方針

まずはガンガン使っていく．必要無さそうなら，その時はその時．

収納時，コードを本体に巻き付けることができる．本体に2か所の凹部があり，コードを固定できる．また端子も本体に差し込んで固定できる．

ただし線が柔らかいため，断線の可能性を常に考えておく必要があると感じた．収納時は，本体凹部に固定はするものの，

  • 大きな曲率（小さな曲率半径）が生じないこと
  • 線に過大な張力が加わらないこと

に留意すべきだろう．

写真レンズ解析プログラム・実装したいものリスト

写真レンズ解析プログラムとは

2020年9月2日に作成を開始した、Pythonによるプログラム。

プログラム作成目的

もともと以下の漠然としたニーズが（自分の中で）あった。

  • 持っているレンズの性能くらいは（できれば他の人よりも）知っておきたい
  • レンズに興味があっても、自分で買ったり作ったりして実験するにはハードルが高すぎる
  • 写真レンズのデータは特許に記されているものの、解析は難しい
  • これを解析するためのツールはピンからキリまであり、その幅が極端に大きい
  • そもそもどんな機能が必要で何が出来れば嬉しいのか

従って勉強も兼ねて作ってしまえ！と思うのは妥当であった。

ここで上記のニーズ、そして素人レベルの技量によって生じた制約は以下のようになった。

  1. 当面は幾何光学による解析のみを行う
  2. 屈折面のみからなる共軸光学系を対象とする
  3. レンズ断面と光路図の描画
  4. 結像性能欠陥（収差）の評価

1つ目は技量による制約で、2つ目は写真用交換レンズの傾向として最低限満たすべき条件である。

当初は3つ目が主題だった。最近は4つ目にシフトしてきたが、3つ目に関係する課題は依然として残されている。

このページの目的

何をすればいいか分からなくなってきたので、

やってきたこと・やりたいことを一覧にし、今後の開発で何を優先すべきかに対する指標とする。

しかし冗長に書いてもうんざりなので、簡単な備忘録的まとめに留める。

いまできること

前提条件

パッケージ、レンズ要素、レンズデータ

パッケージ

  • math
  • numpy
  • matplotlib
  • tkinter
  • openpyxl
  • os
  • threading
  • glob

レンズ要素

屈折のみからなる共軸光学系のみを対象とする。反射や軸外しは想定外。

球面・非球面に対応。ただし非球面は多項式近似のみで、係数のデータを要する。

絞りは1つのみ。手動で元データをいじくり、F値を変える。

初面が凹面だと計算できない可能性がある。魚眼レンズなどの画角が180度を超えるもの（半画角が90度を超えるもの）は対象外。

レンズデータ

xlsx形式のExcel読み込み。規定フォルダ内に、レンズ1つあたり1つのExcelファイルを正しく作れば、プログラムが自動で後述の解析を行える。

必要なデータは以下の通り。

  • 曲率
  • 非球面の場合は各種の係数
  • レンズ間隔
  • 屈折率
  • 分散
  • レンズ有効径
  • 絞り直径
  • 絞り位置
  • バックフォーカス
  • センサー高さ

センサー高さ、バックフォーカスは自分で設定する必要がある。実は古い特許には書いていない場合も多い。

データ上、物界・像界の屈折率は何でも入れられるようになっているが、光線追跡の初期条件を設定するプログラムが空気の1.0以外に対応していないので、これ以外を入れると光線追跡できなくなる可能性がある。従って実のところは1.0に限定。

絞りはレンズの一部ではなく「別の特殊部品」として計算しているため、絞り位置は初面からの距離で設定する。特許データには見られない形式なので、データ作成には注意を要する。

非球面の場合、係数を別のシートに記入する。次数に制限はないが、各次数に対して具体的な数値を入力する必要がある。

解析系

硝材の設定、光路図、収差図、その他の解析、データ保存

硝材の設定

色収差を計算するために、硝材（光学ガラスの種類）を設定する必要がある。

硝材は4社のデータに対応。ただし最近のデータのみで、過去のデータは無い。

特許にある屈折率と分散から、もっともらしいデータを自分で選択する。

光路図

レンズ断面と、センサーの中央・7割高さ・端に至る、主光線・上側マージナル光線・下側マージナル光線の描画。

収差図

  • 中央・7割高さ・端におけるメリジョナル・サジタル横収差の描画
  • 球面収差とOSC、M/Sの各非点収差、中心射影方式による像面湾曲を示す縦収差図
  • 上記のそれぞれに対し、3つの波長をまとめて描画（色収差に相当）

その他の解析

以下は色収差に非対応だが、 * スポットダイアグラム（中央・7割高さ・端） * 幾何光学的MTF（10 LP/mm、30 LP/mmのM/Sそれぞれ）

データ保存

グラフは全てjpegかpngで保存可能。保存したい場所を選び、名前を設定することが可能。

硝材データは保存できず、毎回の計算で設定しなおす。

（特筆すべき）内部の数値計算系

初期条件設定・光路図、スポットダイアグラム、幾何光学的MTF

初期条件設定・光路図

計算したいExcelのデータファイルを選択しボタンを押すと、光路図の計算を開始する。

このとき、まず光線行列解析により、近軸的に以下を求める。

  • 焦点距離
  • 入射瞳半径
  • 射出瞳半径
  • 射出瞳位置

この数値は解析の基本となるので、解析中は保存される。ただし現行のプログラムでは表示されない。

光路図の計算は、センサーの中央・7割高さ・端それぞれにおいて以下のように行われる。

  • 上で計算した射出瞳位置を狙ってセンサー側から主光線の逆光線追跡
  • 主光線データを順方向に変換・保存
  • 主光線データをもとに上下マージナル光線を計算・保存
  • 主光線データをもとに左右マージナル光線を計算・保存

解析中、これら主光線と4本のマージナル光線のデータは常に保存され、収差図の基本として使われる。

スポットダイアグラム

まずレンズ初面の直前において、上で求めた主光線位置を中心に、タテヨコともに入射瞳直径\times 3の正方形領域を直交メッシュに区切る。

各格子点から主光線に平行な光線を出し、それがセンサーに交差する点をまとめて散布図に表示することでスポットダイアグラムを描画する。

幾何光学的MTF

スポットダイアグラムをフーリエ変換する。

今後やりたいこと

プログラム寄りなもの

計算速度

計算速度がボトルネックとなっているのが、スポットダイアグラムとMTF特性図の2つである。

特にMTFについては、現在の計算性能では信頼性のおける結果を得るのは不可能（時間を犠牲にするなら話は別である）。

以下に、最近実装した改善案とその結果を示す。即採用となった結果と、思いのほか性能が上がらなかった結果があり、興味深い内容となった。

最近行った改善

最近行った改善策では、

  1. 今まで1本の光線につき1回光線追跡の関数を使用していた
  2. 関数の呼び出しは内部の変数（特にレンズデータ）呼び出しも含むから、時間かかりそうだぞ
  3. 一気に複数の光線の初期条件を入れられないか？
  4. 行列（numpy配列）で計算できそうだ
  5. 非球面のとき光線とレンズの交点は近似で求める必要があるから、その段階だけは別々に計算しなきゃ

という手順（よりは紆余曲折あったが）を踏んで、確かに場合により高速化できた。以下、スポットダイアグラムの計算において同様の初期条件で計算した結果を示すと、

屈折面数 うち非球面枚数 レンズタイプ 元の計算時間 [s] 改造後の計算時間 [s] 計算時間の比
18 0 望遠 1.99 0.683 0.343
19 0 広角 3.44 0.680 0.198
18 1 広角 4.46 3.76 0.843
27 4 広角 10.9 14.2 1.30

ただし面倒だったので、データは1回の試行による。とはいえ毎回こんな傾向があるのは確かである。比を取って比較するのが適切かは置いておくとして、結果により以下のことが伺える。

  • 球面しかないレンズは早い
  • 改造は球面のみの広角レンズに特に効果的だった
  • 非球面が1面でも入ると改善効果は薄くなる
  • 屈折面も非球面数も多い場合、改造は裏目に出る
  • 改造は屈折面数に関係ないので、時間増は非球面の近似計算に起因すると予想できる
  • 余談だが、計算中に表示させるプログレスバーの動きが格段に良くなった

とすると、次に行うべきは近似計算の改善である。現行のプログラムでは収束の堅牢さと処理の単純さから二分法を選択している。ニュートン法では収束しない場合が散見されたため、Dekker法やBrent法などの堅牢かつ高速なプログラムを用いる必要がある。

Dekker法をまず簡単に組んで用いたところ、計算速度への効果は期待したほどではなかった。

屈折面数 うち非球面枚数 二分法での計算時間 [s] Dekker法での計算時間 [s] 計算時間の比
18 1 3.76 3.37 0.896
27 4 14.2 17.3 1.22

原因として挙げた仮定は、

  • 大して反復回数が減少していないのではないか
  • 処理（特に条件分岐）の複雑化が想定より大きく影響している

まず1つ目を調べたところ、二分法だと反復が最低25回、最大37回だった一方で、Dekker法では最低2～4回、最大11回と大幅に効率化されていることが判明した。

従って今後は条件分岐の効率化を求めることになる。

pythonファイルどうしの関係を見直す

プログラムの規模として、どうしても複数ファイル、複数クラスを横断するほかないが、

  • どの関数が何を必要とし、何をreturnするのか
  • 例外処理はどう行われ何がreturnされるのか
  • 同一ファイル内にどのクラスを入れるか
  • 同一クラス内にどの関数を入れるか
  • つまり、関数をもっと細分化する必要は無いか
  • 逆に一体化して良いものは無いか

といった基本的な部分がまだ固まっていない。より見やすいコードを目指したいし、これをおろそかにすると例外処理で沼にハマりそう。まずは上の2つから見直したい。

硝材の保存

いま、Excelファイルのシート1にレンズの基本データ、シート2に非球面データを入れている。そうすると、硝材データはシート3を新たに作成するのが濃厚な手段だ。

硝材データ案

つまりシートをこんな感じにする。

  1. 基本データ
  2. 非球面係数のデータ
  3. 硝材データ（New！）

ただ、硝材シートを新たに作るとなると問題が生じる。

今のところシート読み込みをシート名ではなくシートの順番（つまり配列のような形）で読み込ませているのだが、現状では球面のみでできたレンズ系に対して非球面のシート（つまりシート2）を要しない。このとき

  1. 基本データ
  2. 硝材データ

の2つで事足りるが、これでは硝材データを非球面データと勘違いして読み込んでしまうだろう。

従って改善策は以下の2つが考えられる。

  • 空のシート2（ダミー非球面シート）を作る
  • プログラム全体をシート名で読み込ませるよう変更する

Excel読み込みで使用しているopenpyxlは内部で勝手にシート番号をつけてくれる他、我々がつけたシート名でも読み込み可能なので、上記はどちらも達成できる。ただし前者だと無用なシートができるし、後者だと命名ミスで動かなくなることが有り得る。

吟味したいが、たぶんシート名読み込みになるだろう。

レンズ系の基礎

パワー配置

特許ではよく「第何レンズ～第何レンズを第1群とし...」といった表記がされる。設計においても、複数のレンズが集まった「群」をひとつの薄肉レンズとみなして初期設計する。

従って群というコンポーネントは非常に重要となる。これを解析するプログラムを含められないか。

コンポーネントの選択さえできれば、あとは光線行列解析で焦点距離や主点位置などが求められるので、計算はそう難しいことではない。

問題はおそらく選択する画面とか選択方法、あとは結果の表示だろう。tkinterの表・チェックボックス・プルダウンなどの機能を見比べる必要がある。

硝材の配置を可視化

硝材の配置がどうなっているかは素人目線で少し気になる。

屈折率の強いレンズとか、分散がおかしいレンズ（異常分散ガラス）がどこに配置されているかが一目で分かれば楽しいだろう。

これを実現する1つの案が、レンズ断面図を光軸上下で分け、上に屈折率、下に分散のカラーバー（数字が大きいほど色が濃くなるやつ）でレンズを塗りつぶす方法である。

あまり一般的ではないが、できれば面白いと思う。見やすくなるかは分からない。

初面凹面の対応

現在、光線の初期位置をレンズ頂点に接する平面としているため、もし初面が凹面だとレンズ内部から発光しているような様子になる。

考えられる改善方法は、初期位置を頂点からオフセットしておくものだ。ただしオフセット量を

  • なにかの定数
  • 凹面が前にせり出している量（サグ量）の最大値

とする2つの方法が考えられる。

実は開発初期、これを考慮して一つ目の方法を取っていた。ただ、定数を足し引きする計算の箇所がうまく管理できず、得られたデータがオフセットを含むか否かが判断しにくかったこともあり、現在は頂点の接平面となっている。

従って、どちらの方法を取るにせよ

  • 欲しいデータはレンズ頂点基準か光線の初期位置基準か
  • どこでオフセットを加えてどこで引くのか

はキッチリと押さえておく必要がありそうだ。ちなみに現行の光線追跡計算上、光線データはレンズ頂点を原点とするのが楽である。

収差関係

収差係数

収差係数といって、設計値と近軸光線追跡データの四則演算だけで収差量が評価できる指標がある。

これを計算して表にするボタンがあると嬉しい。

ただし収差係数の式は多岐に渡る。本国で有名なのはおそらく松居吉哉氏の『レンズ設計法』にある式だと思われる。一方、私が導出できたのは中川治平氏の『レンズ設計工学』にある式で、確か何かしらの定数で除するか否かという点で相違があった。片や3次収差論で世界的に有名なのはBerek氏の『レンズ設計の原理』で、ここでの式がどういうものだったかは覚えていない。

導出できていない式を、自習も兼ねたプログラムに実装するのは忍びない。一方で、敢えてマイノリティーな式で計算させるのも悩ましい。

おそらくカッチリ数式を表記したreadmeの作成は必要になるだろうが、加えて上記2、3パターンの選択形式とするべきか悩む。

歪曲シミュレーション

いま、既に歪曲収差を縦収差図の形で表示できる。しかし実際に格子状の物体がどう歪むかをシミュレーションできれば非常に便利だ。

想定しているのはレンズメーカーのシグマなどがポツポツと出しているような一般的な図で、理想的な格子と歪曲された格子を重ねて描くもの。

基準となる格子を設定する必要があるため、センサーサイズ、特に縦横比を設定する必要がありそうだ。また、格子の描画が簡単にできるものか調査する必要がある。

発展した光学

ここら辺から、聞きかじった知識でしか書けないレベル（難しいの）。ぶっちゃけ計算できなくてもいいかな（諦め）。

波面収差

光線ではなく波として捉え、収差を波面の乱れとして考える方法。

図示はできないが、以下に示すものの基礎となる。

点像分布関数

PSF（point spread function）と呼ばれるもの。計算するには、たぶん瞳関数も必要。

点光源がどのような輝度分布を持ってセンサーで結像するかを表していて、解析できれば嬉しい。

波動光学的MTF

より精密なMTF特性図。

レンズの光線追跡式

光線追跡プログラム開発1周年

先日の2021年9月2日をもって、写真レンズ関係のプログラムが開発1周年を迎えました。

まだまだ進化を続けていこうと思いますので、温かい眼差しで見守っていただけると幸いです。



厳密な光線追跡をしよう

レンズ設計・性能解析の要と言って過言でない光線追跡。設計段階の初期では近軸追跡という三角関数の近似計算がメインとなるものの、いずれにせよ厳密な光線追跡は大事。

簡単に解説すると、平面内（メリジョナル断面内）で考えるだけならn\sin i = n'\sin i'と計算したり、3次元なら外積に相当する計算を繰り返して、厳密な光線追跡を行います。

ズボラなので図を1つしか描いていないけども、この記事を読み込めば2D光線追跡プログラムが組めるレベルに仕立て上げたので、興味がある人はぜひ挑戦してみてください。

2次元に限って言えば、高校1年生の知識だけで読めます。

3次元はまた次の機会に触れようと思います。





基本のき

光線追跡に限らず、レンズには様々な「流儀」があります。まずはプロフェッショナル：レンズの流儀を身につけましょう。今回の記事で、そこらへんの雰囲気を会得していただくと専門書などが読みやすくなるのではないかなと思います。

基本：初級編

「光線」の概念は当然の如く認める。従って波長とレンズ系の規模の間に制約があることを断っておく。

さらに認めてしまいましょうという原理が、スネルの法則の名で有名な

\begin{aligned} n\sin i = n' \sin i' \end{aligned}

である。屈折率と入射角の関係を表す重要な式だが、これを大前提として以下の話を進める。

ちなみに、以下では共軸の屈折光学系だけ考え、軸外し・反射面・トロイダルレンズの存在などは全く考えない（反射球面は案外簡単に議論を拡張できるらしいが踏み込まない）。

うっ、こう考えると制約が多い...

基本：中級編

確か、歪曲収差の記事で「角度は光学屋さんの流儀に則る」という文を複数書いた気がする。これがかなり特殊なので、この段階で明確にしておく。

まず、左に物界、右に像界を持ってくる。もっと広く言えば、だいたい光が左から右へ進むようにする。

レンズと光軸の交点を、レンズ「頂点」とか、屈折面の「頂点」と呼ぶ。

物点・像点はいずれもレンズ頂点から右側（つまり像より）を正に取る。

つまり中学物理の凸レンズは

\dfrac{1}{a} + \dfrac{1}{b} = \dfrac{1}{f}

ではなく

-\dfrac{1}{a} + \dfrac{1}{b} = \dfrac{1}{f}

になる！

まじかよ...（まじです）。

自ずと光線角度の正負に制約が生まれ、右肩下がりを正、右肩上がりを負とする。

例えば、無限に薄い平板上で光軸から距離hだけ離れた点から、平板の後方（つまり右方）にLだけ離れた光軸上の点に向かう光線の場合、角度uを

\tan u = \dfrac{h}{L}~~,\quad\therefore\quad u = \arctan \dfrac{h}{L}

で定義するのが自然な感じがするでしょ？だから右肩下がりが正になる。

具体的な状況は、子午面内の光線追跡を読んでいくと分かるだろう。

基本：上級編

まず直感的に説明すると、レンズの曲率半径は頂点が左に凸の場合を正、逆を負と定めます。

実際には曲率半径が負というのは数学的に誤りなうえ、よく考えれば中級編の基準に沿って頂点の右に曲率中心があるのが正と理解するのが適当でしょう。

いずれにせよ、これら基準に沿って計算していくと、レンズの向き・虚像や実像などでイチイチ場合分けする必要がないことが確認できます（本記事では確認しません）。

素晴らしいな。



「球面」屈折光学系の子午面内（メリジョナル面内）光線追跡

球面だけを含む共軸光学系を例に、光軸を含むレンズ断面内の計算を考えよう。回転対称性から、光線が断面内から外れることがないと示せる。

ただ、最近はこれに特化したコードを組む場面は少ないらしい。大は小を兼ねるので3次元光線追跡が幅を利かせているようだ。

ここでは以下の図をもとに話を進める。詳しくはMax Born, Emil Wolfの『光学の原理（1）』を参照されたい。

f:id:positiveclimb:20210905091420j:plain
2D光線追跡式に利用する記号

まず、レンズ面から左側の領域の屈折率をn、右側をn'とする。

光軸に対して、曲率半径rのレンズ頂点から距離Lの点に角度uで入射する光がある。これとレンズの法線がなす角をiとする。

光線はABに沿って図の左方から入射し、Aで屈折することでCへ向かう。

基本の中級編で話した光線角度（uやu'）は、上図では右肩下がりだから正になる。

ちなみに参考書に則って記号を書いたが、邦書ではLやL'よりsとs'を見かけることが格段に多い。

まず簡単に

\triangle OABに注目する。Oを通るABへの垂線を引くと、その垂線の長さから

\begin{aligned} (\text{垂線の長さ}) = r\sin i &= (L - r)\sin u\\ \therefore\qquad \sin i &= \dfrac{L - r}{r}\sin u \end{aligned}

が導ける。

従ってレンズの設計値rと入射光線の向かう場所Bすなわちr + OB = L、そして光線角度uが分かればレンズ曲面への入射角iが分かるようになった。

次に、どれほど屈折するのか、つまりi'を明らかにしたい。スネルの法則から

\begin{aligned} n\sin i &= n' \sin i'\\ \therefore\qquad \sin i' &= \dfrac{n}{n'}\sin i \end{aligned}

非常に簡単。

また屈折後に光線が光軸となす角u'は、\triangle ABCの外角Cを求めればいいから

u' = u + i - i'

で求まる。

最後に、最初と同様に\triangle OACでOを通るACの垂線の長さから

\begin{aligned} (\text{垂線の長さ}) = r\sin i' &= (L' - r)\sin u'\\ \therefore\qquad L' &= \dfrac{sin i'}{\sin u'} r + r \end{aligned}

これで、晴れて単一屈折面の計算が終了した。

屈折面が複数ある場合、新たにレンズ間隔を定義する必要がある。レンズ間隔は、互いの頂点どうしの距離のことを言い、基本的にdで置く。正負は前のレンズ頂点基準で右を正、左を負とする（物点とかと同様）。ただし反射面を考えないと断った以上、 d \geq 0である。

従って第二面の入力値Lとuを、分かりやすいようにそれぞれL _ 2とu _ 2とすれば、明らかに

\begin{gather} L _ 2 = L' - d\\ u_2 = u' \end{gather}

と引き継げばよいことが分かる。このように複数のレンズがあるときにdを引く行為にも、きちんとLやL'の符号を定める重要性が見いだせる。

お気づきだろうか、除外した場合

ここまで、既に一般的に論じた雰囲気を醸し出しているが、実際には特殊事例を除外している。

早い話、入射・射出角度が0のときに距離LやL'を定義できないのである。

この条件は特に初期条件で与える場合が多いので、しっかり対応しておかないと後で何もできなくなる。

ところで、今まで光軸上の距離だけ考えていて、交点Aの高さを気にしなかった点を解消しよう。これは光線をグラフ化する時、特に重要になる。

法線が光軸となす角がu + iあるいはu' + i'になっていることに気付けば、レンズ頂点と交点の光軸方向距離（以下、サグ量とする*1）は

(\text{距離}) = r - r\cos(u + i) = r - r\cos(u' + i')

で求まるし、光軸からの高さhは

h = r - r\sin(u + i) = r - r\sin(u' + i')

で求まることがすぐ分かる。前者に「レンズ初面から第k面までの距離」を足せば、描画の問題は解決された。

なぜ突然これを計算したかというと、たとえu=0であれ光線高さの概念は存在するため、上式でu = 0とかを代入しても成立するからである。つまり、今度は光線高さから攻めるんだな。

上の式で、u = 0とした場合に

h = r - r\sin i = r(1 - \sin i)~~,\quad\therefore\quad \sin i = 1 - \dfrac{h}{r}

であることがわかる。屈折自体はスネルの法則を使っていいのでi'が分かるし、さっき使った

u' = u + i - i' = i - i'

がそのまま適用できるから、これ以降は先の方法に沿えばいい。

ただし、このときu' = 0が生じた場合は例外で、無理にL'を求める必要は無い。というのも、この場合は次の面でu = 0になって上の計算を適用するが、これ自体が既にLとかを欲しないシステムであるから。

従って先のナンチャッテ一般論に加えて、光線角度での条件分岐と、光線高さhを逐次求めるくらいで対応可能なことが示された。

初期条件の与え方・1回の屈折計算で計算すること

いろいろ方法があると思うが、私の方法を書いておく。

最初の穴だらけの方法だと、屈折前後で

(L, u) \longrightarrow (L', u')

という変換をしていて非常に単純だが、実際にはuのワガママに付き合う必要があることが判明し、光線高さhを導入する流れになった。

結果、私は(L, u) \longrightarrow (L', u')の中で条件分岐させると共に、内部で暗に(h, u) \longrightarrow (h', u')を行うのが楽だろうと思ったので、以下に概略を示す。

初期条件はあとで色々便利なように、レンズ初面の頂点に接する平面上での光線高さをHとして、(H, u)というセットで与える。ちょっと定義が複雑になったが、頭書の図から明らかにu\neq 0 で\dfrac{H}{L} = \tan uであるし、u = 0でH = hである。初面以降はHに相当する概念は必要ない。

これをもとに屈折の計算をし、レンズ上の屈折点が厳密に

(\text{サグ量}~,\quad\text{高さ}) = \left(r - r\cos(u' + i')~,\quad r - r\sin(u' + i')\right)

で求まり、これを用いて光線を描画できる。さらに座標の後者はhとなるから、描画以外にも重要である。

この段階でu' = u_2 \neq 0ならばL'が求まっているから、dを引いてL_2を求め、引き続き計算する。

u' = u_2 = 0ならばh = h_2なので、条件分岐させて特殊事例の計算を適用すれば問題ない。

従って、これでプログラムが組めてしまう。

平面の場合

平面の場合はrの概念が通用しなくなり、i = uかつi' = u'となる。従ってスネルの法則が

n\sin u = n'\sin u'

となり、これ自体でu\longrightarrow u'の変換が成立する。

ただしL\longrightarrow L'の受け渡しは、球面の特殊事例以上にhを活用する必要がある。

i面で\left(\text{高さ},\text{サグ量}\right) = (h _ i, z _ i)とする。レンズ間隔d _ iだけ後ろに平板のi + 1面があるとすれば、その点での高さh _ {i + 1}は

h _ {i + 1} = h _ i - (d _ i - z _ i) \tan u _ i

で求まり、

L' _ {i + 1} = \dfrac{h _ {i + 1}}{\tan u' _ {i + 1}}

である。

このようにして、平板であれ光線追跡を実装することができる。



まとめ

これであなたも光線追跡マスターになれる！と思いきや3次元光線追跡について論じていません。さすがに文量が膨大になるので、この話題は次回に回しましょう。

ただ、物点・像点や角度の流儀は共通ですし、一面ずつ屈折の様子を探っていくという基本のスタンスは全く変わりません。

素人なので計算が冗長かもしれません。「こうすれば簡単なのにな～」と感じる方はぜひコメントください。

では、この辺で。



*1:私は勝手にサグ量と呼んでいるが、厳密に正しい言葉の使い方か判然としない

歪曲収差（ディストーション）が分からない

与太話

つい先日、ようやく子午・球欠像面を計算できるようになりました。主光線が描画できたら

\begin{aligned} -n\frac{\cos^2 i}{t} + n\frac{\cos i}{r} &= -n'\frac{\cos^2 i'}{t'} + n'\frac{\cos i'}{r}\\ \frac{n'}{s'} - \frac{n}{s} &= \frac{1}{r}(n'\cos i' - n\cos i) \end{aligned}

のように、光線に沿ってtやt'、sやs'を測ってやればいいわけです。

ここでtとt'はタンジェンシャル（子午、メリジョナルと同義）光線についての物・像と屈折点の距離、sとs'はサジタル（球欠）光線についてのそれです。そのうち、この導出を載せるかもしれません。

これが計算できた利点は、いわゆる非点収差をグラフ化できたって事です。余談だからこれ以上は踏み込まないで良かろう。

今回のお題

非点収差は解決したとして、話題は歪曲収差。歪曲収差は、四角いものが歪んで見えるレンズ欠陥のこと。つまり物体の高さによって倍率が変わってしまううんだな。こういう像上の倍率を光学屋さんは特に横倍率Mという。

  • 物体の高さが高いほど倍率が小さい方にシフトする → 樽型

  • 物体の高さが高いほど倍率が大きい方にシフトする → 糸巻型

  • 両者が混ざる複雑な形状 → 陣笠型

とかいう分類がある。で、問題はこの倍率についてなんですよ。無限遠の物体の高さってなんだ？それの倍率って？

うーん、考え方が悪い。というか定義式があるはずじゃないか。ググってみると（Zemaxのサイト*1を見ると）

Dist. = 100\times\dfrac{y _ {chief} - y _ {ref}}{y _ {ref}}

である。なるほど、実際の主光線高さy _ {chief}と近軸的に理想な光線高さy _ {ref}のズレを評価するんだな。確かに。で、y _ {ref}ってどうやって求めるの？それが知りたいんだけど。

何も解決しなかった。もうちょうい調べるか。

同サイトによると物体位置によって像位置が変わるから、それに対応したディストーション焦点距離というのを求める必要があるらしい。なにそれおいしいの？

ダメだ、別を当たろう。

文献*2によると、f\thetaレンズについては

Dist. = 100\times\dfrac{f\theta - f\tan\theta}{f\tan\theta}

で歪曲を求めるらしい。

f\thetaレンズっていうのはレジのバーコード読み取り機とかに使われるやつ。レーザーを一定の角速度で回すだけだと、バーコードの中心は遅く、端は速く掃引することになるから、レンズを噛ませてそれを回避している。

あの機械そんなすごいマシンだったのか。

まあ、そのレンズに要求されるのが「像高がf\thetaになる」という条件らしい。確かに入射角\thetaの角速度\dot{\theta}=const.を与えると、焦点距離はf = const.だからバーコードの掃引速度はf\dot{\theta}になって速度一定だな。

じゃあ像高f\theta = y _ {chief}と置いていいか。おお、項の位置がZemaxのy _ {chief}にピタリ当てはまる。

これならy _ {ref}の方もy _ {ref} = f\tan\thetaと読み替えて良さそうだな。

他のサイトでも「一般的なレンズは像高がf\tan\thetaになるからf\thetaレンズは特殊」みたいなことが散見される。ほほう、何にせよ入射角\thetaによって理想像高が決定されるのね。

仮定：f\tan\thetaで理想像高が与えられる

仮にf\tan\thetaで像高が与えられるとしよう。念のため書いておくが、fはレンズ最終面から近軸像点までの距離ではなく、後側主点から近軸像点までの距離を表している。従って自分なりに仮定をまとめると

理想像高f\tan\thetaは、入射角度\thetaで主光線が与えられたとき、後側主点から入射角度\thetaで射出された場合の光線と像面の交点である。

ちょっと用語が多くなったが、この考えはある意味で正当性があるように感じる。面倒な人は次の章まで飛ばしてください。

複雑なレンズ系を等価な1枚の薄肉レンズに置き換えたとき、中学校理科レベルの考えが適用できる。つまり「レンズ中心を通る光線は屈折しない」のである。

実際のレンズ系は厚みを持つから「レンズ中心ってドコ？」となるので、中学生は混乱する。それを回避するのが主点とか節点とかいう光学系主要点の概念である。

主点（後側主点）は、薄肉レンズがそこにあると考えていい場所である。特に平行光線が軸上のある点に集光する場合、実際のレンズと同じパワーを持つ薄肉レンズを主点位置に置いたと考えれば、そこから焦点距離だけ離れたところに集光することが厳密に示せる。薄肉レンズの「平行光線は焦点を通る」条件を満たすわけである。

ただし実際のところ、一般的には「レンズ中心を通る光線は屈折しない」条件を満たすこの「レンズ中心」は主点に一致しない。主点位置に向かって入射させてもレンズ前後で屈折しちゃう場合がある。

主点は基本的に軸上だけで成立するものだと私は理解している（合ってるかは知らない）。

じゃあ角度が保存される点はどこなの？というと、それが節点である。これも近軸的な議論で位置を求められる。

ただし、物空間側と像空間側の屈折率が同じ場合は主点と節点が一致する。

今回初めて知った（ついでに厳密に示せた）。レンズ、難しいな。

まとめると、理想像高y _ {ref}をf \tan\thetaとおく計算は主点と節点が一致した場合を考えて計算していると解釈していいだろう。この点はちょっと危うくないか？

ともかく、この考え方ではf = const.っぽい雰囲気があるし有効焦点距離で事足りる感じがある。結局ここでは、Zemaxみたいに像位置によって歪曲収差量が変わるとかを気にしてディストーション焦点距離を定義するのではなく、焦点位置の時点でどれだけずれているかを考えているんじゃないか？

いざ計算に挑戦

総じてpythonでコードを組んでいる。

適当にdistortionという配列に上記の計算結果を付け足していこう。像高は、割合を意味するheight (0.01\leq height \leq 1) とセンサー対角線長さd_sの半分をかけて求める。-uは入射角で、負号の意味は気にしなくていい（光学屋さんの流儀に則ってる）。本当は今までの\thetaも負にすべきなんだよな。

distortion.append((height[j] * d_s / 2 - f * math.tan(-u)) / (f * math.tan(-u)) * 100)

これをforで回す。計算量多くなりそうな書き方だけど許して。

光線高と入射角が共に0のときに光線追跡できないポンコツ仕様なので（計算するまでもなく光軸に一致するから手を付けていない）、あとでheightに0を追加する。

そのときの歪曲収差量は...0じゃないか？うーん分からん。

ひとまず動かしてみるか。軸ラベルは面倒だから、まだ書かない（不親切）。

f:id:positiveclimb:20210815113334p:plain
なんか歪曲収差が違うんだが。左から球面収差、非点収差、歪曲収差。

...なんかおかしいぞ。

たしかに数式上は、入射角を0に漸近させるとy _ {chief} \longrightarrow 0かつy _ {ref} \longrightarrow 0だから

　\displaystyle \lim_{\theta \rightarrow 0} 100\times\dfrac{y _ {chief} - y _ {ref}}{y _ {ref}} = 100\times\dfrac{0-0}{0}

となるんで、なにか極限値があってもおかしくなさそう。安直に歪曲収差量は\theta = 0で0とかしちゃいかんな。

いやでも普通は0でしょ。

別の文献によると

別の文献*3によると、これも節点（nodal point）を用いた計算をしている。

文献では、節点から焦点までの距離をf、焦点距離をp、主光線の角度をuとし、レンズ系の前後でプライムを付けている。

文献中の式（7）によると、

D _ m = \left(p' \tan u' - f \tan u \right) - \left(p \tan u - f \tan u' \right) M

物点を無限遠に飛ばすと

D _ 0 = p' \tan u' - f \tan u

なるほど、入射角uの方向を保ちながら物点を無限大にするんだな。

このD _ mとかD _ 0は今まで議論してきた式の分子に相当する。

p' \tan u'は今までのy _ {chief}に、uは入射角\thetaに相当するから、上式を書き改めると

D _ 0 = y _ {chief} - f \tan \theta

となる。

あれ、これ前のやつと一緒じゃん。

ひょっとして私が何か間違えてるだけで、式そのものは前ので合ってるのか？

問題発見

さっきの極限の話に戻ると、\theta \longrightarrow 0つまりy _ {ref} \longrightarrow 0でDist. \longrightarrow 0を要求するなら、y _ {ref} = hと置いて、

\displaystyle \begin{aligned} \lim _ {h\rightarrow 0} \dfrac{y _ {chief} - h}{h} &= \lim _ {h\rightarrow 0} \dfrac{\frac{\partial}{\partial h} y _ {chief} - 1}{1}\\ &= \lim _ {h\rightarrow 0} \frac{\partial}{\partial h} y _ {chief} - 1 = 0\\ \therefore \qquad \left. \frac{\partial}{\partial h} \left[y _ {chief}\right]\right| _ {h = 0} &= 1 \end{aligned}

使っていいか分からないロピタルの定理を適用して、上の結果が得られる。

最後の式の書き方も厳密には上側極限が正しい感じがあるので合っているか分からないが、こんなもので良いだろう。

つまり、像中央部では理想像高と実際の像高が傾き1の一次の関係で近似できる必要がある。

まずこの関係をそのままグラフにしてみたんだが、もはや近似関係は誤差レベルでほぼ成立しているように見えたので失敗。条件を変えよう。

この条件によれば、縦軸にy _ {ref}、横軸にy _ {chief} - y _ {ref}を置いてプロットしたとき、傾きが±無限大になる必要がありそうだ。やってみると...

f:id:positiveclimb:20210816140506p:plain
さっきと同じレンズ系に対して、右端のグラフの軸を変更した。原点の傾きが有限値である。

うーん、やっぱり傾きが...。原因は何だ？

解決

この問題は、たとえ像中心部でも実際の像高が理想像高と定数倍異なる（かつそれと一定の差があることも想定できる）ことによる。さっきのグラフでは原点を通っていたから、ひとまず一定の差がある方は無視していいだろう。

理想像高と実際の像高で比例定数≠1の比例関係にあるような欠陥とは...。

あるいは先のグラフのように、収差と理想像高が像中心部で謎の比例関係を持ってしまう欠陥とは...。

思いついたぞ。

私のプログラムでは、像高を「本当の像高」にしているんだ。

像は収差によってボケるため、最良な像の位置は近軸像点とは異なる。先に示したような非点収差に起因する最小錯乱円、あるいは球面収差に起因する火線などなどを吟味して初めて最良像面が決定される。

私のプログラムでは、そうやって決められた近軸像面とは無関係な像面上の主光線高をy _ {chief}としていたんだ！

y _ {ref}がf\tan\thetaのように近軸像点に関係あるfを用いる以上、主光線高は近軸像面基準に直さねばならない。

仮にこのような計算間違いをした場合、最良像面と近軸像面間の距離を\deltaとすると、像中央部では近似的にy _ {chief} - y _ {ref} = \delta\tan\thetaが成立する。従って先のグラフは縦軸がy _ {ref} = f\tan\thetaだったから、傾きf / \deltaが生じる！

さらに、\delta \longrightarrow 0だと傾きが無限大（符号は\deltaによる）に発散するから、さっきの予想と合致する！

ああ、なんということだ、ちまちまとブログにまとめ始めたのに、

こんな間違いは私にしか理解できないじゃないか...

結果と結論

修正した。それで、さっきと同様に縦軸にy _ {ref}、横軸にy _ {chief} - y _ {ref}をとるグラフを描いてみた。

f:id:positiveclimb:20210816144100p:plain
修正後の理想像高-収差グラフ（右端）。原点で傾きが発散している。

うん、大丈夫そうなので、最初のように歪曲収差をグラフ化してみよう。縦軸に理想像高、横軸にDist.（単位は[%]）をプロットすると...

f:id:positiveclimb:20210816144715p:plain
左から順に球面収差とOSC、MとSの像面湾曲、歪曲収差。完璧だ。

いやぁ、素晴らしい。他のレンズ系でも完璧に描画できるのが確認できた。

つまり、

像面湾曲量Dist. [%]は

Dist. = 100\times\dfrac{y _ {chief} - y _ {ref}}{y _ {ref}} = 100\times \dfrac{y _ {chief} - f\tan\theta}{f\tan\theta}

で計算できるが、y _ {chief}は理想像面上の値で、fは物界・像界での屈折率が等しい場合に限り焦点距離長さを適用できるが、その他の場合は節点から理想像点までの距離を適用する

と結論付けられるだろう。Zemaxのディストーション焦点距離は検証できなかったが、ひとまずこれで定義式上の歪曲収差が計算できた！

めでたしめでたし、である。

*1:https://www.zemax.jp/blogs/news/tech-tips-june-2021

*2:中川治平『レンズ設計工学』

*3:Arthur A. Magill, Variation in Distortion with Magnification, Journal of Research of the National Bureau of Standards, Vol. 54, No. 3, March 1955