Maximaでレンズ配置の設計

プログラミング

変数がn個，制約となる式がm個あったら，任意に選んだn-m個の変数を残りのm個の関数として書ける（当たり前だが）．これをやってくれるのがsolve．ここでは以下のような光ファイバ・カップリング用の，ビーム径と広がり角を調整する光学系 ("telescope") の設計に利用している．中学生でも知ってる「レンズの公式」ベースですが．

それにしてもMaximaにはfoldr/foldl, reduce程度のリスト操作関数も用意されていないとは... Common Lispレベルに降りて書くのが普通なのかなぁ．

wxMaxima 0.7.6 http://wxmaxima.sourceforge.net
Maxima 5.16.3 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (aka GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.

(%i1) all_true(l):= not member(false,l)$

(%i2) all_positive(l):=all_true(map(lambda([s], is(rhs(s)>0)), l))$

(%i3) fpprintprec:3$ ratprint:false$

(%i5) filter_positive_sln(l):=sublist(float(l),all_positive)$

(%i6) filter_non_empty_sln(ll):=sublist(flatten1(ll), lambda([l], is(length(l)>2)))$

(%i7) reduce(f,l):=block([acc:first(l)], (for e in rest(l) do acc:f(acc,e), acc))$

(%i8) flatten1(l):=reduce(append,l)$

(%i9) block([a:-2200,D0:1.36,g:400,D3:1.54,l:300,lenses:[25.4,35,40,100,125,150,175,200,300]],filter_non_empty_sln(outermap(lambda([f1,f2],append([f1,f2],
  filter_positive_sln(solve([1/(a+b)+1/c=1/f1, 1/d+1/(e+g)=1/f2, (a+b)/a * d/c * g/(e+g) * D0 = D3, b+c+d+e=l], [b,c,d,e])))),lenses,lenses)));

(%o9) [[25.4,35,[b=141.,c=25.1,d=37.7,e=96.2]],[25.4,40,[b=39.2,c=25.1,d=42.9,e=193.]],[35,40,[b=210.,c=34.4,d=44.3,e=11.7]],[100,
100,[b=9.38,c=95.6,d=127.,e=67.8]]]

(%i15) block([a:-2200,D0:1.36,g:400,D3:1.54,l:320,lenses:[25.4,35,40,100,125,150]],filter_non_empty_sln(outermap(lambda([f1,f2],append([f1,f2],
  filter_positive_sln(solve([1/(a+b)+1/c=1/f1, 1/d+1/(e+g)=1/f2, (a+b)/a * d/c * g/(e+g) * D0 = D3, b+c+d+e=l], [b,c,d,e])))),lenses,lenses)));

(%o15) [[25.4,35,[b=167.,c=25.1,d=37.7,e=90.5]],[25.4,40,[b=65.9,c=25.1,d=42.9,e=186.]],[35,40,[b=234.,c=34.4,d=44.4,e=7.18]],[100,
100,[b=32.9,c=95.6,d=127.,e=64.0]]]

(%i14) block([a:800,D0:0.73,d:400,D2:1.54,lenses:[25.4,35,40,100,125,150,175,200,300]],filter_non_empty_sln(outermap(lambda([f1,f2],append([f1,f2],
  filter_positive_sln(solve([1/(a+b)+1/(c+d)=1/f1, (a+b)/a * D0 = (c+d)/d * D2], [b,c])))),lenses,[0])));

(%o14) []

(%i13) block([a:800,D0:0.73,g:400,D3:1.54,l:260,lenses:[25.4,35,40,100,125,150,175,200,300]],filter_non_empty_sln(outermap(lambda([f1,f2],append([f1,f2],
  filter_positive_sln(solve([1/(a+b)+1/c=1/f1, 1/d+1/(e+g)=1/f2, (a+b)/a * d/c * g/(e+g) * D0 = D3, b+c+d+e=l],[b,c,d,e])))),lenses,lenses)));

(%o13) []