PHPで弾道ミサイルの軌道を計算する

(1/1)
北朝鮮が、相次いで弾道ミサイルを発射している。弾道や着弾地点に関して様々な報道がなされているが、結果の数字ばかりで、それが正しいのか検証が難しい。
そこで今回は、マップ上で発射地点と着弾地点を指定し、弾道の最高点の高さを入力するだけで、軌道要素と弾道を計算する PHP プログラムをつくってみることにする。
ユーザーインターフェースは、「PHP で大圏航路を描く」を流用した。

なお、作者の力不足から、コリオリの力、多段ロケットによる加速、弾頭の形状・材質、地球大気による減速など、物理的な諸条件は含めていない。大まかな検証をすることが目的だと考えていただきたい。

(2021 年 3 月 28 日)PHP8 対応,リファラ・チェック改良。
(2020 年 3 月 29 日)地理院地図・ OSM に対応した。OSM Nominatim Search API を利用できるようにした。

目次

サンプル・プログラムの実行例

弾道ミサイルの軌道計算
弾道高を入力し、発射点と着弾点び住所を入力するか、地図上の地点をクリックする。計算ボタンをクリックすると、2 地点を結ぶ射程(赤い曲線)が描かれ、ミサイルの射角、飛翔経路長、平均速度、飛翔時間などが表示される。

YOLP コンテンツジオコーダ API(Yahoo!JAPAN)は、世界のキーワードから緯度・経度を求める機能が弱い。たとえば「サンティアゴ」でキーワード検索すると、Google Geocoding API チリの首都を第一候補として返すが、YOLP コンテンツジオコーダ APIはドミニカ共和国の県都しか返さない。
そこで、地図上のマーカーをドラッグし、計算ボタンを押すことで射程を再描画する機能を採り入れた。これは Google マップ でも利用できる。

サンプル・プログラム

圧縮ファイルの内容
ballistic.phpサンプル・プログラム本体
pahooGeoCode.php住所・緯度・経度に関わるクラス pahooGeoCode。
使い方は「PHPで住所・ランドマークから最寄り駅を求める」「PHPで住所・ランドマークから緯度・経度を求める」などを参照。include_path が通ったディレクトリに配置すること。

弾道方程式

今回は、報道が流している数字が正しいのか検証することが目的であるので、まず、弾道方程式について、数式を含めた説明をしておく。各々の方程式の証明は、楕円について知っていれば、高校の数学レベルで解法可能である。ただし、動径上の速度については、物理学の知識が求められる。
末尾に、参考サイトや参考書籍を掲げた。ぜひご覧いただきたい。
弾道方程式
弾道ミサイルは人工衛星と異なり地球を周回しないが、地球の中心を焦点の一方とする楕円軌道を描く。これを弾道(ballistic)と呼ぶ。

発射点を a、着弾点を b とすると、地表面に沿った円弧 acb を射程(range)と呼ぶ。acb の左側は地中に潜っていると考えてほしい。
射程を、地球中心から見た角度φを半射程角(semi-range angle)と呼ぶ。射程を S、地球の半径を R(6,378,136m)とすると、φは次の式によって求められる。
 mimetex 


ミサイルの最終段のロケットモーターが燃え切ったときの速度比を voとすると、次の式によって求められる。
 mimetex 

 mimetex 

 mimetex 

 mimetex 


楕円の離心率e、長径 a、短径 b は、次の式によって求められる。
 mimetex 

 mimetex 

 mimetex 


地上から測った軌道の最高点の高さ(軌道高)hm は、次の式によって求められる。
 mimetex 


軌道上の速度は変化する。地球中心に対する軌道上の角度をθ、地上の引力係数をμとすると、θにおける動径 r と、動径上の速度 v は、次の式によって求められる。
 mimetex 

 mimetex 

 mimetex 

準備:pahooGeoCode クラス

