都道府県の長さ

ACTIVE GALATTICさんが，次のように書いていた。

中核都市の金沢と被害が激しい輪島は東京と水戸くらい離れていると石川県の広さを説明する投稿を先日見かけたが、石川県の広さは34~35位とランキング下位なので、むしろ東京や大阪が都道府県としては異様に狭いと認識すべきなのだろう

面積と距離が微妙に絡まっていて，ややモニョる。ある図形が与えられたときに，その面積ではなくて，長さを図形の大きさの指標として与えることはできるだろうか。地図上の領域だと東西南北の端点を考えて，南北間距離，東西間距離を抽出できるが，これでは2つの数字の組になる。あるいは図形の周長も考えられるが，フラクタル的な図形だと面積が等しくても周長はいくらでも長くできる。

そこで考えたのが次のような指標だ。図形を格子状の正方形に分割する。これらの正方形の全ての対を考え，各対の正方形の中心間の距離の平均値を求める。分割する正方形のサイズを小さくした極限として，一つの長さが定まるのではないか。（方法１）

検索してみると，もっとスマートな方法があった。図形内の点をランダムに選ぶ。その際に縦横方向にそれぞれ座標値が一様分布であるとした確率分布を仮定する。こうして選んだ2点間の距離の期待値を求めればよいというものだ。（方法２）

まず，方法１のプログラムをJuliaで作成して計算してみた。

a=zeros(Int64,1000001,2)

function ju(a,m,n)

  k = 0

  for i in 1:m

    for j in 1:n

      k = k + 1

      a[k,1] = i

      a[k,2] = j

    end

  end

end

function su(a,m,n)

  mn = m*n

  sum = 0

  for i = 1:mn

    for j = i:mn

      sum = sum + sqrt((a[i,1]-a[j,1])^2+(a[i,2]-a[j,2])^2)

    end

  end

  return sum/binomial(mn,2)

end

m=800

n=800

ju(a,m,n)

@time su(a,m,n)/m

n=   25    time= 0.016720　0.521829

n=   50　time= 0.018188　0.521513

n= 100　time= 0.065998　0.521433

n= 200　time= 0.788281　0.521412

n= 400　time= 12.1924  　0.521407

n= 800　time= 198.615  　0.521405

nは正方形の1辺の格子点の数，timeは計算時間(秒)であり，最後が計算結果の平均長である。

方法２による理論値が解析的に得られており，( 2+√2+5 log(√2 + 1) )/15 = 0.521405 なので，nを増やすとともに理論値に収束し，n=800では理論値に一致している。

m = 50; n = 50;

a = Flatten[Table[{i, j}, {i, 1, m}, {j, 1, n}], 1];

Timing[Sum[

     Sqrt[(a[[i, 1]] - a[[j, 1]])^2 + (a[[i, 2]] - a[[j, 2]])^2], 

　　　{i, 1, n*n}, {j, i, n*n}]/Binomial[n*n, 2]/n // N]

Mathematicaでは，若干コードが簡単になるけれど，時間は n=50 で 12.4801秒と，600倍以上かかってしまった。

（付）このコードはそのまま乱数による計算に置き換えることもできる。Juliaでは，a[k,1] = rand()，a[k,2] = rand() として最後の結果の/nを除く。Mathematicaでは，Tableの中身を{i,j}から{RandomReal[1], RandomReal[1]}に置き換えて最後の結果の/nを除く。

［１］地域内距離（腰塚武志）

円軌道はむずかしい

万里鏡１号・弾道ミサイルの軌道（２）からの続き

北朝鮮の弾道ミサイルの簡単なシミュレーションコードをMathematicaで作っていた。これを少しアレンジすれば人工衛星を軌道に投入するところまでできそうな気がする。

早速，以前のコードを修正してみた。まずは通常の加速直後に角度方向だけに加速度を加えるようにしたがうまくいなかい。打ち上げ加速は投射角の方向になっているので，動径速度成分が大きく残っているうえ，加速すれば軌道は膨らむ。このため，離心率の大きな長円軌道になって地表にぶつかってしまうか，地球の重力圏から脱出してしまうのだ。

次に，打ち上げ加速の直後に空白時間をおいて，動径速度成分が小さくなったところで角度方向に加速できるようにした。それでもうまくいかない。簡単な試行錯誤では周回軌道にのせるのが難しい。そもそも角度方向に加速するということは面積速度すなわち角運動量をふやし，動径方向の微分方程式で軌道半径を膨らませる方向に作用してしまう。

そこで，後期加速では衛星をその速度ベクトルの方向に加速することにした。$t=0$で速度ベクトルをとりだすところに発散があったので，これを回避するため，地球の自転による2倍面積速度$h(t)$の初期値として，$h(0) = R^2\omega=R^2 \frac{2\pi}{24*3600}=2930$ km$^2/$sを与えた。初期加速度はこれまでの$\ a=0.0446$として$30$秒加速する。その後，800秒程度休止した後に，後期加速度$\ b=0.1445$（ここを微調整した）で$250$秒加速すると，なんとか軌道に投入することができた。投射角は$s=45$度，初期加速における燃料比は$p=0.85$であり，加速方向の角度には$0.3$をかけて動径成分を抑えた。

なかなか難易度の高いゲームである。衛星の軌道高度が1200km程度の準円軌道となっている。これを500kmにしなさいといわれても，こんな単純な2段階制御ではちょっと難しい。なお，プログラムの検証のため，$r=6850$kmの宇宙空間で第一宇宙速度に相当する$v=\sqrt{gr}=8.2$ km/sを角度方向の初速度として与えると，正確に円軌道を描くことが確かめられた。

g = 0.0098; R = 6350; τ = 30; τs = τ*27; τt = 250; p = 0.85;

 a = 0.0446; b = 0.1445 a; s = 45 Degree; T = 15400; 

fs[t_] := 0.3*ArcTan[r[t]*r'[t]/ h[t]]

fr[t_, τ_] :=  a*Sin[s]*HeavisideTheta[τ - t] + 

   b*Sin[fs[t]]*HeavisideTheta[t-τs-τ]

   *HeavisideTheta[τ+τs+τt-t]

ft[t_, τ_] :=  a*Cos[s]*r[t]*HeavisideTheta[τ - t] + 

   b*Cos[fs[t]]*r[t]*HeavisideTheta[t-τs-τ]

   * HeavisideTheta[τ+τs+τt-t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] 

   +h[t]^2/r[t]^3 -g R^2/r[t]^2 +fr[t, τ],

   r[0] == R, r'[0] == 0, 

   h'[t] == -fm[t, τ]*h[t] + ft[t, τ], 

   h[0] == 2930 + 0*Sqrt[g] R^(3/2)}, {r, h}, {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]},

 {t, 0, T}, PlotRange -> {-5, 15}]

 

図：苦労すると有難みがわかる衛星の準円軌道のグラフ

P. S. もう少しがんばると，軌道高度650km（r=6980km）の準円軌道まで達成できた。

g = 0.0098; R = 6350; τ = 25; τs = τ*15.3; τt = 350;  p = 0.85;

a = 0.0446; b = 0.1275 a; s = 45 Degree; T = 15400; 

出生時の体重分布

赤ちゃんの出生時の体重分布がないのか調べてみた。500gきざみはすぐ見つかったのだが，もう少し精度がないかと探したところ，統計でみる日本 e-Stat にあった。人口動態調査 人口動態統計 確定数 保管統計表 出生 ／ 出生数，性・出生時の体重（100g階級）・出生時の身長別 というものだ。

これをExcelでグラフにするが，普通の正規分布と比較しやすくするため，縦軸の100g区間ごとの出生数を全体の出生数に対する割合（%）に直した。なお，これは2021年（最新）のデータであり，男子の（出生数，平均値，標準偏差）は（415,836人，3.047 kg，0.390 kg），女子のそれは（395,687人，2.962 kg，0.358 kg）となっている。

図１：女子（青）と男子（オレンジ）の出生数割合（横軸100gビン数）

