［kc5party］おやきーずと打ち上げ花火のお目出度い関係

お久しぶりです。変態猫もといきまぐれ猫です。
今回はおやきーずの片割れとして書くので、変態猫ですかね。

高専カンファレンスが5歳の誕生日ということで、またまたおやきーずを再起動させて、LTをしてきました。今回は割りとトークはオプションで、メインは自動花火打ち上げ機の方だったんですが、当日に落ちるというまさかの展開が・・・ね・・・・・・。会場で花火の仕組みについても質問されたりしたので、例のごとく反省と仕様を織り交ぜて記事にしてみようと思います。

例のごとくAtom君の後塵を拝したので、彼のエントリーからペタペタしますね。
忘れられないパーティにしようぜ。 #kosenconf via @FromAtom

出来たもの

Github見てください。動かなかったらIssuesに投げてください。（丸投げ
https://github.com/FromAtom/kosen-firework:kosen-firework

何がしたかったの？

Atom君が何か開発したいと主張したので、それに乗ってみました
お祝いということで、花火を上げようと思いました、まる
これまでの開催履歴をバックに、花火がどんどんと上がるイメージです
花火はTwitterのポストに反応して打ちあがります
花火の色は投稿者のアイコンの特徴色っぽいものを使ってます（あまり精度は良くない
ほんとは五分できれいに終わる予定だったんだ！よ！！

何したの？

猫が作ったのは花火のモジュールです。
TwitterAPIの制御とアイコンの丸加工、バックの開催地表示制御（大切！）の辺りはAtom君が作成しました。
あ、開催地情報の打ち込みも僕がしました。カンファWikiの表からデータ引っ張ってきて、Excel先生と戯れながらサクッとデータの用意を完了しました。表計算ソフトは偉大です。マジで。

懺悔を聞こうじゃないか？

特徴色モジュール、時間が無かったので手抜きしましたすいません。あとシミュレーションもオイラー法で超適当です。ホントにすいません。

誇れることは無いのかね？

モジュールの設計は悪くない、はず。
それと、花火の発光パターンを時間tによる関数にしてあります。発光パターンの追加が簡単！！
もし覚えていて、気力があれば、火の玉に対する加速度も時間tの関数にすると、もっと色んな花火が作れるかなぁと狙ってます。

花火の話をしよう

会場で花火の仕組みついて質問された（超嬉しい！！）ので、それについて少しだけ書こうかと思います。
最初にざっくり言ってしまうと、火の玉一つ一つを数値計算でシミュレーションしてるだけだったりします。

花火らしさとは？

花火が花火らしく見えるファクターは何でしょうか？取り敢えず火の玉が中央から周辺に向かって球を形造りながら飛んでいくのは必須ですね。それと重力加速度のシミュレーションも有ったほうが良いでしょう。しかし、それだけでは花火っぽくなりません。具体的には、僕の大好きな垂れ花火が再現出来ませんでした。キーワードは「終端速度」です。
垂れ花火を見てみると分かるんですが、下に垂れ下がった花火の玉は、ある程度まで加速した後、加速が弱まって速度一定で落下している様に思えます。これは多分、空気抵抗により火の玉の落下速度が終端速度に近づいた為だと思われました。はい、これで花火の動きを再現する要素が出揃いましたね。以下の３つです。

球形に見えるような初速度
重力加速度のシミュレーション
空気抵抗のシミュレーション

花火らしい動きを目指して

花火に必要な要素が出揃った所で、それをどう料理するかを考えます。まずは

球形に見えるような初速度

を考えましょう。
これはあまり難しい話ではなく、直感的には、速度と角度の乱数を発生させ、ランダムで発生した速度をランダムな角度で回転させてやればOKです。ただ、よく考えて下さい。花火の広がりは球形なのであって、決して円形ではありません。これは些細ですが、見た目に大きな差を与える違いです。*1
では、球形に広がるような初速度を持たせるにはどうすれば良いでしょうか。色々やり方は有るでしょうが、猫は今回、速度と方向一定の三次元ベクトルをxyz三軸に対してランダムに回転させる事で、球形に広がった初速度を与えています。詳しいアルゴリズムは"Firelayer.pde"のgenereteBall関数をご覧ください。見て分からなかったら質問して下さい。

初速度が与えられた所で、

重力加速度のシミュレーション
空気抵抗のシミュレーション

を考えましょう。正しいかどうかは別にして、今回これらの運動は、微分形を使った運動方程式で下記のように表しました。

$\LARGE{\frac{\text{d}v_x}{\text{d}t}=-\kappa |\vec v|v_x}$
$\LARGE{\frac{\text{d}v_y}{\text{d}t}=-\kappa |\vec v|v_y+g}$
$\LARGE{\frac{\text{d}v_z}{\text{d}t}=-\kappa |\vec v|v_z}$

ここで、 $\kappa$ は空気抵抗係数、 $g$ は重力加速度係数です。それぞれの係数は今回はフィーリングで、もっとも花火らしく見えそうな数値を使っています。一応ハードコーディングは避けて変数にしてあるので、そちらを参照して下さい。

シミュレーションを行なっているのは"FallingFireball.pde"ファイルです。この程度の運動方程式ならシミュレーションしなくても解析的に解けそうな気がしますが、面倒だったので今回は一次のオイラー法でシミュレーションしました。酷い手抜きだ！という声が聞こえますが、その通り、手を抜きました。良いんだよ！それっぽく見えれば！！*2
オイラー法についての説明は、ググればいくらでも情報が出てくるので今回は割愛します。使った猫が言うのも何ですが、出来るだけオイラー法は使わないで下さい。*3猫も時間があればせめて二次のルンゲクッタ法には直すつもりです。

