sikino のすべての投稿

弾道計算, 物理学

球体の抗力係数

1件のコメント

説明

本稿では、完全な球体に働く抗力係数について述べます。
ここでは概要だけを説明します。詳細は以下の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

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

[adsense2]

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

サバイバルゲーム, 弾道計算, 物理学

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

9件のコメント

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

この結果に付随して、
0.25gで0.8Jの軌道は
0.20gで1.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.20gの1.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乗に比例して強くなるため、それに見合うようなの飛距離の伸びは得られません。

G&G バイオBB弾 0.25g 1KG 4000発入り 【G07126】

東京マルイ 0.25g　ベアリングバイオBB弾 生分解バイオ 1300発入り【３袋セット】

補足

「流速チューン」というものがあるそうですね。
流速チューンの問題点と改善案について -週休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弾)のコード

17件のコメント

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

Linux

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弾)の結果

11件のコメント

弾道計算(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弾が描く良い軌道とは

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

軌道であることだ、と考えます。つまり、１つの「最適な軌道」として「上下振れ幅が最小になる軌道」だと言えるでしょう。
上の画像の例でいうならば太線で示した
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はこちら。

No1 H&K MP7A1 (18歳以上電動コンパクトマシンガン)

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

物理学

Schrödinger, Heisenberg, Interaction描像

コメントする

Schrödinger描像、Heisenberg描像、Interaction描像というものがあります。

Schrödinger描像は波動関数による表現方法で、多くの場合で分かりやすい表現なため良く取り入れられます。
ただし、時と場合によってHeisenberg描像の方が分かりやすかったり、Interaction描像の方が最適、というときがあります。

表現方法が違うだけで同じものを表現します。

  1. Schrödinger描像(S描像)

    \(\hat{H}\)が時間依存しない場合、時刻\(t\)での波動関数は初期状態(\(t=0\))の波動関数\(\Psi_S(0)\)を用いて、
    \(
    \displaystyle \Psi_S(t)=e^{-i\frac{\hat{H}}{\hbar}t}{\Psi_s(0)}\ \ \ \ (1)
    \)
    と書くことができます。ここで添え字”\(S\)“はS描像であることを強調しています。
    \(\displaystyle \hat{U}(t)＝e^{-i\frac{\hat{H}}{\hbar}t}\ \ \ (2)\)
    と表記されることが多く、時間発展演算子と呼ばれています。
    今、(2)の両辺を時間微分すると、\(\hat{U}(t)\)の満たす微分方程式
    \(
    \displaystyle i\hbar \frac{d\hat{U}(t)}{dt}=\hat{H}\hat{U}(t), \ \ \ \hat{U}(0)=1\ \ \ (3)
    \)
    が得られます。また、式(3)は実は\(\hat{H}\)が時間依存していても成立します。
    この意味は、式(3)が時間発展演算子のより一般的な記述方法であることを意味しています。

    今、Schrödinger描像で、とある演算子\(\hat{A}\)の期待値を考えます
    (\(\hat{A}_S\)は時間依存してもしなくてもok)。
    (例えば\(\hat{A}_S\)は位置空間で、位置演算子だったら\(\hat{x}=x\)、運動量演算子だったら\(\hat{p}=-i\hbar\frac{d}{dx}\)です。)
    すると、その期待値を表す関数は時間依存し、時刻\(t\)の\(\hat{A}\)の期待値\(\langle\hat{A}\rangle_t\)は
    \(
    \begin{align}
    \langle\hat{A}\rangle_t
    &\displaystyle =\langle{\Psi_S}(t)|\hat{A}_S|\Psi_S(t)\rangle \\
    &\displaystyle =\langle e^{-i\frac{\hat{H}}{\hbar}t} \Psi_S(0)|\hat{A}_S|e^{-i\frac{\hat{H}}{\hbar}t}\Psi_s(0)\rangle
    \end{align}
    \)
    と書けます。これがSchrödinger描像における演算子の期待値を表現しています。

  2. Heisenberg描像(H描像)

    時刻\(t\)でのHeisenberg描像の波動関数\(\Psi_H(t)\)の定義は
    \(
    \begin{align}
    \Psi_H(t)&=e^{i\frac{\hat{H}}{\hbar}t}\Psi_S(t) \\
    &=e^{i\frac{\hat{H}}{\hbar}t}e^{-i\frac{\hat{H}}{\hbar}t}\Psi_S(0) \\
    &=\Psi_S(0)
    \end{align}
    \)
    です。
    Heisenberg描像で、\(\hat{A}\)の期待値\(\langle\hat{A}\rangle_t\)を考えましょう。
    \(
    \begin{align}
    \langle\hat{A}\rangle_t&=\langle\Psi_S(t)|\hat{A}_S|\Psi_S(t)\rangle \\
    &\displaystyle =\langle\Psi_S(t)|e^{-i\frac{\hat{H}}{\hbar}t}e^{i\frac{\hat{H}}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}}{\hbar}t}e^{i\frac{\hat{H}}{\hbar}t}|\Psi_S(t)\rangle \\
    &\displaystyle =\langle\Psi_H(0)|\hat{A}_H(t)|\Psi_H(0)\rangle
    \end{align}
    \)
    すなわち、S描像の演算子\(\hat{A}_S\)との間には、
    \(
    \hat{A}_H(t)=e^{i\frac{\hat{H}}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}}{\hbar}t}
    \)
    という関係があるのです。

    では、H描像で演算子\(\hat{A}_H(t)\)の時間依存性はどうなるのでしょう。
    上の式を両辺微分して\(i\hbar\)を掛けると、
    \(
    \begin{align}
    i\hbar\frac{d}{dt}\hat{A}_H(t)&=
    -\hat{H}e^{i\frac{\hat{H}}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}}{\hbar}t}
    +e^{i\frac{\hat{H}}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}}{\hbar}t}\hat{H} \\
    &=-\hat{H}\hat{A}_H(t)+\hat{A}_H(t)\hat{H} \\
    i\hbar\frac{d}{dt}\hat{A}_H(t)&=\left[\hat{A}_H(t), \hat{H}\right]
    \end{align}
    \)

    もしも、\(\hat{A}_S\)が時間依存する形だったら
    \(
    \displaystyle i\hbar\frac{d}{dt}\hat{A}_H(t)=\left[\hat{A}_H(t), \hat{H}\right]+i\hbar \frac{\partial \hat{A}_H(t)}{\partial t}
    \)
    と書けます。

  3. Interaction描像(I描像)

    前提がありまして、系のハミルトニアン\(\hat{H}\)を\(\hat{H}=\hat{H}_0+\hat{V}\)と分けます。
    そして、\(\hat{H}_0\)の時間発展は分かっているものとします。

    この時、時刻\(t\)でのInteraction描像の波動関数\(\Psi_I(t)\)の定義は
    \(
    \begin{align}
    \Psi_I(t) &= e^{i\frac{\hat{H_0}}{\hbar}t}\Psi_S(t) \\
    &(= e^{i\frac{\hat{H_0}}{\hbar}t}e^{-i\frac{\hat{H}}{\hbar}t}\Psi_S(0) )
    \end{align}
    \)
    です。
    このInteraction描像の波動関数の意味を簡単に説明すると、
    \(\Psi_I(t)\)は、ある時刻で、S描像の波動関数を\(e^{-i\frac{\hat{H}}{\hbar}t}\)に従って時間を進める方向に時間発展させる部分と、\(e^{i\frac{\hat{H}_0}{\hbar}t}\)に従って時間を戻す方向に時間発展させる部分の２つがある、ということです。
    詳しくは後ほど述べます。

    I描像での演算子\(\hat{A}\)の期待値\(\langle\hat{A}\rangle_t\)はどうなるのでしょう。
    \(
    \begin{align}
    \langle\hat{A}\rangle_t&=\langle\Psi_S(t)|\hat{A}_S|\Psi_S(t)\rangle \\
    &\displaystyle =\langle\Psi_S(t)|e^{-i\frac{\hat{H}_0}{\hbar}t}e^{i\frac{\hat{H}_0}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}_0}{\hbar}t}e^{i\frac{\hat{H}_0}{\hbar}t}|\Psi_S(t)\rangle \\
    &\displaystyle =\langle\Psi_I(t)|\hat{A}_I(t)|\Psi_I(t)\rangle
    \end{align}
    \)
    と書けます。

    すなわち、I描像では演算子が時間依存し、波動関数も時間依存します。
    よって、
    Interaction描像での演算子\(\hat{A}_I(t)\)の時間依存を記述する方程式と、
    Interaction描像での波動関数\(\Psi_I(t)\)の時間依存を記述する方程式
    があり、それらを求めてみましょう。

    • 演算子\(\hat{A}_I(t)\)の時間依存

      \(
      \hat{A}_I(t)=e^{i\frac{\hat{H}_0}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}_0}{\hbar}t}
      \)
      演算子の期待値は上式の通りであり、もしも\(\hat{H}_0\)が時間依存しない場合、Heisenberg描像で\(\hat{H}\to \hat{H}_0\)と置き換えたものと同じ形になります。すなわち、
      \(
      \displaystyle i\hbar\frac{d}{dt}\hat{A}_I(t)=\left[\hat{A}_I(t), \hat{H}_0\right]
      \)
      となります。もしも\(\hat{H}_0\)が時間依存する場合は、
      \(
      \displaystyle i\hbar\frac{d}{dt}\hat{A}_I(t)=\left[\hat{A}_I(t), \hat{H}_0\right]+i\hbar \frac{\partial \hat{A}_I(t)}{\partial t}
      \)

    では、I描像で波動関数\(\Psi_I(t)\)が満たす運動方程式はどうなるでしょう。
    I描像の波動関数の定義式に左から\(e^{-i\frac{\hat{H_0}}{\hbar}t}\)を掛けると、
    \(
    \Psi_S(t) = e^{-i\frac{\hat{H_0}}{\hbar}t}\Psi_I(t)
    \)
    であり、これをS描像のSchrödinger方程式に代入します。
    \(
    \begin{align}
    & i\hbar \frac{\partial}{\partial t}\left(e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t)\right)
    =\hat{H}e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t) \\
    & i\hbar \left[-\frac{i}{\hbar}\hat{H}_0 e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t)+e^{-i\frac{\hat{H}_0}{\hbar}t}\frac{\partial \Psi_I(t)}{\partial t}\right]
    =\hat{H}_0e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t)+\hat{V}e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t) \\
    & i\hbar e^{-i\frac{\hat{H}_0}{\hbar}t}\frac{\partial \Psi_I(t)}{\partial t}=\hat{V}e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t) \\
    \end{align}
    \)
    左から\(e^{i\frac{\hat{H}_0}{\hbar}t}\)を作用させると、
    \(
    \displaystyle i\hbar \frac{\partial \Psi_I(t)}{\partial t}=e^{i\frac{\hat{H}_0}{\hbar}t}\hat{V}e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t)
    \)
    ここで、Interaction描像の演算子\(\hat{H}_I(t)\)を
    \(
    \displaystyle \hat{H}_I(t)=e^{i\frac{\hat{H}_0}{\hbar}t}\hat{V}e^{-i\frac{\hat{H}_0}{\hbar}t}
    \)
    と置くと、波動関数\(\Psi_I(t)\)が満たす運動方程式
    \(
    \displaystyle i\hbar \frac{\partial \Psi_I(t)}{\partial t}=\hat{H}_I(t)\Psi_I(t)
    \)
    が得られます。