正規分布の方は，Mathematicaの組み込み関数を使うのだけれど，100g刻みにしているデータと比較するため，μを10倍，σを10倍にしたものを使っている。

In[1]:= f[\[Mu]_, \[Sigma]_] := 

 Integrate[ Exp[-(x - \[Mu])^2/(2*\[Sigma]^2)]/

  Sqrt[2*Pi*\[Sigma]^2], {x, -Infinity, Infinity}]

In[2]:= f[29.62, 3.58]

Out[2]= 1. - 8.59978*10^-15 I

In[3]:= g[\[Mu]_, \[Sigma]_] := 

 Plot[Exp[-(x - \[Mu])^2/(2*\[Sigma]^2)]/Sqrt[2*Pi*\[Sigma]^2],

 {x, 0, 20 \[Sigma]}, PlotRange -> {0, .12}]

In[4]:= g1 = g[29.62, 3.58];

In[5]:= g2 = g[30.47, 3.90];

In[6]:= Show[g1, g2]

Out[6]= -Graphics-

In[7]:= Plot[

 {PDF[NormalDistribution[29.62, 3.58], x] // Evaluate, 

  PDF[NormalDistribution[30.47, 3.90], x] // Evaluate},

 {x, 0, 50}, Filling -> Axis, PlotRange -> {0, 0.13}]

Out[7]= -Graphics-

ほぼそれらしいような気もするが，微妙に正規分布からずれているのかもしれない。例えば，青の4 kg の位置では累積割合で 99.5% になる。この位置は正規分布では，2.86σ相当に対応する。それで積分すると99.8%になってしまうので，少しずれているようだ。

図２：平均値と標準偏差をそろえた正規分布

折り紙で３次方程式（２）

折り紙で３次方程式（１）からの続き

Mathematicaで折り紙の3次方程式問題をplotする関数を書いてみた。

最初に，Plot[{-x - 3, x + 1, 2 x, (x + 2)^2/4,   Sqrt[8 (x)] - 1, -Sqrt[8 (x)] - 1}, {x, -8, 8},  AspectRatio -> Automatic, PlotRange -> {-8, 8}]　こんな風に問題ごとにプロットしてみた。係数を混同していてなかなか正解にたどりつかなかったが，どうにか安定して答えがでてきた。ただし，接線の方程式は目の子であてはめている。

係数，a,b,c,dを与えて，焦点，準線，放物線，接線をすべて描画するように改善したものが次のコードである。最終ステップで，三次方程式を解いた解が3つリストとして出てくるので，これを実数の個数によって場合分けして描画するのはさぞ面倒だろうと心配していたのだが，リストの成分表示からリストに直しただけで求めていた結果が得られてしまった。虚数解が入ってきても何の問題もなくスキップしてくれた。やはりMathematicaは賢い。

そのコードは次のようなものである。

fold3[a_, b_, c_, d_] := 

{g0 = Graphics[{Point[{b, a}], Point[{d, c}], 

      Line[{{-10, -a}, {10, -a}}], Line[{{-d, -10}, {-d, 10}}]}, 

     Axes -> True]; 

 g1 = Plot[{(x - b)^2/(4*a), Sqrt[4 d x] + c, -Sqrt[4 d x] + c},

    {x, -10, 10}, PlotRange -> {-10, 10}, AspectRatio -> Automatic]; 

 sol = Solve[a t^3 + b t^2 + c t + d == 0, t]; 

 t1 = t /. sol; 

 s1 = -( a t1*t1 + b t1);

 g2 = Plot[t1 x + s1, {x, -10, 10}, PlotRange -> {-10, 10}, 

     AspectRatio -> Automatic, PlotStyle -> Red];

 Show[g0, g1, g2]}[[1]]

この関数で，fold3[1, 2, -2, -1]，fold3[1, 1, 1, 1]，fold3[1, 0, -1, 0.0001] を実行した結果を以下に示す。傾きが0の場合は，極限値として得られることになる。

図：fold3の実行結果

ブランコの物理

Physical Reveiw E に，Initial phase and frequency modulations of pumping a playground swing という日本のグループの論文が掲載された。ブランコを漕いで揺らすコツを解明したというものだ，たぶん。アブストラクトをDeepLにかけてGPT-4で整えると次のようになる。

遊具のブランコは、ブランコ本体とそれを揺らす人間からなる動的な結合振動子系である。本研究では、上半身の自然な動きの初期段階がブランコの連続的なポンピングに与える影響を分析するモデルを提案し、3種類の長さのブランコで10人の参加者が実際にポンピングした際の運動データを用いて検証する。

研究の結果から、スイングが垂直位置（中点）にあり、振幅が小さい時に前方へ移動する際、最大リーンバックの初期段階が発生し、最大振幅のポンピングが予測された。振幅が大きくなるにつれて、最適な初期段階はサイクルの初期段階、すなわちスイングの軌道の後端に向かって徐々に移動する。モデルの予測に従って、参加者全員が、スイングの振幅が大きくなるにつれて上半身の動作の初期位相が早くなることが確認された。

この結果から、スイングする人は、上半身の動作の頻度と初期位相を適切に調整することで、ブランコを効果的にポンピングすることができると考えられる。

人間の体を3つの線質量要素で表現し，上半身と下肢の角度をそれぞれ変数として導入したラグランジアンを考えている。このとき，これらの身体の角度については，何通りかのあらかじめ定められたモデルの周期運動だとしている。

現実に近いのだけれど面倒なモデルなので今一つピンと来ない。そこで，人体を振り子につながって距離が固定された２質点で表し，その間で周期的な質量移動をすることでスウィングするというモデルを考えてみた。質量移動は，位相と角振動数の２つの変数で制御できる。速度に比例する減衰振動のある単振子をこの質量移動で励振するというものだ。

Do[

 s0 = Pi/6*0; w0 = Pi; m1 = 1; m2 = 1;

 λ = Pi/12; ν = 0.05/Pi; ε = 0.95 + i*0.01; δ = Pi/8*0;

 m[t_] := (m1 + m2)/20 * Cos[ε w0 t + δ];

 sol = NDSolve[

       {(m1 + m2) s''[t] == -w0^2 * 

        ( (m1 + m[t]) * (Sin[s[t] + λ + ν * s'[t]) 

        + (m2 - m[t]) * (Sin[s[t] - λ + ν * s'[t])), 

        s[0] == s0, s'[0] == 0}, s, {t, 0, 60}];

  s /. sol[[1]];

  f[t_] := s[t] /. sol[[1]];

  g[t_] := s'[t] /. sol[[1]];

  h[i] = Plot[f[t], {t, 0, 30}, AspectRatio -> 1/2, 

   PlotRange -> {-0.6, 0.6}, PlotStyle -> Hue[0.08*i]],

 {i, 0, 7}]

もとの振動子の初期位相は s0 であり，固有角振動数は w0 (=√g/ℓ)，２質点の距離は角度 λ で表現している。速度に比例した減衰係数は ν，質量移動の角振動数と固有角振動数の比が ε，質量移動の位相のズレが δ である。以下の図は，ε = 0.95 + i*0.01 としたもの（上）と，δ = Pi/8*i; としたもの（下）である。固有角振動数比が 0.97のときが最も励振が大きく，2%ふえるだけで減衰に移行してしまう。

図：振幅の質量移動振動への依存性，振動数（上）と位相（下）

凹面鏡

軸対称な凹面鏡の断面の曲線が放物線ならば，鏡の対称軸＝主軸に平行に入射する光線は反射して焦点に集まることが良く知られている。この曲線が一般の形の場合は主軸近傍では2次関数で近似できるので同様に焦点に光が集まる。近傍から外れた場合に生ずる焦点からのズレが収差となる。ところで，焦点を結ぶのは2次関数に限られるのだろうか？

$x$軸上の原点Oに凹面鏡の底を接地させ，$y$軸上にある焦点Fと原点Oの距離を$f$とする。原点近傍の凹面鏡の断面の曲線上の点P$(x, y(x))$について，直線FPの傾きは$-\frac{f}{x}$と近似できる。したがって，点Pにおける鏡面の傾きは，$\frac{x}{2f}$となる。つまり，$\frac{dy}{dx}= \frac{x}{2f}$という微分方程式が成り立つ。つまり，少なくとも原点近傍では $y(x) = \frac{x^2}{4f}$という2次関数でなければならない。

Mathematicaでのシミュレーションコードを書いてみた。

f2[a_] := a^2

t2[a_] := ArcTan[D[f2[x], x]] /. x -> a 

s2[a_] := Tan[Pi/2 + 2*t2[a]]

g2[y_, a_] := s2[a]*(y - a) + f2[a]

TrigExpand[f2[a] - a*s2[a]] // Simplify

Out[-]= 1/4

p0 = Plot[f2[x], {x, 0, 1}, PlotStyle -> {Red}]; 

p2 = Table[Plot[g2[y, 0.1*n], {y, 0, 0.1*n}], {n, 1, 8}];

q2 = Table[b = 0.1*n; 

   Graphics[ Line[{{b, 1}, {b, f2[b]}}, 

   VertexColors -> {Green, Blue}]],{n, 1, 8}];

Show[{p0, p2, q2}, PlotRange -> All, AspectRatio -> Automatic]

図：凹面鏡（放物面鏡）の断面図と焦点への結像

先ほどの議論を一般化する。凹面鏡の断面の曲線を$y(x)$とする。点P$(x,y(x))$における接線の傾きは，$y'(x)$であり，点Pにおける入射光線と法線のなす角度は，$\theta = \arctan{ y'(x)}$。そこで，点Pにおける反射光線の傾きは，$m = \tan(\frac{\pi}{2} +2 \theta) = -\frac{\cos 2 \theta}{\sin 2 \theta}  = \frac{y'(x)^2-1}{2 y'(x)}$となる。

さて，反射光線の方程式は，$Y-y(x)=m(X-x)$である。したがって，焦点の位置 を$(X,Y)=(0,f)$とすると，$f = y(x) - x \cdot  \frac{y'(x)^2-1}{2 y'(x)}$という微分方程式で定まる。先ほどの，$y(x)=\frac{x^2}{4f}$はこの方程式を満足している。

wolframscript

Wolfram Engine からの続き

ChatGPTは，因数分解ができないので，Wolfram Alphaと連携すればよいという話について。

ChatGPT数学が苦手である。Stephen Wolframもそう言っている。BingChatに聞いてみると，WolframAlphaのAPIを使う必要があるようだ…というかそんなことは先刻承知。それにはライセンスが必要なので少々ためらう。何か他の方法はないものかと眺めていると，wolframscriptというWolfram言語のコマンドライン版が見つかった。

とりあえずインストールしてみようと（すでに当初の目的から話がずれつつある），確かめてみれば，自分のMacBook Air M1には既にインストール済みだった。以前，JupyterでWolfram言語を使うためにWolfram Engineを導入した時にインストールしていたらしい。うかつなことだ。

最近，Mathematica本体を立ち上げるのが億劫な場合，コマンドラインでフリーのMathicsを使うことが多かった。ただ機能はかなり制限されている。今回，はじめて wolframscript を実行してみると，Mathicsよりは大分制限が少ないというか，インストール済みのMathematica 12.0 の機能がフルに使えるのかもしれない。そこで次のような単独で実行可能なスクリプトも組めるのだった。拡張子は .wls である。

#!/usr/local/bin/wolframscript

(* ::Package:: *)

n=ToExpression[$ScriptCommandLine[[2]]]

Print[RandomInteger[{1,6},n^2]]

コマンドライン引数（文字列変数）が，$ScriptCommandLine[[2]]以下で取得できる。これを整数型に変換するために，ToExpression[] 関数を使った。これで，プログラム外から与えた変数に対して，Wolfram言語での処理が可能になる。

［１］LangChainのWolfram Alpha連携を試す（npaka）

［２］WolframScript（Wolfram言語＆システム ドキュメントセンター）

弾道ミサイルの軌道（６）

弾道ミサイルの軌道（４）からの続き

H3ロケットの発射がフェイルセーフで中止された次の日，なんだか間の悪いことに北朝鮮のICBM級ミサイルが平壌近傍から発射され，北海道渡島大島西200kmに落下した。

飛行時間は66分，飛行距離は900km，最高高度は5700km，落下速度はマッハ20程度と推測される。ペイロードを調整すれば，射程が14,000kmでアメリカ本土が射程に入ると報道された。

それでは，以前作ったのプログラムで確認してみよう。一生懸命Juliaのフォルダを探していたが，Mathematicaのコードだった。昨年11月のコードのパラメータで，飛行時間と飛行距離を微調整した結果が次の通り。

g = 0.0098; R = 6350; τ= 86; p = 0.75; 

a = 0.0446; s = 86.5 Degree; T = 3960; 

(* g = 0.0098; R = 6350; τ = 87; p = 0.75; a

 = 0.0446; s = 45 Degree; T = 5100; *)

fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ - t]

ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ - t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 - 

 g R^2/r[t]^2 + fr[t, τ], r[0] == R, r'[0] == 0, 

 h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r, h},

 {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, 

{t, 0, T}, PlotRange -> {-4, 8}]

tyx[T_] := {T, f[T] - R, NIntegrate[R d[t]/f[t]^2 , {t, 0, T}]}

v[T_] := Sqrt[(f[T] - f[T - 1])^2 + (R d[T]/f[T]^2)^2]

{tyx[T], v[T]/0.340}

g0 = ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}], 

  f[tt] - R}, {tt, 0, T}]

燃焼時間τが86秒，燃料重量比が0.75，燃焼加速度が 0.0446km/s^2，投射角が86.5度である。これで飛行時間を与えると，飛行距離と最高高度と落下時速度が概ね再現できる。万有引力は距離の2乗に反比例し，コリオリ力や空気抵抗は無視，地球の形を考慮して2次元極座標で質量が変化する1段ロケットの微分方程式を解いている。

燃焼時間τを87秒にして，投射角を45度にすれば，飛行時間が5100秒で飛行距離は14,600km，落下時速度はマッハ24になる。このときペイロードはほとんど変化させていない。燃料重量比は同じで燃焼時間を1秒のばしただけだ。

図：投射角86度のロフテッド軌道と投射角45度の弾道軌道

軸対称電荷分布（５）

軸対称電荷分布（４）からの続き　　

クーロンポテンシャルを場合分けによって書くことができたが，実際に計算する場合には，いちいち分けるのは面倒だ。といっても高々六類型なのでどうってことはないといえばない。

軸対称電荷分布（２）でやったように，ヘヴィサイドの階段関数を使えば，$r_<,\ r_> \ $が表現できるので，クーロンポテンシャルのルジャンドル成分$\ \widetilde{V}_n(r)\ $を一つの関数を使って表わせる。そこで，$R(t)=\dfrac{1}{\sqrt{1+t^2}}, \ a=1.0,\ c= 0.707107$の場合について具体的に計算してみた。

rs[s_, r_] := s* HeavisideTheta[r - s] + r*HeavisideTheta[s - r];

rl[s_, r_] := r* HeavisideTheta[r - s] + s*HeavisideTheta[s - r];

R[t_] := 1/Sqrt[1 + t^2];

v[n_, r_] := 3*NIntegrate[ 1 + 0 * LegendreP[n, t] 

    *rs[s, r]^n * s^2/rl[s, r]^(n + 1), {t, 0, 1}, {s, 0, R[t]}]

v0 = v[0, 0]

2.64412

Plot[{v0, v[0, r]}, {r, 0, 3}]

Plot[{v0, v[2, r]}, {r, 0, 3}, PlotRange -> {2.55, 2.65}]

Plot[{v0, v[4, r]}, {r, 0, 3}, PlotRange -> {2.64, 2.665}]

Plot[{v0, v[6, r]}, {r, 0, 3}, PlotRange -> {2.638, 2.648}]

ここで，被積分関数の中に 1 を加えたのは，積分値が 0 になって計算が不安定になるのを防ぐためである。その部分の定数$\ v_0\ $を取り除くために，ルジャンドル関数の前に 0 をかけて積分することで$\ v_0\ $を求め，実際の計算ではこれを取り除く。

その結果は次のようになった。$\widetilde{V}_4(r),\ \widetilde{V}_6(r)$が，$r<c\ $でゼロになるのは，$\int ds\ $からでてくる$\ R(t)^{n-2}\ $がルジャンドル関数の$n$次未満の線形結合になるためである。

図：クーロンポテンシャルのルジャンドル成分 $\widetilde{V}_n(r)\  (n=0,2,4,6)$

軸対称電荷分布（３）

軸対称電荷分布（２）からの続き　

今から25年前，1997年度の出世・松井・渡会の卒業論文では，変形核の軸対称電荷分布の形として原点からの距離を$R(\cos\theta)$として，2種類の関数形を採用した。なお，$\theta$は対称軸である$z$軸正方向からの方位角を表わし，$x-y$平面に対して対称な分布を考えるので，$R(\cos\theta)$は$\cos\theta$の偶関数に限定される。$t=\cos\theta$とおいて，電荷密度$\rho_0$は，$Ze = 4 \pi \rho_0 \int_0^1 \int_0^{R(t)} s^2 ds\ dt $で与えられる。

(1) 回転楕円体：$R(\cos\theta) = \dfrac{ac}{\sqrt{c^2+(a^2-c^2)\cos^2\theta}}$である。

この場合の電荷密度は，$\rho_0=\dfrac{3Ze}{4\pi a c^2}$となる。

(2) 2次ルジャンドル関数：$R(\cos\theta) = a+(c-a) \cos^2\theta = \alpha + \beta P_2(\cos\theta)$である。

ただし，$\alpha = \dfrac{2a+c}{3c},\quad \beta = \dfrac{2(c-a)}{3}$。

この場合の電荷密度は，$\rho_0=\dfrac{105 Ze}{4 \pi (16 a^3+8a^2c + 6 a c^2 + 5 c^3)}$となる。

ところで，ルジャンドル関数は直交多項式であり，回転楕円体の断面の楕円はルジャンドル関数による展開で近似できる。この近似がどのくらい妥当するかを，長軸と短軸の比が1:2や2:1の場合について確かめてみる。楕円の離心率は，$e=\sqrt{\dfrac{a^2-c^2}{a^2}}$で与えられる，$\varepsilon=\dfrac{e^2 a^2}{c^2}$とおくと，$R(t) = \dfrac{a}{\sqrt{1+ \varepsilon t^2}}$となる。まず，$n=4$までの展開係数を求めておく。

$\dfrac{1}{\sqrt{1+ \varepsilon t^2}} \approx 1 - \dfrac{\varepsilon}{2}t^2 + \dfrac{3}{8}\varepsilon^2 t^4 = \alpha P_0 (t)+ \beta P_2(t) + \gamma P_4(t)$

$= \alpha + \beta \dfrac{3 t^2-1}{2} + \gamma \dfrac{35t^4-30t^2+3}{8}=1-\dfrac{\beta}{2}+\dfrac{3\gamma}{8} +\Bigl( \dfrac{3\beta}{2}-\dfrac{15\gamma}{4}\Bigr) t^2+\dfrac{35\gamma}{8}t^4$

これを解いて，$\alpha = 1 -\dfrac{\varepsilon}{6} + \dfrac{3\varepsilon^2}{40}$，$\beta = -\dfrac{\varepsilon}{3} + \dfrac{3 \varepsilon^2}{14}$，$\gamma = \dfrac{3\varepsilon^2}{35}$

a = 0.5; c = 1.0;

(* a = 1.0; c = 0.5; *)

d[a_,c_,t_] := a/Sqrt[1 + (a^2-c^2)/c^2 * t^2]

4 Pi * Integrate[Integrate[r^2, {r,0,d[t,a,c]}], {t,0,1}]

g0 = Plot[d[t,a,c], {t,0,1}, PlotStyle->Hue[0]]

f[n_,a_,c_] := Integrate[LegendreP[n,t]*d[a,c,t], {t,0,1}]

b = Table[f[n,a,c], {n,0,16,2}]

g[m_] := Plot[Sum[b[[n/2+1]]*LegendreP[n,t]*(2*n+1), {n,0,m,2}],

         {t,0,1}, PlotStyle-> Hue[(m+1)/10]]

Show[Table[{g0, g[k]}, {k, 4, 6, 2}], PlotRange -> {0.4,1}]

図：軸対称電荷分布（左 oblate，右 prolate，横軸 $t=\cos\theta$，縦軸 $R(t)$）

扁平率 $f=\dfrac{a-c}{a}\ $が-1〜1/2の範囲では，4次のルジャンドル展開で回転楕円体をほぼ表現できる。6次ならば十分だ。

軸対称電荷分布（２）

軸対称電荷分布（１）からの続き

前回の$z$軸対称で，$x-y$平面にも対称な電荷分布のつくるクーロンポテンシャルは，

$\displaystyle V(r, s) = \dfrac{-e}{4 \pi \varepsilon_0} \sum_{n=0}^\infty  P_n(s)  \int_{0}^{1}  \int_0^{R(t)} \rho(r', t) P_n (t) \Bigl\{   \dfrac{1}{r_>} \bigl (\dfrac{r_<}{r_>}\bigr)^n \Bigr\}  4\pi r'^2 dr' d t\ $であった（ただし，$s=\cos\theta$，$t=\cos\theta'$）。

最初に球対称の場合で，上記の式を確認する。このとき$R(t) =R$， $\rho(r', t) = \rho(r') = H(R-r') \rho_0$である。$H(x)\ $はヘヴィサイドの階段関数，$\dfrac{4\pi}{3} R^3 \rho_0 = Ze $である。

このとき$\ n=0\ $の項だけが残り，$\displaystyle V(r) =  \dfrac{-e}{4 \pi \varepsilon_0}  \int_0^{R} \rho_0  \dfrac{1}{r_>}  4\pi r'^2 dr' $

(1) $r>R\ $の場合：

$\displaystyle V(r) =  \dfrac{-e}{4 \pi \varepsilon_0}  \int_0^{R} \rho_0  \dfrac{1}{r}  4\pi r'^2 dr'  =  \dfrac{-e}{4 \pi \varepsilon_0}  4\pi  \rho_0  \dfrac{1}{r} \int_0^{R}  r'^2 dr'  =   \dfrac{-Ze^2}{4 \pi \varepsilon_0} \dfrac{1}{r} $

(2) $r<R\ $の場合：

$\displaystyle V(r) =  \dfrac{-e}{4 \pi \varepsilon_0}  4 \pi \rho_0 \Bigl\{ \int_0^r   \dfrac{1}{r}  r'^2 dr'  + \int_r^R  \dfrac{1}{r'} r'^2 dr'  \Bigr\} $

$\displaystyle  =  \dfrac{-e}{4 \pi \varepsilon_0}  4 \pi \rho_0 \Bigl\{ \dfrac{1}{r} \dfrac{r^3}{3}  +   \left [ \dfrac{r'^2}{2} \right ]_r^R  \Bigr\} =  \dfrac{-Ze^2}{4 \pi \varepsilon_0} \Bigl\{\dfrac{r^2}{R^3} + \dfrac{3(R^2-r^2)}{2R^3}\Bigl\}$

$\displaystyle  =  \dfrac{-Ze^2}{4 \pi \varepsilon_0 }\dfrac{1}{R} \Bigl\{\dfrac{3}{2} - \dfrac{r^2}{2R^2}\Bigl\}$

これで，球対称一様電荷分布によるクーロンポテンシャルが再現できた。

Mathematicaでは次のように計算する。

v[r_] := -3＊NIntegrate[s^2/

       (r*HeavisideTheta[r-s] + s*HeavisideTheta[s-r]), {s, 0, 1}]

Plot[v[r], {r, 0, 5}, PlotRange -> {0, -1.6}]

図：球対称一様電荷分布のクーロンポテンシャル(R=1.0, V(R)=-1.0)

Mathics（２）

Mathics（１）からの続き

ちょうど2年前の今ごろ，Mathicsについての記事をかいていたが，完全に記憶の彼方に跳んでいた。別の記事を復習していてたまたま目に付いたのである。

Mathicsは，Pythonベースのフリーでオープンソースの汎用計算機代数システムであり，Mathematica®と互換性のある構文と関数を特徴としている。早速，MacBook Air M1で動かしてみると動かない。pip3 install Mathics で簡単に復活した。

Mathics 5.0.2

on CPython 3.10.8 (main, Oct 13 2022, 09:48:40) [Clang 14.0.0 (clang-1400.0.29.102)]

using SymPy 1.10.1, mpmath 1.2.1, numpy 1.23.4

Copyright (C) 2011-2022 The Mathics Team.

This program comes with ABSOLUTELY NO WARRANTY.

This is free software, and you are welcome to redistribute it

under certain conditions.

See the documentation for the full license.

前回より，かなりバージョンがあがっている。 GUIの画面も実装しているようなので試してみたが，これがなかなか難渋する。どうやら，Three.jsというJavascriptによる三次元グラフックシステムを利用しているらしい。これも7-8年前に流行ったようだが全く知らなかった。MathicsのGUI版の方は，Mathics-Djangoという名前で，PythonのウェブフレームワークのDjangoの上で動いていた。

pip3 install Mathics-Django 

cd mathics-django-master/mathics_django

*You have 21 unapplied migration(s). 

python3 manage.py migrate

python3 manage.py runserver 

*error AttributeError: module 'scipy.linalg' has no attribute 'inv'

pip3 install scipy-stack

python3 manage.py runserver

*error sqlite3.OperationalError: unable to open database file

*error django.db.utils.OperationalError: unable to open database file

mkdir ~/.local/var/mathics 

djangoは，作業ディレクトリ（~/mathics-django-master/mathics_django）で，manage.py runserver を実行すると，ローカルのウェブサーバー（ http://127.0.0.1:8000/）が立ち上がる仕組みになっている（エラーが出たら指示に従ってディレクトリに必要なファイルを持ってくる）。最初のエラーは，scipy-stackのインストールで，次のエラーは，データベースの場所の問題だったが，~/.local/var/mathicsというデフォルトのディレクトリを作成することでようやく解決することができた。

早速実行してみたが，相変わらず，計算量の制限は大きく，1500!（10進 4114桁）くらいまでしか計算できない。本物のMathematicaならば，10000000!（10進 65 657 059桁）でも平気で計算できるのだけれど。

なお，List of computer algebra systems には無事に掲載されていた。

図：MacBook Air M1にインストールしたMathics GUI版

弾道ミサイルの軌道（５）

弾道ミサイルの軌道（４）からの続き

NHKの報道によれば，11月18日（金）午前10時14分，北朝鮮のICBMが発射され，午前11時23分に北海道渡島大島の西200kmの海上に落下した。飛行時間は69分，飛行距離は1,000km，最高高度は6,100km，速度はマッハ22とされている。

また，弾頭質量によっては15,000kmを越える射程となり，米国本土を狙えるとあった。しかし，今回の実験もロフテッド軌道だったので射出角は90度にかなり近いはずだ。これをより低い角度にすれば質量を変化させなくても十分到達できると思うけれどどうなのか。とりあえず，自分のモデルで計算して確かめてみよう。

g = 0.0098; R = 6350; τ = 87; p = 0.75; a = 0.0446; s =  45 Degree; T = 5100; 

g = 0.0098; R = 6350; τ = 87; p = 0.75; a = 0.0446; s =  86.2 Degree; T = 4140; 

fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ - t]

ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ - t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 - g R^2/r[t]^2 + fr[t, τ], 

r[0] == R, r'[0] == 0, h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r, h}, {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, {t, 0, T},  PlotRange -> {-4, 8}]

tyx[T_] := {T, f[T] - R, NIntegrate[R d[t]/f[t]^2 , {t, 0, T}]}

v[T_] := Sqrt[(f[T] - f[T - 1])^2 + (R d[T]/f[T]^2)^2]

{tyx[T], v[T]/0.340}

{{4140, 41.9866, 1036.42}, 22.9628}

g0 = ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}],  f[tt] - R}, {tt, 0, T}]

このモデルの仮定は次のようなものである。(1) 距離の二乗に反比例する万有引力が働く極座標での運動方程式を解く。(2) コリオリの力や空気抵抗は考えない。(3) 推進剤の燃焼噴射によって，投射体の速度ベクトルの方向に等加速度（a=0.0446 km/s^2）を発射時から一定時間加える。これにより投射体の質量は燃料が0になるまで一定の割合で減少する。(4) 多段式ロケットは考えず，一段式として全質量に対する燃料比（p=0.75）を与える。このとき，運動方程式は全質量に依存しない。(5) 飛行時間の観測値を与え，投射角度と推進加速時間を主パラメータとして，飛行時間と最高高度の観測値が再現されるように調整する。

結論としては，弾頭質量（というか全質量に対する燃料比）を変化させなくても，射出角を86.2度から45度に変更すれば，飛行時間5100秒（85分）で到達距離が14,600kmになる。なお，ロフテッド軌道の場合の弾道ミサイルの最高速度はマッハ23となり，報告されたマッハ22をほぼ再現した。

図：同一パラメタで射出角と飛行時間だけを変えた場合の弾道軌道

P. S. 11月19日の北朝鮮の朝鮮中央通信によれば，このロフテッド軌道の弾道ミサイルは，新型大陸間弾道弾「火星17」の試験発射であり，最高高度 6,040.9 km ，飛行距離 999.2 km ，飛行時間 4135 秒とのことだ。なお，同じ火星17は，今年の3月24日には，最高高度 6,248.5 km ，飛行距離 1090 km ，飛行時間 4052 秒で飛行している。

［１］【詳報】北朝鮮ICBM級ミサイル発射 北海道の西 EEZ内落下か（NHK）

［２］北朝鮮のミサイル等関連情報（続報第２報，防衛省・自衛隊）

　　（防衛省からは，わかった部分から順に発表されている。発表時刻を入れるべきでは）

弾道ミサイルの軌道（４）

弾道ミサイルの軌道（３）・ロフテッド軌道（３）からの続き

11月3日文化の日，午前7時50分に例のけたたましいJ-アラートが新潟県・山形県・宮城県に発出された。7時40分にミサイルが発射されたということだったが，続報では，午前7時48分に日本上空を通過したとなり，その後ミサイルは日本海上空で消失したなどと情報は二転三転する。

防衛省が11月3日の落ち着いた段階で発表したものは次の通り。

北朝鮮は本日7時台から8時台にかけ，ICBM級の可能性があるものも含め，少なくとも3発の弾道ミサイルを，東方向に向けて発射しました。詳細については現在分析中ですが，落下したのはいずれも我が国の排他的経済水域（EEZ）外であり，飛翔距離等については以下の通りと推定しています。

① 7時39分頃，北朝鮮西岸付近から発射し，最高高度約2000 km程度で，約750 km程度飛翔し，朝鮮半島東側の日本海に落下。当該ミサイルはICBM級の可能性があります。

② 8時39分頃，北朝鮮内陸部から発射し，最高高度約50 km程度で，約350 km程度飛翔し，朝鮮半島東岸付近に落下。

③ 8時48分頃，北朝鮮内陸部から発射し，最高高度約50 km程度で，約350 km程度飛翔し，朝鮮半島東岸付近に落下。

なお，日本列島を超えて飛翔する可能性があると探知したものについては，その後，当該情報を確認したところ，探知したものは日本列島を超えず，日本海上空にてレーダーから消失したことが確認されました。

航空自衛隊が設置している警戒管制レーダーのJ/FPS-5（ガメラレーダー），J/FPS-7は弾道ミサイルを探知することができるはずだ。

実際に弾道ミサイルの飛来を探知した場合は，まず航空総隊/BMD統合任務部隊にデータが送信され，防衛省を経て内閣総理大臣および各省庁にデータが送られる。総務省は，全国瞬時警報システム（J-アラート）を通じてマスメディアや地方公共団体に警報を伝達する。国民には，報道や防災無線を通じて屋内退避などの自己防衛が指示される（Wikipedia J/FPS-5から引用）。

ということならば，問題は，Jアラート側というより，航空自衛隊の警戒管制レーダーの精度と分析能力にあることになる（あるいは発出の政治的判断を含めて）。こんな調子で，防衛費を今後5年間で48兆円に倍増しても大丈夫なのかな。

なお，失敗ともいわれているICBM級ミサイルの軌道は高高度のロフテッド軌道に準ずるものだ。前回の4600kmとんだものの高度は1000kmだったので，失敗としてもその2倍に達している。先に作ったコードを用いて，加速パラメタを共通（0.0446 km/s^2）にしておけば，投射角度と燃焼時間を調整することで，成功したロフテッド軌道（最大高度6250 km，到達距離1090 km，到達時間 4000 秒）と，今回の（失敗かもしれない）ロフテッド軌道（最大高度 1980 km，到達距離 760 km，到達時間 不明）をそれぞれ再現することはできそうだった。

このとき，11月3日午前7時39分に発射された弾道ミサイルの飛行時間は27分前後=1650秒となるので，落下時間は午前8時6分ごろである。多段式だとすれば，今のモデルではうまく記述できていないのだが，第1段の落下時刻はこれより早かったと思われる。

(* g = 0.0098; R = 6350; τ = 86.2; p = 0.75; a = 0.0446; 

s = 86 Degree; T = 4000; *)

g = 0.0098; R = 6350; τ = 77; p = 0.65; 

a = 0.0446; s = 84 Degree; T = 1650;

fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ - t]

ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ - t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 - 

     g R^2/r[t]^2 + fr[t, τ], r[0] == R, r'[0] == 0, 

   h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r, 

   h}, {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, {t, 0, T}, 

 PlotRange -> {-4, 8}]

tyx[TT_] := {TT, f[TT] - R, NIntegrate[R d[t]/f[t]^2 , {t, 0, TT}]}

tyx[T]

ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}], 

  f[tt] - R}, {tt, 0, T}]

