命中率と、砲外弾道計算...


基本命中計算...
命中率は軍事機密？なのか正確な物が中々無いようなので．．．
軍事雑誌などに載っている「１０００ｍ先の目標に８割命中する」
などといった簡単なデーターから類推して計算する基本プログラムです。
複雑な要素を無視して(^^)．．．
弾丸の分散は正規分布するのを利用して簡単に命中率を計算します。
入力は射距離と命中率一つだけです。
出力は、１００ｍ刻みで１００００ｍまでです…

基本命中率割り出し計算(テキストボックスとボタン式)...

累積命中率計算...
複数の砲で単独目標を射撃する場合の命中率は．．．
命中率３０％で三門なら
三発とも当たる確率は（３０％）×（３０％）×（３０％）＝（２．７％）
三発とも外れる確率は（７０％）×（７０％）×（７０％）＝（３４．３％）
と簡単に暗算できますが．．．
「タンクテクノロジー」にも似たような表が載ってるのですが．．．
門数が増えたり一門づつの命中率が違ってくると暗算できないので作ってみました．．．

簡単に左端の欄に命中％を書き込めば計算してくれます．．．
戦車小隊で集中射撃をする場合にも便利です．．．
小隊内にエースが居たり、不利な位置で命中率が違う車輛が混じっていても、
命中数期待値を計算できます．．．

累積命中率計算(テキストボックスとボタン式)...

命中率が一定の場合の累積命中率計算を計算できるように改造しました。
2000.0902
使い方、確率論的に１％以下の出来事はなかなか発生しない．．．
50％で７発打ち込むと99.2188%
確実に命中させるには６､７発撃ち込む必要があると計算できる。

現代戦用基本命中計算...
どちらも、測遠器修正していない基本式です
ＷＷ２用ではなくて現代戦用になります。
現代兵器の不気味なぐらいの命中精度が実感できます。
（現代戦車のレーザー測遠器は±５ｍ、 初速1500m/sだと2000mで0.008ミル程度ＦＣＳ固有誤差の５％程度です。）
（誤差ミルを全部合計しても0.43ミル程度です）
一次元式の方は幅と高さで計算して積算します。
二次元式の方は基本式なので標的が円になります
長方形の戦車は積算したほうが良いかもしれません。

正規分布表を使った命中計算一次元式

正規分布表を使った命中計算二次元式

「ＰＡＮＺＥＲ」７８年４月号で当時新型だった７４式戦車の （装備開始から二年目です）特集を組んでありました．．．
７４式の特徴であった油圧式懸架装置の解説の中に射撃の左右ズレの最大の原因 砲耳傾斜（カント）の図面が載っていたのでこれを三角関数で解いてみました．．．
現在ではジャイロ、傾斜計で自動測定自動補正するのでまったくの ブラックボックスです.

砲耳傾斜計算

車体の傾斜は、射角を下げ砲身を傾いた方へ向けます．．．
結果、初速（上向き）＝初速×Ｓｉｎ（砲耳傾斜）、減少させ射距離が減り
横ズレ＝射程×Ｓｉｎ（砲耳傾斜）×Ｔａｎ（射角）、を発生させます．．．

計算すると．．．
初速=1500m/s発射角度7ミル_砲耳傾斜ズレ17.7ミル（１度）飛翔距離＝3237ｍ
横ズレ＝0.3865853618947555m射程距離のショート0.48873587650723493m

初速=1500m/s発射角度7ミル_砲耳傾斜ズレ34.4ミル（２度）飛翔距離＝3236ｍ
横ズレ＝0.7512247345384576m射程距離のショート1.8459257810482086m

初速=1500m/s発射角度7ミル_砲耳傾斜ズレ52.1ミル（３度）飛翔距離＝3233.ｍ
横ズレ＝1.1374762399806537m射程距離のショート4.233693477696306m

補正しないと大体車体が３度も傾けば標的に当たらなくなることになります．．．

初速を入力して修正するタイプに改定したのですが(-_-;)．．．
光学測遠器を使う戦闘は朝鮮戦争ぐらいからなのでＷＷ２用ではありません。
ＷＷ２用に使われるシーンとしては砲兵用の１ｍ弱の基線長の光学測遠器と．．．
対戦車砲グループの連携使用です．．．

ボックス式に直しました。
第二次世界大戦用に使うには、光学式測遠器の基線長（0.9m~2m）
人間の目の幅にすると（0.07m×倍率ぐらい）でOKです．．．
初速が下がってくるとかなり命中率が落ちてくるのが解ります。
現代戦車やパンターの1000m/s以上と火力支援用の四号戦車75mm短砲身や
日本の５７ｍｍ短砲身の命中率の違いが実感できます。
（真空弾道計算なので実際より減りが少ないですが．．．）
2001.4.27