さて、定義が終わったところで、Heisenberg描像とInteraction描像の意味を見てみましょう。
簡単にそれぞれの描像を説明すると、

Heisenberg描像：量子-古典対応を成す
Interaction描像：既に分かっている部分を取り除く

という特徴があります。

まずはHeisenberg描像から。
Heisenberg描像でHeisenberg描像の演算子\(\hat{A}_H(t)\)の時間発展を記述する運動方程式は
\(
i\hbar\frac{d}{dt}\hat{A}_H(t)=\left[\hat{A}_H(t), \hat{H}\right] \cdots(A)
\)
です。
今、\(\hat{A}_H(t)\)が位置演算子\(\hat{q}\)と運動量演算子\(\hat{p}\)によって表記されていると考えます。
この時、Schrödinger描像で、それぞれの演算子はハミルトニアンと
\(
\begin{align}
\left[\hat{q}_S,\hat{H}_S\right] &= i\hbar \frac{\hat{p}_S}{m} \\
\left[\hat{p}_S,\hat{H}_S\right] &= -i\hbar \frac{\partial}{\partial \hat{q}_S}\hat{H}_S
\end{align}
\)
という関係があります。
Heisenberg描像でSchrödinger描像のハミルトニアンがどう記述されるのかを考えてみると、
\(
\begin{align}
\hat{H}_S &= \frac{\hat{p}_S^2}{2m}+\hat{V}_S \\
&= \frac{1}{2m}e^{-i\frac{\hat{H}_S}{\hbar}t}\hat{p}_H^2 e^{i\frac{\hat{H}_S}{\hbar}t}+e^{-i\frac{\hat{H}_S}{\hbar}t}\hat{V}_H^2 e^{i\frac{\hat{H}_S}{\hbar}t} \\
&= e^{-i\frac{\hat{H}_S}{\hbar}t}\left( \frac{\hat{p}_H^2}{2m}+\hat{V}_H\right) e^{i\frac{\hat{H}_S}{\hbar}t} \\
\hat{H}_H&= e^{i\frac{\hat{H}_S}{\hbar}t}\hat{H}_S e^{-i\frac{\hat{H}_S}{\hbar}t} \\
&=\hat{H}_S
\end{align}
\)
となり、結局
\(
\hat{H}_S=\hat{H}_H
\)
であることが導かれます。
Heisenberg描像での位置演算子\(\hat{q}\)と運動量演算子\(\hat{p}\)の交換関係は
\(
\begin{align}
\left[\hat{q}_H,\hat{H}_H\right] &= i\hbar \frac{\hat{p}_H}{m} \\
\left[\hat{p}_H,\hat{H}_H\right] &= -i\hbar \frac{\partial}{\partial \hat{q}_H}\hat{H}_H \\
\end{align}
\)
であることが導けるため、\(\hat{q}_H\)と\(\hat{p}_H\)に関して、運動方程式に代入して
\(
\begin{align}
i\hbar\frac{d}{dt}\hat{p}_H(t)&=\left[\hat{p}_H(t), \hat{H}_H\right] \\
&= -i\hbar \frac{\partial}{\partial \hat{q}_H}\hat{H}_H \\
\frac{d}{dt}\hat{p}_H(t)&=-\frac{\partial}{\partial \hat{q}_H}\hat{H}_H
\end{align}
\)
であり、
\(
\begin{align}
i\hbar\frac{d}{dt}\hat{q}_H(t)&=\left[\hat{q}_H(t), \hat{H}_H\right] \\
&= i\hbar \frac{\hat{p}_H}{m} \\
\frac{d}{dt}\hat{q}_H(t)&=\frac{\partial}{\partial \hat{p}_H}\hat{H}_H
\end{align}
\)
となります。
この二つの式は、古典力学のハミルトンの運動方程式に酷似していますよね。ここからHeisenberg描像が量子-古典対応を成す、という所以になるわけです。