図：11月3日朝の弾道ミサイルのロフテッド軌道

P. S. 1  成功したロフテッド軌道（最大高度6250 km，到達距離1090 km，到達時間 4000 秒）の条件で，投射角度を45度にすると，到達時間 4850秒で到達距離 14300 kmとなる。北米大陸は射程内だ。

P. S. 2  アメリカの軍事偵察衛星からでは発射直後の赤外線を感知できるが，日本のレーダーでは，高度100 kmに達しないと感知できないという説明がサンデーモーニングであった。もしそうならば，上記の11月3日朝のロフテッド軌道が感知できたのは，発射後66秒後となる。

P. S. 3  仮に上記弾道ミサイルの第1段が，同じ条件で燃料噴射後60秒で放物運動に至る物体だとすれば（切り離された第1段とは正確には一致しないが），それは（最大高度 1000 km，到達距離 470 km，到達時間 1100 秒）で飛行することになる（高度100 km に達するのは 発射後 64秒後）。

［１］ミサイル防衛について（防衛省・自衛隊）

［２］弾道ミサイル防衛（平成20年3月，防衛省）

［３］国民保護ポータルサイト（内閣官房）

［４］北朝鮮 弾道ミサイル「火星17型が飛行に失敗か」（NHK）

［５］北朝鮮からの射程距離（地図蔵）

