ペルチェ素子霧箱の詳しい説明を書きました。長くなったので別ファイルにしています。
PeltierCloudChamber.html
（追記）以下、コメントにお答えします。
ペルチェ素子霧箱の作り方を投稿しました。
ペルチェ霧箱のかんたんな製作方法
ニコニコ動画に、ペルチェ素子霧箱の作り方を投稿しました。詳しい説明は数日中にこのブログに書きますが、一つだけ。動画で説明したペルチェ素子を2枚重ねる方法は、結構条件がシビアです。確実に現象を見たい場合は、19012-5L31-06CQQという2段重ねた状態で売られているペルチェ素子に、12Vをかけてみてください。
コミックマーケットC81当選しました
コミケですが、当選しました。
3日目 東地区 R – 51 a http://twitcmap.jp/?id=0081-3-RRa-51-a
また、本のタイトルも決まっています。
「最後の歯車式計算機クルタ」クルト・ヘルツシュタルク氏インタビュー
です。
こんな感じの本です。
クルタという歯車式計算機の、開発者に対するインタビューを、許可を得て翻訳しました。
この本を読むと、↓がわかります。
・人類初の表計算機能を実現した、メカニカルなシステムとは？
・「00」や「000」のキーを発明したのは誰？
・足し算だけで引き算を実現する方法とは？
・強制収容所の中でクルタ計算機を設計できたのはなぜ？
・クルタ計算機を、リヒテンシュタインという小国で生産した理由は？
2日間に渡ったインタビューは自伝的要素を多分に含んだ色濃い内容で、19世紀後半～20世紀前半の計算機市場の実情から、計算機設計のアイデア、ユダヤ人差別、ナチス収容所内での計算機設計、アメリカによる開放とドイツからオーストリアへの逃避行、リヒテンシュタイン皇室によるサポート、そして電卓の登場による機械式計算機の衰退まで、読んでいて飽きることはありません。当時のヨーロッパの現実を伝える第一級の史料となっています。計算機に興味がある方、ヨーロッパの歴史に興味がある方、是非お読み下さい。
コミケで本を出します
2年ほど間が空いてしまいました。
突然ですが、サークル「ぴっちぶれんど」としてコミケで本を出します。
「クルタ計算機」という機械式計算機を開発したクルト・ヘルツシュタルク氏へのインタビュー記事(pdf)を許可を得て日本語に翻訳したものです。サークルの紹介ページにも書いていますが、内容はこんな感じです。
2日間に渡ったインタビューは自伝的要素を多分に含んだ色濃い内容で、19世紀後半～20世紀前半の計算機市場の実情から、計算機設計のアイデア、ユダヤ人差別、ナチス収容所内での計算機設計、アメリカによる開放とドイツからオーストリアへの逃避行、リヒテンシュタイン皇室によるサポート、そして電卓の登場による機械式計算機の衰退まで、読んでいて飽きることはありません。当時のヨーロッパの現実を伝える第一級の史料となっています。計算機に興味がある方、ヨーロッパの歴史に興味がある方、是非お読み下さい。
サークルカットはこんな感じで作ってみました。
コミケの出店は抽選だそうで、当選するかどうかわかりませんが、落選した場合もどこかに委託して出したいと思っています。みなさん、よろしくお願いします。
赤外線写真の世界
相変わらず年に1、２回しか更新されないこのブログにお越しいただきありがとうございます。
最近、赤外線写真でいろいろと遊んでいます。普通のカメラには赤外線をカットするフィルターが入っていますが、それを取り外して逆に赤外線しか通さないフィルターを取り付けると赤外線カメラの出来上がりです。
最近は、赤外線だけではなく可視光を少し通すフィルターを取り付けて撮影し、疑似カラーの赤外線写真を撮るのがブームなようで、わたしも630nmフィルターで撮影しています。撮影後のレタッチは必須で、特に赤チャンネルと青チャンネルを入れ替えて空の色を青に合わせたりします。スノー効果で樹木が真っ白に写るのが特徴です。
下は、そうして作った写真です。非日常的な世界をお楽しみください。どれも雪景色ではありませんよ。
「最適化数学で新年の挨拶」の解説
みなさん、あけましておめでとうございます。
新年開始とともに、ニコニコに動画を投稿しました。
おおむね好評だったようで嬉しい限りです。コメントをくださった方、ありがとうございます。
もっと詳しく説明しろとのコメントがいくつか寄せられていたので、解説します。2と9の文字は、両方とも4つの点とその周りを公転する4点、更にそれを回る4点を元に作られます。以下、それぞれを恒星、惑星、衛星と呼ぶことにします。4衛星を2組にわけ2本の直線を作り、その交点に文字を描かせます。
それぞれの点の座標は、時刻をtとして次のように表しました。
p1=(a1*cos(2(2t-e1)),a1*sin(2(2t-e1)))+p5
p2=(a2*cos(2(2t-e2)),a2*sin(2(2t-e2)))+p6
p3=(a3*cos(2(2t-e3)),a3*sin(2(2t-e3)))+p7
p4=(a4*cos(2(2t-e4)),a4*sin(2(2t-e4)))+p8
p5=(a5*cos(3(t-e5)),a5*sin(3(t-e5)))+(c5,d5)
p6=(a6*cos(3(t-e6)),a6*sin(3(t-e6)))+(c6,d6)
p7=(a7*cos(t-e7),a7*sin(t-e7))+(c7,d7)
p8=(a8*cos(t-e8),a8*sin(t-e8))+(c8,d8)
(c5,d5)～(c8,d8)が4恒星の座標、p5～p8が4惑星の座標、p1～p4が4
衛星の座標となります。a1～a8は公転半径の大きさ、e1～e8は公転での位相を決定します。またp1とp2を通る直線、p3とp4を通る直線を結びその交点を求め、これを点pとします。pの座標は、複雑な形になりますが方程式を解けば式として表すことが可能です。a*,c*,d*,e*の合計24変数を調整し、pが文字を描くようにすれば良いわけです。
なおそれぞれの式のtには1,2,3の係数が掛かっていますが、これは各点が回転する速度を表します。t=2Piで全ての点が元に戻るようにしたかったので、係数は整数である必要があります。この値は、実は手動で調整しました。手で描くとわかるのですが、「2」も「9」も一度普通になぞって逆方向もなぞり元の位置に戻ると、X軸方向は3往復、Y軸方向は1往復します。このことから縦に動く線とその3倍の速さで横に動く線ができることを期待してp5～p8でのtの係数、すなわちは惑星の公転速度は1および3と決定しました。これは上手く働いたようです。p1～p4でのtの係数、これは衛星の公転速度になりますが、いろいろと試行した末に上の形になっています。今見ると変な係数の付け方で計算していますが、どうせ最適化すれば同じ式になる（ことが期待できる）ので気にしなくて良いかと思います。
さて、pにどのような軌跡を描かせたいのかを決めます。これはフリーハンドで描いてその座標を読みました。まず「2」はこんな感じです。
t=0Pi (-0.437739, 0.64179)
t=1Pi/28 (-0.286784, 0.88061)
t=2Pi/28 (-0.0754772, 0.940305)
t=3Pi/28 (0.166043, 1.)
t=4Pi/28 (0.437739, 0.940305)
t=5Pi/28 (0.588657, 0.910438)
t=6Pi/28 (0.679222, 0.761181)
t=7Pi/28 (0.739611, 0.611923)
t=8Pi/28 (0.709435, 0.432838)
t=9Pi/28 (0.588657, 0.223886)
t=10Pi/28 (0.528304, 0.164191)
t=11Pi/28 (0.407526, -0.0149336)
t=12Pi/28 (0.226395, -0.164191)
t=13Pi/28 (0.0150882, -0.283581)
t=14Pi/28 (-0.166043, -0.432838)
t=15Pi/28 (-0.316997, -0.582095)
t=16Pi/28 (-0.467915, -0.731353)
t=17Pi/28 (-0.649046, -0.850743)
t=18Pi/28 (-0.8, -0.970133)
t=19Pi/28 (-0.679258, -1.)
t=20Pi/28 (-0.467915, -1.)
t=21Pi/28 (-0.256608, -0.970133)
t=22Pi/28 (-0.0754772, -0.970133)
t=23Pi/28 (0.13583, -0.940305)
t=24Pi/28 (0.286784, -0.940305)
t=25Pi/28 (0.437739, -0.970133)
t=26Pi/28 (0.618869, -0.970133)
t=27Pi/28 (0.739611, -0.940305)
t=28Pi/28 (0.8, -0.910438)
t=29Pi/28 (0.739611, -0.940305)
t=30Pi/28 (0.618869, -0.970133)
t=31Pi/28 (0.437739, -0.970133)
t=32Pi/28 (0.286784, -0.940305)
t=33Pi/28 (0.13583, -0.940305)
t=34Pi/28 (-0.0754772, -0.970133)
t=35Pi/28 (-0.256608, -0.970133)
t=36Pi/28 (-0.467915, -1.)
t=37Pi/28 (-0.679258, -1.)
t=38Pi/28 (-0.8, -0.970133)
t=39Pi/28 (-0.649046, -0.850743)
t=40Pi/28 (-0.467915, -0.731353)
t=41Pi/28 (-0.316997, -0.582095)
t=42Pi/28 (-0.166043, -0.432838)
t=43Pi/28 (0.0150882, -0.283581)
t=44Pi/28 (0.226395, -0.164191)
t=45Pi/28 (0.407526, -0.0149336)
t=46Pi/28 (0.528304, 0.164191)
t=47Pi/28 (0.588657, 0.223886)
t=48Pi/28 (0.709435, 0.432838)
t=49Pi/28 (0.739611, 0.611923)
t=50Pi/28 (0.679222, 0.761181)
t=51Pi/28 (0.588657, 0.910438)
t=52Pi/28 (0.437739, 0.940305)
t=53Pi/28 (0.166043, 1.)
t=54Pi/28 (-0.0754772, 0.940305)
t=55Pi/28 (-0.286784, 0.88061)
また「9」はこうなりました。
t=0Pi (0.65, 0.830183)
t=Pi/21 (0.273674, 0.924517)
t=2Pi/21 (0.034217, 1.)
t=3Pi/21 (-0.273674, 1.)
t=4Pi/21 (-0.444739, 0.924517)
t=5Pi/21 (-0.65, 0.716971)
t=6Pi/21 (-0.615804, 0.547154)
t=7Pi/21 (-0.513173, 0.396214)
t=8Pi/21 (-0.205282, 0.377363)
t=9Pi/21 (0.034217, 0.433969)
t=10Pi/21 (0.239478, 0.547154)
t=11Pi/21 (0.444739, 0.667552)
t=12Pi/21 (0.513132, 0.735849)
t=13Pi/21 (0.615762, 0.84906)
t=14Pi/21 (0.513132, 0.547154)
t=15Pi/21 (0.478935, 0.45282)
t=16Pi/21 (0.478935, 0.30188)
t=17Pi/21 (0.410501, 0.0188773)
t=18Pi/21 (0.273674, -0.226423)
t=19Pi/21 (0.136848, -0.509426)
t=20Pi/21 (-0.0000205261, -0.811306)
t=21Pi/21 (-0.136848, -1.)
t=22Pi/21 (-0.136848, -1.)
t=23Pi/21 (-0.0000205261, -0.811306)
t=24Pi/21 (0.136848, -0.509426)
t=25Pi/21 (0.273674, -0.226423)
t=26Pi/21 (0.410501, 0.0188773)
t=27Pi/21 (0.478935, 0.30188)
t=28Pi/21 (0.478935, 0.45282)
t=29Pi/21 (0.513132, 0.547154)
t=30Pi/21 (0.615762, 0.84906)
t=31Pi/21 (0.513132, 0.735849)
t=32Pi/21 (0.444739, 0.667552)
t=33Pi/21 (0.239478, 0.547154)
t=34Pi/21 (0.034217, 0.433969)
t=35Pi/21 (-0.205282, 0.377363)
t=36Pi/21 (-0.513173, 0.396214)
t=37Pi/21 (-0.615804, 0.547154)
t=38Pi/21 (-0.65, 0.716971)
t=39Pi/21 (-0.444739, 0.924517)
t=40Pi/21 (-0.273674, 1.)
t=41Pi/21 (0.034217, 1.)
改めてみると字が下手すぎですね。フリーハンドで描いたのは、別に味のある文字を使いたかったとかいうのではなく、単に時間が無かったからに他なりません。
時刻tにおけるpの座標をp(t)とします。「2」の場合はt=0～55Pi/28、「9」の場合はt=0～41Pi/21まで変化させ、上で求めた各点とp(t)との距離の2乗をそれぞれの時刻で計算し、合計します。この式は24変数を含むスカラーになることがわかります。この巨大な式を評価式として最小化すれば良いわけです。ただ評価式が複雑すぎて解析的な解は求まらないかと思います。そこで最小化には遺伝的アルゴリズムと滑降シンプレックス法を組み合わせた自作のプログラムを使用しました。なお滑降シンプレックス法は、名前こそ似ていますが線形計画法で登場するシンプレックス法とは直接の関係は無いので注意が必要です。
結果はこのようになりました。
「2」の場合
a1=-0.281248,a2=-0.217597,a3=-0.133595,a4=0.135594,
a5=0.233587,a6=-0.533821,a7=-0.865466,a8=-0.774918,c5=0.296162,
c6=0.128995,c7=-0.0170449,c8=-0.310357,d5=1.06265,d6=-0.197444,
d7=-0.199307,d8=-0.308207,e1=0.311494,e2=0.218138,e3=0.0167031,
e4=1.52846,e5=0.952566,e6=0.0411405,e7=1.25651,e8=1.3886
「9」の場合
a1=-0.197302,a2=0.145299,a3=-0.197725,a4=-0.200862,
a5=-0.524991,a6=0.779246,a7=-0.472475,a8=1.29185,c5=0.117048,
c6=0.560647,c7=0.611663,c8=-0.875045,d5=0.415859,d6=-0.461184,
d7=1.32743,d8=-0.340324,e1=1.35564,e2=-0.552802,e3=-1.25281,
e4=0.787604,e5=0.866171,e6=-0.150847,e7=-0.496076,e8=4.18745
円の半径a*の値が負になったりしていますが、p1～p8の式のまま素直に描画すれば良いので特に問題はありません。最適化の限界で、またそもそも変数の数の限界から、元の点とぴったり同じ軌跡を描くことはできません。できあがった動画を見ても、往復で軌跡がずれているのがわかります。なお、円の数を増やせばもっと綺麗な文字の形になるのかとのコメントがありましたが、確かに綺麗になる可能性はありますが変数の数が増えるので探索空間が爆発的に増大し最適化が困難になるかもしれません。この辺りのさじ加減にはなかなか難しいものがあります。
以上が「2」と「9」の描き方です。「0」については楕円で代用しています。あのように棒を動かすとその棒上の点が楕円を描くことは簡単に証明可能です。
好評だったので来年の元旦にもまた何か投稿しようと思います。今度は遺伝的プログラミングなどを使って、どの円とどの円を組にするかなども含めて最適化すると面白いかな、と思っています。あるいは全く違う生物的なネタも考えています。ただ、「2010」という文字が簡単すぎて今ひとつ創作意欲をかき立てられないですね。まだ1年間あるので、いろいろと考えて見ることにします。
ご無沙汰しております
ご無沙汰にもほどがありますが、とりあえず生存報告ということで。
いろんなことを放ったらかしにしたままですみません。
最近機械式計算機に凝っていて、リヒテンシュタイン製のクルタ計算機をいくつか手に入れました。分解してみるとその巧みな機構に感嘆するばかりです。この計算機が生まれた過程もとても興味深い。
そのうち、ニコニコにクルタ関係の動画を投稿するかもしれません。もしかすると冬コミで同人誌作ろうか、という話もあります。
期待せずに待っていてください。
2007 WD5 は本当に火星にぶつからないのか その3（再々計算）
これの更に続きです。WD5の位置計算は、太陽とWD5の2体問題でやっていたのですが、今回は火星がWD5に及ぼす引力の影響も考慮した微分方程式で計算してみました。
20000モンテカルロですが、うーん、結果は変わりませんね。平均、標準偏差もほぼ同等です。火星に接近するケースが増えるかと思ったのですが。ちなみに正規分布を仮定して平均・標準偏差から衝突確率を求めると、約0.0031%となります（前回も）。
なぜ音階は「ドレミファソラシド」なのか
ピアノには白黒合わせて1オクターブに12の鍵盤がありますが、小学校で最初に習う曲は白鍵の「ドレミファソラシ」だけを使って演奏できるものばかりで、黒鍵はあまり使われません。なぜ12音中この7音が選ばれたのか、私にとってはずっと謎でしたが、最近ブルーバックスの「音律と音階の科学」（小方厚著）を読んでその理由がわかったのでまとめてみようと思います。要はこの本の紹介をやりたいのですが、できるだけ丸写しにならないように配慮はしている積もりです。ちゃんと勉強したい人は、是非ブルーバックスを読んでください。多数の図や表が使われていて理系の人間にとってもわかりやすい解説ですし、音律、音階の歴史的な経緯や今後を見据えた発展などが面白く、とても楽しめました。
結論から言うと、綺麗なハーモニー（和音）を出すためにこの音階が決められているようです。これを解説してみます。
まず、臨界帯域幅という用語があるようです。二つの純音を同時に鳴らして、その周波数差を徐々に広げていきます。それらが独立した2音として聞こえ、うなりも、ゴロゴロするような不快な音も聞こえなくなった時の周波数差を臨界帯域幅と言うようです。
臨界帯域幅をグラフにするとこうなります。元のグラフは、J.G.Roederer “The Physics and Psychophysics of Music” に載っているようです。
これは、多くの被験者に心理的な実験を行って得られたグラフなのでしょう。見ての通り、線形ではありません。近似式を作ってみました。数値計算上、高速で精度良く計算できるように式を変換しています。
(97.27500666976535 + x*(-0.10438654962335651 + x*(0.00013795361645053944 + x*(-7.50643237852489e-8 + (3.864047759875421e-11 – 2.715433778943346e-15*x)*x))))/
(1 + x*(-0.0009358884960842764 + x*(4.5256275668615123e-7 + (3.595960994872951e-11 – 2.546811069035789e-15*x)*x)))
次に、2つの純音を鳴らした時の不協和度のグラフがあります。これも心理学的実験から得られたもので、被験者がどれぐらいの違和感を感じるのかをテストして作られています。
近似式はこうなります。
(-0.00020410890928780975+x*(11.293912104998519 + x*(-64.80401525171715 + x*(142.52964847415188 + x*(-128.65558599526295 + 40.28982394173413*x)))))/
(1 + x*(-1.1833111299032766 + x*(-4.205411534000957 + x*(0.8554237901431951 + 21.877444099945798*x))))
横軸は純音の周波数差なのですが、先ほどの臨界帯域幅に対する比となっています。従って、周波数によって横軸のスケールが変わることになります。こちらのグラフは、W.A.Sethares, “Tuning, Timbre, Spectrum, Scale” 2nd. ed に元図があります。
実際にC4（ド）音付近の2つの純音を聴いてみましょう。
・周波数差が臨界帯域幅の0.275倍
グラフでは、最も不協和度が高いところです。うなりがはっきり聞こえて違和感がありますね。
・周波数差が臨界帯域幅の0.6倍
うなりが消えて、2音の平均周波数の単音のみが聞こえます。ただし、ゴロゴロ感があります。
・周波数差が臨界帯域幅の0.8倍
そろそろ音が2つ聞こえるようになります。ゴロゴロ感が少しだけ残っています。
・周波数差が臨界帯域幅の1倍
より強く、2音と認識できるのではないでしょうか。違和感がかなり少なくなりました。
・周波数差が臨界帯域幅の1.2倍
ここまで来ると、和音の違和感が全くないはずなのですが、どうでしょう。
もちろん、人によって感じ方に差はあるかと思います。
以上は純音の話でした。これだけなら話は簡単なのですが、実際のピアノの音は2倍音、3倍音といった倍音が豊富に含まれています。ピアノだけではなく他の楽器や人の歌声も同様です。この倍音を含めた不協和度を考える必要があります。2つの鍵盤から出る倍音を計算し、それらの2音の組み合わせ全てを求めてその不協和度を足し算すれば良いわけです。S1(x),S2(x)をそれぞれ鍵盤1,鍵盤2から出る音の振幅（周波数xの関数）、U(x,h)を周波数xでh離れた音との不協和度とすると、全ての倍音を含めた不協和度は次に示すような畳み込み積分となります。ブルーバックスには振幅をどう使うのか明確な記述が無かったので、2つの音の振幅をかけ算してみました。
実際には、馬鹿正直に積分する必要はなくて、離散的に足せば良いでしょう。ド(C)の音を基準とした不協和度のグラフはこうなります。P109図37です。ピアノの倍音の大きさに関しては、ブルーバックスに載っている記述と同じ設定を使っています。
不協和度が小さいほど協和している音なので、注目すべきは下に凸の点です。このボトムの位置が、平均律ではおおよそミのフラット(E♭)、ミ(E)、ファ(F)、ソ(G)、ラ(A)の周波数となります。この「おおよそ」というのが重要で、ピッタリの位置ではないわけです。正確にボトムの周波数で音階を作るのが純正律なわけですが、後に述べるようにいろいろと不都合があって、現在はややずれがある平均律がメジャーに使われています。平均律は、図で言えば横軸が100で割り切れる位置にある音階です。
いよいよ問題の核心部分です。ミのフラットはすぐ隣にミがあるのでちょっと横に置いて（短音階で使われます）、ボトムにミ、ファ、ソ、ラが出てきました。グラフを見ると、ドと最も協和する音はソのようです。ではそのソと協和する音は何かを探します。グラフを平行移動すれば良いわけです（厳密には臨界帯域幅が変わるので少し形が変わります）。すると、シ(B)、ド(C)、レ(D)、ミ(E)が出てきて、これでドレミファソラシドの全てが揃いました。これが、音階がドレミファソラシドである理由です。
次に、グラフでボトムの位置を数値計算で求めてみました。横軸はセントと言って100セントが平均律の半音階に当たる単位です。これを使って書くと、315.64, 386.31, 498,04, 701.96, 884.36にボトムがあります。これを1200で割って2の肩に置くとド音との周波数比がでてきて、それぞれ1.2000=6/5, 1.2500=5/4, 1.3333=4/3, 1.5000=3/2, 1.6667=5/3となり、見事に純正律で使われる音階と一致します。当然、この音階を使った方が良く協和するハーモニーが奏でられるのですが、それぞれの間隔が均等でないのがネックになります。特に移調や転調を行った時に、曲調ががらっと変わってしまい、あまり使われなくなってしまったようです。もっとも、どれぐらい酷いことになるのか実際の音を私は聞いたことが無いので、何とも言えませんが。
ちなみに不協和度のグラフですが、ドを基音として3和音のグラフを描くとこうなります。P117の図41を再現してみました。やはり、グラフのボトムの位置に現れる和音が、良く使われるようです。
2007 WD5 は本当に火星にぶつからないのか その2（再計算）
前回の日記の続きです。
すみません、計算の途中で足し算と引き算を間違えていました。
改めて、正しい方法で再計算しました。モンテカルロ部分のみの修正なので、中心ケースの最接近時刻・距離は変わりません。
CGはこんな感じです。誤差を考慮した100シミュレーションの軌跡を細い線で表示しています。
モンテカルロは、50000回に増やしてみました。結果はこうなります。
1000kmごとのヒストグラムも作ってみました。
正規分布っぽい結果で良い感じです。このうち火星に衝突するのは2ケースでした。2/50000=0.004%ということでNASAの発表0.01%と比べると少ないのですが、ケースが少ないので何とも言えません。マルコフ連鎖モンテカルロ法とかありますが、知りたいケース付近を集中的に試行するのでこういう場合に使われるのでしょうね。いずれにせよ、NASAの発表と矛盾は無いと思われます。ただ、2体問題の計算しかやっていないので、若干距離が遠くなる方向に計算が偏っているはずです。次に同じような天文ショーがあるときは、N体問題で解きたいと思います。
というわけで、NASA発表の追試は終了です。NASAの発表文を読むと、なんとなくモンテカルロじゃなくて解析的に確率を求めているような気がします。どんな式を使っているのでしょうか。
結論
WD5 2007 が火星にぶつかる可能性はほとんどありません。
最近のコメント