注意として、Heisenberg描像は、量子-古典対応を成しますが、演算子それ自身を測定することなどできないのです。あくまで期待値と対応するだけです。
また、Heisenberg描像は古典力学への回帰を見出し、エーレンファストの定理と関連します。

続いて相互作用描像です。
系のハミルトニアン\(\hat{H}\)を\(\hat{H}=\hat{H}_0+\hat{V}\)と分けます。
そして、\(\hat{H}_0\)の時間発展は分かっている時、Interaction描像の波動関数\(\Psi_I(t)\)は
\(
\begin{align}
\Psi_I(t) = e^{i\frac{\hat{H_0}}{\hbar}t}e^{-i\frac{\hat{H}}{\hbar}t}\Psi_S(0)
\end{align}
\)
でした。じっくり見ていきましょう。
例を考えます。
\(
\displaystyle \hat{H}=\hat{H}_0+V(x)=-\frac{\hbar^2}{2m}\nabla^2+V(x)
\)
で、ポテンシャル\(V(x)\)が単調減少で図のような山なりの波形をしている場合を考えます。
Interaction描像イラスト_c
初期状態\(\Psi_S(0)\)があったとすると時間経過すると、その波動関数の形は広がりながら右へ移動します。
Interaction描像では、この、広がる、という動きと、右へ移動する、という動きが合わさったと考えます。
広がる動きは、ポテンシャルが無いときの波束の動きです。ポテンシャルがない時、右に動くことはしません。
右への動きは、ポテンシャルしか無いときの波束の動きです。質点を考えれば、ポテンシャルの山を転がり、右へ動きます。

Interaction描像の波動関数\(\Psi_I(t)\)は既に分かっている広がる動き(\(\hat{H}_0=-\frac{\hbar^2}{2m}\nabla^2\)の部分)を取り除き、波束が右に動く描像(\(V(x)\)の部分)だけを切り抜くのです。

言葉で言えば、Interaction描像は、
\(t’\)秒後の\(I\)描像における波動関数は、\(t’\)秒後の\(S\)描像における波動関数に対して、\(t’\)秒間の\(\hat{H}_0\)による逆方向の時間発展を行えばよい、ということです。

あくまでもイメージなので、細かな議論はご了承を。

参考文献

David J. Tannor著、山下晃一ほか訳『入門 量子ダイナミクス（上）』(2011), p.231~236

入門 量子ダイナミクス(上) 時間依存の量子力学を中心に

gnuplot

gnuplotのカラーマップ

コメントする

配色(pm3dの色)について
gnuplot ver4.6で確かめています。

一番良いと思ったカラーマップは、matlabで使われている”jet”だと感じました。