ドーナツ地球（３）

ドーナツ地球（２）からの続き 

12個の球で近似するくらいならば，いっそのことドーナツ型の内部に一様分布する無数の質点モデルをとったほうがよいかもしれない。$\bm{r}$の位置における単位質量に対する万有引力のポテンシャルは，$V(\bm{r}) = \sum_{i=1}^{N} \dfrac{G m_i }{\sqrt{(\bm{r}-\bm{r}_i )^2}}$である。$N$ 個の質点は，質量 $m_i = \frac{M}{N}$，位置ベクトル$\ \bm{r}_i\ $で与えられ，空間に一様分布している。

トーラス（ドーナツ）は，円環の半径が$R$で，断面が半径$a$の円形であるとする。$\ \bm{r}=(x,y,z)\ $とすると，$x=\cos\phi (R + r \cos\theta)$，$ y=\sin\phi (R + r \sin\theta)$，$z = r \sin\theta )$ と表わすことができる。ただし，$\phi$はトーラス中心に対し円環上で$x$軸から測った方位角，$\theta$はトーラス断面円の中心に対し$x-y$平面から測った仰角で，トーラス中心に遠い方から測ったものである。

トーラス断面円の中心からの距離が$r$であり，$(x-R \cos\phi)^2+(y-R \sin\phi)^2 +z^2 = r^2\ $である。そこで，トーラス面上の点は，$(x-R \cos\phi)^2+(y-R \sin\phi)^2 +z^2 = a^2\ $を，トーラス内部の点は不等式$\ (\sqrt{x^2+y^2}-R)^2 + z^2 < a^2\ $を満足する。