空気抵抗についてもう少し（2013/6/19追記）

空気抵抗の式が納得いかん！！ってお言葉を頂いたので、ちょろっと書き足しますね。
今回の大元のアイディアは、「空気抵抗は速度の二乗に比例する」というものです。それを元に、空気抵抗力 $\vec D$ について、空気抵抗係数を $\kappa$ として以下の様な式を考えます。
$\LARGE{|\vec D|=\kappa|\vec v|^2}$
もちろんこれはスカラー値なので、速度ベクトルとは逆方向に、空気抵抗力分の大きさを持つベクトルが実際の空気抵抗 $\vec D$ になります。では、速度ベクトルとは逆向きの単位ベクトルをスカラー値に掛け合わせましょう。速度ベクトルとは逆の単位ベクトルは $-\frac{\vec v}{|\vec v|}$ で表すことができます。これを掛け合わせると、
$\LARGE{\vec D=-\kappa|\vec v|^2\frac{\vec v}{|\vec v|}}$
です。分母分子に $|\vec v|$ が共通して出ているのでこれを取り除きましょう。
$\LARGE{\vec D=-\kappa|\vec v|\vec v}$
あとはこれをxyz軸の各成分に分解してやるだけです。
繰り返しになりますが、猫自身はあまりこれが正しい式だとは思っていません。ただ、今回は花火っぽく見えればそれで良かったのと、この式を突っ込んでみたらそれっぽく見えたので、この式を採用しました。

当日のアレコレ

そんなこんなでなんちゃってシミュレーションを持ち込んでLTをかまして来ました。

とりあえずLTはどうだったのか

あ…ありのまま　あの時　起こった事を話すぜ！
「おれは　奴の隣で出番を待っていたと思ったら　いつのまにかマイクを握っていた」
な…　何を言っているのか　わからねーと思うが、
おれも　何がおきたのか　わからなかった…
頭がどうにかなりそうだった…　誤操作だとかメモリリークだとか
そんなチャチなもんじゃあ　断じてねえ
もっと恐ろしいバグの片鱗を　味わったぜ…

もっと端的に結論をいうと

プログラム止まる
アトムが修正作業に入る
マイク投げ出される
プログラム再起動する
最後までしゃべり続ける
穴が有ったら入りたかった

じゃあ反省を聞こうか

例外処理不足です。
大抵の場合、例外が起きた時に備えて色々予防線を張るんですが、コーディングがギリギリになってしまって、例外に対する準備が圧倒的に不足していました。長時間ランニングテストはしたんですけど、何で僕の端末でしてたんでしょうね、Atomの端末でしろって話ですよね。ネットが突如切れたらどうなるかとか、何でテストしなかったんでしょうね。例外が起きても取り敢えず死なないように、何で出来なかったんでしょうね。人が集中するのが分かってる回線に頼らないで、何で自前の回線でやんなかったんでしょうね*4。
後わりと、開発に掛り切りになって、トーク内容の検分が不十分だったなぁって気はしてます。物作るだけなら引き篭もってりゃイイんだよ。タダの実装屋にはならないと、あれほど誓ったというのに。
次への検討課題満載で終わったLTでした。

ところでパーティーは楽しかったですか？

楽しくないパーティーとは何なのか。

大学院入学以降高専界隈の人間とどっぷり話す機会が減った私ですが、やっぱり高専生は楽しかったです。髪型のせいで認識されないことも多々ありましたが、分かってたよ多分そうなるって。

パーティーと言うことも有り、概念的なLTが多かったですが、流石の流暢なトークだったように思います。それに飛び入り発表の質の高さ。絶句しながらも、「うん、高専カンファっぽい」って思った自分が居ました。北海道か仙台には参加したいなぁ・・・お金有るかなぁ・・・・・・。

まとめ

何はともあれ、僕の誘いに乗ってくれたAtom君、急に投げたお仕事を瞬殺仕事人して下さったイトウカイトさん、お付き合い頂きありがとうございました！
それと、楽しいパーティーを形作って下さった発表者の皆様、参加者の皆様、何より運営のじゅーんさん、いがいがさん、会場提供のVOYAGE様、ありがとうございました！！

*1:上記の方法で円形に広げた火の玉は、ある円形領域を均等に埋め尽くす様な分布を取ります。一方、三次元的な球形の広がりを二次元に投影した場合、円形領域の周辺部分は、中心部分よりも火の玉が厚く存在するはずですよね？これは多分、景色を透明なビニールボール越しに見た時をイメージすると理解してもらえると思います。中心部分はビニールの層と視線が直角に交差するため、向こうの景色がよく見えます。しかし、周辺に行くにつれてビニールの層は視線と平行に近づいていき、見た目の厚みが増した結果、向こうの景色が見え難くなるはずです。

*2:良い子は真似せずせめて二次のルンゲクッタ法を使いましょう。可能なら四次のルンゲクッタを。計算スピードが心配なら、もう解析的に解いてから使って下さい。

*3:オイラー法は所謂一次元近似で、非常に精度が悪いです。

*4:自前のWiMAX回線は凄く安定してました。負荷テストもそっちでやってて、発表前から動かしてて殺し忘れてたランニングテストのプログラムは発表後も動いてたんだ。それ見て絶句したよホントに。