sikino のすべての投稿

非線形Schrödinger方程式のソリトン解

非線形シュレディンガー方程式
\(
\displaystyle i\frac{\partial \Psi}{\partial t}=-\frac{1}{2}\frac{\partial^2 \Psi}{\partial x^2}-g|\Psi|^2\Psi
\)

にはある解析解が存在します。それがソリトン(soliton)解と呼ばれるもので,上式のソリトン解は
\(
\displaystyle \Psi(x,t)=\sqrt{\Omega}\ {\rm sech}\left\{\sqrt{\Omega}\left(x\sqrt{g}-\frac{gV}{2}t\right)\right\}\cdot \exp\left\{i\frac{V\sqrt{g}}{2}x-i\frac{g}{2}\left(\frac{V^2}{4}-\Omega\right)t\right\}
\)

です。(\(g>0\),\(\Omega\):ソリトンの振幅、\(V\):ソリトンの速度に関するパラメータ、ソリトン自体の速度は\(V\sqrt{g}/2\))

[adsense1]

ソリトンの歴史的背景


「非線形」とは重ね合わせの原理が成り立たない系です。

1844年、スコットランドのJ.Scott-Russellによって孤立した波(solitary wave)を観測した事が報告されました J.Scott-Russellによる報告”Report on Waves”(リンク先のSR44.pdf, 16.3MB))。
当時の認識では、波は波動方程式で記述され、その波の速度\(v\)は\(v=f\lambda\)の元、一定である。だからパルス状の波は異なる波長の波の重ね合わせで書けているはずで、時間と共に分散していくはず。なのになぜ時間が経過しても孤立した波が存在できるのか?という事で大きな論争となりました。

60年後の1895年、オランダのKortewegとde Vriesによって”浅い水の波”を記述する非線形偏微分方程式(KdV方程式)が提出され、この方程式の特解として孤立波が与えられました。
孤立波は、

  1. 空間的に局在した波が、その性質(速さや形)を変えずに伝搬する
  2. 孤立波は互いの衝突に対して安定であり、各々の個別性を保持する

という性質を持つ非線形波動と定義されます[1]。
2番目の、粒子のような性質を持つことから、solitary に接頭語-on をつけ、soliton(ソリトン)と名づけられました。

その後、1981年に佐藤幹夫がソリトンの統一理論(佐藤理論やKP理論)を発表しました。
これによりソリトン方程式(ソリトンを記述し,かつ厳密に解ける方程式)に決着が付きました。
ソリトン方程式は非線形なのに厳密に解ける、可積分系である。

ソリトン方程式を解く方法は([4]を引用しますが)

上で指摘したように,logistic方程式が解けるからくりとソリトン方程式が解けるからくりはよく似ています.違いは,logistic方程式が変数変換一発で線形常微分方程式になってしまったのに対し,ソリトン方程式の場合は変数変換で双線形形式になり,双線形形式の解として行列式が現れ,行列式の中身に簡単な線形方程式が現れるというところです.しかし,離散化で保存するべき構造は明らかです.まず,解の中身の線形方程式を離散化し,行列式の構造をそのまま使って双線形形式を作る.最後に変数変換して非線形のレベルに戻ればよい.

となるそうです。

また、ソリトン方程式の特徴である、無限個の対称性(無限個の保存量)は、Gardner変換という変換をすることで証明できるそうです[5]。
これ以上はこの分野の専門家ではないので話せません。

ちなみに津波もソリトンの一つとみなせます。

ソリトン解が生まれるイメージ


なぜソリトン解が生まれるのでしょうか。
今、孤立した波(空間的に凸)を考えます。この時、

エネルギー的に安定になろうとして密度を均一にするために広がろうとする効果

粒子間を結び付ける引力相互作用(例えば水面だったら水と水との分子間力等)のため集まろうとする効果

のつり合いによって、丁度均衡が保たれるとき、このソリトン解が生まれます。

・・・実は、ソリトン解には2種類あります。
それは明るい(Bright)ソリトン解と暗い(Dark)ソリトン解です。
今まで話していたのは全て明るいソリトン解です。
暗いソリトン解とはどういったものでしょう。
暗いソリトン解とは、ある部分が空間的に凹んでいる、孤立した解です。

エネルギー的に安定になろうとして密度を均一にするためにその凹みを埋めようとする効果

粒子間の斥力相互作用のために粒子間を避けようとする効果

のつり合いによって、丁度均衡が保たれるとき、この暗いソリトン解が生まれます。
暗いソリトン解が生まれるのは斥力相互作用の時で、斥力相互作用を持つ系というのは、調べた限りでは量子力学のボーズ・アインシュタイン凝縮体で、暗いソリトンは渦ソリトンという形で現れるそうです。これ以上の具体例は分かりませんでした。もしも具体例を知っているという方は教えていただければ幸いです。
暗いソリトンの解析解は参考文献[1]の本に紹介されているので、それをご参考にしてください。

非線形シュレディンガー方程式におけるソリトン解


では本題の、非線形シュレディンガー方程式における(明るい)ソリトン解を考えましょう。
下の形の非線形シュレディンガー方程式を考えます。
\(
\displaystyle i\frac{\partial \Psi}{\partial t}=-\frac{1}{2}\frac{\partial^2 \Psi}{\partial x^2}-g|\Psi|^2\Psi
\)

ここで、\(\Psi=\Psi(x,t)\)で、\(g\)はの値で相互作用の強さ(この場合、引力相互作用)を表します。

この非線形シュレディンガー方程式のソリトン解\(\Psi(x,t)\)は、
\(
\displaystyle \Psi(x,t)=\sqrt{\Omega}\ {\rm sech}\left\{\sqrt{\Omega}\left(x\sqrt{g}-\frac{gV}{2}t\right)\right\}\cdot \exp\left\{i\frac{V\sqrt{g}}{2}x-i\frac{g}{2}\left(\frac{V^2}{4}-\Omega\right)t\right\}
\)

であり、\(\Omega\)はソリトンの振幅の大きさ、\(V\)はソリトンの速度を決めるパラメータを表します。ソリトン自体の速度は\(V\sqrt{g}/2\)となります([2]を参考)。
また、\({\rm sech}(x)\)は双曲線関数の一種(双曲線正割と呼ばれる)であり、
\(
\displaystyle {\rm sech}(x)=\frac{1}{\rm cosh(x)}=\frac{2}{e^{x}+e^{-x}}
\)
を表します。

解析解のプロット

解析解をプロットします。gnuplotコードは下のほうに載せておきます。
\(g=2, V=1, \Omega=1\)とすると、以下の振る舞いが観測されます。
ここで紫はソリトン解の実部、緑は虚部、青は絶対値2乗を表します。
動画は1枚当たり原子単位系で0.1秒、合計で10秒間のシミュレーションです。
また、このソリトンの速度は\(V\sqrt{g}/2\sim 0.7071\)です。
soliton1

[adsense2]

gnuplotコード


非線形シュレディンガー方程式のソリトン解(解析解)を出力します。
gnuplot上で下のスクリプトを実行してください。
(ただし、gnuplot ver4.6以降に限ります。)

omega=1e0
V=1e0
g=2e0
x0=-5e0 # initial position

set xr[-10:10]
set yr[-1.5:1.5]
set samples 1000
set xl "x[a.u.]"

sech(x)=2e0/(exp(x)+exp(-x))
amp(x,t)=sqrt(omega)*sech(sqrt(omega)*((x-x0)*sqrt(g)-g*V*t*0.5e0))
phase(x,t)=V*sqrt(g)*x*0.5e0-g*t*0.5e0*(V*V*0.25e0-omega)
soliton(x,t)=amp(x,t)*exp({0e0,1e0}*phase(x,t))


#set term gif animate delay 10 optimize size 960,720
#set output 'movie.gif'
do for[i=0:100:2]{
   t=i*0.1e0
   plot abs(soliton(x,t))**2 lw 3 lc 1 lt 1 ti sprintf("|\psi|^2, t=%2.1f",t),\
        real(soliton(x,t)) lw 3 lc 2 lt 2 ti "Real",\
        imag(soliton(x,t)) lw 3 lc 3 lt 3 ti "imag"
}
#set out
#set terminal wxt enhanced

もっとソリトンについて知りたい方はまず参考文献[3]を読むことをお勧めします。
その後、[4]を読み、[5]を読み、[1]の本を読むのが良いと思われます。
[3],[4]は簡単な表現を用いてソリトンとその後の発展について記述されています。

参考文献


[1]和達三樹著 『非線形波動 (現代物理学叢書) 』岩波書店 (2000年) p.7
[2]和達三樹著 『非線形波動 (現代物理学叢書) 』岩波書店 (2000年) p.29

[3]ソリトンの数学 – Researchmap
[4]ソリトン ~ 不思議な波が運んできた,古くて新しい数学の物語 ~
[5]〔連載〕非線形波動―ソリトンを中心として―第5章 逆散乱法


↑画像のフォントはキユマヤ園様による数式フォント -びゅんびゅん→SSSです!


モスクワでSIMカードの購入

まとめ

  1. 2015年10月、モスクワにあるМегаФонでプリペイドsimカードを買った時のお話
  2. 店員にまずは英語で”I want a prepaid sim card.”と話す。
    ダメなら”Предоплаченные Cим-карты, Пожалуйста.“と言う。
  3. パスポートを渡し、サインして、お金払ってお仕舞い。
  4. ※標準, micro, nanoの違いはありません。simを切り離すことでどれでも対応できます。

モスクワ市内であればそこそこ英語は通じます。
可能な限り空港で買うのが良いでしょう。
シェレメーチエヴォ国際空港(SVO)の場合、一度空港の出発ロビーへ行くのが良いでしょう。確か到着ロビーは1階で、出発ロビーは3階か4階です。そこまで行けば歩いていれば確実に見つかります。


ロシアのモスクワでsimカードを購入してインターネットに繋げるまでの体験です。
この情報は2015年10月のお話で、あくまで僕が体験したことを思い出しながら書いています。
当たり前ですがsimフリーじゃなきゃできません。

英語は申し訳程度の保険的な扱いです。
モスクワでは英語は通じません。英語の案内版すらありません。

僕は
МегаФон
というお店で買いました(英語名Megafon,メガフォン)。

まずお店に入ります。
まずは店員にダメ元で英語で話しかけるのがいいでしょう。
通じれば何でもいいですが、例えば
“I want a prepaid sim card.”
とかでいいでしょう。
これでダメでしたらロシア語で、
プリペイドsimカードを下さい。

Предоплаченные Cим-карты, Пожалуйста.

と言うか、この文を見せればいいでしょう。
※スマートフォン等にgoogle翻訳のオフライン版を事前に入れておくことをお勧めします。
僕の場合はたまたま英語が通じる店員が店内にいたので困らずに変えました。

これで通じたら、どのプランがいいか、を聞かれると同時にパスポートの提示が求められます。
店員にパスポートを渡します。
その後、simカードが入ったカードの束が渡され、
「好きな電話番号を選んでくれ」
と言われます。
このプリペイドsimに付属する電話番号は(多分)ロシア国内でМегаФонの電話へ掛ける場合に使えるものです。
まぁ、掛けるつもりが無いのならば何でもいいでしょう。