0037: class pahooGeoCode {
0038:     var $items;      //検索結果格納用
0039:     var $error;      //エラーフラグ
0040:     var $hits;       //検索ヒット件数
0041:     var $webapi; //直前に呼び出したWebAPI URL
0042: 
0043:     //Google Cloud Platform APIキー
0044:     //https://cloud.google.com/maps-platform/
0045:     //※Google Maps APIを利用しないのなら登録不要
0046:     var $GOOGLE_API_KEY_1 = '**************************';   //HTTPリファラ用
0047:     var $GOOGLE_API_KEY_2 = '**************************';   //IP制限用
0048: 
0049:     //Yahoo! JAPAN Webサービス アプリケーションID
0050:     //https://e.developer.yahoo.co.jp/register
0051:     //※Yahoo! JAPAN Webサービスを利用しないのなら登録不要
0052:     var $YAHOO_APPLICATION_ID = '*****************************';

地図描画や住所検索を行うために、クラスファイル "pahooGeoCode.php" を使用する。組み込み関数  require_once  を使って読めるディレクトリに配置する。ディレクトリは、設定ファイル php.ini に記述されているオプション設定 include_path に設定しておく。
クラスについては「PHP でクラスを使ってテキストの読みやすさを調べる」を参照されたい。

地図や住所検索として Google を利用するのであれば、Google Cloud Platform API キー が必要で、その入手方法は「Google Cloud Platform - WebAPI の登録方法」を、Yahoo!JAPAN を利用するのであれば、Yahoo! JAPAN Web サービス アプリケーション IDが必要で、その入手方法は「Yahoo!JAPAN デベロッパーネットワーク - WebAPI の登録方法」を、それぞれ参照されたい。

準備:地図サービスの選択

0034: //地図描画サービスの選択
0035: //    0:Google
0036: //    2:地理院地図・OSM
0037: define('MAPSERVICE', 0);
0038: 
0039: //住所検索サービスの選択
0040: //    0:Google
0041: //    1:Yahoo!JAPAN
0042: //   11:HeartRails Geo API
0043: //   12:OSM Nominatim Search API
0044: define('GEOSERVICE', 12);
0045: 
0046: //逆ジオコーディングサービスの選択
0047: //    0:Google
0048: //    1:Yahoo!JAPAN
0049: //   11:HeartRails Geo API
0050: //   21:簡易ジオコーディングサービス
0051: define('REVGEOSERVICE', 1);

表示する地図は、Google マップ、[Yahoo!マップ]blue]、地理院地図・オープンストリートマップ(OSM)から選べる。あらかじめ、定数 MAPSERVIC に値を設定すること。
住所検索サービスは、GoogleYahoo!JAPANHeartRails Geo APIOSM Nominatim Search API から選べる。あらかじめ、定数 GEOSERVICE に値を設定すること。
逆ジオコーディングサービスは、GoogleYahoo!JAPANHeartRails Geo API簡易ジオコーディングサービスから選べる。あらかじめ、定数 REVGEOSERVICE に値を設定すること。

解説:初期値

0053: //マップの表示サイズ(単位:ピクセル)
0054: define('MAP_WIDTH',  600);
0055: define('MAP_HEIGHT', 400);
0056: //マップID
0057: define('MAPID', 'map_id');
0058: 
0059: //初期値
0060: define('DEF_LONGITUDE0', 165.0);       //中心座標(経度)
0061: define('DEF_LATITUDE0',  30.5);       //         (緯度)
0062: define('DEF_TYPE',     'roadmap');        //マップタイプ
0063: define('DEF_ZOOM',     4);              //ズーム
0064: define('DEF_CATEGORY', 'world');        //カテゴリ
0065: 
0066: //発射点
0067: define('DEF_FROM',     '東倉里ミサイル発射場');        //地名
0068: define('DEF_CAT_FROM', 'world');                        //カテゴリ
0069: define('DEF_LATITUDE1',  39.660075);                  //緯度
0070: define('DEF_LONGITUDE1', 124.705294);                  //経度
0071: 
0072: //着弾点
0073: define('DEF_TO',     'ホノルル');                    //地名
0074: define('DEF_CAT_TO', 'world');                        //カテゴリ
0075: define('DEF_LATITUDE2',  21.308889);                  //経度
0076: define('DEF_LONGITUDE2', -157.826111);                  //経度
0077: 
0078: define('DEF_HM', 1500);               //弾道高(km

軌道計算プログラムの本体 "ballistic.php" の冒頭に、初期値を定義している。
初回起動時に表示する発射点は、北朝鮮の東倉里 (トンチャンリ) ミサイル発射場、着弾地点はハワイ・ホノルルとしてある。弾道高は 1500km である。

解説:マーカー設定など

ユーザー関数 jsDragMarker_YOLPmapjsDrawLine_YOLPmapjsDragMarker_GmapjsDrawLine_Gmap は、マーカー設定などを行うための KavaScript で、「PHP で大圏航路を描く」で作成したものと同じである。

解説:射程,弾道高から軌道要素を求める