正規分布表（光学測遠器用）を使った命中計算
計算です...
入力ボックスは、誤差ミル、初速m/s、標的車体高さｍ、標的までの距離m、基線長mです。
半角で入力してください．．．
１００ｍ〜３０００ｍまでの命中確率も１００ｍきざみで出力されます。
誤差は、機械的要素だけではなく、射手の技量によっても変わります。

正規分布表（光学測遠器用）を使った命中計算
二次元計算です...

修正した正規分布表（光学測遠器用）を使った命中計算ですが縦方向はこれで良いようですが．．．
偏流を考えていなかったので横方向はまだ修正しないといけないようです

せっかく計算しても一次式の方は表を使用しているので精度がなくなってしまって ます(T_T)．．．
それから縦誤差としては．．．
空気抵抗を計算した場合射角が大きくなるのと
地球が丸いため地平線修正（約4.4kmで1.8m地面が下がるため）
4km近くを狙う現代戦車用はこの修正地も必要になります．．．
後．．．
日本特有の問題ですが基本値の重力加速度が他の先進国より低緯度にあるので小さくなります...

７．
コリオリの力...
コリオリの力ベクトル表示...
（注、暫定版です、砲弾の飛翔データ-から一回だけ計算したものです）
（実際はｚ軸の余りの部分を再計算しないといけないのですが）
（解析していくと元々の方程式に精度が無いので厳密に出しても）
（無駄なような気がしてきてそのままになってます．．．）
（ほんとは、怠けてるだけだったりする）
コリオリの力角度限定...
三次元座標でXYZともプラスになる東〜南の方角のみ計算したもの
（注、式を立てる時に基点からのずれ角度と落下点からのずれ角度の）
（差が微小なので同一として計算しています．．．）
自転速度差

ベクトル計算で自転速度差だけを取り出して計算したもの

落下物体計算

垂直投射計算

定点でのコリオリの力の計算です．．．

簡単に説明すると、上からの落下物体は元々地上より大きなスピードで

回転運動しているため地上との速度差分東にずれる力がプラスされますが

垂直投射の方は丸い地面が回転して進んでしまうため西にずれます．．．


名前の由来
Coliolis、G.Coliolis(1792-1843)フランスの軍人
コリオリの力は、大砲の弾道解析に関連してはじめて見い出された．．．

現象の説明
回転運動差

丸い地球の回転速度は緯度により違います。
発射地点と着弾地点の速度差が弾道の曲がりになります。
着弾地点が同じ緯度になる真東真西へ撃った場合は曲がりません。

遠心力

自転角速度により砲弾は地面に縛られている戦車の砲口を離れると
まっすぐに飛び逆に回転運動する地球に乗っている戦車から見ると
弾道が曲がっているように見えます。
力の成分には、上向きの力がかなり有り．．．
船の場合東に進めば遠心力が増加して軽くなり、西へ進めば遠心力が弱まって
重くなる、この現象はエトヴェス効果と呼ばれています。
（船の上で重力測定をするためには船の速度を正確に測る必要があった）
（エトヴェス＝重力偏差計の発明者の名前）
（重力計は、油田探査に使われて発展してきたためアメリカ製が主流）
エトヴェス補正
船用ですが、最近の戦車は砲安定装置付きで走りながら射撃するので造ってみました．．

力の働く方向を細かく見ると（北半球）．．．

速度の方向、それる方向
北方向＝＝＞　東へ
東方向＝＝＞　南と上へ
南方向＝＝＞　西へ
西方向＝＝＞　北と下へ
上方向＝＝＞　西へ
下方向＝＝＞　東へ

常に右へ曲がって行きます（南半球では左に曲がる）

プログラム
射撃方向を入力するとベクトル分解して答えを出しますが
北半球専用なので南半球は数値（左右前後）のプラスマイナスを
逆にしてください、上にかかる力はそのままです。

ちょっと表示が変ですが...
ｘｙ座標より着弾の射程が増えたり減ったりするのは精密計算を見てもらえば
解りますが、ｚ座標が上下して飛翔時間が増減しているためです。
ｘｙ座標で新しい方位を出し＝＞ｚ座標の上下で増減した数値を縦ずれ（射程ずれ）
に加えているためです...

惑星角速度、赤道における正規重力、惑星赤道半径は、地球以外の
火星や金星などのデータ-を入れるためのものです．．．

惑星角速度は３６０度÷自転時間（秒）÷（１８０／π）
赤道地表面自転速度＝２×半径（ｍ）×π÷自転時間（秒）−（360÷公転周期day）
地球、0.00007292115rad/s、6378.13659km、9.7803278m/s2
火星、0.00007087919rad/s、3397.00000km、3.716524564m/s2
金星、0.00000029924rad/s、6052.00000km、8.900098298m/s2
水星、0.00000123993rad/s、2440.00000km、3.716524564m/s2
月、   0.00000266169rad/s、1738.00000km、1.662655726m/s2

