弾道計算(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)←今ここ
バレル内部での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
を参照しました。

SurfacePro3について

Surface Pro 3 corei5 128GB の情報について。


調べ方
フリーソフト→Speccyを使いました。
自分のパソコンのCPUやHDDの温度・各ハードウェアの詳細などをわかりやすくまとめて表示するフリーソフト「Speccy」 -GIGAZINEより

ソフトウェアについて


VMwareによる仮想環境
少しもたつきますが動きます(VMwarePlayer+Linux Mint cinnamon17.1)。
Linux Mint cinnamon17.1ではなくて、軽いディストリビューションであるLinux Mint LXDEを使えば割とサクサクです。
また、CentOSはサクサク動きます。
一応タッチパネルも動作します。ただし、過信はいけません。時々変な振る舞いをします。

※仮想環境を使用している最中にキーボードを畳んだりするとその後使おうとしたとき、私だけかもしれませんが、キーボードが認識されなくなることがあります。この時には、
カーソルをwindowsにするキーボードを外して10秒位待つキーボードをつけるwindowsで認識していることを確認→カーソルをVMwareに戻す
の手順が確実です。

ちなみに私はLinuxに慣れていないのでVMwarePlayer+Linux Mint cinnamon17.1で動かしています。

ゲーム
PSO2→動きます。フル画面ではカクつきます。
設定を最低にし、windowモードで半分くらいの画面サイズにすれば割とサクサクです。

スペック


  • オペレーティングシステム

    Windows 8.1 Pro 64-bit

  • CPU

    コア数 2
    スレッド 4
    名称 Intel Core i3/i5/i7 4xxx
    コードネーム Haswell ULT
    ソケット形式 Socket 1168 BGA
    製造プロセス 22nm
    CPUの完全な名称 Intel Core i5-4300U CPU @ 1.90GHz

  • メモリ

    4.00GB デュアル-Channel DDR3 @ 798 MHz (12-15-15-34)

  • マザーボード

    Microsoft Corporation Surface Pro 3 (SOCKET 0)

  • グラフィック

    Generic PnP Monitor (2160×1440@59Hz)
    Intel HD Graphics Family (Microsoft)

  • ストレージ

    119GB SAMSUNG MZMTE128HMGR-000MV (SSD)

  • 光学ドライブ

    DTSOFT Virtual CdRom Device
    ※仮想ドライブなので、物理的にあるCDやDVD、blu-rayを読み込ませるには外付けのドライブが必要になります。

  • オーディオ

    Realtek High Definition Audio

その他、
・Wi-Fi 802.11ac/802.11 a/b/g/n
・Bluetooth® 4.0
・5.0メガpxの背面カメラ
・5.0メガpxの前面カメラ
(使用した感じ、背面カメラのほうが画質が良かったよ?)
・マイク付き
・イヤホンジャック×1
・microSDポート×1
・USB3.0ポート×1
・MiniDisplayPort×1(外部ディスプレイに繋げたい時は、このポートを使います。Macと同じものです。変換コネクタがなければ買うしかありません。)
僕はこれを買いました。MiniDisplayPort⇔D-sub15ピン変換用のものです。

また、

バッテリー

バッテリー名…X863568, LG化学製
定格電圧…7.4V
バッテリー容量42157mWh, \(\sim\)5500mAh
でした。(調べ方はここ、または下に書いたので参考にどうぞ。)
どうやら違う人もいるみたいです。

追記 2015/04/14)
8ヶ月くらい使った体感ですが、満充電にしておいてから、バッテリーのみで駆動させた場合、
WindowsMediaPlayerで音楽を流しながら、ちょこちょこパワポを使ったりとしてsurface pro3を使ったとき、バッテリーの持ちは大体4~5時間くらいです。
CPU100%でずっと動かすなら、1時間程度です。
モバイルバッテリーで充電する場合は充電口が特殊なため、ACアダプターを別に購入し、工作しないといけません。
純正のACアダプターは1万円くらいします。
また、
モバイルバッテリー→USBからSurfacePro3へ給電
はできません。


スポンサーリンク



バッテリーの調べ方

コマンドプロンプトにて

powercfg /batteryreport

と入力するとhtml形式で結果がファイルに出力され、その内容を見ることができます。
(Surface Pro 3 i5 battery capacity -Microsoft Communityより。)
それによると、僕の場合は

Installed batteries