2015年10月時点ではМегаФонのプリペイドsimのデータプランは5つ存在し、
199pub, 0.5GBまで, ロシアの番号でМегаФонの番号なら通じる, 通話可能時間300分
390pub, 3GBまで, ロシアの番号でМегаФонの番号と家への電話なら通じる, 通話可能時間400分
590pub, 4GBまで, ロシアの番号でМегаФонの番号と家への電話なら通じる, 通話可能時間600分
1290pub, 8GBまで, ロシアの番号全てに通じる, 通話可能時間1500分
2700pub, 10GBまで, ロシアの番号全てに通じる, 通話可能時間5000分
です。

で、サインを求められるので、サインしてお金を払って終了です。
このサインは漢字でも平仮名でもアルファベットでも何でも構いません。

最後にこんなsimカードがもらえます。
ロシアモスクワsimカード_c

時間としては10~15分程度だったでしょうか。早かった印象があります。

– 1 = 1 ?

が発端です。間違いであることは直観で分かります。
問題はどこに間違いがあるか?です。

改めて問題を書けば、
\(
\begin{align}
-1&=i\times i \\
&=\sqrt{(-1)}\times\sqrt{(-1)} \\
&=\sqrt{(-1)\times(-1)} \\
&=\sqrt{1} \\
&=1
\end{align}
\)
です。

結論を先に言えば、2行目から3行目にかけての変形がダメです。
定義域の考慮が抜けています。
なぜダメなんでしょうか。詳しく見ていきましょう。

間違いが生じるのはなぜ?


\(i\)を極座標形式で考えます。
すると、
\(
\displaystyle i=0+i=e^{i\frac{\pi}{2}}
\)

ですよね。これは正しいです。
そして次に2行目から3行目を極形式でゆっくりと書いていけば、
\(
\begin{align}
-1&=i\times i \\
&=e^{i\frac{\pi}{2}}\times e^{i\frac{\pi}{2}} \\
&=\left(e^{i\pi}\right)^{1/2} \times\left(e^{i\pi}\right)^{1/2} \cdots \mbox{2行目に対応} \\
&=\left(e^{i2\pi}\right)^{1/2}\cdots \mbox{3行目に対応}
\end{align}
\)
です。

もしも、\(e^{i\theta}\)の\(\theta\)の範囲を勝手に\(0 \le \theta \lt 2\pi\)にしていたら。
この場合、角度\(2\pi\)はこの定義域内に含まれていません。
\(e^{i\theta}\)は\(2\pi\)の周期性を持つはずだ、だから
\(
e^{i2\pi}=e^{i0}=1
\)

であるはずだ、と思ってしまうんです。
これを行ってしまうと、最後の式は、
\(
\begin{align}
\left(e^{i2\pi}\right)^{1/2}&=\left(e^{i0}\right)^{1/2} \\
&=1^{1/2} \\
&=1
\end{align}
\)

となります。

よって間違った答え\(-1=1\)が生まれます。

間違いを起こさないためには


この問題を解決するためには定義域をきちっと示してあげればいいのです。
別の複素数を掛ける場合、通常、\(\theta\)の定義域も考慮しなければなりません。
故に、厳密に書けば、
\(
\displaystyle i=e^{i\frac{\pi}{2}}\ \ \ (0\le\theta \lt 2\pi)
\)

であり、
\(
\begin{align}
-1&=i\times i \\
&=e^{i\frac{\pi}{2}}\times e^{i\frac{\pi}{2}}\ \ \ \ (0\le\theta \lt 4\pi)\\
&=\left(e^{i2\pi}\right)^{1/2}\ \ \ (0\le\theta \lt 4\pi)
\end{align}
\)
なのです。\(e^{i\theta}\)は本来, 周期関数であって多価関数であることを忘れてはいけません。
故に、定義域内に収まっていないからと言って\(2\pi\)を加える必要はないのです。
だから、形はそのままであり、計算は
\(
\begin{align}
&\left(e^{i2\pi}\right)^{1/2}\ \ \ (0\le \theta\lt 4\pi)\\
&=e^{i\pi}\ \ \ (0\le \theta\lt 2\pi)\\
&=-1
\end{align}
\)
だから、やはり\(-1=-1\)なのです。

※上の式\(-1=i\times i\cdots\)は長いので、
\(
\sqrt{(-1)^2}
\)

の値は何ですか?という問いかけでもいいと思います。
この場合、この答えは\(\pm 1\)です。どちらでも正解です。
なぜならば、\(X=\sqrt{(-1)^2}\)とおいて、
\(
\begin{align}
X&=\sqrt{(-1)^2} \\
X^2&=1 \\
&(X+1)(X-1)=0 \\
\end{align}
\)
だからです。


この問題が顕著に表れるのは数値計算です。コンピュータは定義域の考慮はできません。人為的に入れるしか手段は無く、必ず\(2\pi\)の範囲に収まっていなければならないのです。なので、数値計算で複素数を扱う場合、この問題が発生していないか確かめましょう。符号が知らない間に変わっている、軽視いているととんでもない問題になります。

離散フーリエ変換と高速フーリエ変換(fortran90)

こちらのページは後程消去いたします。

以下のページに統合しますので、ご参照ください。
https://slpr.sakura.ne.jp/qp/dft-for-numerical-calculation/

本ページにはミスもあり、上の方が正しいですのでご参照ください。


数値計算で離散フーリエ変換を行う方法です。
言語は
fortran90

intelのMKL(マス・カーネル・ライブラリー)
を用いて離散フーリエ変換を行いたいと考えます。


高速フーリエ変換は離散フーリエ変換を高速に行う手法のことです。

通常の離散フーリエの計算量のオーダーは、\(O(N^2)\)であり、
基数2の高速フーリエの計算量のオーダーは、\(O(Nlog_2 N)\)です。

基数2は良く言われる通常の配列の個数\(2^m\)個のことを表しています。

離散フーリエ変換を行う際には周期境界条件が課されます。
分点の端で周期境界条件が満たされていない場合、それはその点でステップ関数が含まれることを意味します。
依って予期しない高周波成分が出現するので、注意してください。
一番良い判定方法は関数の対数微分が一致するかどうかを判定するのが良いでしょう。

MKLの離散フーリエ変換ルーチンの概要1


MKLの離散フーリエ変換ルーチンは非常に優秀です。
渡す配列サイズによってプログラムが自動的に判断し、最適な手法で離散フーリエ変換を行います。
MKLのマニュアルによると、

注 : DFT関数は任意の長さをサポートしている。
これらのルーチンは、基数 2 だけでなく、基数 3、5、7、11 に対しても高
性能で広範な機能性を提供する。対応する基数の一覧は、使用しているラ
イブラリー・バージョンの「インテルMKLテクニカル・ユーザー・ノー
ト」を参照のこと。

とあります。
高速フーリエ変換の考え方は\(2^m\)の時だけでなく、\(3^m, 5^m, 7^m, 11^m, \cdots\)の時にも同様に考えることができます。
この時、何の整数乗にするかを区別するために “基数” という言葉を使います。
プログラムを作る場合、通常は一番シンプルな基数”2″が選ばれます。

具体的に配列サイズ\(N\)が
\(N= 2^m\)、ただしmは正の整数(配列をAとするとサイズ\(A(0:2^m-1)\))の時に基数2の高速フーリエ変換が行われます。
それ以外の時は通常の離散フーリエ変換が行われているはずです。

余り自信が無いのではっきりとしたことは言えませんが、分点量(2~20000程度まで)に対するMKLのフーリエ変換の計算速度はほとんどわからないくらいでした。特定の基数の時だけ早いんだろうと予想していましたが、そんなことはなく、それだけ優秀であるという事でしょう。

もう6年前に書かれている記事ですので今はどうなっているかわかりませんが、FFTWと計算速度を比較したページがありましたので載せておきます。
一言でいえば、MKLによる離散フーリエ変換はFFTWよりも早いよ、ということです。
インテル謹製の数値演算ライブラリ「MKL」を使ってプログラムを高速化 3ページ -OSDN, (2009)

MKLの離散フーリエ変換ルーチンの概要2


MKLでの1次元、長さnの順方向離散フーリエ変換は、
\(
\displaystyle z_k=\sum_{j=0}^{n-1}w_je^{-i2\pi jk/n}
\)

に従い行われます(MKLリファレンスマニュアル、11-26)。
数値計算の世界では、余分な係数\(\frac{1}{2\pi}\)をなくすため、指数関数の肩に\(2\pi\)が掛けられたものが良く使用されます。

※ここでは、フーリエ変換前の空間を”位置空間“と呼び、フーリエ変換後の空間を”波数空間“と呼ぶことにします。
位置波数ではなく、時間周波数としても同じことなのでどちらでもいいのですが、ここでは、位置波数の呼び方に統一いたします。

DFTexampleNeven
上図のように位置空間上で区間\(L\)の領域で\(N\)個に分割します。
上図はNが偶数の場合です。
MKLの離散フーリエ変換ルーチンはいかなるNに対してもフーリエ変換が可能ですが、奇数の場合、以下のように波数成分の格納が変わります。
周波数格納順

\(
f(x)=\sin(\alpha x)
\)

を離散フーリエ変換した場合、
波数空間でのフーリエ変換後のピークは\(k=\frac{\alpha}{2\pi}\)に現れます。
この時、f(x)のフーリエ変換後の関数\(g(k)=F_{[f(x)]}\)は複素数として格納され、上記の場合では
\(
\displaystyle g(k)=i\left\{\frac{N}{2}\delta(k-\frac{\alpha}{2\pi})+\frac{N}{2}\delta(k+\frac{\alpha}{2\pi})\right\}
\)

となります。ここで\(i\)は虚数単位を表し、\(\delta(k)\)はデルタ関数を表します。
波数空間での最大の周波数値\(f_c\)はナイキストの定理より、
\(
\begin{align}
f_c = \frac{1}{2\Delta}, \Delta=\frac{L}{N}
\end{align}
\)
と表されます。これは、波数空間では\(\pm f_c\)までサンプリング可能であることを示しています。
フーリエ変換後の値\(g(k)\)の分点上の和\(S\)、すなわち
\(
\displaystyle S=\sum_{j=0}^{N-1}g(k_j)
\)

は離散フーリエ変換の性質により、\(S\)は分点の総数\(N\)に等しくなります。

逆離散フーリエ変換によって元の空間に戻る場合、規格化因子\(N\)で割らなければなりません。

MKLの1次元離散フーリエ変換ルーチンのプログラム例


使い方はFFT Code Examples -Intel®Developer Zone を参考としました。

例として
関数
\(
f(x)=\sin(4x)
\)
をフーリエ変換し、元に戻すプログラムを書きましょう。
プログラムの流れとしては、

関数\(f(x)\)の位置空間での書き出し(ファイル”before.d”に書き出し、横軸:位置\(x\), 縦軸\(f(x)\))
↓ 順方向フーリエ変換(規格化因子で割られない)
関数\(g(k)=F_{[f(x)]}\)の書き出し(ファイル”k.d”に書き出し、横軸:\(k/(2\pi)\), 縦軸\(g(k)\))
↓ 逆方向フーリエ変換(規格化因子Nで割られる)
関数\(f(x)=F^{-1}_{[g(k)]}\)の書き出し(ファイル”after.d”に書き出し、横軸:x, 縦軸f(x))
になっています。

コンパイルを行うにあたり重要な事が2点あります。
1, MKLを使う事
2, fftサブルーチンのにinclude文を付け加える事(かも?)

このincludeは、MKLにある離散フーリエ変換ルーチンを使いますよーというサインであり、パスはMKLのバージョンや、インストールした状況によって変わります。
うまく、ファイル”mkl_dfti.f90″を見つけてください。パスを上から辿っていくのがいいと思います。