注、地球以外の有効桁数は４桁ぐらいです．．．

方程式は、元々ω2Rcosψなどの微小値を無視して作られているので
（＝3.37m/s2）あまり下まで出しても無意味です．．．
無理して計算するもの無駄で回転系の基本方程式を作る時のこのくらいの誤差が元からある
そうです．．．
（もし、作るなら最初に作りかけたベクトル型で作るのが良さそうです．．．）

月面射撃の場合の注意事項、燃焼のための酸素は自前で準備するとして．．．
ニュートンの第二法則は質量なので、砲弾の初速は変わりません．．．
（発射された後の重力は1/6なので真空計算だとそのまま六倍の射程距離になります。）

地球の自転速度．．．
始め、角速度は力学本によると
ω＝０．００００７２９２rad/sec
ω＝０．００００７３rad/sec
ω＝０．００００７rad/sec
などの簡単な数値で地球の公転周期誤差も引かれいないそうなので
理科年表から正確な数値を取ると
ω＝０．００００７２９２１１５rad/sec

この数値は引かれているかどうか計算すると．．．
1日＝２４時間だと
３６０度÷（２４×６０×６０）÷（１８０／π）＝0.000072722052

1日＝２３時間５６分４．０９０５秒だと
360度÷(23×60×60+56×60+4.0905)÷(180/π)=0.00007292115858

どうやら理科年表はちゃんと修正された数値のようです...

赤道面速度を計算してみると...

(ω×(180/π)/360)×(6378136.59×2π)＝465.1011097m/s

ωに24時間で割った(0.000072722052)を入れてみると...
463.8311808m/s
差は僅かですか1.2699m/sでした。

遠心力
地球の重力を万有引力定数と地球の半径から導くと...
mg0=G(Mm)/R2

g0=G(M)/R2
398600441800000×6378136.592=9.79828673m/s2

g0=赤道の重力（北極のように遠心力の無い場合）
G=万有引力定数、M=地球の全質量、R=地球の半径

遠心力込みの重力は、g≒g0-(ω2R/g0)cos2α
α=緯度
 0=g-0.0344246m/s2
15=g-0.0211860m/s2
30=g-0.02581846m/s2
45=g-0.172123067m/s2
60=g-0.08606153353m/s2
75=g-0.02306011841m/s2
89=g-0.0000010485272m/s2
0=g-0m/s
正規重力計算は
γ＝97803.85(1+0.0053024sin2α-0.0000059sin22α)
δ＝0.87-0.0000956H（大気補正、天体計算用の数値なので大気の底で使えるように大気引力を引く）
λ＝0.3067H（標高補正、高さは3センチ精度で8桁目が机の上と床で7桁目が動く）
α＝月と太陽の引力影響（7桁めが1増えるか減るか程度）
凾＝g0+α-γ-δ+λ

H＝高度（単位ｍ）、α緯度（地理緯度）、

注「γ計算」理科年表に載っている方法はもう少し複雑

高度補正はフリーエアと呼ばれて岩石の無い空中の勾配を現している 岩石が余分に詰まっている分質量が増し重力が増えそうですが、 山も平地も大体同じになります。
これは巨大な岩石を支えるため地下の密度が下がって 釣り合うというアイソシスターという構造のためで、 よって数値はそのまま使えます...

（隕石の落下後やホットプルームなどが有ると重力異常が起こりますが通常はそのまま使っても かまわないでしょう
油田探査や、地殻調査をするわけでは無いのですから．．．）

正規重力と遠心力補正
大体、緯度で5ガル、高さ３０００ｍで１ガル違ってくる．．．(９．８m/s2＝９８０ガル)
この方法だと四桁ぐらいの精度です．．．
月と太陽の引力影響は標高補正１ｍ未満ですこのボックスは地形補正などを入力して 使ってください．．．
海溝、ホットプルーム、隕石落下後の油田地帯などなど．．．

検証用のデータ-...
１．
w=0.0000729(rad/s)、g=9.8m/s3、V=400m/s、緯度=45度、仰角=45度、 投射方向は真南
凾=64.8m、ずれ角度=0.23度（西へ）

２．
北半球では毎秒ωsinα度で右にずれ南半球では左にずれる。 （注、ｔ（飛翔時間が小さい場合のみ）
緯度３５度
１秒＝0.00239575度=0.0425911mil
２秒＝0.004791499度=0.851822mil
３秒＝0.071872494度=0.127773mil

３．
緯度40度、前上に100m/sで投げ上げ
（投げ上げだけど速度なので９０度(1800mil)で．．．）
西へ1.7cmずれる