m = 1000000; R = 2.0; r = 1.0;

t = Table[{RandomReal[{-3, 3}], RandomReal[{-3, 3}], 

    RandomReal[{-3, 3}]}, m];

(* s = Select[t,(#[[1]]^2+#[[2]]^2+#[[3]]^2)<1 &]; *)

s = Select[t, ((Sqrt[#[[1]]^2 + #[[2]]^2] - R)^2 + #[[3]]^2) < 1 &]; 

smax = Length[s]

n[w_] := 1/Sqrt[w[[1]]^2 + w[[2]]^2 + w[[3]]^2 + 1/smax]

ListPointPlot3D[s, BoxRatios -> Automatic]

図１：トーラス内部に一様分布する質点

p = {4 - z, 0, 0};

d = Table[p - s[[i]], {i, 1, smax}];

w = 1/smax*Sum[n[d[[k]]], {k, 1, smax}];

v[z_] := w

g[1] = Plot[v[z], {z, 0, 8}]

p = {0, 0, z - 4};

d = Table[p - s[[i]], {i, 1, smax}];

w = 1/smax*Sum[n[d[[k]]], {k, 1, smax}];

v[z_] := w

g[2] = Plot[v[z], {z, 0, 8}]

図２：万有引力ポテンシャルの値（左x軸上，右z軸上）

不連続な質点モデルで計算しているため，たまたま観測点が質点と重なると大きなスパイクが現れる。この現象を緩和するため，分母にくる距離の項 n[d[[k]]] において，1/smaxという微小項を加えている。

ドーナツ地球（２）

ドーナツ地球（１）からの続き

前回のモデル式をMathematicaで書くと次のようになる。地球の半径を$R \rightarrow r$としており，重力加速度の大きさは，${\rm [g]}=\dfrac{GM}{R^2}$を単位として得られたものになる。変数$t$が，重力の大きさを検討している大円上の角度になる。

v = Table[{2 r Cos[Pi/6 i], 2 r Sin[Pi/6 i], 0}, {i, 1, 12}];

d = Table[{2 r + r Cos[t], 0, r Sin[t]} - v[[i]], {i, 1, 12}];

n[w_] := w/Sqrt[w[[1]]^2 + w[[2]]^2 + w[[3]]^2]

m[i_] := (1 + Mod[i + 1, 2])/2

f = Sum[n[d[[i]]]*m[i], {i, 1, 12}];

g[t_] := f /. r -> 1.0

Plot[g[t], {t, 0, Pi}]

Plot[{Sqrt[g[t][[1]]^2 + g[t][[3]]^2], 

    180/Pi*ArcTan[g[t][[3]]/g[t][[1]]]}, {t, 0, Pi}];

vp = Table[{{Cos[t], Sin[t]}, {-g[t][[1]], -g[t][[3]]}}, 

    {t, 0., Pi, Pi/30}];

g1 = ListVectorPlot[vp, AspectRatio -> 0.5];

g2 = Graphics[{White, Disk[{0, -0.05}, 0.95, {0, Pi}]}];

Show[g1, g2] 

 

図１：重力加速度の方向依存性（大きさ，x成分，z成分 vs 大円上の角度 t ）

図２：断面における重力加速度ベクトル場の様子（x = -2.0 がドーナツの中心）

ドーナツ地球の外側で8G，内側で2G（ただしドーナツの中心を向く）となった。ということは，断面図が円形の状態では安定な物質分布とならないわけか。

弾道ミサイルの軌道（３）

弾道ミサイルの軌道（２）からの続き

10月4日の弾道ミサイルは，本来グアムを標的とする射程5,000kmの火星12型相当の中距離弾道弾なので，日本を対象とした射程1,500kmのノドン改良型とは違う。が，面倒なので，前回求めた解の推進加速度を7割程度に落とし加速角度を30°にして，日本を狙ったときの到達時間が適当な値となる解があるかを確認してみた。

g = 0.0098; R = 6350; τ = 60; p = 0.75; a = 0.032; s =  30 Degree;

fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ-t]

ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ-t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ-t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 - 

     g R^2/r[t]^2 + fr[t, τ], r[0] == R, r'[0] == 0, 

   h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r, 

   h}, {t, 0, 360}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, 360}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, {t, 0, 299},  PlotRange -> {-4, 7}]

(f[t] - R) /. FindRoot[D[f[t], t] == 0, {t, 200}]

　86.6527

FindRoot[D[f[t], t] == 0, {t, 300}]

　　{t -> 199.993}

{d[τ]/f[τ], f[τ+1] - f[τ]}

　　{4.14703, 0.971161}

{Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2], 

 Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2]/.34}

　　{4.25923, 12.5271}

{NIntegrate[R d[t]/f[t]^2 , {t, 0, 200}], 

 NIntegrate[R d[t]/f[t]^2 , {t, 0, 360}]}

　　{651.9, 1305.61}

ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, T}],  f[T] - R}, {T, 0, 360}]