※2016/02/27時点での最新バージョン(ifort –version で確認できます)
ifort (IFORT) 16.0.2 20160204
において、mkl_dfti.fのデフォルトでの場所は、
“/opt/intel/compilers_and_libraries_2016/linux/mkl/include/mkl_dfti.f90”
です。また、このバージョンだけかどうかわかりませんが、include文が無くてもどうやら動くようです。

プログラムのコンパイルと実行は

$ ifort -mkl main.f90
$ ./a.out

で良いでしょう。

サブルーチンのコンセプトは、
ただ順/逆方向離散フーリエ変換をしてくれるサブルーチンです。速度は追求していません。

速度を追及する書き方をしているのは下のベンチマークの項目です。プログラムの可読性が悪くなるので、中身を知らないと調整が出来ません。
ベンチマークに書いてあるように、サブルーチンdftの中身の

  Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,1,n)
  Status = DftiCommitDescriptor(hand)
  Status = DftiFreeDescriptor(hand)

はFFTを使うための準備で、これは配列の大きさが変わらない限り1度実行するだけで良いのです。
これを排除すると3-4倍ほど速くなります。

上で用いているサブルーチン
dft,dfts,dftfはそれぞれフーリエ変換、フーリエ変換後の空間の値、ソート用のルーチンであることを意味し、その中身は以下のようなものです。例のごとく、MKLのルーチンはそのままでは使いにくいので1クッション挟みます。

ここ↑に書いてあるサブルーチン

  1. dft(N,z,FB)
  2. dftf(N,f,h)
  3. dfts(N,f,z,mf,mz)

を説明します。

  1. dft(N,z,FB)

    このルーチンによって順/逆方向フーリエ変換が行われます。
    引数はそれぞれ

    • (入力) N : integer, データ配列のサイズ
    • (入出力) z : complex(kind(0d0)), データ配列z(0:N-1)
    • (入力) FB : character, 順方向(“forward”)、逆方向(“backward”)の指定

    です。データ配列に上書きして値が返ってきます。
    サイズNのデータ配列z,順方向の離散フーリエ変換をしたい場合は、

    call dft(size(z,1),z,"forward")

    とするのがいいでしょう。
    規格化は順方向の時には行わず、逆方向の時にだけ\(1/N\)を掛けています。

  2. dftf(N,f,h)

    このルーチンは波数空間での横軸の値f(0:N-1)を与えます。
    引数はそれぞれ

    • (入力) N : integer, データ配列のサイズ
    • (出力) f : double precision, 波数空間の横軸の値f(0:N-1)
    • (入力) h : double precision, 位置空間でのサンプルレート(刻み幅)

    です。
    波数fの順番は少しややこしく、
    周波数格納順
    となっています。
    また、ここで出力されるfは波数を\(2\pi\)で割ったものが出力されます。
    具体的には、\(\sin(x)\)をMKLのフーリエ変換を行いグラフを書くとフーリエ変換後の空間(波数空間)で、ピークの位置は\(\pm 0.159115494309\cdots (\sim \frac{1}{2\pi})\)に表れることを意味します。
    また、変換先の空間をグラフで描画したい時に値を小さい順に並べ替えたい必要が出てきたときは並べ替え用のルーチンdftsを使用してください。

  3. dfts(N,f,z,mf,mz)

    このルーチンは変換先の波数空間でグラフを描画したい時に、順番を値を小さい順に並べ替えます。
    引数はそれぞれ

    • (入力) N : integer, データ配列のサイズ
    • (入力) f : double precision, 波数空間の横軸の値f(0:N-1)
    • (入力) z : complex(kind(0d0)), 波数空間でのデータ配列z(0:N-1)
    • (出力) mf : double precision, 並び替えられた波数空間の横軸の値mf(0:N-1)
    • (出力) mz : complex(kind(0d0)), 並び替えられた波数空間でのデータ配列z(0:N-1)

    です。上書きして引数を返す必要性がないだろう、と判断して別の配列に書き出させています。
    もしもメモリが足らないとか、気に食わない人は書き換えてください。

[adsense1]

出力ファイル”before.d”, “k.d”, “after.d”を書いてみます。
位置空間での区間\([-3\pi\sim3\pi]\)で、分点数\(N=288\)で行った場合、
位置空間での波形、この場合は\(\sin(4x)\)が描かれます。
位置空間

波数空間では
\(
\displaystyle g(k)=i\left\{\frac{288}{2}\delta(k-\frac{4}{2\pi})+\frac{288}{2}\delta(k+\frac{4}{2\pi})\right\}
\)

が得られるはずで、実際にプロットして確かめてみると
N288_c
となります。
※Nは2以上のどんな整数でもokです。偶数、奇数、素数であっても。

※実部に現れる小さなピークに関して
原因はよくわかりません。\(N=2^m\)ではない時になんかこれが出てきます。
不安な人は\(N=2^m\)に固定して使ってください。
N=1024で行うとこうなります。
N1024

ベンチマーク

同じ計算機上でのベンチマーク結果を載せます。
コンパイラ:ifort (IFORT) 16.0.2 20160204
MKL:/opt/intel/compilers_and_libraries_2016.2.181
計算スレッド数1
ベンチマークdft_mkl_c

ベンチマーク用プログラムはこちら↓

並列計算時の計算速度

call mkl_set_num_threads(Ncpu)

を用いてFFTを並列計算で行った時、計算時間のCPU個数(8コア16スレッド)と配列の要素数依存性を調べます。
フーリエ変換するのはガウス関数とします。
結果はこのようになりました。

複数の線は試行回数を表しています。プログラムは全く同じで、実行した時刻だけが違います。
”計算1回”とは(順方向→逆方向フーリエ変換)の1セットを”1回”としています。
8コア16スレッドの環境で実行したので、物理的なCPUの個数8個で計算速度はおおよそ打ち止めになっています。

計算に用いたコードはこちら

2次元のフーリエ変換について


ルーチン等々はほぼ同じです。
下のプログラムは
\(
f(x,y)=\sin(4x)+\sin(3y)
\)
をフーリエ変換するものです。
フーリエ変換実行前の,位置空間での関数の形(“before.d”, x軸:\(x\), y軸:\(y\), z軸\(f(x,y)\)),
フーリエ変換し、波数空間で見た関数の形(“k.d”, x軸:\(k_x/(2\pi)\), y軸:\(k_y/(2\pi)\), z軸\(Imag\{g(k_x,k_y)\}\)),
逆フーリエ変換した後の関数の形(“after.d”, x軸:\(x\), y軸:\(y\), z軸\(f(x,y)\))は
before_k_after_2D_c
となります。
波数空間で\(k_x\)だけを見た時,\(k_x=0\)にあたかもピークが見えます。
これはy成分があるために出てくるもので、証明は以下のようにできます。
関数\(f(x,y)\)は
\(
\begin{align}
f(x,y)&=\sin(4x)+\sin(3y)\\
&=i\left\{-\frac{1}{2}e^{i4x}+\frac{1}{2}e^{-i4x}-\frac{1}{2}e^{i3y}+\frac{1}{2}e^{-i3y}\right\}
\end{align}
\)
なので、連続のフーリエ変換を考えるとすると、
\(
\begin{align}
&\int f(x,y)e^{-ik_x x}e^{-ik_y y}dxdy\\
&=i\left\{-\frac{1}{2}\int e^{i(4-k_x)x} dx \int e^{-ik_y y} dy+\frac{1}{2}\int e^{-i(4+k_x)x} dx \int e^{-ik_y y} dy \right. \\
&\;\;\;\;\;\;\;\;\;\;\;\; \left.-\frac{1}{2}\int e^{-ik_x x} dx \int e^{i(3-k_y)y} dy+\frac{1}{2}\int e^{-ik_x x} dx \int e^{-i(3+k_y)y} dy\right\} \\
&=i\pi\left\{-\delta(k_x-4)\delta(k_y)+\delta(k_x+4)\delta(k_y)-\delta(k_x)\delta(k_y-3)+\delta(k_x)\delta(k_y+3)\right\}
\end{align}
\)
となります。なので、フーリエ変換によってピークが出てくる位置というのは、
\((k_x,k_y)=(4,0), (-4,0), (0,3), (0,-3)\)
の4点となります。なので、y成分がいるために\(k_x=0\)にあたかもピークがあるように見えてしまうんですね。(証明終わり)

2次元のフーリエ変換はそのまま2次元データの配列を渡すことは出来ないため、一度、1次元配列への変換(サブルーチンswap1d2d)を呼び出しています。本来は2次元データと1次元データの間の関係ではEquivalenceを使うべきなのですが、何をしているのかの(個人的な)見やすさを考慮してわざわざ別の配列を用意して格納しています。

include "/opt/intel/mkl/include/mkl_dfti.f90"
program main
  implicit none
  integer::i,j,Nx,Ny
  double precision,allocatable::x(:),kx(:),mkx(:)
  double precision,allocatable::y(:),ky(:),mky(:)
  double precision::hx,xmin,xmax
  double precision::hy,ymin,ymax
  complex(kind(0d0)),allocatable::z(:,:),mz(:,:)

double precision,parameter::pi=dacos(-1.d0)
  complex(kind(0d0))::func
  external::func

Nx=100; Ny=100
  allocate(x(0:Nx-1),kx(0:Nx-1),mkx(0:Nx-1))
  x=0d0; kx=0d0; mkx=0d0
  allocate(y(0:Ny-1),ky(0:Ny-1),mky(0:Ny-1))
  y=0d0; ky=0d0; mky=0d0
  allocate(z(0:Nx-1,0:Ny-1),mz(0:Nx-1,0:Ny-1))
  z=dcmplx(0d0,0d0); mz=dcmplx(0d0,0d0)

xmin=-3d0<em>pi
  xmax=3d0</em>pi
  hx=(xmax-xmin)/dble(Nx)
  ymin=-3d0<em>pi
  ymax=3d0</em>pi
  hy=(ymax-ymin)/dble(Ny)

do i=0,Nx-1
     x(i)=xmin+dble(i)<em>hx
  enddo
  do j=0,Ny-1
     y(j)=ymin+dble(j)</em>hy
  enddo
  do i=0,Nx-1
     do j=0,Ny-1
        z(i,j)=func(x(i),y(j))<br />
     enddo
  enddo

open(21,file="before.d")
  do i=0,Nx-1
     do j=0,Ny-1
        write(21,'(4e20.10e3)')x(i),y(j),dble(z(i,j)),dimag(z(i,j))
     enddo
     write(21,*)
  enddo
  close(21)

call dftf(Nx,kx,hx); call dftf(Ny,ky,hy) ! get frequency

call dft2d(Nx,Ny,z,"forward") ! forward dft

call dfts2d(Nx,Ny,kx,ky,z,mkx,mky,mz) !sort frequency
  open(21,file="k.d")
  do i=0,Nx-1
     do j=0,Ny-1
        write(21,'(4e20.10e3)')mkx(i),mky(j),dble(mz(i,j)),dimag(mz(i,j))
     enddo
     write(21,*)
  enddo
  close(21)

call dft2d(Nx,Ny,z,"backward") ! backward dft
  open(21,file="after.d")
  do i=0,Nx-1
     do j=0,Ny-1
        write(21,'(4e20.10e3)')x(i),y(j),dble(z(i,j)),dimag(z(i,j))
     enddo
     write(21,*)
  enddo
  close(21)