ps.
方程式だけ作っても現象の正確な理解と正誤ができないので...
結局力学参考書を図書館で乱読して探すことになってしまいました。
書かれていることは、一応力学書３０数冊分の成果です...
（コンピューターで検索したら...）
（ただの力学じゃなくて構造力学、熱力学、から政治力学まで出てきてしまって...）
（５００件の中から選り分ける羽目になりました...）
（結局「当たり」が無かったということですが...）
（注、この愚痴は地球物理で検索しなかったと言う馬鹿なことをしてしまった結果生じた。）


追加
コリオリや弾道計算は要人暗殺やその他の過激な活動に結びつくと世間では考えられているらしく（地域別の重力数値を見るのに本人確認が必要なのを見ても解る様に）あまり正確に計算できるものはご法度のようです。コリオリの計算も一応無害で実用に使えない真空弾道とセットにしておきます。しかし、飛翔距離と飛翔時間を実際の砲弾や弾丸と合わせてしまうと正確な数値が出てしまいますが．．．一応ワンクッションおいてあるということです．．．

砲外弾道計算
砲外弾道各部の名称解説
原点=発射砲の砲口中心（厳密には発射瞬間の砲口中心）
落点=弾道計算し易くするため原点を地表面として計算した着弾点
弾道基線=原点と落点を結ぶ弾道計算用仮想地表面に引かれた水平直線。
射程=原点から落点までの直線距離
最大射程=弾丸の到達できる最大の射程
（真空弾道で砲口での姿勢変化を考えない場合、発射角度４５度が最大射程になる。）

発射射線=実際の砲身軸の延長線、
射線=砲外弾道計算上の砲身軸の延長線（弾丸の運動方向）
発射角度=実際の砲身軸の角度（跳起角＋射角）
射角=砲外弾道計算上の砲身軸角度（弾丸の飛翔角度）
跳起角=砲口で弾丸が軸不一致で飛び出してくるのを補正するための角度
（原因は、砲身内壁と砲弾直径との誤差、射撃時の砲口移動、非対称砲口爆風など）

落角=砲弾の原点着弾時の仮想地面との角度（真空弾道の場合発射角と同じになる）
落速=砲弾の原点着弾時の速度
落線=砲弾の原点着弾時の砲弾姿勢の延長線

在速=弾道上のある点における弾丸の速度

弾道最高点=弾道上で一番高い場所（弾道の垂直方向の運動が重力に引かれて０になった地点）
最大弾道高=弾道最高点の弾道計算用仮想地表面からの高さ

偏流=ライフル弾の場合の回転によるカーブやシュートする水平角度。

目標の地点が計算用仮想地表面と違っている場合（つまり実戦では）．．．
弾着角＝落角
弾着線＝落線
目標面＝弾道基線
目標高＝弾道基線からの高さ
などの用語を使う、一々分けずに混同される場合も多い．．．
Galileo,「ガリレオ」
古代では、弓矢や砲弾は直線に飛んで行ってある地点で方向を換えて落下してくる と信じられていました(-_-;)...
軍事学上の必要からイタリアルネッサンス期に飛翔物体に作用する力が 科学的に計算され...
水平運動（等速運動）と垂直運動（重力による落下）の合成で曲線になると いう計算式が初めて？作られました...
その１、真空弾道計算（１７世紀的に空気抵抗を無視したもの）
単純計算で空気抵抗を計算しない場合は初速と発射角度だけで弾道が出せます...
空気抵抗を計算するには定数を導入しても個々の砲弾の直径と重量が必要になるので まず簡単な真空弾道式から．．．
真空弾道の計算(テキストボックスとボタン式).
入力は初速ｍ／ｓと発射角度（ミル）だけです．．．

式は、（初速度Ｖ、発射角度θ、経過時間ｔ、重力加速度ｇ
速度成分を分割すると
ｘ成分（水平方向）、Ｖｘ＝Ｖ×Cosθ..........真空弾道の場合水平速度は常に一定です
ｙ成分（垂直方向）、Ｖｙ＝Ｖ×Sinθ−ｇ×ｔ..........重力加速度分減らされていきます。

位置は単純に速度成分に経過時間ｔを掛ければ出てきます
(Ｖｘ)×ｔ＝(Ｖ×Cosθ)×ｔ
(Ｖｙ)×ｔ＝(Ｖ×Sinθ−ｇ×ｔ)×ｔ

飛翔距離は
L＝(Ｖｘ)×(Ｖｙ)÷ｇ ＝＞ V2×sin(2×θ)÷ｇ
（速度式Vyを0と置いてｔが0となる場所を計算し速度式Vxに代入します）
（加法定理で見易く整理）

