リッチー・クレチァン望遠鏡の理論（『吉田正太郎,新版 反射望遠鏡光学入門』（他））

[「M」]

はらぱん さん 　みなさん　こんにちは

まずは訂正。

数式の清書はLaTexを使っているのですが、(2026 1月中という自分で決めた
締切り間際に)書式を変えた際の手作業でy^8の項をy^6の項,y^10の項をy^8と
ズラす間違いをし、そのまま投稿してしまいました（トホホ。）
余計な手間をおかけしました。

Python(/FORTRAN?)形式では

コード: 全て選択

Ds= -(12*m**3+(12*e-36)*m**2+(36-12*e)*m+e**2+e-12)/(384*e**3*m**5)

これですと
ｍ＝0.326667; e＝0.2525　で
Ds=174.534461 となります。

はらぱん>数式処理ソフトＭａｘｉｍａで厳密解を展開してみました。
「M」も(しばらく前から)Maxima(https://ja.wikipedia.org/wiki/Maxima)を
使っています!　多倍長整数/多倍長浮動小数点数の計算ができるのも便利ですね。その一方、
なかなか「使える」ようにならない。　強力な機能が「満載」で1000ページにおよぶ
Maxima Manualで迷子になってしまいます。
※たとえば「Taylor展開」は見つかるのですが「媒介変数表現の展開」が見つからない...。
近頃は「自分なりのidiomを集められればよい」と考えています。まあ自然言語と同じで
「習うより慣れろ」なのですかね。

はらぱん>これをｔにそのまま代入するのではなく、テイラー展開したものを代入して全体をテイラー展開するのだそうです。
※「そのまま代入」しても、ちゃんと結果が出ると思いますが？

はらぱん>mの値を負にしないと正しい答えが出ません。
yの0次項が「m+e」となっていますが「e-m」となる筈ですね(主鏡厳密式でt=0とする。）
主鏡厳密式の入力で「m」の符号が逆になっていませんか？

はらぱん>さて副鏡に関してですが、同様にＭａｘｉｍａで解いているのですが上手くいきません。
はらぱん>極座標などというのは全くわからず詰まってしまいました。

副鏡面の極座標が(ρ, u)のとき、直交座標は(ρ*cos(u), ρ*sin(u))。
副鏡厳密式の左辺は1/ρで右辺は(uでなく)tの式で面倒ですが(主鏡厳密式と同様に）
「倍角公式」で　ρ = ＜uの式＞と変形できる。さて、
※主鏡厳密式の場合は、（tを消去して）うまくxがyの関数(陽形式):
※　　＜x = H(y)＞のかたち
※に出来た。そこで Ｍａｘｉｍａの「taylor()」一発でMaclaurin展開できたが
副鏡厳密式の場合は、
x = ρ(u(t))*cos(u(t))
y = ρ(u(t))*sin(u(t))
を 「陽形式」にできるか？　「媒介変数tを消去できるか？」とも言えますね。
(でも... ... デキソウニ アリマセン！) 　そこで方向転換して...

「M」>「我流」で媒介変数表現の展開法をひねりだしました

...ということなのですが、当方のやり方を、いきなり紹介してしまうことは差し控えようと思います。
(Ｍａｘｉｍａを楽しんでつかっていらっしゃる様子ですので:-)

※Ｍａｘｉｍａでの「媒介変数表現のMaclaurin展開」も（うまく探せば）見つかるような
※気がしてきました。　理由は...あらためて。

※なお、
※『媒介変数表示関数の微分法　dy/dx=dy/dt/dx/dt』:-
※https://examist.jp/mathematics/derivation/baikaihensuu-bibunhou/
※などを参考にしました。

[はらぱん]

みなさま　こんにちは。

mの符号、ｔへの代入、Ｍさんのご指摘のとおりでした。
ありがとうございます。
戴いたヒントをもとに、もう少し勉強します。

[はらぱん]

みなさまこんにちは。

Ｍさん、あるいは何方かお助けをお願いします。
RC鏡の副鏡式の展開。
数ⅡＢで赤点をとって数学人生が終わった者には難しい問題です。
Ｍさんに教えていただいたサイトの問題をＭａｘｉｍａで解いて準備体操をし、以下ようなコードを書きました。


/* 1/p の定義 */
inv_p : t/e + (1/m) * ((1 - t/e)^(1/(1-e))) / ((1 - t)^(e/(1-e)));
p:1/inv_p;

/*pのテーラー展開*/
p_t_y:taylor(subst(t=(1 - sqrt(1 - y^2)),p),y,0,8);

/*直交座標 x =*cos(U), y = p*sin(U) に変換 */
/* cos(u) = 1-2t, sin(u) =y */
x_c : p_t_y * (1 - 2*t);
y_c:p_t_y* y;

/*微分*/
xd_dt: diff(x_c,t);
yd_dt: diff(y_c,y);

/*dy/dx*/
dy_dx: yd_dt/xd_dt;

/*d2y/dx2*/
d_dt: diff(dy_dx,y);
inv_xd_dt: 1/xd_dt;
d2y_dx2: d_dt * inv_xd_dt;


コピペで動くと思います。
直交座標に変換したのちの微分で接線が、も一度微分すると接点が求まるという考えでよいのでしょうか。
しかしＹの奇数乗になってしまいます・・・

[「M」]

はらぱん さん 　みなさん　こんにちは

はらぱん>コピペで動くと思います。　...
はらぱん>しかしＹの奇数乗になってしまいます・・・
うまく行っていれば、dx/dy、d2x/dy2,...では(xもyも消えて)変数はtだけになる筈。
副鏡厳密式のyと主鏡厳密式のyが「混線」している感じですね。
また、 （副鏡の）xを（副鏡の）yの関数として展開したい:-
　　 x(y) = As*y^0 + Bs*y^2 + Cs*y^4 + Ds*y^6 + ...
ので(dy/dx,d2y/dx2,...でなく)dx/dy、d2x/dy2,...を求めることになりますね...
そこで

コード: 全て選択

/* 1/p の定義 */
inv_p : t/e + (1/m) * ((1 - t/e)^(1/(1-e))) / ((1 - t)^(e/(1-e)));
p:1/inv_p;

/*直交座標 x = p*cos(U), y = p*sin(U) に変換 */
/* cos(u) = sqrt(1-4*(1-t)*t), sin(u) = sqrt(4*(1-t)*t) */

x_c: p * sqrt(1-4*(1-t)*t);
y_c: p * sqrt(4*(1-t)*t);

/*微分*/
xd_dt: diff(x_c,t);
yd_dt: diff(y_c,t);

/*dx/dy*/
dx_dy: xd_dt/yd_dt;

/*d2x/dy2*/
d_dt: diff(dx_dy,t);
inv_yd_dt: 1/yd_dt;
d2x_dy2: d_dt * inv_yd_dt;

factor(subst(t=0, ratsimp(dx_dy)))/1!;  /*1次の項の係数*/	
factor(subst(t=0, ratsimp(d2x_dy2)))/2!;/*2次の項の係数 = Bs*/

とすれば 「yの1次の項の係数」(=0)、「yの2次の項の係数」が出るようになるはずです。
※なお、pを先にテーラー展開してから dx/dy、d2x/dy2,...を
※組み立てるやり方では、「無視した高次項が結果に影響しないこと」の吟味が必要だと思います。

※※ここで使っているMaximaの機能の(おおざっぱな)説明:-
※※コメント　　/* ... */　　
※※代入　　　<左辺>:<右辺>
※※書き換え　subst(<変数>=<式2>,<式1>)
※※微分　　　diff(<式>,<変数>)
※※簡約　　　ratsimp(<式>)
※※因数分解　factor(<式>)
※※※(その他) 四則「+ - * /」,累乗「^」,平方根「sqrt()」,階乗「!」

「M」は、こんなぐあいに書いています:- ※冗長な記述が多いです。

コード: 全て選択

rho(t):= 1 / (  t/e + (1/m) * ((1 - t/e)^(1/(1-e))) / ((1 - t)^(e/(1-e)))
 		); 

/* t = sin(u/2)**2 */
/* cos(u) = sqrt(1-4*(1-t)*t) */
/* sin(u) = sqrt(4*(1-t)*t) */
				
/*rho(t)*cos(u) */
r_param_xx(t):=''rho(t)*sqrt(1-4*(1-t)*t);

/*rho(t)*sin(u) */
r_param_yy(t):=''rho(t)*sqrt(4*(1-t)*t);

/*rho(t)*sin(u)**2 */
/*
r_param_yy(t):=''rho(t)*sqrt(4*(1-t)*t)**2;
*/

r_param_xxd_1(t):=''(diff(r_param_xx(t),t))$
r_param_yyd_1(t):=''(diff(r_param_yy(t),t))$

r_param_xx_yy_d1(t):= ''(ratsimp(r_param_xxd_1(t)/r_param_yyd_1(t)))$
r_param_xx_yy_d2(t):= ''(ratsimp(diff(r_param_xx_yy_d1(t),t)/r_param_yyd_1(t)))$

r_param_xx_yy_d3(t):= ''(ratsimp(diff(r_param_xx_yy_d2(t),t)/r_param_yyd_1(t)))$

factor(ratsimp(r_param_xx_yy_d1(0)));
factor(ratsimp(r_param_xx_yy_d2(0)));
factor(ratsimp(r_param_xx_yy_d3(0)));

※※Maximaの機能の(おおざっぱな)説明(追加):-
※※関数定義　　<関数名>(<引数名>,...):=<右辺>
※※評価の抑止　「’」


※はらぱん>直交座標に変換したのちの微分で接線が、も一度微分すると接点が求まるという考えでよいのでしょうか。
※とりあえずの答えは
※「ここでは『形式的に式を操作』しているので、図形的な意味は考えていません...」でしょうか？
※※とはいっても......宿題とさせてください...(「M」は「ななめ下」に考え込んじゃう人なので:-)


※※※>...赤点...
※※※「逆マジック点数」というのもいただきました（次回の試験でxx点越えないと「可」に届かない。）
※※※※※「可山優三」は、とてもムリ。セイゼイいって「不可山良二」（=それだれ？）

閑話休題。
　副鏡厳密式の展開では（主鏡の場合と様子が違って）次数を増やした時の
作業メモリーの容量＋計算量の増加がいちじるしいです。「M」の普段づかいの計算機(*1)では
上のやり方でyの次数を0,1,2,...とふやしていくと6次の係数項まで届きませんでした。
「クレチァンもこれ以上進むことができなかったのか...」という状況が、すこし分かった次第です。
さてそこで、また「考え込んじゃった」わけですが...（つづく）

(*1)ほぼ20年前の計算機　＝　CPU: 2.2GHz, RAM: 4GByte。

※なお、次回投稿は週末以降になると思います。

[はらぱん]

みなさま　こんにちは。

Ｍさん、添削ありがとうございます。
早速試してみたところ、正しい式が出力されました。
おっしゃるように計算負荷はかなりのもので、ＣＰＵ-i９、メモリー１６Ｇでｙ⁶まで７分ほどかかります。
ｙ⁷までだと放っておいたところ約１時間。
なのでｙ⁸は諦めました。
しかし、式の誤植と合わせて長年気になっていたことが一応解決され、何とも言えないすっきりとした気持ちです。
Ｍａｘｉｍａも使い始めて３週間ほどなので思うようになりません。
ゆっくりとコードを読んで慣れたいと思います。
本当にありがとうございます。
不可山良二さんは救世主です。

[「M」]

はらぱん さん 　みなさん　こんにちは

はらぱん>ｙ⁶まで７分ほどかかります。
さすがi９は速いですね。　single-threadかつ(ほとんど)整数系演算のみの場合は
クロックの差だけ...ではないようです。
とはいえ、当方には、新しい計算機のための予算がない。　当初
普段はつかわないRAMの多い機械を持ち出し、Lisp処理系を速いという評判のSBCLに切り替え
（MaximaはCommon Lispで書いてある）...明け方までかかってｙ⁶の項を計算しました。
　計算量が、問題の大きさ（この場合は、yの次数を上げていく段数）につき指数関数になっている感じです。
昔、ある先生に「カシコイ算法の工夫は大事だが、問題の大きさを抑える工夫も大事だよ」とおそわり
ましたが、この級数展開の問題では、yの偶数次だけを計算する手があることに気付くのに
時間がかかりました。　その一方『カシコイ算法』も気になる。Schwarzschild,Chrétienは
どうやったのでしょう。　『新版 反射望遠鏡光学入門』 の(9.14)を手がかりにして、いろいろ
やってみて...格段に速く計算できるようになった...と思ったのですが、どうもおかしい。
Maximaをよく理解できていないことを（またまた）思い知らされました。　
この問題ちょっと寝かせて別の方向から取り組もうとおもいます。

　なお、以下の版だと、さほど待たなくても答えがでる筈です。

コード: 全て選択

rho(t):= 1 / (  t/e + (1/m) * ((1 - t/e)^(1/(1-e))) / ((1 - t)^(e/(1-e)))
                );

/* t = sin(u/2)**2 */
/* cos(u) = sqrt(1-4*(1-t)*t) */
/* sin(u) = sqrt(4*(1-t)*t) */

/*rho(t)*cos(u) */
r_param_xx(t):=rho(t)*sqrt(1-4*(1-t)*t);

/*rho(t)*sin(u)**2 */
r_param_yy(t):=rho(t)*sqrt(4*(1-t)*t)**2;

r_param_xxd_1(t):=(diff(r_param_xx(t),t))$
r_param_yyd_1(t):=(diff(r_param_yy(t),t))$

r_param_xx_yy_d1(t):= ((r_param_xxd_1(t)/r_param_yyd_1(t)))$
r_param_xx_yy_d2(t):= ((diff(r_param_xx_yy_d1(t),t)/r_param_yyd_1(t)))$
r_param_xx_yy_d3(t):= ((diff(r_param_xx_yy_d2(t),t)/r_param_yyd_1(t)))$
r_param_xx_yy_d4(t):= ((diff(r_param_xx_yy_d3(t),t)/r_param_yyd_1(t)))$
r_param_xx_yy_d5(t):= ((diff(r_param_xx_yy_d4(t),t)/r_param_yyd_1(t)))$
r_param_xx_yy_d6(t):= ((diff(r_param_xx_yy_d5(t),t)/r_param_yyd_1(t)))$

define(xx_yy_d1(t), ((r_param_xx_yy_d1(t))))$
define(xx_yy_d2(t), ((r_param_xx_yy_d2(t))))$
define(xx_yy_d3(t), ((r_param_xx_yy_d3(t))))$
define(xx_yy_d4(t), ((r_param_xx_yy_d4(t))))$
define(xx_yy_d5(t), ((r_param_xx_yy_d5(t))))$
define(xx_yy_d6(t), ((r_param_xx_yy_d6(t))))$

ratsimp(rho(0));
ratsimp(xx_yy_d1(0));
ratsimp(xx_yy_d2(0));
ratsimp(xx_yy_d3(0));
ratsimp(xx_yy_d4(0));
ratsimp(xx_yy_d5(0));
ratsimp(xx_yy_d6(0));

※比べると、前の版での「悪あがき」が目立ちますね...

[はらぱん]

Ｍさん　みなさん　こんにちは。

Ｍさんに教えていただいた、ＲＣ鏡副鏡の厳密式を繰り返し微分するのは高階微分というものだと知りました。
厳密式直接の微分は計算負荷が大きくｙ＾6までしか計算できませんでしたが、厳密式をテーラー展開したものを微分したところ、割と短時間でもう少し高次の係数を求めることが出来ました。
厳密式を１２次まで展開した式は最高３０回程度微分が出来ました。
誤差が気になります。
以下のように厳密式と高階微分による展開式に数値を入れて確かめてみました。

コード: 全て選択

/*RC副鏡高階微分による展開式の誤差検証*/
e:0.2525;   /*面間隔*/
m:0.326667;  /*副鏡から焦点位置*/
U:20*(%pi/180);  /*ヴァージェンスアングル(ラジアンに変換)*/
/*t:sin(U/2)^2;*/
y:sin(U);

inv_p : t/e + (1/m) * ((1 - t/e)^(1/(1-e))) / ((1 - t)^(e/(1-e)));    /*厳密解*/
p:1/inv_p;
/*直交座標 x = p*cos(U), y = p*sin(U) に変換 */
/* cos(u) = sqrt(1-4*(1-t)*t), sin(u) = sqrt(4*(1-t)*t) */
x_c: p * sqrt(1-4*(1-t)*t)$
y_c: p * sqrt(4*(1-t)*t)$
X_E:float(ev(x_c,t:sin(U/2)^2))$     /*角度Ｕの時の厳密解によるｘ*/
Y_E:float(ev(y_c,t:sin(U/2)^2))$      /*角度Ｕの時の厳密解によるy*/

/*テーラー展開*/
x_c_ty: taylor( (p * sqrt(1-4*(1-t)*t)),t,0,12 ) $    
y_c_ty: taylor(  (p*sqrt(4*(1-t)*t)),t,0,12)   $
X_T:float(ev(x_c_ty,t:sin(U/2)^2))$  /*角度Ｕの時のテーラー展開式によるｘ*/
Y_T:float(ev(y_c_ty,t:sin(U/2)^2))$  /*角度Ｕの時のテーラー展開式によるy*/

/*高階微分による展開式、有効と思われる次数まで*/
X: m+1.2755078922378549*Y^2
-11.518947586175601*Y^4
+174.53446111285976*Y^6
-3262.2463269834407*Y^8
+67892.03149531294*Y^10
-1509805.1530738475*Y^12
+3.512282857830617·10^7*Y^14
-8.441722091157774·10^8*Y^16
+2.079776770709023·10^10*Y^18
-5.22432006704646·10^11*Y^20
+1.3330023635674469·10^13*Y^22
-3.4452229299281925·10^14*Y^24$
/*+9.00078062129218·10^15*Y^26
-2.3730995418972886·10^17*Y^28
+6.306217105212293·10^18*Y^30
-1.6873002802071778·10^20*Y^32
+4.5417362626753976·10^21*Y^34
-1.2290125299802198·10^23*Y^36
+3.3415176620517266·10^24*Y^38$*/

/*float(y);*/
float(Y_E);
float(Y_T);
print("ヴァージェンスアングル=",float(U/(%pi/180)))$
print("厳密解　　　　=",float(X_E) )$    /*角度Ｕの時の厳密解によるｘ*/
print("テーラー展開式　 =",float(X_T) )$   /*角度Ｕの時のテーラー展開式によるｘ*/
/*角度Uの時の厳密解のyの値を高階微分による展開式に入れる*/
y_h:float(Y_E)$
print("高階微分展開式=",float(ev(X,Y:y_h)) )$

ヴァージェンスアングル（Ｕ)を順次書き換えてみると、２０°程度が限界で、以降急速に誤差が大きくなります。
F値でいうとF1.6あたりです。
厳密式からの展開次数を６０次にしてから微分してみましたが精度は上がりませんでした。
同様に主鏡の展開式の精度も確かめて見ました。
ｙ＾60までの展開式で、ヴァージェンスアングル５０°くらいが限界でした。

モデルケースにしているｅ=0.2525,ｍ=0.326667で有効と思われる係数を入力、口径を少し大きくしては、係数を１項ずつ増やして収差を確認するシミュレーションをしました。
初めは口径を増やす（Fを明るくする）と増える収差は、係数の項を増やせば補正されていましたが、ｙ＾22の項以降は変化しなくなりました。
結局Ｆ1.6程度が限界でした。
テーラー展開した式の高階微分によって得られた展開式の係数で、ｙ＾22以上は意味のない数値のようです。

以下光路図と縦収差曲線です。
焦点距離１では分かりづらいので２０００にスケーリングしています。

主鏡径：1200
副鏡径：500
面間隔：504.9970
焦点距離：1999.9999
バックフォーカス：653.3343
Ｆ：1.6666
光路図の最も内側の光線がＦ10に相当します。

RC鏡を初めて光線追跡したとき、Ｆを明るくすると収差が急速に増えるのがわかりました。
円錐係数による曲面のせいだと思いましたが、クレチァンの面形状の展開式がｙ＾４までしかわからなかったので、どうすることもできませんでした。
今回ようやく高次の項の算出法がわかり、明るい口径比での光線追跡が出来ました。
しかし展開式次数が全然足らないようです。
『写真レンズの科学』ではF0.54の非球面レンズは190次の非球面と書いてありました。
RC鏡もそのくらい必要なのかもしれません。
『光学機器大全』の中の非球面レンズの章で、吉田先生が級数展開ではせいぜいＦ1程度までと書かれているのはこういうことなのでしょうか。
面形状を展開式表す方法では、限界があるようです。
今のところこんな具合で暫く停滞しそうです。
「アマチュアは無駄なことをするものだと呆れるばかりだ」といったところです。