stop
end program main

function func(x,y)
  implicit none
  double precision::x,y
  complex(kind(0d0))::func

func=dcmplx(dsin(4d0<em>y)+dsin(3d0</em>x),0d0)

return
end function func

がメインのプログラムで、以下が上で呼び出しているサブルーチンの中身です。

subroutine dftf(N,f,h)
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::h
  double precision,intent(out)::f(0:N-1)

integer::i
  double precision::mf(0:N-1)

if(mod(N,2).eq.0)then
     do i=0,N-1
        mf(i)=(dble(i+1)-dble(N)<em>0.5d0)/(dble(N)</em>h)
     enddo

<pre><code> do i=0,N-1
    if(i.le.N/2-2)then
       f(i+N/2+1)=mf(i)
    else
       f(i-N/2+1)=mf(i)
    endif
 enddo

else
do i=0,N-1
mf(i)=(dble(i)-dble(N-1)0.5d0)/(dble(N)h)
enddo


 do i=0,N-1
    if(i.le.(N-1)/2-1)then
       f(i+(N-1)/2+1)=mf(i)
    else
       f(i-(N-1)/2)=mf(i)
    endif
 enddo

endif

return
end subroutine dftf
!——————————–
subroutine dft2d(Nx,Ny,z,FB)
!sikinote
!date : 2015/12/29
!developer : sikino
!CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/deed.ja)
use MKL_DFTI
implicit none
integer,intent(in)::Nx,Ny
complex(kind(0d0)),intent(inout)::z(0:Nx-1,0:Ny-1)
character(*),intent(in)::FB

complex(kind(0d0))::z1(0:Nx*Ny-1)
integer::Status
TYPE(DFTI_DESCRIPTOR),POINTER::hand
integer::id(1:2)

!DFT : Discrete Fourier Transform
!
!n –> number of data.
!z(i,j) –> value of data at x and y
!FB –> “forward” : Forward DFT
! –> “backward” : Backward DFT

call swap1d2d(nx,ny,z,z1,1)
id(1)=nx; id(2)=ny

Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,2,id)
Status = DftiCommitDescriptor(hand)

if(trim(FB).eq.”forward”)then
Status = DftiComputeForward(hand,z1)
elseif(trim(FB).eq.”backward”)then
Status = DftiComputeBackward(hand,z1)
else
write(6,*)”DFT string different”
stop
endif
Status = DftiFreeDescriptor(hand)

call swap1d2d(nx,ny,z,z1,-1)

if(trim(FB).eq.”backward”)then
z=z/dble(nx*ny)
end if

return
end subroutine dft2d
!———————————————
subroutine swap1d2d(nx,ny,z2,z1,isign)
!change array matrix between 1D and 2D
! if isign == 1 –> from 2D to 1D
! if isign == -1 –> from 1D to 2D
implicit none
integer,intent(in)::nx,ny,isign
complex(kind(0d0)),dimension(0:nx-1,0:ny-1),intent(inout)::z2
complex(kind(0d0)),dimension(0:nx*ny-1),intent(inout)::z1
integer::j1,j2,k

if(isign.eq.1)then
do j2=0,ny-1
do j1=0,nx-1
k=j2nx+j1
z1(k)=z2(j1,j2)
enddo
enddo
elseif(isign.eq.-1)then
do j2=0,ny-1
do j1=0,nx-1
k=j2
nx+j1
z2(j1,j2)=z1(k)
enddo
enddo
endif

return
end subroutine swap1d2d
!—————————-
subroutine dfts2d(Nx,Ny,fx,fy,z,mfx,mfy,mz)
implicit none
integer,intent(in)::Nx,Ny
double precision,intent(in)::fx(0:Nx-1),fy(0:Ny-1)
complex(kind(0d0))::z(0:Nx-1,0:Ny-1)
double precision,intent(out)::mfx(0:Nx-1),mfy(0:Ny-1)
complex(kind(0d0)),intent(out)::mz(0:Nx-1,0:Ny-1)
complex(kind(0d0))::mmz(0:Nx-1,0:Ny-1)
integer::i,j,k,l

if(mod(Ny,2).eq.0)then
do i=0,Ny-1
if(i.le.Ny/2)then
j=i+Ny/2-1
mfy(j)=fy(i)
mz(0:Nx-1,j)=z(0:Nx-1,i)
else
j=i-Ny/2-1
mfy(j)=fy(i)
mz(0:Nx-1,j)=z(0:Nx-1,i)
endif
enddo
else
do i=0,Ny-1
if(i.le.(Ny-1)/2)then
j=i+(Ny-1)/2
mfy(j)=fy(i)
mz(0:Nx-1,j)=z(0:Nx-1,i)
else
j=i-(Ny-1)/2-1
mfy(j)=fy(i)
mz(0:Nx-1,j)=z(0:Nx-1,i)
endif
enddo
endif

mmz=mz

if(mod(Nx,2).eq.0)then
do k=0,Nx-1
if(k.le.Nx/2)then
l=k+Nx/2-1
mfx(l)=fx(k)
mz(l,0:Ny-1)=mmz(k,0:Ny-1)
else
l=k-Nx/2-1
mfx(l)=fx(k)
mz(l,0:Ny-1)=mmz(k,0:Ny-1)
endif
enddo
else
do k=0,Nx-1
if(k.le.(Nx-1)/2)then
l=k+(Nx-1)/2
mfx(l)=fx(k)
mz(l,0:Ny-1)=mmz(k,0:Ny-1)
else
l=k-(Nx-1)/2-1
mfx(l)=fx(k)
mz(l,0:Ny-1)=mmz(k,0:Ny-1)
endif
enddo
endif

return
end subroutine dfts2d

昔のMKLの離散フーリエ変換ルーチンについて


MKLは汎用性を高くするために仕様変更が行われ、
以前は

call zfft1d(...)

とするだけで順/逆方向フーリエ変換は使えていたようですが、いつからか(2004~2006年辺り?)廃止され、使えなくなっています。

間違い、勘違いがあると思います。
いかなる問題が発生しても私は責任を一切負いません。
それを念頭に置いたうえ使用してください。

参考文献


インテル®マス・カーネル・ライブラリーリファレンスマニュアル(2006年)
ニューメリカルレシピinC

[adsense2]

球体の抗力係数

説明


本稿では、完全な球体に働く抗力係数について述べます。
ここでは概要だけを説明します。詳細は以下のPDFをご覧ください。

https://slpr.sakura.ne.jp/qp/supplement_data/drag_coefficient_air/drag_coefficient.pdf

概要

速度\(v(=|\mathbf{v}|)\)で動く半径\(R\)の完全な球体に働く空気抵抗力の大きさ\(F_d(=|\mathbf{F}_d|)\)は,
\begin{align}
F_d =\frac{1}{2} C_d(R_e) \rho \pi R^2 v^2
\end{align}
と書けます. 空気抵抗力\(\mathbf{F}_d\)は, 球体の速度方向と反対に働くので, 方向を含めれば
\begin{align}
\mathbf{F}_d =-\frac{1}{2} C_d(R_e) \rho \pi R^2 v^2\frac{\mathbf{v}}{v}
\end{align}
と書けます.

つまり、重力加速度\(g\)の重力下における質量\(m\)の完全な球体に対する運動方程式は
\begin{align}
m\frac{d^2}{dt^2}\mathbf{r} =-mg\mathbf{k}-\frac{1}{2} C_d(R_e) \rho \pi R^2 v^2\frac{\mathbf{v}}{v}
\end{align}
となります。ここで、\(\mathbf{k}\)は重力の方向の単位ベクトル(多くの場合では鉛直上向きを\(z\)軸の正の方向、単位ベクトル\(\mathbf{k}\)とするので、負号が付いています)です。
速度が非常に遅い極限\((R_e\to 0)\)においては、\( C_d(R_e)|_{R_e\to 0}\to \frac{24}{R_e}\)となります。