図：弾道ミサイルの軌道（横軸 水平距離 km，縦軸 高度 km）

60秒加速後の速度がマッハ 12.5となり，360秒（これは初期値として与えた）で1300 km離れた日本に到達する。最高高度は90 km弱となる。このあたりは推進加速度とその角度を調整すればなんとでもなる。いずれにせよ，カップ麺ができる3分以内に対応する必要があるのだけれど，Ｊアラートにインプットする情報分析や対応措置は時間内に間に合うのだろうか。

弾道ミサイルの軌道（２）

弾道ミサイルの軌道（１）からの続き

今回のモデルによるシミュレーションで，軌道を再現する加速時間は180秒=3分だったが，これを少し変えるだけで，軌道は大きく変わってしまう。つまり，ミサイル発射直後の軌道推定は困難であり，少なくとも初期加速が終る時点までの数分間は待つ必要がある。と思ったのだが，コロラド先生は燃焼時間が1分だといっていた。

そこで，各方向の加速度を0.069 km/s^2 (7G)，加速時間を60 sにしたところ，vt = 4.07 km/s，vr = 3.58 km/s，v0 = 5.42 km/s（マッハ15.9）が得られた。このときの，最高高度が1160 km，到達距離が4270 km である。これでよいかと思ったが，調べてみると，そもそも推進剤の質量が全質量に対して非常に大きな割合を占めることがわかった。やり直し。