弾道最高点の高さ
ｈ＝(1/2)×Vx2÷ｇ ＝＞ V2×sin2θ÷2ｇ
弾道最高点到達までの時間
ｔa＝Vsinθ÷ｇ
飛翔時間
ｔ＝２×(Vsinθ÷ｇ)
.
射角計算(テキストボックスとボタン式)

真空弾道の場合は計算で射角が出せます．．．
sin2θ＝2×((g)×(x)÷Ｖ2)

旧射角計算(テキストボックスとボタン式)
式を何度も計算させて角度を出す空気抵抗有タイプと同じ方法で解いたもの
-------------------------------------------------------------------
放物線弾道式（ｔ消去）
ｙ＝ｘ*tanθ−（ｇ／（２Ｖ2）*（Cos2θ））*ｘ2
-------------------------------------------------------------------
障害物計算

障害物計算用.
入力は、初速と．．．
(x)障害物までの距離＝L
(y)障害物の高さ＝ｈ
障害物までの距離を仰角で現すと、tanθ＝ｈ÷L
（三つの内二つを入力）
空きを自動計算します。
注、空白を検知して計算するので必ず自動計算させるボックスは
空にしておいてください．．．

木や建物など障害物の上を超える発射角度計算は．．．
解いてみるとxy座標が定まった場合の角度解は二つ有り
計算式は...
障害物の視差を傾斜面として下の（傾斜面用）に代入して答えをまず一つ出します．．．
発射角度をα1α2、障害物の視差をβ

α1-β＝９０−α2

色々やってみたのですが．．．
この式が一番簡単なようです．．．

真空弾道の計算(傾斜面用）.
入力は初速ｍ／ｓと発射角度（ミル）、地面の傾斜角度、の三つです．．．

傾斜面で射撃した場合の最大射程を計算できるボタンを付けました...
式は、（初速度Ｖ、発射角度α、傾斜角度β、経過時間ｔ、重力加速度ｇ
方程式は...
射手自身も傾いているため座標を傾斜面上に設定し重力成分の方を傾けます。
ｘ軸方向（sinβ）、ｙ軸方向（cosβ）
速度成分を分割すると

ｘ成分（水平方向）、Ｖｘ＝Ｖ×Cosα−g×t×sinβ
ｙ成分（垂直方向）、Ｖｙ＝Ｖ×Sinα−ｇ×ｔ×cosβ
（下り傾斜、傾斜がマイナスの場合）
（ｘ成分（水平方向）、Ｖｘ＝Ｖ×Cosα＋g×t×sinβ）

位置は速度式を積分して．．．
（水平位置）X＝(Ｖ×Cosα×ｔ)−(1/2)×ｇ×ｔ2×sinβ
（垂直位置）Y＝(Ｖ×Sinα×ｔ)−(1/2)×ｇ×ｔ2×cosβ
（下り傾斜、傾斜がマイナスの場合）
（（垂直位置）Y＝(Ｖ×Sinα×ｔ)＋(1/2)×ｇ×ｔ2×cosβ）

このままだと傾いているので．．．
斜面上、 x＝x-ytanβ、ｙ＝y/cosβ)
弾道基線上、x＝(x-ytanβ)cosβ、ｙ＝y/cosβ)

飛翔距離（斜面上）は...
L＝V2.×(sin(2α+β)-sinβ)÷(g×cos2β)
（速度式Vyを0と置いてｔが0となる場所を計算し速度式Vxに代入します）
弾道基線上は、Ｌ×cosβ

弾道最高点到達までの時間（この式は傾斜無しで計算）
ｔa＝Vsin(α+β)÷ｇ

弾道最高点の高さ
(弾道基線上仮想ｈ＝V2×sin2(α+β)÷2ｇ)
斜面上ｈ＝到達時間ｔaからｘｙを出し...
h-(ｘ-y×tanβ)×sinβ

飛翔時間は速度Y式のy速度を0と置いて．．．
t＝(2×V×Sinα)÷(ｇ×cosβ)
.
最大射程は．．．
仰角α＝(90-β)÷2、下り傾斜の場合（仰角α＝(90+β)÷2、）
最大射程Ｌ＝(V2)÷(g+g×sinβ)
弾道基線上の射程距離Ｌ×cosβ

弾道計算（比例抵抗）2001.5.27...
単位質量当たりの比例係数ｋの値は空白にしておかないと自動計算せずに 書かれてある値を使用します．．．
レイノルズ数が小さい場合用ですＦ∝ηｖ （微生物サイズ用？）
　
レイノルズ数Ｒ＝ρｖａ／η
ρ＝流体の密度
ｖ＝流体の速度
ａ＝物体の寸法（断面積）
η＝流体の粘性
（流体力学で、空気抵抗を考えると速度や形状によって大きく違ってきて）
（一般に使う速度の二乗に比例するというのとは違ってきます．．．）
（人間サイズで考えた場合は乱流、渦などが発生し二乗でOKですが）
（小さな物体の場合は流れが乱れないので比例抵抗だそうです。）