set palette defined ( 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

colorjet
参考：Matlab colorbar with Gnuplot

※一番下に黒または白を追加した以下の物も良いと思います。

set palette defined ( -1 '#000030', 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

set palette defined ( -1 '#ffffff', 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

※正負をはっきりと見せたい時。

set palette defined ( 0 '#0fffee',1 '#0090ff', 2 '#000fff',3 '#000090',4 '#ffffff',5 '#7f0000', 6 '#ee0000', 7 '#ff7000', 8 '#ffee00')

set palette defined ( 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#ffffff',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

追記）2023/01/13
matlabでは、最近parulaというカラーマップが採用されているようです。

set palette defined ( 0 '#352a87',1 '#2053d4',2 '#0d75dc',3 '#0c93d2',4 '#07a9c2',5 '#38b99e',6 '#7cbf7b',7 '#b7bd64',8 '#f1b94a',9 '#fad32a',10 '#f9fb0e')

グレースケールにした時にも画像を変更する必要もなく、このままで綺麗に出るようです。
また、\(f(x,y)=\sin(xy)\)をjetと比較してみましょう。

グレースケールとした場合は、以下のようになります。

確かにparulaの方が高い低いがはっきりと分かります。

カラーの場合は偏りが生じないjetの方が綺麗ですかね？
しかし、parulaの方が柔らかな印象です。
グレースケールで印刷されることが予想されるならば、parulaの方がいいですね。
印象が大きく変わるかと思います。好きな方をご使用ください。

配色について。
gnuplotでは色の指定方法として3種類あるそうです((2),(3),(4))。

お勧めのカラーの指定順に,
(1)外部サイトgnuplottingにあるカラーマップを利用する
(2)cubehelixで指定する
(3)definedで指定する
(4)rgbformulaeで指定する

方法を紹介します。

(1)gnuplottingにあるカラーマップの利用

一番早く、綺麗なカラーマップを得るためには”gnuplotting(http://www.gnuplotting.org/)”という外部サイトに紹介されているgnuplotのカラーマップを使うことです。
ここに紹介されているカラーマップは、https://github.com/Gnuplotting/gnuplot-palettesより画像を引用して、
overview((c)gnuplotting)です。

使い方
  1. Gnuplotting/gnuplot-palettes -githubに飛びます。
  2. 右上にある “Clone or download” をクリックしてすべてのファイルをダウンロード
  3. カラーマップを見て使いたい名前の “***.pal” ファイルをgnuplot上でload
    で使えます。線の色なども上書きされるので、もしも線の色だけを変えたければ、linestyleなどで変更する必要があります。

(2)cubehelixで指定する

gnuplot version4.6以上

set palette cubehelix start 0.5 cycles -1.5 saturation 1

この色の指定方法はD.A.Greenによって発表された、宇宙の、強度イメージを色付けするために作られた色の指定方法です[1]。

それぞれのパラメータが何を表すかは折りたたんでおきますので、気になる方は下のリンクを開いてください。

[adsense1]

cubehelixでの具体例

対象とするデータは、水素原子の(n,l,m)=(4,2,0)の解析解です。
※このデータが欲しい方はこちらHydrogen_xz_nlm420.txt(1.5MB)に置きましたので、コピペしてください。
Ctrl+aで全選択できます。

3Dの表示ではこういうグラフです。
hydrogennlm420

特に記述が無いものはgammaの値はデフォルトの1.5です。

1.hydro420_s0.5c-1.5s1g1.5_c

set palette cubehelix start 0.5 cycles -1.5 saturation 1

きちんとした裏付けの下カラーリングされているので、論文等々にお勧めです。

2.hydro420_s0.5c-1s3g1.5_c

set palette cubehelix start 0.5 cycles -1 saturation 3

もしくは

 set palette defined(0"#000000",1"#5b13b8",2"#1c60ff",3"#00b7d8",4"#01eb75",5"#55f02e",6"#ccd73b",7"#ffc68c",8"#ffd5e2",9"#ffffff")

3.hydro420_s0.5c-1s3g3_c

set palette cubehelix start 0.5 cycles -1 saturation 3
set palette gamma 3

もしくは↑に近い色で

set palette defined(0"#000000",0.1"#8e007a",2.1"#579dff",3.1"#2de4ff",4.1"#4bffab",5.1"#98fe7d",6.1"#e4eb8c",7.1"#ffe1c0",8.1"#ffe9f0",9.1"#ffffff")

4.hydro420_s1c0s2g1.5_c

set palette cubehelix start 1 cycles 0 saturation 2

5.hydro420_s2c0s3g1.5_c

set palette cubehelix start 2 cycles 0 saturation 3

6.hydro420_s2c1s2g1.5_c

set palette cubehelix start 2 cycles 1 saturation 2

7.hydro420_s3c0.5s3g1.5_c

set palette cubehelix start 3 cycles 0.5 saturation 3

8.hydro420_s1c5s3g1.5_c

set palette cubehelix start 1 cycles 5 saturation 3

9.hydro420_s1c0s2g1.5neg_c

set palette cubehelix start 1 cycles 0 saturation 2 negative

10.hydro420_s1c-1s2g1.5neg_c

set palette cubehelix start 1 cycles -1 saturation 2 negative

11.hydro420_s3c0.5s1g1.5neg_c

set palette cubehelix start 3 cycles 0.5 saturation 1 negative

12.hydro420_s3c-2s2g1.5neg_c

set palette cubehelix start 3 cycles -2 saturation 2 negative

13.hydro420_s1c5s3g1.5neg_c

set palette cubehelix start 1 cycles 5 saturation 3 negative

82.54.4 Cubehelix

(3)definedで指定する

set palette defined(-3"blue",0"white",5"red",10"yellow")

z軸が-3になったとき青色(blue)、0になったとき白(white),5になったとき赤(red)、10になったとき黄色 (yellow・になるような配色です。
zの範囲が-50 < z < 100だった場合は-3,0,5,10の比で表されます。
すなわち、z=-50で青、z=-16付近で白、z=42付近で赤、z=100で黄色となるように自動的に調整されます。
色は有名な色は上の例のように言葉で記述できますが、16進数のカラーコードによる表示もできます。

カラーマップ例

カラーマップ例を載せます。

set palette defined ( 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

matlabで使われている”jet”のカラーマップです。
colorjet
参考：Matlab colorbar with Gnuplot

set palette defined(0"#ff0000",1"#ff8000",2"#ffff00",3"#80ff00",4"#00ff00",5"#00ff80",6"#00ffff",7"#0080ff",8"#0000ff",9"#8000ff",10"#ff00ff",11"#ff0080",12"#ff0000")

cyclic_c
色は色相環を元に作っています。
周期的な関数を出力するときに役に立ちます。例えば、リーマン面を角度に依存させて表示する際などに利用すると良いです。こんなように。
Riemann_angle

set palette defined(1"#0000ff",2"#0080ff",3"#00ffff",4"#00ff80",5"#00ff00",6"#80ff00",7"#ffff00",8"#ff8000",9"#ff0000")

cl2_c
色相環を元に作成したものです。

set palette defined (-9 "purple",-6 "blue",-3 "turquoise", 0 "#f5f5f5", 3 "gold", 6 "orange", 9 "red")

cwhite
※匿名希望さん提供

set palette defined(0"#ffffff",0.8"#00008b",1.8"#2ca9e1",3"#008000",4.2"#ffff00",5"#eb6101",5.5"#8b0000")

カラフルtr_c
白⇒青⇒水色⇒黄色⇒オレンジ⇒赤
細かい変化を見たい時にいい感じになります。

set palette defined(0"#000000",1"#0000ff",2"#1e90ff",3"#00ffff")

黒から青tr_c
黒⇒青⇒明るい青
不連続な関数を表示させる時にこのマップを使うといい感じになります。

set palette defined(0"#ffffff",1"#0000ff",2"#ff0000",3"#ffff00")

白青黄tr_c
白⇒青⇒赤⇒黄
上の２つが気に入らなかったりした時に。

set palette defined(0"#00008b",1"#2ca9e1",2"#38b48b",3.5"#ffff00",5"#eb6101",5.3"#c9171e")

カラフルatr_c
青⇒水色⇒黄色⇒オレンジ⇒赤
一番下が白じゃないほうがいいときに。

(4) rgbformulaeで指定する

set palette rgbformulae 33,13,10

”set palette rgbformulae” のあとの３つの整数で赤(R)/緑(G)/青(B)に対応しています。
整数で赤、緑、青をどのような変化具合をさせるのかということを指定します。
直線的に色を変化させるのか、2次、3次関数的なのかはたまたsin,cos的なのかを。どういった整数を使えばどういう関数に相当するのか？
これを知るには gnuplot上で
”show palette rgbformulae”
と打ってください。gnuplot ver.4.6では、以下の画像のように出力されます。
rgbformulae_c
0~36まで出てきます。指定は-36~36までできます。負の値は反転することを意味しています。

また、いくつかのrgbformulaeを使った例が ”help palette rgbformulae” で見れます。実際に入力するとこんな感じ。
helprgb_c
helpで現れるカラーマップはこんな感じです。
rgbformulae_help_c

また、今使っているカラーマップをrgbformulaeで記述をしたい場合、

show palette fit2rgbformulae

とやれば、今使っているカラーマップをrgbformulaeで表示する場合の３つの数字を出してくれます。
あくまで似た図です。

デフォルトのカラーマップを変更する

デフォルトのカラーマップを変更するには、ホームディレクトリにて、
.gnuplot
という名前のファイルを作り、その中に

set palette defined(0"#ffffff",0.8"#00008b",1.8"#2ca9e1",3"#008000",4.2"#ffff00",5"#eb6101",5.5"#8b0000")

と記述します。これでデフォルトのカラーマップが変更されます。

ちなみに、僕が使っている.gnuplotファイルは以下のものです。

set terminal wxt dashed noraise enhanced font 'Times New Roman,20'

#matlab jet
#set palette defined ( 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

set palette defined ( -1 '#ffffff', 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

set linetype  1 lc rgb "black" lw 1
set linetype  2 lc rgb "red" lw 1
set linetype  3 lc rgb "blue" lw 1
set linetype  4 lc rgb "forest-green"
set linetype  5 lc rgb "magenta" lw 1
set linetype  6 lc rgb "#FFD900" lw 1
set linetype  7 lc rgb "cyan" lw 1
set linetype  8 lc rgb "#7F00FF" lw 1
set linetype  9 lc rgb "#FF7F00" lw 1
set linetype  10 lc rgb "#00FF7F" lw 1
set linetype  11 lc rgb "gray" lw 1
set linetype  12 lc rgb "gray40" lw 1

上の設定を適応すると、カラーマップはjetに黒を追加したもの、線は以下のようになります。

配色のあれこれ、gnuplot demoより。
カラフルに塗りたい http://gnuplot.sourceforge.net/demo/pm3dcolors.html
モノクロに塗りたい http://gnuplot.sourceforge.net/demo/pm3dgamma.html
点に色をつけたい http://gnuplot.sourceforge.net/demo/rgb_variable.html
線に色をつけたい http://gnuplot.sourceforge.net/demo/rainbow.html

http://ayapin-film.sakura.ne.jp/Gnuplot/Primer/Misc/colormap.html

[1]D. A. Green(2011), ‘A colour scheme for the display of astronomical intensity images’, Bulletin of the Astronomical Society of India, 39, 289.(2011BASI…39..289G at ADS.)

またDave Green’s `cubehelix’ colour schemeをご覧下さい。

[adsense2]

Linux

LinuxMintでアップデート失敗

コメントする

僕の環境は、
VMware(R) Player バージョン7.1.0 build-2496824
に、LinuxMint 17.1 Rebeccaをインストールしている環境です。

sudo apt-get update
sudo apt-get dist-upgrade

をした後、再起動したあと、
LinuxMintの起動後にモードを選択する画面が表示され、

udevadm trigger is not permitted while udev is unconfigured.
udevadm settle is not permitted...
...
ALERT! /dev/disk/by-uuid/....-....-....-.... -.... does not exist.
Dropping to shell!
BusyBox v1.21.1 (Ubuntu 1:1.21.0-1ubuntu1) built-in shell (ash)
Enter 'help' for a list of built-in commands.

(initramfs) _

画像ではこんな感じになりました。
20150511-51639_b_c1
※uuidを晒しても安全か分からなかったため、念のため隠しています。

こうなってしまった後の場合の復帰方法を書きます。手順を書けば、

  1. liveCD からLinuxを起動
  2. “sudo fdisk -l” を入力し、ブートディスクを確認する, 僕の場合は /dev/sda1 でした。
  3. “sudo mkdir /media/newroot” を入力
  4. “sudo mount /dev/sda1 /media/newroot”, を入力, ブートディスクを作成したディレクトリに変更
  5. “sudo chroot /media/newroot” を入力
  6. “sudo apt-get update”
  7. “sudo apt-get dist-upgrade”
  8. 再起動で幸せに!

で治りました。
Thread: “udevadm trigger is not permitted while udev is unconfigured.”より。
ただし、6,7の手順ではすんなりアップデートはできず、apt-get update入力後失敗し、
>代わりにconfigure…を入力してください
と出て、こちらを改めて入力しました。
また、dist-upgradeもすんなりはできず、オプション追加して再入力してください、とあったのでその通りにしました。

詳しく書いていきます。

  1. liveCD からLinuxを起動

    必要なのは使っているLinuxの.isoファイルです。
    LiveCDはLinuxをインストールしたときに使用した,LinuxOSが入っているCD/DVD/USBのことです。
    ですが、物理的なCDドライブは特に必要ありません。.isoイメージファイルさえあればokです。
    LinuxMint17.1の場合はisoイメージファイルを
    Download Linux Mint 17.1 Rebecca
    よりダウンロードできます。日本の場合はno-codecs版です。
    ※ちなみに、VMwareでUSBからのブートは出来ません。Biosに項目があるかもしれませんが、できません。
    ダウンロードしたら、VMwareの設定で、起動時に仮想ドライブを起動するようにします。
    仮想マシンの設定を選択して、設定します。
    20150512-164029a_ca_c
    ↑を選択してから
    20150512-213544_c_c
    を選択してokを。

    そして、問題になっているLinuxを起動し、Bios画面を表示させ、CDromから起動するように変更します。
    起動の途中でEscを押せばBios画面に入れるのですが、一瞬で読み込みが終わり、起動画面は消えるので結構シビアです。
    確実にBios画面に入るためにはVMwareをCDブート、USBブートする/VMwareのBIOS画面を出すより、このようにすればいいそうです。

    仮想OSを作成したフォルダ内の「*.vmx」ファイルに、以下の1行を追加する事で、ESCキーを受け付ける時間を延ばす事が可能です。
    BIOS.BootDelay = “5000”
    単位はミリ秒。上記の場合、5秒間の猶予が出来ます。
    2～3秒に設定していると、意外にあっという間に通過します。お勧めは4～6秒くらいです。
    余談ですが、以下の一行を追加すると、毎回BIOS画面に入れます
    bios.forceSetupOnce = “TRUE”

    そうして起動します。するとGUIでLinuxMintが起動するかと思います。そうしたらターミナルを開いて次のステップに進みます。

  2. “sudo fdisk -l” を入力し、ブートディスクを確認する

    次に、今使われているブートディスクは何なのか、を調べます。

    sudo fdisk -l

    を入力するとこの画像のように表示されることでしょう。
    20150512-162004_c
    ここで、左から２列目、”Boot” とある欄を見ますとアスタリスクが付いていることが分かります。
    なので、これが自分のBootディスク、ということになります。
    上の画面の場合は/dev/sda1, ということです。

  3. “sudo mkdir /media/newroot” を入力

    説明はありません。

  4. “sudo mount /dev/sda1 /media/newroot” を入力

    説明はありません。

  5. “sudo chroot /media/newroot” を入力

    説明はありません。

  6. “sudo apt-get update”

    ここでちょっと躓きました。僕の場合だけかもしれませんが、エラーが出されました。
    > Err http://…
    みたいな文がずらずらと。ただ、よくよくエラー文を読むと、
    代わりにconfigure…を入力してください
    とあったのでこれを入力しました。そうしたら、とりあえずokとなりました。

  7. “sudo apt-get dist-upgrade”

    次にこれを入力します。ここでもエラーが出ましたが、エラー文に
    オプションを追加して入力してください、
    とあったのでその通りにしたらokとなりました。

  8. 再起動

    再起動すれば元通り！
    この時、Bios画面か、仮想CD/DVDドライブから読み込ませないようにしておかないと意味がないので外しておきましょう。

途中、sudo に関して、

sudo: unable to resolve host

が出てきました。一応解決するためにはこのリンク先通りにやれば解決します。
これが必要な手順か、は知りません。

カーネルが、という場合もあるようです。その時は手順6,7が違うものになるそうです。

参考先

Thread: “udevadm trigger is not permitted while udev is unconfigured.”　このエラーのフォーラム
Download Linux Mint 17.1 RebeccaLinux Mint 17.1 Rebeccaのダウンロード
VMwareをCDブート、USBブートする/VMwareのBIOS画面を出すBIOS画面を確実に出すためには。
sudo: unable to resolve host が表示されたらsudoの解決
と
とある詳しい方。

物理学

演算子の種類と説明

3件のコメント

\(\hat{A}\)がエルミート演算子ならば、\(e^{i\hat{A}}\)はユニタリー演算子である。

この文の意味が分かる人はこのページはいらないと思います。

前提として、一度勉強した人が思い出す、という体を想定しています。
詳しくブラケット表記であるとか、正しく知りたい人は、
J. J. Sakurai著 桜井 明夫訳『現代の量子力学〈上〉』 (物理学叢書) (1989)

現代の量子力学(上) 第2版 (物理学叢書)

や
David J.Tannor著 山下晃一訳 『入門 量子ダイナミクス 時間依存の量子力学を中心に（上）』 化学同人
などを参考にしてください。このページの参考先もこの２つです。

[adsense1]

エルミート演算子(Hermite Operator)

エルミート演算子は自己随伴演算子とも呼ばれます。

  • ディラックのブラ・ケット表記で考えます。
    演算子\(\hat{A}\)の状態\(u\)と状態vによる内積を\(\langle u|\hat{A}|v\rangle \)と、表現します。
    関数による表現では、
    \(
    \displaystyle \langle u|\hat{A}|v\rangle =\int u^* \hat{A} v dx
    \)
    となるわけです。この時、
    \(
    \langle u|\hat{A}|v\rangle^*=\langle v|\hat{A}^{\dagger}|u\rangle
    \)
    を関係を満たす演算子を随伴演算子と呼び、\(\hat{A}^{\dagger}\)と表現します(\(\hat{A}\)のエルミート共役を取る、とも言います)。

    状態による表記では、
    \(
    \displaystyle \hat{A}|v\rangle =\lambda |v\rangle
    \)
    に対してエルミート共役を取ると、
    \(
    \displaystyle \langle v|\hat{A}^{\dagger} =\lambda^* \langle v|
    \)
    と書けます。

    そして、たまたま\(\hat{A}^{\dagger}\)が\(\hat{A}\)に等しい場合、すなわち、
    \(
    \hat{A}^{\dagger}=\hat{A}
    \)
    を満たすとき、演算子\(\hat{A}\)はエルミート演算子(自己随伴演算子)だ、と呼びます。

  • エルミート演算子の持つ性質
    1. エルミート演算子の固有値は実数である。
      \(
      \hat{A}|v\rangle=\lambda |v\rangle
      \)
      左から\(v\)を掛けて、内積を取り、式変形します。式変形は２通り考えられて、

      \(
      \begin{eqnarray}
      \left\{
      \begin{aligned}
      \langle v|\hat{A}|v\rangle&=\langle v|\lambda |v\rangle=\lambda\langle v|v\rangle \\
      \langle v|\hat{A}|v\rangle&=\langle v|\hat{A}^{\dagger}|v\rangle=
      \langle v|\lambda^*|v\rangle=\lambda^*\langle v|v\rangle
      \end{aligned}
      \right .
      \end{eqnarray}
      \)

      となります。同じものから出発したので値は同じものになるはずです。なので
      \(\lambda=\lambda^*\)
      これを満たす\(\lambda\)は虚数部がゼロでなければなりません。
      よってエルミート演算子の固有値は実数である、となります。

    2. あるエルミート演算子が異なる2つの固有値を持つ場合、これらの固有値に対応する固有関数は互いに直交する。
      ２つの固有値を\(\lambda_1, \lambda_2\)と書いて、
      それぞれの固有値に属する固有ベクトルを\(|v_1\rangle, |v_2\rangle\)と書くことにします。すなわち、
      \(
      \hat{A}|v_1\rangle=\lambda_1 |v_1\rangle \ ,\ \ \ \ \hat{A}|v_2\rangle=\lambda_2 |v_2\rangle
      \)
      であるとします。ここで\(\lambda_1\neq \lambda_2\)です。
      1番目の式に左から左から\(\langle v_2|\)を作用させて式変形します。式変形は２通り考えられて、
      \(
      \begin{eqnarray}
      \left\{
      \begin{aligned}
      \langle v_2|\hat{A} |v_1\rangle&=\langle v_2|\lambda|v_1\rangle=\lambda_1\langle v_2|v_1\rangle \\
      \langle v_2|\hat{A}|v_1\rangle&=\langle v_2|\hat{A}^{\dagger}|v_1\rangle=
      \lambda_2^*\langle v_2|v_1\rangle=\lambda_2\langle v_2|v_1\rangle
      \end{aligned}
      \right .
      \end{eqnarray}
      \)
      同じものから出発したので、
      \(
      \begin{align}
      \lambda_1\langle v_2|v_1\rangle=\lambda_2\langle v_2|v_1\rangle \\
      \rightarrow (\lambda_1-\lambda_2)\langle v_2|v_1\rangle=0
      \end{align}
      \)
      仮定より、\(\lambda_1\neq \lambda_2\)なので、\(\langle v_2|v_1\rangle=0\)になるほかありません。
      \(\langle v_2|v_1\rangle=0\)は内積がゼロ、すなわち直交である、と言っているので仮定は示されました。

    • ある状態\(\phi\)が、演算子\(A\)の固有値\(\lambda_i\)に属する固有状態\({v_i}\)の組で書かれるとき、
    すなわち、
    \(
    \displaystyle |\phi\rangle=\sum_ia_i|v_i\rangle
    \)
    で書かれるとき、
    \(
    \displaystyle \frac{\langle\phi|\hat{A}|\phi\rangle}{\langle\phi|\phi\rangle}
    =\sum_i p_i \lambda_i,\ \ \ p_i\equiv \frac{|a_i|^2}{\sum_i|a_i|^2}\ \ \cdots (a)
    \)

  • 量子力学の中心的教義は以下の3つを主張しています。
      (i)     全ての観測量にはエルミート演算子\(\hat{A}\)が結び付けられる。
      (ii)   \(\hat{A}\)に属する観測量の測定について、起こりうる結果は\(A\)の固有値\(\{A_i\}\)のみ。
      (iii)  系が状態\(|\phi\rangle=\sum_ia_i|v_i\rangle\)にあれば、
      \(\lambda_i\)の値を得る確率は\((a)\)で与えた\(p_i\)で与えられる。
      また、\(\displaystyle \sum_ip_i\lambda_i\)は\(\hat{A}\)を測定した平均値、もしくは期待値である。
  • おまけ
    随伴行列\(A^{\dagger}\)は、\(\hat{A}\)の転置行列の複素共役で与えられます。
    \(
    (A_{ij})^{\dagger}=A^*_{ji}
    \)
    特に、エルミート行列の場合、\(A^{\dagger}=A^*\)なので、
    \(
    A_{ij}=A^*_{ji}
    \)
    対角要素については\(A_{ii}=A^*_{ii}\)が成り立つので、エルミート行列の対角要素は実数でなければならないことがわかります。

逆演算子

\(\hat{A}\)の逆\(\hat{A}^{-1}\)を意味します。

  • 定義

    \(\hat{A}|u\rangle=|v\rangle\ \ \ \cdots (1)\)
    のとき、逆演算子\(\hat{A}^{-1}\)は
    \(|u\rangle=\hat{A}^{-1}|v\rangle\ \ \ \cdots (2)\)
    と定義されます。

  • 性質

    (1)の左から\(\hat{A}^{-1}\)を作用させると、
    \(\hat{A}^{-1}\hat{A}|u\rangle=\hat{A}^{-1}|v\rangle\ \ \ \cdots (1-1)\)
    (2)の左から\(\hat{A}\)を作用させると、
    \(\hat{A}\hat{A}^{-1}|u\rangle=\hat{A}|v\rangle\ \ \ \cdots (2-1)\)

    (1)と(2-1), (2)と(1-1)を見比べれば、明らかに
    \(\hat{A}\hat{A}^{-1}=\hat{A}^{-1}\hat{A}=\mathbf{1}\)
    となります。

  • 逆演算子の存在
    全ての演算子に逆演算子があるわけではありません。以下の通り、逆が存在しない場合があることを示せます。
    演算子\(\hat{A}\)が2つの異なる初期ベクトルを同じ終ベクトルに写す状況を考えます。式で表せば、
    \(
    \begin{eqnarray}
    \left\{
    \begin{aligned}
    \hat{A}|u_1\rangle&=|v\rangle \ \ \ (3)\\
    \hat{A}|u_2\rangle&=|v\rangle \ \ \ (4)
    \end{aligned}
    \right.
    \end{eqnarray}
    \)
    ここで、
    \(\hat{A}^{-1}|v\rangle\)を考えた時、それを\(|u_1\rangle\)か、\(|u_2\rangle\)かを決める術はありません。
    よって逆が存在しないことになります。
    また、 \((4)-(3)\)を行うと、
    \(
    (4)-(3)=\hat{A}(|u_1\rangle-|u_2\rangle)=\mathbf{0}
    \)
    となります。この意味は、\(\hat{A}\)は\(\mathbf{0}\)ではない、いくつかのベクトルを消去する、ということを表現しています。

[adsense2]

ユニタリー演算子

ユニタリー演算子は演算子\(\hat{A}\)の逆演算子\(\hat{A}^{-1}\)が
\(\hat{A}\)のエルミート共役\(\hat{A}^{\dagger}\)に等しいとき、その演算子はユニタリー演算子だ、と定義されます。すなわち、
\(
\hat{A}^{-1}=\hat{A}^{\dagger}
\)
を満たすとき、と定義されます。
また、別の表現では
\(
\hat{A}^{\dagger}\hat{A}=\hat{A}\hat{A}^{\dagger}=\mathbf{1}
\)
という場合もありますが同じことです。
ユニタリー演算子はよく、\(\hat{U}\)という表記がされます。

  • ユニタリー演算子の性質
    1. ユニタリー演算子はノルムを保存する。
      ユニタリー演算子が状態\(|u\rangle\)に作用した場合を考えます。この時、ノルムは自身の内積を取ればいいので、
      \(\langle u|\hat{A}^{\dagger}\)と\(\hat{A}|u\rangle\)を作用させればノルムになります。故に、
      \(
      \langle u|\hat{A}^{\dagger}\hat{A}|u\rangle=\langle u|u\rangle
      \)
      であり、ノルムは変わりません。
    2. ユニタリー演算子の固有値は絶対値\(\mathbf{1}\)を必ず含む。
      あるユニタリー演算子を\(\hat{U}\)と書き、固有値問題を
      \(
      \hat{U}|v\rangle=\lambda |v\rangle\ \ \ (5)
      \)
      と書きます。

      両辺のエルミート共役をとって、
      \(
      \displaystyle \langle v|\hat{U}^{\dagger} =\lambda^* \langle v|
      \)
      ユニタリー演算子の性質を使うと、
      \(
      \displaystyle \langle v|\hat{U}^{-1} =\lambda^* \langle v|
      \)
      右から\(\hat{U}\)を作用させて、
      \(
      \begin{align}
      \displaystyle \langle v|\hat{U}^{-1}\hat{U} &=\lambda^* \langle v|\hat{U} \\
      \displaystyle \langle v|\hat{U} &=\frac{1}{\lambda^*}\langle v| \ \ \ \ (6)
      \end{align}
      \)

      (5)に左から\(\langle v|\)を作用させると、
      \(
      \langle v|\hat{U}|v\rangle=\langle v|\lambda |v\rangle=\lambda\langle v |v\rangle\ \ \ (7)
      \)
      であり、また、内積を以下のように変形し、(6)を使うと、
      \(
      \displaystyle \langle v|\hat{U}|v\rangle=\frac{1}{\lambda^*}\langle v|v\rangle\ \ \ (8)
      \)
      従って、式(7),(8)は同じものから出発したので等しいはずで、
      \(
      \begin{align}
      \lambda\langle v |v\rangle&=\frac{1}{\lambda^*}\langle v|v\rangle \\
      &\rightarrow |\lambda|^2=1
      &\rightarrow |\lambda|=1
      \end{align}
      \)
      となるため、ユニタリー演算子の固有値は必ず絶対値1を含みます。
      ※これは、\(e^{i\theta},\ \ \ \theta\)は実数、であればいいと言っています。1である必要性はありません。

    3. 異なる固有値を持つユニタリー演算子の固有状態は直交する
      異なる２つの固有値を\(\lambda_1, \lambda_2\)と書き、
      それぞれの固有値に属する固有ベクトルを\(|v_1\rangle, |v_2\rangle\)と書くことにします。
      \(|v_2\rangle\)と\(|v_1\rangle\)による内積をそれぞれ考えると、
      \(
      \langle v_2|\hat{U}|v_1\rangle=\langle v_2|\lambda |v_1\rangle=\lambda_1\langle v_2 |v_1\rangle\ \ \ (9)
      \)
      と変形できるし、
      \(
      \langle v_2|\hat{U}|v_1\rangle=\frac{1}{\lambda_2^*}\langle v_2|v_1\rangle=\lambda_2\langle v_2|v_1\rangle \ \ \ (10)
      \)
      とも変形できます。最後の変形は2.の性質、絶対値1を持つことを利用しています。
      2つの固有値\(\lambda_1, \lambda_2\)は違う値を仮定したので、(9)=(10)が成り立つためには
      \(\langle v_2|v_1\rangle=0\)
      でなければなりません。よって、ユニタリー演算子の異なる固有値に属する固有状態は直交していなければなりません。

ユニタリー演算子とエルミート演算子

エルミート演算子からユニタリー演算子を作る方法を示します。
この方法は運動量演算子の導出でも用いるので、覚えておいて損はないかと思います。

\(\hat{A}\)がエルミート演算子ならば、\(e^{i\hat{A}}\)はユニタリー演算子である。
演算子\(e^{i\hat{A}}\)の逆行列を考えて、それがエルミート共役に等しくなるか、見てみます。
定数に対するエルミート共役は単なる複素共役、\(\hat{U}=\hat{U}^{\dagger}\)であることを利用すると、
\(
\begin{align}
\left(e^{i\hat{A}}\right)^{-1}&=e^{-i\hat{A}} \\
&=1+(-i\hat{A})+\frac{(-i\hat{A})^2}{2!}+\cdots \\
&=1+(i\hat{A})^{\dagger}+\frac{\left\{(i\hat{A})^{\dagger}\right\}^2}{2!}+\cdots \\
&=\left(e^{i\hat{A}}\right)^{\dagger}
\end{align}
\)
となり、ユニタリー演算子の定義\(\hat{U}^{-1}=\hat{U}^{\dagger}\)を満たしていることがわかります。

逆行列とユニタリー行列

行列表記と演算子表記は同じものです。
量子力学では、ユニタリー行列の逆行列に関心があります。数式ならば、
\(\hat{U}^{-1}_{ij}=\hat{U}^{\dagger}_{ij}=\hat{U}^{*}_{ji}\)
ということです。

ユニタリー行列の性質の一つに、ユニタリー行列の列ベクトルは正規直交ベクトルである、ことを示しましょう。
逆行列との積は\(\mathbf{1}\)に等しいので
\(
(\hat{U}^{-1}\hat{U})_{ik}=\delta_{ik}
\)
です。これを踏まえ、ユニタリー演算子の性質を使って、
\(
\begin{align}
(\hat{U}^{\dagger}\hat{U})_{ik}&=\sum_j\hat{U}_{ij}^{\dagger}\hat{U}_{jk} \\
&=\sum_j\hat{U}_{ji}^{*}\hat{U}_{jk}=\delta_{ik}
\end{align}
\)
という結果が得られます。最後の式から、行列\(U\)の2つの列ベクトルの内積、これが\(\delta_{ik}\)に等しいことがわかります。
故にユニタリー行列の列ベクトルは正規直交ベクトルであることが示されました。

ちなみに、逆行列を得るための一般的な手続きは、
①  行列を対角形に変形
↓
②  対角要素を逆数に
↓
③  逆変換
という流れで行われます。