完全な球体の抗力係数\(C_d\)における近似式は、レイノルズ数を\(R_e\)として以下の通り与えられます。

  • \(R_e\lt 2\times 10^5\)層流領域のみ
    \begin{align}
    C_d(R_e)&=\left[\frac{24}{R_e}(1+0.15 R_e^{0.687})\right]+\frac{0.42}{1+(\frac{42500}{R_e^{1.16}})},~~~\mbox{式(1)}
    \end{align}
    ただし、\(R_e \lt 2\times 10^5~~\mbox{and}~~K_n \lt 0.01 ~~\mbox{and}~~ M_a\to 0 \)
  • \(R_e\lt 2\times 10^5\)の層流領域のみ、マッハ数\(M_a\)依存性込み
    \begin{align}
    C_d(R_e)
    =\frac{24}{R_e}(1+0.15 R_e^{0.687})H_M
    +\frac{0.42 C_M}{1+(\frac{42500}{R_e^{1.16C_M}})+\frac{G_M}{R_e^{0.5}}} ,~~~\mbox{式(2)}
    \end{align}
    ここで、\(C_M, G_M, H_M\)は以下の通り与えられます。
    \(\displaystyle
    \begin{eqnarray}
    C_M=\left\{~~
    \begin{aligned}
    &1.65+0.65\tanh(4M_a-3.4)~~~~~\mbox{for} ~~~M_a\lt 1.5\\
    &2.18-0.13\tanh(0.9M_a-2.7)~~~\mbox{for} ~~~M_a\gt 1.5
    \end{aligned}
    \right.
    \end{eqnarray}
    \)

  • \(\displaystyle
    \begin{eqnarray}
    G_M=\left\{~~
    \begin{aligned}
    &166M_a^3+3.29M_a^2-10.9M_a+20~~~\mbox{for} ~~~M_a\lt 0.8\\
    &5+40M_a^{-3}~~~\hspace{9.4em}\mbox{for} ~~~ M_a\gt 0.8
    \end{aligned}
    \right.
    \end{eqnarray}
    \)

    \(\displaystyle
    \begin{eqnarray}
    H_M=\left\{~~
    \begin{aligned}
    &0.0239M_a^3+0.212M_a^2-0.074M_a+1~~~\mbox{for} ~~~ M_a\lt 1\\
    &0.93+\frac{1}{3.5+M_a^5}~~~\hspace{8.8em}\mbox{for} ~~~ M_a\gt 1
    \end{aligned}
    \right.
    \end{eqnarray}
    \)

  • 定性的、広範囲
    \begin{align}
    C_d
    =\frac{24}{R_e}
    +\frac{2.6\left(\frac{R_e}{5.0}\right)}{1+\left(\frac{R_e}{5.0}\right)^{1.52}}
    + \frac{0.411\left(\frac{R_e}{263000}\right)^{-7.94}}{1+\left(\frac{R_e}{263000}\right)^{-8.00}}
    + \left(\frac{0.25\frac{R_e}{10^6}}{1+\frac{R_e}{10^6}}\right),~~~\mbox{式(3)}
    \end{align}

対象物が流体中を高速で動かなければ、層流領域で実験値を良く合う式(1)(図中の赤線), 式(2)で十分でしょう。
温度なども含めたいならば、PDFで説明しているレイノルズ数、マッハ数を置き換えて利用することを勧めます。
また精度はあまりいらない状況で、定性的な理解だけを対象とするならば、乱流領域の振る舞いを含む式(3)(図中の水色線)を用いると良いです。

以降、詳細は
https://slpr.sakura.ne.jp/qp/supplement_data/drag_coefficient_air/drag_coefficient.pdf
をご覧ください。

fortran90でファイルの存在を確かめる

fortran90です。
ファイルが存在するかしないかを確かめるためのプログラムです。
組み込み関数”access”を用いると簡単に確認できます。

program main
  implicit none
  integer::access
 
  write(6,*)access( "./temp.d", " ")
 
  stop
end program main

上の文は、もしもファイル “./temp.d” が今いるディレクトリに存在するかどうかを返すものです。

また、access関数はファイルの有無(引数は空白文字” “)のみならず、そのファイルの読み取り権限(引数は”r”)、書き込み権限(引数は”w”)、実行権限(引数は”x”)を調べることができます。

ゼロが返ってきたらその権限(ファイルの有無)があるということです。
実行例は

 $ ls
main.f90  temp.d
 $ gfortran main.f90
 $ ./a.out
           0
 $

となります。intel®fortranコンパイラでももちろんできます。

[adsense1]

ここより以下のものはaccess関数を用いない方法です。


ファイルをopenできたら存在する、errorが返ってきたら存在しない、と判断します。
そのプログラムとサブルーチンはこんな感じ。

program main
  implicit none
  integer::ox
 
  call checkfile("temp.d",ox)
 
  write(6,*)ox
 
  stop
end program main

subroutine checkfile(fname,ox)
  implicit none
  character(*),intent(in)::fname
  integer,intent(out)::ox
 
  open(100,file=trim(fname),status='old',err=999)
  close(100)
  write(6,'(3A)')"file '",fname,"' exist"
  ox=1
  return
 
999 continue
  close(100)
  write(6,'(3A)')"file '",fname,"' don't exist"
  ox=0

  return
end subroutine checkfile

”ox”は〇×のつもりで使っています。
ファイル名”temp.d”をinputパラメータとして渡し、そのファイルが存在したらox=1を、存在しなかったらox=0を返します。

実行例は、

 $ ls
   temp.d main.f90
 $ gfortran main.f90
 $ ./a.out
   file 'temp.d' exist
           1
 $ gfortran main.f90
 $ ./a.out
   file 'temp1.d' don't exist
           0

となります。2回目の実行ではtemp.dをtemp1.dに変えて存在しない場合を表示しています。

[adsense2]

参考
ファイルの有無を調べる -知識の箱

弾道計算(BB弾)の結果2、違う重さでの比較

本稿の主要な結果は、様々なパラメータでの、BB弾の軌道の詳細なデータです。

この結果に付随して、
0.25g0.8Jの軌道は
0.20g1.2Jの軌道と同じであることが分かりました。

上下振れ幅が最小になる軌道における、重さとエネルギー、ゼロイン位置を様々にとった時のそれぞれのデータです。
弾道計算を数値計算によって行い、結果を考察することが本稿の目的となります。
※本稿では規定のエネルギーを超える場合のデータがありますが、このデータは全て数値シミュレーションであるため、こういったエアガンを持っている訳ではありません。

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ←今ここ
弾道計算(BB弾)のコード(fortran90)
バレル内部でのBB弾の方程式
水中下でのBB弾の弾道計算


様々なパラメータの最適な軌道


ここで載せるデータは、以下の組み合わせです。

BB弾の重さ: 0.20g , 0.25g
エネルギー:0.60J 0.65J 0.70J 0.75J 0.80J 0.85J 0.90J
ゼロイン位置:25m 30m 35m 40m 45m 50m 55m 60m(0.25gのみ)
各々でどのような軌道を描くかの計算結果を載せます。
破線はホップ無しの軌道に対応しています。

↓これは低画質です。高画質版はこのリンク(4.3MB)を踏んでください。
弾道詳細データ_高画質2 - コピー_c


最適な軌道に合わせるためには下向きに初速を与えます。
その時、下向きに何度傾けて撃てばいいか、のデータはこちらです。

ゼロインと下向き角度


前提


まず、規定速度うんぬんは別にしまして、事実を述べます。

  1. BB弾は重いほど弾はまっすぐ飛ぶ
  2. 射出速度が早いほど弾はまっすぐ飛ぶ
  3. ホップを掛けるほど弾は上下に揺れる

これらは紛れもない事実です。
よって、BB弾を飛ばす最高の条件とは、
「射出速度を早くし、ホップは余りかけないで重いBB弾を使うこと」
となります。
だからこそ、射出速度を速めようとして規定速度の話になります。

最適な軌道とは?


BB弾を飛ばすうえで最適な軌道とはどんな軌道でしょう?
それはホップによる上下方向の揺れを最小限に抑える軌道です。
この最適な軌道とは、BB弾の重さとエネルギーとゼロイン位置を決めた時の理想的な軌道、ということです。

ゼロイン位置と上下の振れ幅とは何かは、下の画像をご覧ください。
0.25,0.8_compressed

重さとエネルギーを決めても、どの角度で射出すればいいか、回転数はいくつか、など他のパラメータが残ります。
それを決めるため、ホップによる上下方向の揺れを最小限に抑える軌道を定めます。

この軌道を ”最適な軌道” と呼ぶことにします。

重さの違いによる軌道の具体的な影響


今度は、ゼロイン位置を固定します。
そして、BB弾の重さとエネルギーを変化させたとき、どのような軌道をたどるか見てみましょう。
zeroin50_2_cc

この結果から分かることは2点あります。

1つ目は、エネルギーを上げていくと軌道の変化が小さくなる
2つ目は、重いBB弾の軌道ほどまっすぐ飛ぶ

ということです。

エネルギーを上げると軌道の変化が小さくなる

0.20gの場合でゼロインを50mに合わせ、エネルギーを0.1Jずつ増やしていったとき、上下振れ幅がどのくらい減少していくのか見てみましょう。

0.60J→0.70J … 14.2cm
0.70J→0.80J … 10.4cm
0.80J→0.90J … 8.0cm
0.90J→1.00J … 6.2cm
1.00J→1.10J … 5.0cm
1.10J→1.20J … 4.2cm
1.20J→1.30J … 3.6cm
1.30J→1.40J … 3.0cm

となります。エネルギーを変えるよりも重さを変えるべきです。
それだけで軌道は大きく変わります。

例えば、0.80Jで射出できるエアガンで、使うBB弾を0.25gで射出すると、この時、軌道は0.20g1.2Jの軌道にほぼ一致します。
m020J12m025J08_2_c

まぁ、軌道自体は一致するのですが、残念ながら到達時間は一致しません。
約0.1秒(まばたき程度)、0.20g1.2Jのほうが早くなります。
zeroin50_2_time_c

追記)
0.25gで0.90Jは、0.20gの1.35Jの軌道におおよそ相当します。
0.30gで0.80Jは、0.20gの1.60J, 0.25gの1.05Jの軌道におおよそ相当します。
0.30gで0.90Jは、0.20gの1.80J, 0.25gの1.20Jの軌道におおよそ相当します。

0.25gで0.8Jでうまく合わせられたら50m飛ばしても上下の振れ幅40cmです。
上に20cm、下に20cmしかずれません。
また、エネルギーを上げても空気抵抗の強さは速度の2乗に比例して強くなるため、それに見合うようなの飛距離の伸びは得られません。


補足


「流速チューン」というものがあるそうですね。
流速チューンの問題点と改善案について -週休5日

流速チューン比較テスト その1「初速変化」 -Gunsmithバトン
にあるように。

本稿ではこれらのことについては全く触れていません。
というのも、銃口から射出された直後の初速と回転数の2つだけが軌道に影響するためであり、発射するまでの過程に何があったか?などどうでもいいのです。
数値計算の観点から言いますと、通常の軌道よりも「まっすぐ遠くまで飛ぶ」場合、初速が上がっていること以外には考えられません
BB弾を飛んでいる最中で加速させる要因など無いのです。


※1
なぜBB弾を使う限り初速を稼ぐことが無駄なのか?
これは空気抵抗が関係します。
空気抵抗の大きさは弾道計算(BB弾)の理論で書いたように、BB弾の半径\(R\)と速度\(v\)にのみ依存します。
空気抵抗力を半径\(R\)と速度\(v\)の関数として\(F_d(R,v)\)、質量\(m\)とすると、その運動方程式は
\(
m\frac{d^2\vec{r}}{dt^2}=F_d(R,v)
\)
であり、両辺をmで割れば
\(
\frac{d^2\vec{r}}{dt^2}=F_d(R,v)/m
\)
となります。BB弾の半径や、温度、空気は変えられないので、可変なパラメータは質量と速度のみです。
そして右辺だけに注目すれば、質量が大きいほど空気抵抗力があたかも小さくなるのです。
すなわち、空気抵抗力を減らそうとするにはBB弾の重さを変えるほかない。
また、初速を変えても空気抵抗は速度の2乗で効いてくるため、皮肉なことに速度が早いほど空気抵抗力が強くなっていくのです。
よって、高いエネルギーでは初速を上げてもそれに見合っただけの結果は得ることはできないのです。

弾道計算(BB弾)のコード

2016/03/13更新
2016/03/18更新
2016/03/20更新
2016/07/23更新
2016/07/25更新
弾道計算(BB弾)のfortran90での計算コードです。


2016/03/18 以降のコードに関して
コードは改変、複製、配布、営利、非営利目的に使用していただいて構いません。
ただし著作権を放棄している訳ではありません。
また、このコードを利用した場合、クレジットの表記をお願いします。


Linuxで、OSがUbuntuやLinuxMintの場合はfortranコンパイラは

sudo apt-get install gfortran

でインストールできます。
intel®fortranコンパイラとgfortranコンパイラで正常に動くことは確かめています。

もろもろのファイルをmain.f90に記述した後、
gfortranコンパイラなら

gfortran main.f90
./a.out

intel®fortranコンパイラなら

ifort main.f90
./a.out

で実行、ファイル”orbit.txt”が生成されます。

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)←今ここ
弾道計算のコード(Excel)
バレル内部でのBB弾の方程式
水中下でのBB弾の弾道計算

できること


・BB弾の弾道計算(重力、慣性抵抗、粘性抵抗、ホップ有、回転の減衰を考慮)ができます。
式で表すと
BB弾の運動方程式(式1)
\(
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=
-mg\vec{k}+\left\{-6\pi \eta R |\vec{V}|-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|^2\right\}\frac{\vec{V}}{|\vec{V}|}
-\frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{V}|\frac{\vec{V}\times\vec{\omega}}{|\vec{V}\times\vec{\omega}|}
\)

とBB弾の回転角速度の減衰を記述する方程式(式2)
\(
\begin{align}
\frac{d\omega_z}{dt}&=N_z/I \\
N_z&=\frac{\rho C_f}{2}R^3\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin^2\theta
\end{align}
\)
の2つを考慮したものになります。

・ゼロインを決めた時の最適な軌道
・上下方向が最小の振れ幅になる軌道
も出力することができます。

※BB弾の回転の減衰を考えるか考えないかは、inputパラメータ”omgdecay”により制御できます。
回転の減衰に疑問がある場合はこれを利用してください。

注意事項


・とんでもないパラメータでは計算できない場合がありますので、そこら辺はご了承ください。
・僕のプログラミング知識ではこれが限界です。
僕はこのプログラムを使用したことによって生じた、いかなる場合でも一切の責任を負いません。

主なルーチンと説明

計算理論


式1の時間発展はルンゲ-クッタ-フェールベルグ法による適応刻み幅制御を行っています。デフォルトでは8桁の精度で(※1)一致するようになっています。モジュールGBL内のtol=1d-8を変更することで精度を制御することができます。
式2の右辺の積分は、デフォルトでは15次のガウス=クロンロッド求積法を2回用いています。今回の問題に対してはあまり良い方法とは言えません。ですが、5-6桁はロンバーグ積分法による結果と一致するので、実用上問題ないと判断し、計算速度を優先しました。
式2の回転減衰の効果はルンゲクッタ法を行う際の連立方程式に組み込んでいません。これだけはオイラー法で行っています。なので、1ステップ目の刻み幅をどうするかで、本来はずれないはずですが、変更するとそれが伝搬してゼロイン位置がずれます。デフォルトでは全ての1ステップ目の刻み幅は統一してあるので上記問題は起こりません。本当に厳密にやるならば、コードを変更して減衰を連立方程式内に組み込んでください。
2016/03/13更新
コードを変更して減衰を連立方程式内に組み込みました。計算は若干遅くなりますが、上記問題はもう起こり得ません。

※1 相対誤差もしくは絶対誤差で、です。この理由は、相対誤差に統一してしまうと値が非常に小さいとき、機械誤差によっていつまでたっても収束しないためです。なので、元になる値が1より小さいときは絶対誤差を、1より大きいときは相対誤差を取るようにしています。

2016/03/13更新
2016/03/18更新
2016/03/20更新
2016/07/23更新
2016/07/25更新
2017/10/07更新

inputファイルとコード


inputというファイル名のファイルを作り、中に

&input
m=0.25d-3,          ! Weight of bullet [kg]
a=3d-3,             ! Radious of bullet [m]
energy=0.90d0,      ! Energy of bullet [J]
g=9.80665d0,        ! Gravity [m s^-2]
temperature=20d0,   ! Temperature [degree of Celsius]
pressure=101325d0,  ! Atmospheric pressure [Pa]
moisture=0.6d0,     ! Value of moisture 0~1
omgdecay="yes",     ! Introduce omega decay. "yes" or "no"
search_zeroin="yes", ! Search good vz and omgy. "yes" or "no"
zeroinx=50d0,       !   +-- if yes, set zeroin length [m]
search_updown="yes",! Search good vz and omgy. "yes" or "no"
updownz=0.20d0,     !   +-- if yes, set up-down henght [m]
theta=-1.7d0,
rx=0.d0,
ry=0.d0,  
rz=1.0d0,
vy=0.d0,
ux=0.d0,
uy=0.d0,
uz=0.d0,
omgx=0.d0, !*2pi
omgy=-340d0, !*2pi
omgz=0.d0, !*2pi
outputfile="orbit", ! Output file name prefix
stept=0.01d0,       ! Output interval [s]
ik=0, ! ik=0 exact, ik=1 approximation
&end

と記述してください。このinputファイルの例では、

重さ(m)0.25g
エネルギー(energy)0.90J
重力加速度(g)9.80665m/(s^2)
温度(temperature)20度
大気圧(pressure)101325Pa
湿度(moisture)60%
回転角速度の減衰(omgdecay)有り
ゼロインによる最適な弾道の探索(search_zeroin)しない
↑がyesなら、ゼロインの位置(zeroinx)50m
上下幅による最適な弾道の探索(search)する
↑がyesなら、上下幅(updownz)0.1m
銃口の位置(rz)1m
y方向の速度(vy)0m/s
射出角度(theta)-0.5度
風速(ux=uy=uz=0)0m/s
ホップ(omgy)-210回転/s
出力ファイル名の接頭語(outputfile)orbit
出力の時間間隔(stept)0.01秒ごとに出力
回転の減衰を表す積分の近似方法(ik=0 : 積分の直接計算、ik=1 : 積分を近似)

となっています。
回転の減衰を考えたくない時は、omgdecay=”no”としてください。

そして下のコードを展開し、コピー&ペーストして一つのファイル(例えばmain.f90という名前)に入れてください。
(カーソルをmodule GBL…のmに合わせた後、一番下までスクロールしてshift+クリックで楽に選択できます。)
個人的にファイルをダウンロードして…は嫌いなのでこういう形で公開しています。

2016/03/13 更新
2016/03/18 更新
2016/03/20 更新
2016/07/08 更新
2016/07/23 更新
2016/07/25 更新
2016/08/06 更新
ソースコードは1523行あります。


使い方


・fortranコンパイラ(ifortやgfortran)
・gnuplot
を想定しています。

2017/02/04追)
Bernadotte66さんの報告より、windows上でgfortranコンパイラを用いる際、gfortranコンパイラがver5.1以前のものである場合エラーが出てしまうようです。ver6.2以上のコンパイラを推奨します。