0502: /**
0503:  * 射程,射角から軌道要素を求める
0504:  * @param double $dt   射程(メートル)
0505:  * @param double $aof  射角(度)
0506:  * @param array  $elements  軌道要素を格納する配列
0507:  * @return なし
0508: */
0509: function getOrbitalElementSub($dt$aof, &$elements) {
0510:     static $r = 6378136;     //地球半径(メートル)
0511: 
0512:     //ラジアン変換
0513:     $aof = deg2rad(90 - $aof);
0514: 
0515:     $sra = $dt / (2 * $r);       //半射程角
0516: 
0517:     $ax = pow(sin($aof), 2) * (pow(cos($sra), 2) - pow(sin($aof), 2));
0518:     $bx = pow(sin($aof), 2) * (1 - pow(cos($sra), 2));
0519:     $cx = pow(cos($sra), 2) - 1;
0520: 
0521:     $vo = sqrt((-$bx + sqrt(pow($bx, 2) - $ax * $cx)) / $ax);   //燃え切り速度比
0522:     $e = sqrt(1 - pow($vo, 2) * (2 - pow($vo, 2)) * pow(sin($aof), 2));    //離心率
0523:     $a = $r / (2 - pow($vo, 2));      //長径
0524:     $b = $a * sqrt(1 - pow($e, 2));       //短径
0525:     $hmm = $a * (1 + $e) - $r;         //弾道高
0526: 
0527:     $elements['sra'] = $sra;
0528:     $elements['vo']  = $vo;
0529:     $elements['e']   = $e;
0530:     $elements['a']   = $a;
0531:     $elements['b']   = $b;
0532:     $elements['hmm'] = $hmm;
0533: }

0535: /**
0536:  * 射程,弾道高から軌道要素を求める
0537:  * @param double $dt  射程(メートル)
0538:  * @param double $hmm 弾道高(メートル)
0539:  * @param array  $elements  軌道要素を格納する配列
0540:  * @return なし
0541: */
0542: function getOrbitalElement($dt$hmm, &$elements) {
0543:     $flag = TRUE;
0544:     $aof = 0.0;       //射角(度)
0545: 
0546:     //誤差0.1度で射角を求める
0547:     while ($flag && ($aof < 90.0)) {
0548:         $aof += 1.0;
0549:         getOrbitalElementSub($dt$aof$elements);
0550:         if ($elements['hmm'] == $hmm) {
0551:             $flag = FALSE;
0552:             break;
0553:         } else if ($elements['hmm'] > $hmm) {
0554:             $aof -= 1.0;
0555:             while ($flag && ($aof < 90.0)) {
0556:                 $aof += 0.1;
0557:                 getOrbitalElementSub($dt$aof$elements);
0558:                 if ($elements['hmm'] >= $hmm) {
0559:                     $flag = FALSE;
0560:                     break;
0561:                 }
0562:             }
0563:         }
0564:     }
0565:     $elements['dt'] = $dt;
0566:     $elements['aof'] = $aof;
0567: }

北朝鮮の弾道ミサイルは、迎撃を回避するため超高高度に達するロフテッド軌道をとるという報道もある。
そこで本プログラムでは、射角の代わりに弾道高を入力して、軌道要素を求める関数 getOrbitalElement を用意した。

弾道高から射角を逆算することは難しいので、射角を 0.1 度ずつ増やし、ユーザー関数 getOrbitalElementSub を呼び出し、与えられた弾道高にほぼ近づいたら計算を打ち切るという方法で軌道要素を求めている。

計算回数を減らすため、初段は 1 度ずつ射角を増やし、近づいてきたら 0.1 度ずつ増やすようにしている。

解説:軌道要素から、飛翔経路長、速度を求める