教科書によく入門用に載っている比例抵抗式です．．．
しかし、この段階ですでに飛翔距離や飛行時間を計算しようとすると
展開式になってしまいます。

飛翔距離
(2*Vy/g)×{1-(1/3)*k*vy/g+(2/9)*(k*vy/g)2.........}

飛翔時間
(2*Vy*Vx/g)×{1-(4/3)*k*vy/g+(14/9)*(k*vy/g)2.........}

Vx＝水平初速
Vy＝垂直初速
ｇ＝重力加速
k＝単位質量当たりの比例抵抗係数
*＝コンピューター用掛け算記号

基本方程式は、
運動方程式
ｍｘ＝−ｍｋｘ、　ｍｙ＝−ｍｇ−ｍｋｙ
加速度
ｘ＝−ｋｘ、　ｙ＝−ｇ−ｋｙ
速度
ｘ＝Ｖx0ｅ-kt、　ｙ＝−ｇ／ｋ＋（Ｖy0＋ｇ／ｋ）ｅ-kt）

位置
ｘ＝−（Ｖx0／ｋ）×（１−ｅ-kt）

　ｙ＝−（ｇ／ｋ）ｔ＋（１／ｋ）×（Ｖx0＋ｇ／ｋ）（１−ｅ-kt））

Ｖx0＝水平初速
Ｖx0＝垂直初速
ｔ＝飛翔時間

多分今まで作った中で一番役に立たないもの．．．
方程式の説明用です。
比較用
弾道計算（比例抵抗､ルンゲ・クッタ法）2001.5.30...
上の比例式をルンゲ・クッタ法で解いたものです。
精度は、
速度は、小数点以下６桁．．．
距離は、小数点以下１桁．．．
でした．．．
計算間隔を0.1秒から0.01秒にすると小数点以下四桁まで合いますが
計算時間がかなりかかります...
そこで、
以前図書館でプログラムの本をめくっていた時に見たステップを減らす方法を使おうと
探したのですが、今回発見できず...（ちゃんとコピー取っておけば良かった）

ルンゲ・クッタ法は、運動方程式を作って放り込むだけで答えが出るので簡単です。

速度
Ｖx＝Ｖx−ｋＶxｔ......Ｖy＝Ｖy−ｇｔ−ｋＶyｔ
距離
Ｘ＝Ｖxｔ−(1/2)ｋＶxｔ2..... Ｙ＝Ｖy−(1/2)ｇｔ2−(1/2)ｋＶyｔ2

Vx＝水平初速...Vy＝垂直初速...ｇ＝重力加速...t＝分割時間
x＝水平距離...y＝垂直距離...k＝単位質量当たりの比例抵抗係数

ルンゲ・クッタ法
ｙn+1＝ｙn＋（ｋ1＋２ｋ2＋２ｋ3＋ｋ4)／６

ｋ1＝ｈ×∫(ｘn,ｙn)
ｋ2＝ｈ×∫(ｘn＋ｈ／２,ｙn＋ｋ1／２)
ｋ3＝ｈ×∫(ｘn＋ｈ／２,ｙn＋ｋ2／２)
ｋ3＝ｈ×∫(ｘn＋ｈ,ｙn＋ｋ3)

オイラー法とか修正オイラー法とかもボタンを作って比べられるようにしようか...
と思ったのですが精度がルンゲ・クッタ法以下なのでやっぱ辞めました。
2001.5.30...
空気抵抗の二乗に比例する弾道計算
弾道の計算（空気抵抗定数計算）

.
(du/dt)=-Cd*(q/2)*A*V*Vcosθ
.
(du/dt)=-Cd*(q/2)*A*V*Vsinθ-W
.
W=弾丸重量(Kg)
Cd=空気抵抗係数
q=空気密度(1.226kg/m3)
A=弾丸断面積(m3)
d=弾丸直径(m)
g=重力加速度(9.8m/s2)
(*=コンピューター記号、四則演算の「×」と同じ)
.
どうやら現代兵器はKdでは無くてCd表示が多いようです。
発射角度は、高速弾5mil〜10mi、低速弾20mil〜30millぐらいで試してみてください。
105mm戦車砲弾の弾道特性
砲弾種類重量kg直径mm初速m/sマッハ空気抵抗係数Cd距離当たりの速度低下/km
APCBC17.410510203.00.2887m/s
APDS4.016114754.340.1599m/s
APFSDS3.572615004.410.3548m/s
HEAT10.1710511753.460.46282m/s
HEAT5.1510511203.290.22231m/s
HESH11.261057302.150.64220m/s
HE15.51056501.910.3169m/s
HE+FS7.21058002.350.27159m/s
Panzer95/12_P59「戦車砲と戦車砲弾の発達」より