質量が変化する場合の極座標の運動方程式では，$m \ddot{\bm{r}}\ $の項に$\ \dot{m} \dot{\bm{r}}\ $が加わることになる。これを極座標にすれば，$\dot{m} (\dot{r} \bm{e}_r + r \dot{\theta} \bm{e}_\theta)\ $である。そこで，運動方程式の各成分は次の通り。

$\ddot{r} - \dfrac{h^2}{r^3} + \dfrac{\dot{m}}{m} \dot{r} = \frac{1}{m} F_r = -g \dfrac{R^2}{r^2} + \alpha H(\tau- t)$

$\dfrac{1}{r }(\dot{h} +  \dfrac{\dot{m}}{m} h) = \frac{1}{m} F_\theta = 0 + \beta H(\tau-t) $

なお，$H(x) = 1 (x>0) ; =0 (x<0)$はヘヴィサイドの階段関数である。$\alpha, \beta$は，加速開始時刻 $t=0$から加速終了時刻 $t=\tau$ まで動径方向と角度方向に加わる加速度を表わす。

出発前の弾道ミサイルの全質量を$m_0$，全質量に対する推進剤の割合を $p$とすると，加速中（$0 \le t \le \tau$）の単位時間当たりの噴射質量は，$\dfrac{p m_0}{\tau}$となる。そこで，時刻 $t$における弾道ミサイルの質量は，$m(t) = m_0 - \dfrac{p m_0}{\tau} t  = m_0 (1 - p\dfrac{t}{\tau})$，$\dot{m}(t) = - m_0 \dfrac{p}{\tau}$。

したがって，$\dfrac{\dot{m}}{m} = - \dfrac{p}{\tau - p t} \quad (0 \le t \le \tau)$ であり，$t>\tau$ では $\ \dfrac{\dot{m}}{m} = 0\ $となる。

なお，最高高度は，$\dot{r}(t)=0$となる時刻 $t_p$に対応する$r(t_p)$で与えられ，到達距離は，$\int_0^T \dfrac{R h(t)}{r(t)^2}dt$で得られる。ただし，$T$は到達時間である。

準備ができたので再計算してみる。全重量に対する推進剤の比率は，p=0.75（3/4）と仮定した。加速時間 60s，加速度 0.0446 km/s^2（4.6 G）で到達時間1320 sが再現される。最高高度は，t=682sのとき 1060 km。加速終了時の速度は，vt = 4.69 km/s，vr = 3.32 km/s，v0 = 5.75 km/s （マッハ16.9）である。また，到達距離は4930 kmであり，6-7%の誤差で報告値に一致した。

g = 0.0098; R = 6350; τ = 60; p = 0.75; a= 0.0446; s = Pi/4

fr[t_,τ_] := a * Sin[s] * HeavisideTheta[τ - t]

ft[t_,τ_] := a * Cos[s] * r[t] * HeavisideTheta[τ- t]

fm[t_,τ_] := -p / (τ- p * t) * HeavisideTheta[τ - t]

sol = 

 NDSolve[{r''[t] == -fm[t,τ] * r'[t] + h[t]^2/r[t]^3 - 

   g R^2/r[t]^2 + fr[t,τ], r[0] == R, r'[0] == 0, 

   h'[t] == -fm[t,τ] * h[t] + ft[t,τ], h[0] == 0},

   {r, h}, {t, 0, 1320}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, 1320}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, {t, 0, 1319}, PlotRange -> {-4, 5}]

(f[t] - R) /. FindRoot[D[f[t], t] == 0, {t, 660}]

{d[τ]/f[τ], f[τ+1] - f[τ]}

{Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2], 

 Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2]/.34}

{NIntegrate[R d[t]/f[t]^2 , {t, 0, τ}], 

 NIntegrate[R d[t]/f[t]^2 , {t, 0, 1320}]} 

図１：軌道半径の時間変化

図２：速度vr，速度vtの大圏射影，速度vt の時間変化

推進時間 60s ，加速度 4.55G，加速角 45° で，高度 1059 km，到達距離 4933 km

推進時間 90s ，加速度 3.22G，加速角 49° で，高度 1038 km，到達距離 4813 km