0470: /**
0471:  * 軌道要素から、飛翔経路長、速度を求める
0472:  * @param array  $elements  軌道要素
0473:  * @return なし
0474: */
0475: function getBallistic(&$elements) {
0476:     static $u = 3.98602e+14;       //地心引力係数
0477:     $n = 100;                        //分割数
0478: 
0479:     $dw = $elements['sra'] / $n;
0480:     $ds = 0.0;                        //飛翔経路長
0481:     $vs = 0.0;
0482:     $vmax = 0.0;                  //最大速度
0483:     for ($i = 0; $i <= 100; $i++) {
0484:         $r = $elements['a'] * (1.0 - pow($elements['e'], 2)) / (1.0 - $elements['e'] * cos($dw * $i));
0485:         $v = sqrt($u * ((2.0 / $r) - (1.0 / $elements['a'])));
0486:         if ($v > $vmax$vmax = $v;
0487:         $vs += $v;
0488:         if ($i == 0)    $tmp1 = $r;
0489:         else            $tmp2 = $r;
0490:         if ($i > 1) {
0491:             $dx = $tmp2 * sin($dw);
0492:             $dy = $tmp1 - $tmp2 * cos($dw);
0493:             $ds += sqrt(pow($dx, 2) + pow($dy, 2));
0494:             $tmp1 = $tmp2;
0495:         }
0496:     }
0497:     $elements['ds'] = $ds * 2;
0498:     $elements['vm'] = $vs / ($n + 1);
0499:     $elements['vmax'] = $vmax;
0500: }

ミサイルの飛翔経路長は、楕円の弧の一部なのだが、これを計算するには無限級数を扱わなければならず、計算時間がかかる。
そこで、「弾道方程式」で紹介した動径を使って、動径角を 100 分割して、それを足していくことで飛翔経路長を近似している。
同時に速度を 100 回求め、平均速度と最高速度を計算している。

解説:軌道要素の計算

0865:     //軌道要素の計算
0866:     if (($hm > 0) && ($items['longitude1'] != '') && ($items['latitude1'] != '') && ($items['longitude2'] != '') && ($items['latitude2'] != '')) {
0867:         $dt = $pgc1->distance($items['longitude1'], $items['latitude1'], $items['longitude2'], $items['latitude2']);
0868:         $items['hmm'] = $hm * 1000;
0869:         getOrbitalElement($dt$hm * 1000, $items);
0870:         getBallistic($items);
0871:     }
0872: 
0873:     //地図描画サービス
0874:     switch (MAPSERVICE) {
0875:         //Google
0876:         case 0:
0877:             $call  = jsDragMarker_Gmap($items2[$i]['latitude'], $items2[$i]['longitude']);
0878:             $call .= jsDrawLine_Gmap(1, $items1, 'rgb(255, 0, 0)', 1);
0879:             $call .= jsDrawLine_Gmap(2, $items2, 'rgb(0, 0, 0)', 1);
0880:             $zd = 1;
0881:             break;
0882:         //Yahoo!JAPAN
0883:         case 1:
0884:             $call  = jsDragMarker_YOLPmap($items2[$i]['latitude'], $items2[$i]['longitude']);
0885:             $call .= jsDrawLine_YOLPmap(1, $items1, 'FF0000', 1);
0886:             $call .= jsDrawLine_YOLPmap(2, $items2, '000000', 1);
0887:             $zd = 0;
0888:             break;
0889:         //地理院地図・OSM
0890:         case 2:
0891:             $call  = jsDragMarker_Leaflet($items2[$i]['latitude'], $items2[$i]['longitude']);
0892:             $call .= jsDrawLine_Leaflet(1, $items1, 'FF0000', 1);
0893:             $call .= jsDrawLine_Leaflet(2, $items2, '000000', 1);
0894:             $zd = 0;
0895:             break;
0896:         default:
0897:             $call = '';
0898:             $zd = 0;
0899:             break;
0900:     }
0901:     $lat = $items2[$i]['latitude'];
0902:     $lng = $items2[$i]['longitude'];
0903: 
0904: //エラー発生時
0905: else {

マップから得た 2 地点の緯度・経度を射程に変換するには、メソッド distance を使った。このメソッドは、「PHP で 2 地点間の直線距離を求める」で紹介している。

結語

東倉里から東京を狙う場合、弾道高100km なら、発射から 5 分半で着弾する。
一方、弾道高1000km のロフテッド軌道だと、着弾まで 16 分半かかる。J アラートが出てから避難するには十分な時間だ。

さて、東倉里からホノルルを狙う弾道高1500km の軌道では、ミサイルの平均速度はマッハ 16、最高速度はマッハ 20 に達する。これは最終段のロケットモーターが燃え尽きた後の速度であるから、打ち上げ直後の速度はもっと遅いし、着弾時には空気抵抗の影響で減速するだろう。
ただ、着弾前に攻撃するミサイル迎撃システム「PAC3」が演習で撃ち落とした飛翔体はマッハ 4 である。

参考サイト

参考書籍

表紙 惑星探査機の軌道計算入門
著者 半揚稔雄
出版社 日本評論社
サイズ 単行本
発売日 2017年09月19日頃
価格 2,420円(税込)
rakuten
ISBN 9784535788459
数値計算をしながら宇宙飛行の楽しさを味わう!人工衛星や惑星探査機における軌道計算と軌道決定のカラクリを、高校数学、物理の知識をもとに分かりやすく紹介!
 
表紙 天文計算入門
著者 長谷川一郎
出版社 恒星社厚生閣
サイズ 単行本
発売日 1996年01月25日頃
価格 2,750円(税込)
rakuten
ISBN 9784769908180
 
(この項おわり)
header