Linuxでは、gfortranコンパイラver4.8.4で問題ないことは確かめています。

それでコンパイルして実行してください。
一連の流れとしては、以下の通りです。

$ gfortran main.f90
$ ./a.out
&INPUT
 M=  2.5000000000000001E-004,
 A=  3.0000000000000001E-003,
 ENERGY= 0.90000000000000002     ,
 G=  9.8066499999999994     ,
 TEMPERATURE=  20.000000000000000     ,
 PRESSURE=  101325.00000000000     ,
 MOISTURE= 0.59999999999999998     ,
 OMGDECAY="yes                                             ",
 SEARCH_ZEROIN="yes                                             ",
 ZEROINX=  50.000000000000000     ,
 SEARCH_UPDOWN="yes                                             ",
 UPDOWNZ= 0.20000000000000001     ,
 THETA=-0.50000000000000000     ,
 RX=  0.0000000000000000     ,
 RY=  0.0000000000000000     ,
 RZ=  1.0000000000000000     ,
 VY=  0.0000000000000000     ,
 UX=  0.0000000000000000     ,
 UY=  0.0000000000000000     ,
 UZ=  0.0000000000000000     ,
 OMGX=  0.0000000000000000     ,
 OMGY= -210.00000000000000     ,
 OMGZ=  0.0000000000000000     ,
 STEPT=  1.0000000000000000E-002,
 OUTPUTFILE="orbit                                           ",
 IK=          0,
 /
Reach ground at  0.980E+00
===============================
 search vz and wy by zeroin -->   50.00000
 ---- false_position
 Sequence 1 of 25 : omgy/2pi  -767.96000000 : zeroin   100.00000000
 Sequence 2 of 25 : omgy/2pi  -735.92000000 : zeroin   100.00000000
 Sequence 3 of 25 : omgy/2pi  -703.88000000 : zeroin   100.00000000
 Sequence 4 of 25 : omgy/2pi  -671.84000000 : zeroin   100.00000000
 Sequence 5 of 25 : omgy/2pi  -639.80000000 : zeroin   100.00000000
 Sequence 6 of 25 : omgy/2pi  -607.76000000 : zeroin   100.00000000
 Sequence 7 of 25 : omgy/2pi  -575.72000000 : zeroin    17.22935680
 Sequence 8 of 25 : omgy/2pi  -543.68000000 : zeroin    15.16813299
 Sequence 9 of 25 : omgy/2pi  -511.64000000 : zeroin    12.90696782
 Sequence 10 of 25 : omgy/2pi  -479.60000000 : zeroin    10.42356698
 Sequence 11 of 25 : omgy/2pi  -447.56000000 : zeroin     7.68923851
 Sequence 12 of 25 : omgy/2pi  -415.52000000 : zeroin     4.66717441
 Sequence 13 of 25 : omgy/2pi  -383.48000000 : zeroin     1.30910154
 Sequence 14 of 25 : omgy/2pi  -351.44000000 : zeroin    -2.44916818
  +--- Root find, go to conversion phase --+  
Conversion =>       0.00033361%
Conversion =>       0.00902931%
Conversion =>       0.24389215%
Conversion =>       6.77148580%
Conversion =>     104.78145034%
 ---- false_position
 Sequence 1 of 1 : omgy/2pi  -366.89192543 : zeroin    -0.57095843
  +--- Root find, go to conversion phase --+  
Conversion =>      0.00075051%
Conversion =>      0.08150935%
Conversion =>      8.77790138%
Conversion =>    579.19923386%
  -371.80132940800553       -2.5094934272766087    
Reach ground at  0.170E+01
===============================
 search vz and wy by up-down height -->    0.20000
 ---- false_position
 Sequence 1 of 25 : omgy/2pi  -959.96000000 : zeroin    19.21215678
 Sequence 2 of 25 : omgy/2pi  -919.92000000 : zeroin    17.82618788
 Sequence 3 of 25 : omgy/2pi  -879.88000000 : zeroin    16.40387340
 Sequence 4 of 25 : omgy/2pi  -839.84000000 : zeroin    14.95914779
 Sequence 5 of 25 : omgy/2pi  -799.80000000 : zeroin     9.80000000
 Sequence 6 of 25 : omgy/2pi  -759.76000000 : zeroin     9.80000000
 Sequence 7 of 25 : omgy/2pi  -719.72000000 : zeroin     9.80000000
 Sequence 8 of 25 : omgy/2pi  -679.68000000 : zeroin     9.80000000
 Sequence 9 of 25 : omgy/2pi  -639.64000000 : zeroin     9.80000000
 Sequence 10 of 25 : omgy/2pi  -599.60000000 : zeroin     9.80000000
 Sequence 11 of 25 : omgy/2pi  -559.56000000 : zeroin     1.45287636
 Sequence 12 of 25 : omgy/2pi  -519.52000000 : zeroin     1.10739833
 Sequence 13 of 25 : omgy/2pi  -479.48000000 : zeroin     0.80320434
 Sequence 14 of 25 : omgy/2pi  -439.44000000 : zeroin     0.54089117
 Sequence 15 of 25 : omgy/2pi  -399.40000000 : zeroin     0.32062896
 Sequence 16 of 25 : omgy/2pi  -359.36000000 : zeroin     0.14208536
 Sequence 17 of 25 : omgy/2pi  -319.32000000 : zeroin     0.00427472
 Sequence 18 of 25 : omgy/2pi  -279.28000000 : zeroin    -0.09464709
  +--- Root find, go to conversion phase --+  
Conversion =>       0.00183551%
Conversion =>       0.00989716%
Conversion =>       0.05297720%
Conversion =>       0.28399231%
Conversion =>       1.52644398%
Conversion =>       8.26417579%
Conversion =>      39.67268127%
Conversion =>     179.58494188%
 ---- false_position
 Sequence 1 of 1 : omgy/2pi  -312.86039405 : zeroin    -0.01405310
  +--- Root find, go to conversion phase --+  
Conversion =>      0.00065868%
Conversion =>      0.03115036%
Conversion =>      1.49063057%
Conversion =>    102.25144566%
  -317.78761751909860       -1.5473682498931858    
Reach ground at  0.149E+01
     2.681[CPU sec]
$ gnuplot

    G N U P L O T
    Version 4.6 patchlevel 4    last modified 2013-10-02
    Build System: Linux x86_64

    Copyright (C) 1986-1993, 1998, 2004, 2007-2013
    Thomas Williams, Colin Kelley and many others

    gnuplot home:     http://www.gnuplot.info
    faq, bugs, etc:   type "help FAQ"
    immediate help:   type "help"  (plot window: hit 'h')

Terminal type set to 'wxt'
gnuplot> plot "orbit.txt" u 2:4 w l
gnuplot> replot "orbit_h.txt" u 2:4 w l
gnuplot> replot "orbit_opt.txt" u 2:4 w l

gnuplot上での出力が
orbitGK
のようになったら成功です。