推進時間120s ，加速度 2.57G，加速角 52.5° で，高度 1015 km，到達距離 4697 km

というわけで，報告値の 1-2%の範囲に収めることも可能な定性的モデルができた。

（付録）最後のパラメタでは，津軽海峡に到達する発射後420 s，距離1368km，高度811kmの速度は3.9 km/s であり，30秒ほどで津軽海峡（幅130 km）を通過することになる。

弾道ミサイルの軌道（３）に続く

弾道ミサイルの軌道（１）

Ｊアラートからの続き

このたびの弾道ミサイルが，大圏コースで4600kmを飛行したということは，角度にして40度強だ。また，高度1000kmというのは地球半径の15%にあたる。さすがに地表面を平面として一様重力場で考えるというのではちょっとマズイ気がする。

そこで，地球の中心を原点とする極座標系での運動方程式を考える。運動する物体の座標を$( r(t), \theta(t) )$とする。運動方程式は，$ m ( \ddot{r} - r \dot{\theta} ) = F_r, \quad \dfrac{m}{r} \dfrac{d}{dt} (r^2 \theta) = F_\theta $ となる。ここで面積速度の2倍を$h(t) = r^2 \dot{\theta}$と定義すると，$ m ( \ddot{r} - \dfrac{h^2}{r^3} ) = F_r, \quad \dfrac{m}{r} \dot{h} = F_\theta $ となる。

(1) 万有引力だけが働く場合：　$F_r = \dfrac{GMm}{r^2} = mg \dfrac{R^2}{r^2}, \quad F_\theta =0 $となる。与えられた条件は，到達距離 L=4600km，到達時間 T=1320 s，最高高度 H=1000 km，平均水平速度 vt =3.55 km である。また，地表重力加速度 g=9.8 m/s^2，地球半径 R=6350 kmとする。Mathematicaのコードで調整するパラメータは，鉛直方向の初速度だけである。到達距離・時間の条件を満たすものとして vr = 4.15  km/s が得られる。このときの速度はv0 = √(vt^2+vr^2) = 5.46 km/s（マッハ16）となる。しかし，最高高度が，1300 kmとなってうまく合わない。

g = 0.0098; R = 6350; vr = 4.15; vt = 3.55; h = R vt

sol = NDSolve[{r''[t] == h^2/r[t]^3 - g R^2/r[t]^2, 

      r[0] == R, r'[0] == vr}, r, {t, 0, 1320}]

f[t_] := r[t] /. sol[[1]]; f[660] - R

Plot[{6350, f[t]}, {t, 0, 1320}]

図１：初速度のみ与えたモデル（横軸: t ，縦軸: 軌道半径 r ）

 (2) 初期加速度が一定時間働く場合：　加速度の値（簡単のため動径方向と角度方向は等しいと仮定），加速時間の2つをパラメタとする。先ほどのように到達距離・時間の条件を満たすようにパラメタを探すと，加速度 0.025 km/s^2 （2.5G）と加速時間 180 s の値が得られた。このときの vt = d[τ]/f[τ] = 4.38 km/s，vr = f[τ+1] - f[τ] = 2.96 km/s，v0 = 5.29 km/s（マッハ15.5）となる。また，最高高度は1020 kmとなり，この場合は全体として辻褄があうことになる。

g = 0.0098; R = 6350; τ= 180;

fr[t_,τ_] := 0.025 * HeavisideTheta[τ- t]

ft[t_,τ_] := 0.025 * r[t] * HeavisideTheta[τ- t]

sol = NDSolve[{r''[t] == h[t]^2/r[t]^3 - g R^2/r[t]^2 + fr[t,τ], r[0] == R, r'[0] == 0, h'[t] == ft[t,τ], h[0] == 0},

{r, h}, {t, 0, 1320}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, 1320}]

図２：初速度0から一定の加速をする場合の軌道半径（横軸: t ，縦軸: 軌道半径 r ）

なお，加速終了時の射影水平速度は，vt = R d[180]/f[180]^2 = 4.16 km/s である。そこで180秒のあいだに進む水平距離は，vt τ/2 = 370kmとなる。残りは，1030km/4.16km/s = 248 s なので，計 428秒で1400 km（津軽海峡上空）に達する。

図の印象で騙されていたが，到達距離4600 kmの確認が済んでいなかった。解けているのは 角速度 $\dot{\theta} = \dfrac{h(t)}{r(t)^2}$なので，これを積分した $R \theta(T)$ が必要なのだ。下図より角速度の平均値が 0.00059だと仮定すると，370 km + 0.00059*R*1140 s = 370 + 4270 = 4640 kmとなる。1%の誤差でOKだった。

図３：初速度0から一定の加速をする場合の角速度（横軸: t ，縦軸: 角速度$\dot{\theta}$  ）

最高高度は，t=679s で1050 kmとなった。加速終了時の速度は，vt = d[τ] / f[τ] = 4.69 km/s，vr = f[τ+1] - f[τ] = 3.31 km/s，v0 = √(vt^2+vr^2) = 5.74 km/s （マッハ16.9）であり，ほぼ報道結果が再現された（P. S. と思ったが・・・）。

弾道ミサイルの軌道（２）に続く。

プリズムによる光の分散

授業の課題で次のような問題を出した。「三角プリズムに太陽光を入射させたとき，プリズムから2mの距離で光のスペクトルはどのくらいの幅に広がるか」。 もちろん，入射角やプリズムの種類で結果は変わるのだけれど，そのあたりも含めて考察せよというものだ。

自分でも考えてみる。可視光におけるプリズムの屈折率$n$の波長依存性は，[1]の上越教育大学の紀要を参考にすると，400nm（紫）から700nm（赤）の範囲で，フリントガラス：$n$=1.65〜1.61，クラウンガラス：$n$=1.535〜1.515くらいになる。なお，紫（380-430 nm），青（430-490 nm），緑（490-550 nm），黄（550-590 nm），橙（590-640 nm），赤（640-770 nm）である。

三角プリズムの底面に垂直に$x$軸をとり，図のように入射角，屈折角，出射角などを定義すると，出射する光線の$x$軸に対する傾き$a$が求まる。スリットを通してプリズムに入射する光の場合，再びプリズムから空気中にでる出射点の位置は分散で数mm程度ズレるかもしれないが，2m先のスペクトルの位置は基本的に光線の傾き$a$で決まる。

図：プリズムによる光の分散

図より，$a = \tan \bigl[\dfrac{\pi}{3} - \sin^{-1}\bigr\{ n \sin \bigl(\dfrac{\pi}{3} - \sin^{-1}\dfrac{\sin \alpha}{n}\bigr)\bigl\}\bigr]$である。これを使って，Mathematicaでスペクトルの分散幅を計算する。ただし，入射角は0.9 rad とした。

f[n_, t_] := Tan[Pi/3 - ArcSin[n Sin[Pi/3 - ArcSin[Sin[t]/n]]]];
g1 = Plot[Table[f[1.605 + 0.01 k, t], {k, 1, 5}], {t, Pi/4, Pi/3}, PlotRange -> {-0.2, 0.3}];
g2 = Plot[Table[f[1.51 + 0.005 k, t], {k, 1, 4}], {t, Pi/4, Pi/3}, PlotRange -> {-0.2, 0.3}];
Show[g1, g2]
200*Table[f[1.605 + 0.01 k, 0.9], {k, 1, 5}]
{13.1215, 9.49956, 5.78634, 1.96916, -1.96692}
200*Table[f[1.51 + 0.005 k, 0.9], {k, 1, 4}]
{46.1667, 44.6014, 43.0304, 41.4533}

図：出射方向aの入射角依存性と波長依存性

結局，入射角が 51.6度（0.9ラジアン）として，プリズムから2m先の位置では，フリントガラスで15cm程度，クラウンガラスで5cm程度に可視光スペクトルが広がる。

［１］屈折率の波長依存性の簡易測定（西山保子・上田淳一）