弾道計算(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
を参照しました。


「弾道計算(BB弾)のコード」への17件のフィードバック

  1. 初めまして
    とても興味深い研究をされていますね、難しくて理解が中々追いつきませんが・・・

    私もパラメータを変えて色々調べてみたくなり
    本日gfortranとgnuplotで上記コードと実行結果のinputファイルを合わせたのですが
    どうにも計算結果が合わないのです・・・
    inputファイル読み込み後のReach ground at 0.490E+00あたりからずれてきます

    無知ですみません、どこか問題あるのでしょうか・・・

    1. 初めまして。ご連絡ありがとうございます。

      恐らく、その原因は私が1週間ほど前にプログラムを変更したからだと思います。
      その際にプログラムだけを変更し、画像やその他データは変更以前のままにしてしまいました。
      その変更に伴って現在のコードでは、最適な軌道を探すことはできなくなっているようです。
      今週末か来週末くらいには訂正いたしますので、それまでお待ちください。
      私の間違いです。お手数おかけしてすみません。

      1. sikinoさん:
        ご返答ありがとうございます
        そういう事だったんですね、調査頂き少し安心致しました

        研究の妨げにならない程度に進めていただければ結構ですので
        気長にお待ちしております

        よろしくお願いいたします

        1. purinさん

          プログラムの変更を行いました。
          問題なく動くと思うので試してみてください。

          もしも不具合がありましたらご連絡いただけると幸いです。

          1. こんにちは

            試したところ、うまく動作致しました
            ありがとうございます!

            私の銃は0.28gでは弾道を見る限り全然ホップが足りていないようです・・・調整の目指す方向が数値から見えるのはやっぱり素晴らしいです。
            また弾道データにスコープの照準ラインを追加して、ゼロイン時のホップ調整の参考に致します。

          2. purinさん

            動作してよかったです。
            スコープのラインを入れれば確かに見やすそうですね。
            恐らく線に色を付けたほうがスコープで覗くことを考えた際良いと思います。
            例えばgnuplot上では

            splot "orbit_opt.txt" u 2:3:4:2 w l lw 5 lc palette

            とすることで
            color_bullet
            というような画像が得られます。ご参考ください。

  2. 初めまして、飛翔する物体について調べたいと思い、何か良いサイトはないかと思い検索をかけていたところこのサイトに辿り着きました。
    球体に対する力に関する方程式や、RK法の解説など詳しく解説がされており、この分野についての教養を深めることができたと思います。

    弾道計算のプログラムについても公開をされていたので、試しに私も計算をシミュレーションしてみようと考えたのですが、gfortranにおいてコンパイルを行った後に、./a.outにてプログラムを実行しようとした際に

    At line 1457 of file main.f90 (unit=10,file=’./input’)
    Fortran rutime error: End of file

    と表示がされ、エラーにてストップしてしまいました。
    最初は作成をしたinputファイルに実行の許可がないために起きているのかと思い、chmodにて実行許可を出してみたりもしたのですが、実行することができませんでした。
    このエラーについて何かご存知であれば教えて頂ければと思い、コメントさせて頂きました。
    お手数ですが、お返事頂ければと思います。

    1. 初めまして。
      方程式、数値計算法を含めご覧いただき嬉しい限りです。

      gfortranでコンパイルした際、確かに
      At line 1457 of file main.f90 (unit=10,file=’./input’)
      Fortran rutime error: End of file
      というエラーがでました。現在のところ詳しい理由は分かりません。
      inputの中の変数 “ik” に原因があるようです。とりあえずの処置としては、inputファイルの中のikを消し、プログラム本文にて直接ik=0と書いてください。こうすればエラーが出ません。
      幾日かの間に原因を見つけたいと思います。

      ifortでは問題なく動作したため、gfortranの問題だと考えています。

      1. >>sikinoさん

        お返事いただきありがとうございます。
        ご教授頂いたように、inputファイル中の変数ikを消し、一先ず本文プログラム中のprogram main内における1461行目に「ik=0」と追加をし、gfortranにてコンパイルを行い再実行を行ったのですが、先程投稿した内容と同様のエラーが出てきてしまいました。
        此方あまりfortranには詳しくないのですが、この場合のinputファイルはlinuxにおいて拡張子なしのファイルで「input」と命名する形でよろしいのでしょうか。

        ifortにおいては問題なく実行できるとのことですので、入手ができましたらそちらの方でも試してみたいと思います。

        1. adeluguさん

          inputファイルはおっしゃる通りの意味です。
          あるディレクトリ”dir”の中に
          dir
          ├── input
          └── main.f90
          というファイルが作られることを想定しています。
          「inputファイル」は拡張なしのファイル名”input”というファイルのことです。
          “dir”, “main.f90″の名前は何でも構いませんが、”input”はinputという名前でなければなりません。

          試していただきたいのが以下の2つの手順です。
          ・内容は変わっていませんが、新たなディレクトリを用意し、新たに上記プログラムと”input”をコピペしてください。
          ・inputファイルの最後、&end の後に改行を入れてください。
          —————————
          &input
          ….
          &end
          —————————
          ではなく、
          —————————
          &input
          ….
          &end

          —————————
          としてみてください。恐らくgfortranの仕様です。

          1. >>sikinoさん

            ご教授頂いた内容の通り、&end後に1行改行を入れた結果コードを無事に動かすことができました。
            私の知識不足により貴重なお時間を割いて返信をいただけたこと、お礼申し上げます。
            gfortranに限らず、様々なプログラム、コンパイラを使用する上での仕様を理解しておくことの重要さを今回のことで理解できた気がします・・・。

          2. adeluguさん

            動作したようで安心しました。
            計算をどうぞお楽しみください。

  3. こんにちは。
    ナポレオン時代の砲弾の弾道計算をしたいと思い、色々と調べて
    このホームページを見つけました。プログラムを使わせて頂き
    試算した所、精度の良い結果を得る事ができました。
    大変ありがとうございます。

    なお、私の環境(OS : Windows10 64bit)では、gfortranの古いVer5.1を
    使うとinputのデータを読込む際に(write文で)Segmentation faultエラーが
    発生しました。新しいVer6.2では発生しませんので、コンパイラーの問題
    のようです。Ver5.1以前のものは使わないように注記されては如何でしょうか。

    1. Bernadotte66さん
      こんにちは。
      ご報告ありがとうございます。
      砲弾に対しても良い結果が得られたとのことで、嬉しく思います。

      若干詳細に入るのですが、砲弾…おそらく球形だと思いますが、これは物理的に本当に下回転がかかる仕組みになっているのか、と疑問に思いました。
      そうでなければ、別の要因を考える必要が出てきます。

      gfortranコンパイラに問題があることは感じていました。
      それも指摘して頂いた通り、inputファイルに関連する問題です。
      私はlinux mint, gfortranのver4.8.4を使用しており、
      ”inputファイルの最後に空行を入れる事”
      で何とか回避ができる、ということが分かっています。
      ただし、これはreadの部分であって、writeの部分ではありません。
      もしかしたら、上記問題とは別の問題がwindows版のgfortranコンパイラでは起こっているのかもしれません。
      windows版についてgfortranのバージョンに関して注記いたしました。ご連絡ありがとうございました。

  4. >>sikinoさん

     弾道計算した砲弾は鉄の球体です。また、ナポレオン時代の大砲の砲身は、
    内部に砲弾の回転を掛ける仕組みはありません。ライフルのような溝を
    切っている訳ではなく、単純な筒状のものです。

     従って、揚力が効いている根拠となる物理的なデータは、ありません・・・。
    しかし、次の理由で揚力か、それに相当するものが効いていると思います。
    1)砲弾が地面に跳ね返って進む跳弾の計算をした場合に、地面に初めて
    着弾する距離はパラメータの調整で合わせる事はできる。
    2)しかし、同じ計算式を使って2回目と3回目の着弾距離が実測値(本に
    記載されたもの)と同様な値になるのは、計算式が正しいからである。

    計算例)1つの砲弾が2回地面に跳ね返って3回着弾する距離の計算結果
        ⇒大砲の射出口を基準0mとして水平方向に測った距離。

          1回目  2回目  3回目
    実測値 365m  548m  640m
    計算値 367m  563m  641m

    1. Bernadotte66さん

      ご返信ありがとうございます。
      現在のところ、1),2)に関しましては以下のように考えています。

      1)に関してはまさにその通りで、パラメータは調整するべきだと思います。
      2)に関しては状況を反映していないので賛成はできません。

      状況を把握するためにBernadotte66さんのページを拝見させていただきました。

      余計なおせっかいかもしれませんが、私なりの意見を述べさせていただきたいと思います。
      私の結論としては「初速が異なっている」です。

      初速が異なる、を考えた理由としまして
      着弾距離:Artillery Equipments of the Napoleonic Wars p28
      初速:Manuel d’Artillerie p212の8ポンド砲
      と別の文献から得られていることから、両者で使用している大砲が同一のものであるか不明であった点です。
      前者の本を閲覧できる状況にないので、私が直接確認することが出来ませんでした。
      なので、もし異なる大砲であり、着弾距離を基準にする場合初速や仰角、砲台の高さも若干変わるかもしれません。

      また、砲弾の描像としましては、以下の事が起きていると考えます。
      [発射~1回目の着弾]
      私は発射直後、1回目の着弾までは回転はほとんど起きていない、と考えます。
      もしも発射時の高さが決まっていれば、初速は1回目の着弾位置を回転無しで逆算するのが良いと思います。
      これが余りにも的外れな初速でしたら別の要因を考えなければいけません。

      [1回目の着弾~2回目の着弾]
      1回目の着弾後は地面と接触し、前方向への回転が発生するはずです。
      この場合、1回目と2回目の着弾までの間はトップスピンになるわけですから、早くに着弾すると思います。
      また無視できるかもしれませんが、砲弾の運動エネルギーは回転のエネルギーへいくらか変換されるはずです。
      もちろん、Bernadotte66さんの考察通り反発係数も考える必要が出てきます。

      [2回目の着弾~3回目の着弾]
      その後の着弾では回転は弱まり、3回目の着弾を迎える、というのが正しい描像だと考えています。

      (もしかして
      1364 Pied/second = 416 m/s
      ではないでしょうか。
      1364 Pied/secondはManuel d’Artillerie p212を参考にしましたので、間違いはないと思います。)

  5. こんにちは。
    私のブログまで見て頂き、ご意見ありがとうございます。

    いつかの点について、補足致します。
    1) 着弾距離と初速を記述している大砲が同一かどうか
     個体としては同一ではないと思います。
    但し、フランス軍の8ポンド砲は、グリボーバル・システムと呼ばれる大砲の
    標準化を基に大量に製造されたものです。個体差のバラつきはありますが、
    平均的な性能を検討するので同じで良いと考えました。

    2) Piedの換算
     国際フィートでは0.3048mですので、1364 Pied/second = 416 m/sになります。
    しかし、Manuel d’Artillerieは1836年に出版されたもので、国際フィートが
    制定された1958年よりも前です。下記のホームページに記載されたフランスで
    1799年まで使われていた0.324839mの方が適切と思って採用しました。
    https://fr.wikipedia.org/wiki/Pied_(unit%C3%A9)#Le_pied_entre_1668_et_1799

    3) ご指摘の回転なしで逆算した初速
     850m/sとなり、Manuel d’Artillerieの2倍近くです。
    フランス軍8ポンド砲は、砲弾と火薬が一体化された弾薬を使うので、
    初速がこれほど違うのは考えにくいです。

    4) 個体差のバラつきなどを考慮
     10%以上の差はないと思いますが、以下のようにデータを10%変更した場合の
    第1回着弾距離は275mです。本に記載されている365mと余りにも違うので、
    個体差のバラつきではないと思います。

    半径 0.0515m ⇒ 0.0464m
    初速 443m/s ⇒ 487m/s
    砲身の高さ 1.11m ⇒ 1.22m
    仰角 0° ⇒ 0.1°
    温度 20℃  ⇒ 0℃ (バラつきでなく、条件を良くする為)
    湿度 40% ⇒ 20% (バラつきでなく、条件を良くする為)

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です