Information about each currently installed battery

BATTERY1
NAME(バッテリーの名前) X863568
MANUFACTURER(工場) LGC_LGC(LG Chem Ltd., LG化学)
SERIAL NUMBER 1F11
CHEMISTRY LION
DESIGN CAPACITY(バッテリーの理論値) 42,157 mWh
FULL CHARGE CAPACITY(実際にバッテリーをフル充電したら) 40,250 mWh
CYCLE COUNT(何回0→100%を繰り返したか(?)バッテリーの全消耗量だと…) 138

と出ました。()内は僕の付け足しです。実際には出ません。

また、上の場合のままでは表記がmWhなのでちょっとイメージしづらいのでmAh表示に直します。
バッテリーの型番”X863568”で検索したところ、ebayでこのバッテリーを販売している人が何人かいました。
もしかしたらリンクが切れるかもしれませんが張っておきます。(2015/01/11確認)
MICROSOFT SURFACE PRO 2 GENUINE BATTERY 7.4V 5500MAH X863568-006 P21GU9 (BAT7)
MICROSOFT SURFACE PRO 2 128GB VERSION BATTERY X863568-006 7.4V 5500mAh
ここでは定格電圧7.4V,5500mAhのバッテリーということでした。

ebayで出ていた表示は本当なのか?念のため、42157mWhから(正しいであろう)定格電圧7.4Vを使ってmAhでの表示を求めてみると、
\(
\begin{align}
X[\mbox{mAh}]&=\frac{42157[\mbox{mWh}]}{7.4[\mbox{V}]} \\
X[\mbox{mAh}]&=5696[\mbox{mAh}]
\end{align}
\)
なのでいい値です。大体5500mAhなのでしょう。
第7回:スマホ時代の新必需品「モバイルバッテリー」はどれを選べば良い?とかより。

また、”CYCLE COUNT“はバッテリーの寿命の指標になるものなので、例えば中古で買ったり、友人から買わない?みたいに言われたときに出来ればチェックしましょう。
Surface Pro のバッテリーおよび電源を見ると、このバッテリーはリチウムイオン電池を使用している、とのことでした。
ほとんどのリチウムイオン電池の充電のサイクル回数は500回程度だ、と言われ、約500回でリチウムイオン電池の性能は60%に落ちるそうです。(リチウムイオン電池の概要 -電気設備の知識と技術)

僕の場合は購入してから5~6ヶ月経過して138回充電を繰り返した、ということなので、このペースだと残り寿命はおおよそ15ヶ月,1年ちょっとです。

スポンサーリンク

対角化

fortran90でLapackを用います。

倍精度ルーチン

複素一般行列の対角化, diag(N,A,Ev)
複素エルミート行列の対角化, zdiag(N,A,Ev)
実エルミート行列の対角化, herdiag(N,A,Ev), diag_slct(N,A,Nq,Ev)
実対称三重対角行列の対角化, tridiag(N,A,Ev)

4倍精度ルーチン

複素一般行列の対角化, qdiag(N,A,Ev)

複素一般行列         :行列要素は複素数で、値は適当に入っている。
複素エルミート行列      :行列要素は複素数で、\(A_{i,j}=A^*_{j,i}\)の関係がある。
実エルミート行列(実対称行列):行列要素は実数で、\(A_{i,j}=A_{j,i}\)の関係がある。
実三重対角行列        :行列要素は実数で、\(A_{i,i}\ne 0, ~~A_{i,i-1}=A_{i,i+1}\)、それ以外の要素は0。

(エルミートな行列と言われたら、その固有値は必ず実数です。)

複素一般行列の場合


wolframのページで、
Eigenvectors[{-2,2,4},{-2,4,2},{-2,1,4}]
と入力することで固有値、固有ベクトルを出してくれます。これで確かめを行いました。

・この行列を入れた時、正しい固有値と固有ベクトルを出力してくれる事を確認しています。
・40×40程度の行列でも確認しています。
・複素数の固有値、固有ベクトルでも確認しています。

Lapackを使います。

ifort -mkl main.f90

か、g95, もしくはgfortranでlapackにリンク

gfortran -llapack main.f90

して、コンパイルしてください。

