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 クラス

  37: class pahooGeoCode {
  38:     var $items;     //検索結果格納用
  39:     var $error;     //エラー・フラグ
  40:     var $errmsg;    //エラー・メッセージ
  41:     var $hits;      //検索ヒット件数
  42:     var $webapi;    //直前に呼び出したWebAPI URL
  43: 
  44:     //Google Cloud Platform APIキー
  45:     //https://cloud.google.com/maps-platform/
  46:     //※Google Maps APIを利用しないのなら登録不要
  47:     var $GOOGLE_API_KEY_1 = '**************************';   //HTTPリファラ用
  48:     var $GOOGLE_API_KEY_2 = '**************************';   //IP制限用
  49: 
  50:     //Yahoo! JAPAN Webサービス アプリケーションID
  51:     //https://e.developer.yahoo.co.jp/register
  52:     //※Yahoo! JAPAN Webサービスを利用しないのなら登録不要
  53:     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の登録方法」を、それぞれ参照されたい。

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

  34: //地図描画サービスの選択
  35: //    0:Google
  36: //    2:地理院地図・OSM
  37: define('MAPSERVICE', 0);
  38: 
  39: //住所検索サービスの選択
  40: //    0:Google
  41: //    1:Yahoo!JAPAN
  42: //   11:HeartRails Geo API
  43: //   12:OSM Nominatim Search API
  44: define('GEOSERVICE', 11);
  45: 
  46: //逆ジオコーディングサービスの選択
  47: //    0:Google
  48: //    1:Yahoo!JAPAN
  49: //   11:HeartRails Geo API
  50: //   21:簡易ジオコーディングサービス
  51: define('REVGEOSERVICE', 1);

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

解説:初期値

  53: //マップの表示サイズ(単位:ピクセル)
  54: define('MAP_WIDTH',  600);
  55: define('MAP_HEIGHT', 400);
  56: //マップID
  57: define('MAPID', 'map_id');
  58: 
  59: //初期値
  60: define('DEF_LONGITUDE0', 165.0);        //中心座標(経度)
  61: define('DEF_LATITUDE0',  30.5);     //         (緯度)
  62: define('DEF_TYPE',      'roadmap');     //マップタイプ
  63: define('DEF_ZOOM',      4);             //ズーム
  64: define('DEF_CATEGORY', 'world');        //カテゴリ
  65: 
  66: //発射点
  67: define('DEF_FROM',      '東倉里ミサイル発射場');        //地名
  68: define('DEF_CAT_FROM',  'world');                       //カテゴリ
  69: define('DEF_LATITUDE1',  39.660075);                    //緯度
  70: define('DEF_LONGITUDE1', 124.705294);                   //経度
  71: 
  72: //着弾点
  73: define('DEF_TO',        'ホノルル');                    //地名
  74: define('DEF_CAT_TO',    'world');                       //カテゴリ
  75: define('DEF_LATITUDE2',  21.308889);                    //経度
  76: define('DEF_LONGITUDE2', -157.826111);                  //経度
  77: 
  78: define('DEF_HM', 1500);             //弾道高(km)

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

解説:マーカー設定など

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

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

 502: /**
 503:  * 射程,射角から軌道要素を求める
 504:  * @param   double $dt   射程(メートル)
 505:  * @param   double $aof  射角(度)
 506:  * @param   array  $elements  軌道要素を格納する配列
 507:  * @return  なし
 508: */
 509: function getOrbitalElementSub($dt, $aof, &$elements) {
 510:     static $r = 6378136;        //地球半径(メートル)
 511: 
 512:     //ラジアン変換
 513:     $aof = deg2rad(90 - $aof);
 514: 
 515:     $sra = $dt / (2 * $r);      //半射程角
 516: 
 517:     $ax = pow(sin($aof), 2* (pow(cos($sra), 2- pow(sin($aof), 2));
 518:     $bx = pow(sin($aof), 2* (1 - pow(cos($sra), 2));
 519:     $cx = pow(cos($sra), 2- 1;
 520: 
 521:     $vo = sqrt((-$bx + sqrt(pow($bx, 2- $ax * $cx)) / $ax);   //燃え切り速度比
 522:     $e = sqrt(1 - pow($vo, 2* (2 - pow($vo, 2)) * pow(sin($aof), 2)); //離心率
 523:     $a = $r / (2 - pow($vo, 2));        //長径
 524:     $b = $a * sqrt(1 - pow($e, 2));     //短径
 525:     $hmm = $a * (1 + $e- $r;          //弾道高
 526: 
 527:     $elements['sra'] = $sra;
 528:     $elements['vo']  = $vo;
 529:     $elements['e']   = $e;
 530:     $elements['a']   = $a;
 531:     $elements['b']   = $b;
 532:     $elements['hmm'] = $hmm;
 533: }

 535: /**
 536:  * 射程,弾道高から軌道要素を求める
 537:  * @param   double $dt  射程(メートル)
 538:  * @param   double $hmm 弾道高(メートル)
 539:  * @param   array  $elements  軌道要素を格納する配列
 540:  * @return  なし
 541: */
 542: function getOrbitalElement($dt, $hmm, &$elements) {
 543:     $flag = TRUE;
 544:     $aof = 0.0;     //射角(度)
 545: 
 546:     //誤差0.1度で射角を求める
 547:     while ($flag && ($aof < 90.0)) {
 548:         $aof +1.0;
 549:         getOrbitalElementSub($dt, $aof, $elements);
 550:         if ($elements['hmm'] == $hmm) {
 551:             $flag = FALSE;
 552:             break;
 553:         } else if ($elements['hmm'> $hmm) {
 554:             $aof -1.0;
 555:             while ($flag && ($aof < 90.0)) {
 556:                 $aof +0.1;
 557:                 getOrbitalElementSub($dt, $aof, $elements);
 558:                 if ($elements['hmm'>$hmm) {
 559:                     $flag = FALSE;
 560:                     break;
 561:                 }
 562:             }
 563:         }
 564:     }
 565:     $elements['dt'] = $dt;
 566:     $elements['aof'] = $aof;
 567: }

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

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

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

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

 470: /**
 471:  * 軌道要素から、飛翔経路長、速度を求める
 472:  * @param   array  $elements  軌道要素
 473:  * @return  なし
 474: */
 475: function getBallistic(&$elements) {
 476:     static $u = 3.98602e+14;        //地心引力係数
 477:     $n = 100;                       //分割数
 478: 
 479:     $dw = $elements['sra'] / $n;
 480:     $ds = 0.0;                      //飛翔経路長
 481:     $vs = 0.0;
 482:     $vmax = 0.0;                    //最大速度
 483:     for ($i = 0$i <100$i++) {
 484:         $r = $elements['a'* (1.0 - pow($elements['e'], 2)) / (1.0 - $elements['e'* cos($dw * $i));
 485:         $v = sqrt($u * ((2.0 / $r- (1.0 / $elements['a'])));
 486:         if ($v > $vmax$vmax = $v;
 487:         $vs +$v;
 488:         if ($i == 0)    $tmp1 = $r;
 489:         else            $tmp2 = $r;
 490:         if ($i > 1) {
 491:             $dx = $tmp2 * sin($dw);
 492:             $dy = $tmp1 - $tmp2 * cos($dw);
 493:             $ds +sqrt(pow($dx, 2+ pow($dy, 2));
 494:             $tmp1 = $tmp2;
 495:         }
 496:     }
 497:     $elements['ds'] = $ds * 2;
 498:     $elements['vm'] = $vs / ($n + 1);
 499:     $elements['vmax'] = $vmax;
 500: }

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

解説:軌道要素の計算

 865:     //軌道要素の計算
 866:     if (($hm > 0&& ($items['longitude1'!''&& ($items['latitude1'!''&& ($items['longitude2'!''&& ($items['latitude2'!'')) {
 867:         $dt = $pgc1->distance($items['longitude1'], $items['latitude1'], $items['longitude2'], $items['latitude2']);
 868:         $items['hmm'] = $hm * 1000;
 869:         getOrbitalElement($dt, $hm * 1000, $items);
 870:         getBallistic($items);
 871:     }
 872: 
 873:     //地図描画サービス
 874:     switch (MAPSERVICE) {
 875:         //Google
 876:         case 0:
 877:             $call  = jsDragMarker_Gmap($items2[$i]['latitude'], $items2[$i]['longitude']);
 878:             $call .jsDrawLine_Gmap(1, $items1, 'rgb(255, 0, 0)', 1);
 879:             $call .jsDrawLine_Gmap(2, $items2, 'rgb(0, 0, 0)', 1);
 880:             $zd = 1;
 881:             break;
 882:         //Yahoo!JAPAN
 883:         case 1:
 884:             $call  = jsDragMarker_YOLPmap($items2[$i]['latitude'], $items2[$i]['longitude']);
 885:             $call .jsDrawLine_YOLPmap(1, $items1, 'FF0000', 1);
 886:             $call .jsDrawLine_YOLPmap(2, $items2, '000000', 1);
 887:             $zd = 0;
 888:             break;
 889:         //地理院地図・OSM
 890:         case 2:
 891:             $call  = jsDragMarker_Leaflet($items2[$i]['latitude'], $items2[$i]['longitude']);
 892:             $call .jsDrawLine_Leaflet(1, $items1, 'FF0000', 1);
 893:             $call .jsDrawLine_Leaflet(2, $items2, '000000', 1);
 894:             $zd = 0;
 895:             break;
 896:         default:
 897:             $call = '';
 898:             $zd = 0;
 899:             break;
 900:     }
 901:     $lat = $items2[$i]['latitude'];
 902:     $lng = $items2[$i]['longitude'];
 903: 
 904: //エラー発生時
 905: else {

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

結語

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

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

参考サイト

参考書籍

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