出力ファイルは最大で3つ生成(inputパラメータ “search_??????” )に依存します。
orbit.txtは、inputによって与えられたパラメータでそのまま計算したもの
orbit_opt.txtはそのエネルギー、重さで上下振れ幅が最小でゼロイン位置がinputのzeroinxのもの
orbit_h.txtはそのエネルギー、重さで上下振れ幅が最小で上下振れ幅がinputのupdownzのもの
となります。

出力ファイルはそれぞれ14列あり、
時刻t[秒],位置x,y,z[m]:x,y,z方向の速度[m/s]:BB弾の回転数x,y,z[回転/s]:x,y,z方向の風速:エネルギー[J]
になっています。
出力ファイルの中身は下のように書いてあるので、見てください。

# m:    0.25000E-3[kg]
# g:       9.80665[m s^{-2}]
# a:    0.30000E-2[m]
# temperature:      20.00000[degree Celsius]
# moisture:       0.60000[no-dimension]
# pressure:  101325.00000[Pa]
# eta:  0.18197E-4[kg m^{-1} s^{-1}]
# rho:     1.20123[kg m^{-3}]
# tol:  0.10000E-7
# eps:  0.10000E-8
#         t[s]         x[m]         y[m]         z[m]      vx[m/s]      vy[m/s]      vz[m/s]    wx[rot/s]    wy[rot/s]    wz[rot/s]   windx[m/s]   windy[m/s]   windz[m/s]    Energy[J]
     0.000000     0.000000     0.000000     1.000000    84.849583     0.000000    -0.740471     0.000000  -210.000000     0.000000     0.000000     0.000000     0.000000     0.900000
     0.010000     0.838863     0.000000     0.992911    82.937695     0.000000    -0.678265     0.000000  -209.261019     0.000000     0.000000     0.000000     0.000000     0.859890
    ...

実銃のパラメータでの計算


実銃は特に想定していませんが、計算ができると言えばできます。
実銃と比較して考慮されていない点は3つあり、

  1. 滑らかな球として考えているので抗力係数\(C_d\)が異なる
  2. 銃弾の歳差運動が考慮されていない(球ではなく細長いので)
  3. ライフリングによる回転が異なる
  4. コリオリ力を考慮に入れていない

1番目はかなり重大です。
私のプログラムを使うと1000メートル地点でのエネルギーの値が2倍ほど大きく出ます。
ライフリングによる横回転だけなら問題なく考慮できます。バックスピンで入っていたものを変えればいいだけです。
ただし、風が無い、またはコリオリ力を考慮しない状況では計算上は横回転があっても無くても変わりません。
風とコリオリ力がある時にライフリングによる弾道のずれが影響します。
2番目はどんな力が働いたりするか私には見当がつきません。専門の方じゃないと分かるまでに時間がかかると思います。
コリオリ力はプログラムで入れればいいです。これはそこまで辛くないでしょう。
ライフリングの効果は入れられるには入れられますが、ちょっと変な結果になるので使わないほうが良いでしょう。
なのでこの効果は入れていません。
パッと思いつく限り、違う点は以上です。

よって、本プログラムによる実銃を想定した計算とは、
風が無く、赤道上で実銃と同じ口径の滑らかな球が実銃のエネルギーで射出されたときの弾道
となります。

ちなみに、.338ラプアマグナム弾の場合を想定したinputパラメータはこちら
角度を使って計算しているのでそこのとこのプログラムは各自直してください。

&input
! For .338 lapua magnum( but treated as sphire)
m=16.2d-3,          ! Weight of bullet [kg]
a=4.29d-3,      ! Radious of bullet [m]
energy=6562.16d0,      ! Energy of bullet [J]
g=9.80665d0,        ! Acceleration of Gravity [m s^-2]
temperature=20d0,   ! Temperature [degree of Celsius]
pressure=101325d0,  ! Atmospheric pressure [Pa]
moisture=0.6d0,     ! Value of moisture 0~1
omgdecay="no",     ! Introduce omega decay. "yes" or "no"
search_zeroin="no", ! Search good vz and omgy. "yes" or "no"
zeroinx=50d0,       !   +-- if yes, set zeroin length [m]
search_updown="no",! Search good vz and omgy. "yes" or "no"
updownz=0.10d0,     !   +-- if yes, set up-down henght [m]
theta=0.6d0,        ! default, theta is neglected. change program main.
rx=0.d0,
ry=0.d0,  
rz=1.0d0,
vy=0.d0,
vz=0.d0,
ux=0.d0,
uy=0.d0,
uz=0.d0,
omgx=0.d0, !*2pi [omgx*2*pi rotate/s]
omgy=0.d0, !*2pi [omgy*2*pi rotate/s]
omgz=0.d0, !*2pi
outputfile="orbit", ! Output file name prefix
stept=0.01d0,       ! Output interval [s]
&end

上記パラメータで計算するとこんな軌道が得られます(点は0.05秒ごとの軌跡)
lapua338_c

.338ラプアマグナム弾のデータは↓
.338 ラプアマグナム弾 / .338 Lapua Magnum -MEDIAGUN DATABASE
を参照しました。

改良版


これから最適な回転量、角度を探すプログラムに更新したいと思いますが、とりあえず
改良したプログラムを載せておきます。
input
main.f90

ffmpegでavi→gif変換

LinuxOSの場合です。

ffmpegを使ってaviからgifアニメーションを作る場合のコマンドです。

ffmpegのインストール手順はCompile FFmpeg on Ubuntu, Debian, or Mintここ。

オリジナルのまま、aviからgifに変換したい場合

ffmpeg -i video.avi -t 10 out.gif

サイズを800×400にしてgifで出力させたい場合

ffmpeg -i video.avi -t 10 -s 800x400 out.gif

参考


How to generate gif from avi using ffmpeg? [closed]

弾道計算(BB弾)の結果

弾道計算(BB弾)の数値計算結果を中心にまとめます。

目次

  1. 目的
  2. 結論
  3. 設定した条件
  4. 実測データとの比較
  5. BB弾の様々な軌道
  6. 「良い軌道」の詳細なデータ
  7. 動画と静止画
  8. 「良い軌道」へ調節するには
  9. 温度の違いについて
  10. 100m飛ぶ?
  11. 僕が実際に試した時の報告
  12. お礼
  13. アイディア募集
  14. 参考文献

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果←今ここ
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)
弾道計算のコード(Excel)
バレル内部でのBB弾の方程式
水中下でのBB弾の弾道計算

目的


屋外で行われるサバイバルゲームで優位に立つために、BB弾は重い球が良いのか、軽い球が良いのかを数値計算によって,良い軌道を描く弾道の軌道のパラメータを探索して、現実で調節するための方法を模索します。

結論


サバゲーを行うにあたり、風があったり、遠くを狙う必要のある屋外では

BB弾は重い方が良い

という結論が得られました。
BB弾は重ければ重いほどまっすぐ、遠くまで飛ぶ。また、当たった時に判定がされやすい事になります。
屋内などの風が無く、接近戦(30m以内)が主になる場合ではほぼ差は見られないため、お値段的に軽い球を使うのがいいと思います。

重さの違いにありがちな話は本当?


「0.20gのほうが着弾が早い」
   →ほとんど嘘
    初速は確かに0.20gのほうが早いですが、最大でも15mの地点で0.01秒程度の違いしかありません。30m以上では0.25gのほうが早く着弾します。

「0.25gはまっすぐ飛ぶ」
   →本当
    →風の影響を受けにくいため、横にぶれにくいだけでなく、縦にもぶれが少くなります。

「0.20gのほうが遠くに飛ぶ」
   →
    →0.25gのほうが遠くまで飛びます。

「0.25gのほうがヒットを取りやすい」
   →本当
    →弾の持つエネルギーは0.25gの方がどの地点でも大きくなるため、当たった時に判定されやすいことになります。

まとめると、
0.20gが優れているのは値段だけ
ということになります。

[adsense1]
以下、詳細なデータを載せます。

設定した条件


ニュートン力学の範囲内で計算します。
考慮した力と条件は、

  1. 重力
  2. 空気の慣性抵抗力、粘性抵抗力
  3. 回転による揚力(Magnus力)
  4. BB弾の回転の減衰

です。コリオリ力は重力に比べ3桁程度小さい力なので考慮に入れていません。
空気の粘性率は、空気の密度は気温20度、乾燥した空気中のものを使用しています。
僕が調べた範囲で他のページではBB弾の回転の減衰が考慮されていませんでした。本稿ではそれを考慮したことが一つの特徴です。
運動方程式の導出や、理論に関してはこちらの弾道計算(BB弾)の理論をご覧ください。

実測データとの比較


計算が実測のデータと合っているのかをまず初めに検証します。
東京マルイ G36CにおけるノーマルホップとG-HOPの弾道特性実測データと、本稿の数値計算結果との比較を行います。

青線、実測結果は黒の点と線です。
※実測データにおいて、BB弾の回転数は測定されていなかったため(難しいため)、最高到達高さが2.4mになるように数値計算の回転数を決めました。

全体的に非常に良い一致を見せていることが分かるかと思います。50m,60m地点においては10~20cm程度の差しかないのではないでしょうか。この実測データの再現、という点においては十分に信頼できる数値計算結果になっていると思います。
ただし、打ち出し直後の浮き上がり具合に若干の差があります。
このずれの原因として、実測データの高低さの加味が十分ではないか、射出角度が若干違う、ということが考えられます。
実測時ではなく、改めて高低差を加味した、と実測データのページにあるので、場所が完全に同じではないかもしれないと考えられます。
風で揺らいでいる、ということは無いでしょう。もしも風で揺らいでいるとしたら、遠くで差が出るはずです。射出直後でこれだけ揺らいで遠くで揺らがないというのは納得ができません。

本稿の数値計算は実測データと十分に一致している、と確認したうえで、いろいろとパラメータを変えて計算してみました。

BB弾の様々な軌道


まず、一番気になるBB弾の軌道を見ます。
エネルギーを0.90Jで固定して、パラメータを様々に変え0.20gと0.25gのBB弾で比較してみます。
その結果がこちらです。赤文字が0.20g, 青文字が0.25gであることを意味しています。

上記グラフは、角度0, 1.2, 2.4°、回転数240, 280, 320回転/秒の組み合わせによる軌道を計算したものです。

BB弾が描く良い軌道とは

  • まっすぐ飛ぶ
  • 遠くまで飛ぶ

軌道であることだ、と考えます。つまり、1つの「最適な軌道」として「上下振れ幅が最小になる軌道」だと言えるでしょう。
上の画像の例でいうならば太線で示した
0.20g → BB弾が280回転/秒で回転し、下向きに1.2°方向に射出されたとき
0.25g → BB弾が320回転/sで回転し、下向きに1.2°方向に射出されたとき

がそれに近い軌道でしょう。

0.25gでこの「良い軌道」に合わせることができたなら、50m先の敵までは、喉~胸あたりに照準を合わせることで確実にヒットが得られるでしょう。
もう少し、この「良い軌道」について調べてみます。

「良い軌道」の詳細なデータ


ここでは、重さの違う0.20gと0.25gとを比較するために、ゼロインを50[m]に合わせ、上下方向の振れ幅が最小になる時で比較します。
ゼロインを50mに合わせた時、0.20gでは60cmの振れ幅、0.25gでは40cm程度の振れ幅になります。
0.20gよりも0.25gの方が小さい振れ幅です。予想通りといえば予想通りの結果ですね。