データ-を入れて「距離当たりの速度低化/km」（つまり１０００ｍごとにこれだけ速度が低化します）
を比べてみてくださいぴったり合います．．．

弾道の計算（二乗抵抗､傾斜面用）
方程式は...
射手自身も傾いているため座標を傾斜面上に設定し重力成分の方を傾けます。
（初速度Ｖ、発射角度α、傾斜角度β、経過時間ｔ、重力加速度ｇ）
ｘ軸方向（sinβ）、ｙ軸方向（cosβ）
速度成分を分割すると

ｘ成分（水平方向）、Ｖｘ＝Ｖ×Cosα−g×t×sinβ-kv2cosα×t
ｙ成分（垂直方向）、Ｖｙ＝Ｖ×Sinα−ｇ×ｔ×cosβ-kv2sinα×t
（下り傾斜、傾斜がマイナスの場合）
（ｘ成分（水平方向）、Ｖｘ＝Ｖ×Cosα＋g×t×sinβ-kv2cosα×t）
表示する時は方位を直して、
登り坂、Ｖx＝Ｖx×cosβ、Ｖy＝Ｖy×cosβ
下り坂、Ｖx＝Ｖx÷cosβ、Ｖy＝Ｖy÷cosβ

位置は速度式を積分して．．．
（水平位置）X＝(Ｖ×Cosα×ｔ)−(1/2)×ｇ×ｔ2×sinβ -(1/2)×kv2cosα×t2
（垂直位置）Y＝(Ｖ×Sinα×ｔ)−(1/2)×ｇ×ｔ2×cosβ -(1/2)×kv2sinα×t2

（下り傾斜、傾斜がマイナスの場合）
（（垂直位置）Y＝(Ｖ×Sinα×ｔ)＋(1/2)×ｇ×ｔ2×cosβ） -(1/2)×kv2sinα×t2


このままだと傾いているので．．．
斜面上、 x＝x-ytanβ、ｙ＝y/cosβ)
弾道基線上、x＝(x-ytanβ)cosβ、ｙ＝y/cosβ)

弾道の計算（空気抵抗定数計算）を利用して...
jacob de Marre
1/2*P*Vr2=c1*(d1.5)*S1.4         S∝((p/d1.5)0.714)*Vr1.43
Krupp
1/2*P*Vr2=c1*(d(5/3))*S1.3     S∝((p/d1.67)0.75)*Vr1.5
Moisson
1/2*P*Vr2=(c1/Vr)*d*S2            S∝((p/d)0.5)*Vr1.5
などの式を代入して計算させています...
この方程式達の基本原理は、

弾丸が装甲板を貫通した穴の体積・・・
（π／４）×穴（砲弾）の直径の二乗×貫通穴の深さ

砲弾の運動エネルギー
（１／２×重力定数）×砲弾重量×砲弾速度の二乗

この二つの間には相関関係があり比例式が成り立つはずなのですが
実験しても正確に比例式にならないので色々修正を加えた式が
各国で作られて使われています...
（注意、どれも古典式です現代の弓矢型砲弾用は別の式が存在します）

注意、初速は．．．
砲弾の速度は砲口爆風のため砲口を離れてもまだ加速を続けます。
そのため、初速は砲口爆風の無くなった場所の速度から砲外弾道式を使って
逆算して計算さます。
初速とは、弾道計算用にもとめられるもので実際の砲口の速度はもっと低い．．．
（つまり０ｍでの貫通力表示は正しく無い．．．）
＃
＃
＃


貫通力計算（空気抵抗グラフ入力計算）
空気抵抗はマッハ１の手前で急速に大きくなり
0.14MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
0.15MMMMMMMMMMMMMMSSSMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
0.14MMMMMMMMMMMMMSMMMSSMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
0.13MMMMMMMMMMMMSMMMMMMSSMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
0.12MMMMMMMMMMMMSMMMMMMMMSSMMMMMMMMMMMMMMMMMMMMMMMMMMM
0.11MMMMMMMMMMMSMMMMMMMMMMMSSSMMMMMMMMMMMMMMMMMMMMMMMM
0.10MMMMMMMMMMMSMMMMMMMMMMMMMMMSSSMMMMMMMMMMMMMMMMMMMM
0.09MMMMMMMMMMSMMMMMMMMMMMMMMMMMMMSSSSMMMMMMMMMMMMMMMM
0.08MMMMMSSSSSMMMMMMMMMMMMMMMMMMMMMMMMSSSSSMMMMMMMMMMM
0.07SSSSSMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMSSSSSSMMMMM
//////////////////////////////////////////////////////
0.02MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
0.01MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
0.00_______0.8__1_______2_______3___________4_________
２〜３と増えるにしたがって小さくなっていきます。
「何故マッハ近くで抵抗が増すかというとマッハ0.8~9辺りから砲弾側面の速度が
音速を超え（圧縮されるため）不規則な垂直の衝撃波が発生し砲弾周囲に乱気流を
発生させるためです。
マッハ1.9辺りまで来ると砲弾の舳先から発生する衝撃波がコーン状になって気流が
安定しだします...」