大きさN×Nの正方行列Aを入力すると, 固有値Ev(1:N)とそれぞれの固有値の要素に属する固有ベクトルを上書きして返します。(1列目に固有値Ev(1)に属する固有ベクトルが、2列目に固有値Ev(2)に属する固有ベクトルが・・・という感じに。)
また固有値は実部の大小を比較して小さい順に並べ替えられ、固有ベクトルもそれに対応した列に格納されます。
計算速度は\(O(N^3)\)です。

複素エルミート行列の対角化


対称行列に限る場合は計算を早くすることができます。
ここでは複素エルミート行列を対角化して固有値を出力するケースを考えましょう。

そのまま使うのはやはり面倒です。なぜなら、このルーチンを用いるためにはワーク配列を何個か宣言して用いる必要があるからです。
このワーク配列は計算に直接関係ありません。
僕らが欲しいのは行列を入れたら固有値と固有ベクトルを返してくれさえすればいいものです。
下はサブルーチンzdiagは僕が作ったもので、複素エルミート行列\(A(1:N,1:N)\)を入れると固有値の配列\(Ev(1:N)\)と\(A\)に上書きして固有ベクトル\(A(1:N,1:N)\)を入れて返します。

エルミート行列なので、実は全ての行列要素を求める必要ありません。
下の例では、
エルミート行列\(A(i,j)\)は,\(j\ge i\)の上三角の要素だけ分かっていればいいです。
なので、コメントアウトしている行列でも同じ固有値と固有ベクトルが得られます。
以下のプログラムをlapackと共にコンパイルし、実行すると、

>./a.out   
(-0.42419E+01)
( 0.21580E+01)
( 0.80839E+01)
(-0.88841E+00,-0.92008E-01)  (-0.15371E+00,-0.13692E+00)  (-0.39557E+00, 0.58357E-01)  
( 0.17660E+00, 0.14837E-01)  ( 0.55488E+00,-0.56209E+00)  (-0.35177E+00, 0.47012E+00)  
( 0.41334E+00, 0.00000E+00)  (-0.57775E+00, 0.00000E+00)  (-0.70382E+00, 0.00000E+00)  
>

という結果が得られます。
1番目の固有値(-0.42419E+01)に対応する固有ベクトルは1列目に、
2番目の固有値( 0.21580E+01)に対応する固有ベクトルは2列目…というように代入されていきます。

プログラム例)
確かめ用のリンクEigenvectors[{-2,2+i,4},{2-i,4,1-2i},{4,1+2i,4}] -wolfram alpha

行と列を入れ替えたければ、以下のルーチンを使用してください。

subroutine swap2d(N,A)
  implicit none
  integer::N,i,j
  double precision::temp,A(1:N,1:N)
 
  do i=1,N
     do j=i+1,N
        temp=A(i,j)
        A(i,j)=A(j,i)
        A(j,i)=temp
     enddo
  enddo
 
  return
end subroutine swap2d

実エルミート行列の場合


以下のサブルーチンを使うと良いでしょう。

\(N\times N\)の実対称行列で、固有値の小さい方から順に\(Nq\)個だけ欲しい場合


実エルミート行列の、
完全対角化(herdiag)と、
選択的に固有値を計算する(diag_slct)
の比較(実時間計測)を行いました。


横軸が選択的に何個固有値を計算するか、を表し、縦軸は実際に掛かった時間をあらわします。
色は行列サイズによって変化し、黄色,紫、緑の順に行列サイズが\(2^16\times 2^16, 2^14\times 2^14, 2^12\times 2^12\)を表します。
破線はherdiagによる対角化時間、実線はdiag_slctによる計算時間を表します。
また、図中の値2,4,8は計算に用いたCPUの数です。

ここから分かるのは、欲しい固有値の数が行列サイズ\(N\times N\)の\(N/4\)より小さい場合は選択的に計算した方が早い、ということです。
選択的に計算する方は計算時に要求されるメモリも少なくて済みます。


(2016/08/27)

4倍精度の一般行列の対角化


4倍精度用の一般行列(複素非エルミートを含む)を対角化し、その固有値と固有ベクトルを出力するfortranコードが以下のページ
4倍精度対角化 -ただの備忘録
で公開されていました。
ここで公開されているサブルーチン”qeispack.f90
を使うことで対角化できます。
gfortranとifortで動くことは確かめました。
コンパイルは例えば

ifort qeispack.f90 main.f90

のようにしてください。
本ページと同じように書き加えて実装すると以下のようになります。

実対称三重対角行列の場合


実対称三角行列の場合も置いておきます。