次に着弾時刻を見てみましょう。30m地点までは0.20gの方が若干早く着弾し、それ以降は0.25gのほうが早く着弾します。50mの遠距離で比較しますと、0.20gと0.25gとの着弾時刻の差は約0.13秒。まばたき程度の差があります。
ちなみに15mでの着弾の時間差は0.01秒です。とても微々たるもので、この差は認識できないと思われます。よって着弾時刻も上下の振れ幅も少ない、遠距離にも対応できる0.25gを使うのがいいでしょう。

また、エネルギーは発射した瞬間から着弾するまで、どこでも0.25gの方が高いため、遠くでも判定がされやすいことになります。

遠距離を狙いたい場合は着弾速度、上下振れ幅に影響するだけでなく、判定もされやすい重い球を使いましょう。


2016/03/13
計算コードを変更し、動画化しました。点の間隔は0.04秒で、gifアニメの速度の時間間隔と実際の時間間隔は同じにしてあります(すなわち、gifアニメの1秒と現実の1秒が同じ)。
bullet_orbit
こちらは静止画。
orbit_st_c

ゼロイン距離と上下方向振れ幅との関係


さて、上下方向の振れ幅を最小限に抑えながらゼロインを適当な位置に合わせたい場合について考えます。
遠くにゼロインを合わせたい時、飛ばせばと飛ばすほどホップを強くかけなければなりません。その結果、上にも下にもブレることになります。
この関係性がどうなっているかを調べますと、以下のようになります。
zeroin_width2_c

エネルギー0.9Jの場合で、上下の振れ幅を最小に抑えるパラメータで考えますと、段々と上下方向の振れ幅が大きくなっていることがわかります。
これは空気抵抗によって弾速が減少するためにその分、ホップによって時間を稼がなければならないことに由来します。
例えば、
0.20gのBB弾で50mにゼロインを合わせた時は上下振れ幅が60cm(射出位置から30cm下がり、ホップによって射出位置から30cm上昇し、ゼロインを迎える)ですし、
0.25gのBB弾で50mにゼロインを合わせた時は上下振れ幅が40cm(射出位置から20cm下がり、ホップによって射出位置から20cm上昇し、ゼロインを迎える)です。
この重さの違いによる上下方向の振れ幅は顕著に表れます
60mなんかにゼロインを合わせようとすると、0.20gでは上下方向に1.6mずれることになり、0.25gでは上下方向に1mずれることになります。
もしかしたら、0.20gで60mに合わせた場合、フィールドで立って撃たない限り、先に地面についてしまってダメになるでしょう。

これらから、遠距離ではますます重いBB弾を使ったほうが良いという結果が得られるわけです。

※ただし、中距離に当てるつもりがない場合はホップをある程度かけて上向きに射出したほうがいいでしょう。

「良い軌道」へ調節するには


さて、サバゲーに行き、シューティングレンジで調節をすることを考えます。
BB弾の回転数など知っているわけがありません。また、角度1.2°などそんなこと分かるわけがありません。
さらに、BB弾がどこまで飛んだかは大抵の場合は10m単位でしか分かるはずがないのです。
せめてジュール数(弾速)は知っているとします。この状況で良い軌道へ調節することを考えましょう。

0.25g, 0.9J(初速84.85m/s)の場合

ここでは、0.25g, 0.9J(初速84.85[m/s])の場合を考えます。

50mにゼロインを合わせることを考えます。
調節方法_c

この場合の手順は、
手順1,    射出する銃口の高さ50m地点にある同じ高さにある的に照準を合わせる。

手順2,    ホップ無しで、50m先の的に照準に合わせたままBB弾が地面に落ちた時の位置が、おおよそ
    \(10\times(\mbox{射出高さ[m]}) +10[\mathrm{m}]\)
    になるように調節する。

手順3,    その照準に合わせたまま、徐々にホップを掛けていき、50[m]の的に当たるまでホップを掛けて調節する。

です。

詳しく説明すると、
手順1では角度を決めるための基準点を探しています。
手順2では角度をホップ無しの場合で、到達距離を見ることによって決定します。何点か計算したら、射出する高さと良い角度の時の落下地点は、\(10\times(\mbox{射出高さ[m]}) +10[\mathrm{m}]\)という関係でした。
なので、地面からの銃口の位置が分かれば、良い角度での射出は、到達距離を見れば分かるのです。
※空気抵抗のない場合、\((\mbox{落下距離})\propto \sqrt{射出高さ}\)になります。ただし、短い区間であれば線形補間で良いでしょう。ちゃんと線形補間すると、\(10.40\times(\mbox{射出高さ}[\mathrm{m}])+9.87[\mathrm{m}]\)となるのですが、実用上は簡単な\(10\times(\mbox{射出高さ[m]}) +10[\mathrm{m}]\)で十分でしょう。

手順3では、その照準のまま、だんだんホップを掛けていき、50mに合わせればいいのです。

0.20g, 0.533J(初速73m/s)の場合

もう一つ具体例、0.20g, 73m/sの場合を考えます。
これは、僕が持っている電動のmp7a1のなにもカスタムしていない時の値です。
調節方法m0200533J_c
50mに合わせると上に40cm、下に40cmも上下するので、無理に狙わず、40mで合わせています。
(※上のグラフとはエネルギーが異なるため、”ゼロイン距離と上下方向振れ幅との関係”で紹介したグラフとは上下振れ幅が異なっています。)
この時のパラメータは、エネルギー0.533J, 下向きに角度1.81度、回転数42回転/sです。
mp7a1はこちら。

2丁持っている場合


2丁銃を持っている場合、精密射撃用と遠距離用と分けることができます。
30m内では上下方向のブレが6cm程度になるので、近接用と割り切りってしまい、
遠距離では上のパラメータで調節しておけば、どんな距離でも気持ち良く狙えます。
例えば、基本的に上のパラメータで調節した銃を使っておき、相手が隠れた時など、壁の小さい隙間を狙って近距離用で精密射撃を行う、なんてことができるかと思います。

温度の違いについて


夏場と冬場ではどのくらい弾道が変わるのでしょうか?
0.25g, 0.9Jで、温度10℃、20℃、30℃で撃った時の軌道の違いを考えます。
青が10℃、黒が20℃、赤が30℃です。
基準は20℃です。20℃でゼロイン50mの時の最小の振れ幅に合わせた時に他の温度ではどうなるかを考えます。

結論として、10℃変わるごとに約1m飛距離が変わります。
なので夏場の方が冬場よりも2m程度よく飛ぶことになります。

2019/01/23追記)

100m飛ぶ?


弾道計算本の公開から1年経過しました。
発射に至るまでの計算結果に引き続き、少しだけ中身を公開します。
ここで公開するのはBB弾は100m飛ぶのか?という点です。
結論は、
0.90Jで0.30gまでのBB弾では100m飛ぶことは無い
です。
どんなに回転を掛けても飛びません。
下のグラフは、0.90Jで、高さ1.0mから射出した時の着弾点のグラフで、
横軸:回転数、縦軸:射出角度
で表したものです。

枠内の数字は地面に到達するまでの飛距離を等高線で表しています。
現実にエアガンを撃った時に掛けられる回転数は、大きく見積もっても各画像左下の白い点線内です。

ちなみに、最大飛距離を出す射出角度で回転を掛けた場合の弾道は全く現実的ではなく、0.20g(上の画像左上)では、

の軌道を描きます。(r/s は 回転/秒を表します。)

僕が実際に試した時の感想


2015年12月、内部をカスタムしていないmp7a1でopsに行き、実際に撃ってきて理論通りの軌道になるだろうか、と試してきました。
レーザーサイト、その他もろもろは持っていないので考察ではなく、感想としてみてください。
僕個人の感想をまとめると、

  • 0.25gと0.20gの差は明らかに感じられる。特に集弾性が素晴らしい。
  • 最小の振れ幅になる軌道は確かにその通りになるが、少なくとも起伏の激しいフィールドでは使いにくい。

でした。
0.20gと0.25gの差は明らかだと感じました。軌道に対しては体感では少しだけ分かる程度でしたが、集弾性がとにかくいいのです。着弾点でばらけることがほぼ無いのです。実際にどの程度影響していたか分かりませんが、0.20gから0.25gに変えた時、キル数が上がりました。
とても気持ちよく戦えるので、初心者の方を誘ったりする際は是非0.25gを使わせて、サバゲーマーの仲間入りをさせてください。

一度下がってから上がる、これは理論上では一番良い軌道であることは確かなのですが、起伏の激しいフィールドでは、この”一度下がる”が非常に厄介です。一度下がるため、そこに障害物(例えば斜面等越しに当てるとき)があったら当たらないのです。少し上向きに撃とうものなら強ホップのため、敵の頭上を通り越します。
結果、敵から反撃を受け、こちらからは当てられない状況が生まれました。
起伏の激しいフィールドでは、30~35m程度まで数cmの振れ幅でまっすぐに進む軌道が一番良いと思います。
それよりも遠い場合は山なりに何発か撃てばいいので、非常に適していると感じました。

今度は平坦なフィールドで試してみたいと思います。

[adsense2]

お礼


@hthsi アイディア提供。実際にサバゲーをやっている上で気になること、○○の違いが知りたい!といったアイディアを弾道計算を行う初期の頃からいっぱいいただきました。重さの違いにありがちな話はほんと?の部分があるといいよね、等の分かりやすい説明は@hthsiの提案です。ありがとうございます!
また、サバゲーるにて告知をしていただきました。
0.25g弾の優位性について

募集


2015/07/06
現在も計算を行い、まとめている最中です。気になること、面白いアイディア、計算に関する質問等がありましたらコメントでもtwitter(@sig_colon)でも、メールでも構いませんのでお知らせください。
計算するまでもない質問の場合はコメントで返信、特に面白そうな場合は、綺麗にまとめて本稿に追加いたします。

2016/02/07
Tech Report  ハイスピードカメラによるエアソフトガンの考察
より、ガスガンと電動ガンの銃口付近での空気の流れが見て取れます。
非常に素晴らしく、実際の状況なので信頼できます。15000fpsの高速度カメラだそうです。

ひとつ、気になるのは報告されているBB弾の回転数が私の導き出した値と2倍近く違うことです。
回転数が違うことによって、数値計算結果の回転数以外に影響はありません。ただのホップの強さの目安と見ているので、もしも間違いがあったとしてもホップの強さと回転数との換算が違うだけです。結果には影響しません。
管理人様にお尋ねしたところ、回転数は管理人様自身が導き出したものではなく、考察程度としてみるのがよい、と返答いただきました。残念ながら確証が取れないので回転数の議論はどなたかやっていただけるまで待つことにします。
動画があり、非常に説得力のあるサイトなので、是非ともご覧いただくことをお勧めします。

参考文献


東京マルイ G36CにおけるノーマルホップとG-HOPの弾道特性

2016/03/22 追)
ガンジニア様のハイスピードカメラによる映像


より、本サイト内で表記している回転数にミスがある疑いが濃厚になってきました。
サイト内に「○○回転している」や「○○回転/s」と言う表記があるのはたぶんです。
これから検証していきますのでお待ちください。

2016/07/09追)この問題は解決しました。詳細はBB弾の回転量について(実験との比較)をご覧ください。

コメント返信用、画像

2015/07/27

ウエダさん宛て
縦軸:高さ[m]
横軸:到達距離[m]
ウエダさんパラメータ_c