砲弾の形状によっても抵抗係数は違ってきます。
データーは弾丸の「長さＬ／直径Ｄ」＝3.5と5の二つです。
（ほんとは砲弾の形ごとにグラフを書かなければいけないのですが_(._.)_...）

あまり精度が無いのでデーターは折れ線グラフのようなデーターに
なってしまうと思います。（マッハ４までを0.1刻みで４０データー）


貫通力計算､空気抵抗グラフ入力計算
（全部まとめて一つにしました）

注、普通コンピューターで簡単に弾道計算ができますといってプログラムが載っているのは...
単純に分割時間ごとに空気抵抗による減速を引いただけのものが多いようです。
このプログラムも同じような簡易版で、解答は、分割時間点で折れ曲がる
折れ線グラフ（多角形）になります...
つまり分割時間を小さくすればするほど弾道軌道に近づいていくことになりますが。
実用的じゃないので、サンプルは０．１秒で計算してあります。

使用する場合は、最低０．１秒以下の分割時間で使ってください...
（分割時間を変更すれば当然計算結果はずれます...）
それから、このプログラムは、戦車砲用に作ったので、
高度による大気密度の補正をしていません...
（そのため、最大仰角の射程が実際より短くなります）
（パンタークラスで18~20度の仰角をかけると高度1000mを簡単に越える）
（地表1062.2mb、1000mで898.7mb、気温-9.1度）
2000.10.28...


貫通力用の変数は、「実験定数n」と実験定数２の二つです．．．
「実験定数（貫通力計算用）基本「1.00」」は微調整用に使ってください。


空気抵抗には色々な種類があってこのグラフも雑多な抵抗の合力です...

造波抗力
弾丸頭部、安定羽部分で発生する垂直空気圧力．．．
＃
＃
（ボートテイル型砲弾の場合は弾丸後部でも発生する）

摩擦抗力
弾丸表面と空気の摩擦による抗力（回転による摩擦もある）
翼安定弾と回転安定弾では、係数が違い
回転体の場合は１５％増しになる。

弾底抗力
弾底に発生する真空状態の吸引力による抗力

＃
＃
＃

ジャイロ効果と空気抵抗．．．
ＷＷ２の弾丸は現代のＡＰＦＳＤＳ（有羽弾）と違い回転により安定させています。
しかし、弾道は曲線を描きジャイロ安定しようとする砲弾とは違った姿勢を
強制します。空気の流れが砲弾の向きを変えます。
つまり、弾道に対して斜めに砲弾が飛んでいくことになり
そのぶん空気抵抗が増します。
＃
＃
＃


つい貫通力計算式に走ってしまいましたが(^^)...
その３、砲外弾道積分計算式

世界初の計算機ＥＮＩＡＣで計算された積分方式と思ったのですが．．．
右辺にｙを含むつまり答の増減に連動して値が変わる式は解くのが難しく
速度Ｖの二乗に連動している場合式は解けないようです(-_-;)．．．
つまりコンピューターで値を順番に入力していって答えの近似値を出すという
作業が必要で積分式を解くではなく無理やり計算式を回して答えの近似値を
吐かせる．．．
ということらしいです(^^)．．．

二乗でない場合は

Ｖｘ＝Ｖcosθ*ｅ(Ｃd*t)

Ｖｙ＝（Ｖsinθ+（g/Ｃd））*（ｅ(Ｃd*t)）−(g/Ｃd)

となります。
＃
＃
＃
基本的な解法は．．．
オイラー法
ｙn+1＝ｙn＋ｈ×∫(ｘn,ｙn)

修正オイラー法
ｙn+1＝ｙn＋（ｋ1＋ｋ2)／２

ｋ1＝ｈ×∫(ｘn,ｙn)
ｋ2＝ｈ×∫(ｘn+1,ｙn＋ｋ1)

ルンゲ・クッタ法
ｙn+1＝ｙn＋（ｋ1＋２ｋ2＋２ｋ3＋ｋ4)／６

ｋ1＝ｈ×∫(ｘn,ｙn)
ｋ2＝ｈ×∫(ｘn＋ｈ／２,ｙn＋ｋ1／２)
ｋ3＝ｈ×∫(ｘn＋ｈ／２,ｙn＋ｋ2／２)
ｋ3＝ｈ×∫(ｘn＋ｈ,ｙn＋ｋ3)

などなのです...


戦史研究＆戦車戦へ戻る