ドーナツ地球（３）

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

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（ただしドーナツの中心を向く）となった。ということは，断面図が円形の状態では安定な物質分布とならないわけか。

ドーナツ地球（１）

ドーナツ地球というのが目に入った。次のようなイメージである。 

図：ドーナツ地球のイメージ（どこぞのYouTubeから引用）

詳しい話が，Gigazineの2014年の「もしもドーナツ状の地球が存在したらそこはどんな世界なのか」にあった。その元記事は，Anders Sandbert のdistributed brainにある Torus-Earth である。重力の大きさ（等ポテンシャル面）の分布図があったけれど，いまいちピンと来ない。

自分で計算してみようとしたけれど，どうやって積分したらよいかで行き詰まる。そもそも，これを太さのない線密度だけのリングだとしても，楕円積分が出てきてどうしましょう状態になりそうだ。円環状の電荷分布が作る電場や電位の演習問題と等価だが，円環の中心を通る対称軸上の値を求める問題解答しかみあたらない。

そこで，このドーナツを球に分割して，その和で重力を近似してみることにする。ドーナツの穴が1地球分空いていて，まわりに6地球を配置すると収まりがよい。取り巻く地球間の6つの隙間は，それぞれ約1/2地球分の体積となる。そこで，その分の小地球を6地球の間に重ねて配置し，重力はこれらの12個の大小の地球を各中心に質量が集中した質点だと近似して計算する。

6地球の1つ（中心の座標が$\bm{R_N}$）を取り出して，トーラスの円環軸に垂直な断面の大円をとる。大円上の点Pの座標を$\bm{r}=  \bm{R_N} + (R \cos \theta, R \sin \theta ) \quad 0 \le \theta \le \pi $とする。P上の質量$m$の物体に働く重力加速度は $\bm{F}/m = \sum_{i=1}^{12} G M_i \dfrac{\bm{r}-\bm{R}_i }{|\bm{r}-\bm{R}_i |^2}$となる。

図：ドーナツ地球を球体の集合で近似する

ウィーンの変位則

非常勤の授業準備をしていると，自分の理解の穴がいたるところに待ち受けている。ウィーンの変位則は，プランクによる量子の発見にいたる重要なマイルストーンであるが，ここにも落とし穴をひとつ発見。

ウィーンの変位則は，黒体放射のエネルギーがピークとなる電磁波（光）の波長$\ \lambda_m \ $と黒体の温度$\ T \ $との間に $\lambda_m = \dfrac{b}{T}$という反比例の関係があるというものだ。例えば，高温の物体からの電磁波は短波長側にシフトすることになり，溶鉱炉中の鉄の色は温度が上昇するとともに，赤から黄に変化する。温度の低い星は赤く，温度の高い星は青いのもこのためだ。

古典電磁気学と熱力学だけから導いた黒体放射のスペクトルは，レイリー・ジーンズの法則を満たすが，これはウィーンの変位則を満足しない。黒体放射スペクトルの短波長側によく当てはまり，ウィーンの変位則が成り立つのはヴィーンの放射法則だ（慣用的に同じ人名をウィーンとヴィーンで使い分けるのはなんとかならんのか）。これら両者を統合して電磁波の全波長領域のスペクトルを再現したのがプランクの法則で，量子の発見につながった。

授業の演習課題で，プランクの法則からウィーンの変位則を導く問題をつくったところ，面倒なことに，指数関数を含んだ非線形の方程式がでてきた。そういえば，ランベルトのW関数というのは見覚えがあった。しかたがないので，ヴィーンの法則に置き換えたが，こんどは，振動数$\ \nu \ $の関数にするか，波長$\ \lambda \ $の関数にするかで，ピークは変わってくる。$\nu_{\rm max} \cdot \lambda_{\rm max}=c\ $が満たされないのを見落としていた。

そういえば，昔，岩波の「現代物理学の基礎」で量子力学の勉強を始めようとしたときに，黒体放射の式がよくわからずにいきなり挫折したことを思い出した。$d\nu = -\frac{c}{\lambda^2}\ d\lambda \ $を忘れてはいけない。

(1) ウィーンの放射法則

$u(\nu) = \dfrac{8\pi h}{c^3}\nu^3 e^{-h\nu / kT}$より，$\dfrac{du(\nu)}{d\nu} =0 \rightarrow \nu = \frac{3kT}{h}$

$w(\lambda) = \dfrac{8\pi hc}{\lambda^5} e^{-h / kT \lambda}$より，$\dfrac{dw(\lambda)}{d\lambda} =0 \rightarrow \lambda = \frac{hc}{5kT}$

(2) プランクの法則

$u(\nu) = \dfrac{8\pi h}{c^3}\nu^3 \dfrac{1}{e^{h\nu / kT}-1}$より，$\dfrac{du(\nu)}{d\nu} =0 \rightarrow (3-x) e^x = 3 \quad (x=\frac{h\nu}{kT}=2.821= 3. + {\rm ProductLog[-3 E\hat{}-3]})$

$w(\lambda) = \dfrac{8\pi hc}{\lambda^5} \dfrac{1}{e^{h / kT \lambda}-1}$より，$\dfrac{dw(\lambda)}{d\lambda} =0 \rightarrow (5-x) e^x=5 \quad (x=\frac{h}{kT\ \lambda}=4.965=5. + {\rm ProductLog[-5 E\hat{}-5]})$

図：ランベルトW関数の例，$(n-x)e^x=n$の逆数の図（$n=3,5$）

［１］物理授業におけるランベルトのＷ関数の活用法（伊東正人）

宗教法人法

統一教会（世界平和統一家庭連合）の宗教法人としての解散命令は，宗教法人法に規定されている。宗教法人法は，宗教団体に宗教法人という種類の法人格を与えることが目的となっている。裁判所が宗教法人に対する解散命令を出せば，その宗教法人は再び法人格を持たない宗教団体に回帰する。以下，関連事項の抜粋を引用する。

第一章　総則

（この法律の目的）

第一条　この法律は、宗教団体が、礼拝の施設その他の財産を所有し、これを維持運用し、その他その目的達成のための業務及び事業を運営することに資するため、宗教団体に法律上の能力を与えることを目的とする。

２　憲法で保障された信教の自由は、すべての国政において尊重されなければならない。従つて、この法律のいかなる規定も、個人、集団又は団体が、その保障された自由に基いて、教義をひろめ、儀式行事を行い、その他宗教上の行為を行うことを制限するものと解釈してはならない。

（宗教団体の定義）

第二条　この法律において「宗教団体」とは、宗教の教義をひろめ、儀式行事を行い、及び信者を教化育成することを主たる目的とする左に掲げる団体をいう。

一　礼拝の施設を備える神社、寺院、教会、修道院その他これらに類する団体

二　前号に掲げる団体を包括する教派、宗派、教団、教会、修道会、司教区その他これらに類する団体

（法人格）

第四条　宗教団体は、この法律により、法人となることができる。

２　この法律において「宗教法人」とは、この法律により法人となつた宗教団体をいう。

（公益事業その他の事業）

第六条　宗教法人は、公益事業を行うことができる。

２　宗教法人は、その目的に反しない限り、公益事業以外の事業を行うことができる。この場合において、収益を生じたときは、これを当該宗教法人、当該宗教法人を包括する宗教団体又は当該宗教法人が援助する宗教法人若しくは公益事業のために使用しなければならない。

第六章　解散

（解散の事由）

第四十三条　宗教法人は、任意に解散することができる。

２　宗教法人は、前項の場合のほか、次に掲げる事由によつて解散する。 

五　第八十一条第一項の規定による裁判所の解散命令

第八章　宗教法人審議会

（設置及び所掌事務）

第七十一条　文部科学省に宗教法人審議会を置く。

２　宗教法人審議会は、この法律の規定によりその権限に属させられた事項を処理する。

３　宗教法人審議会は、所轄庁がこの法律の規定による権限（前項に規定する事項に係るものに限る。）を行使するに際し留意すべき事項に関し、文部科学大臣に意見を述べることができる。

４　宗教法人審議会は、宗教団体における信仰、規律、慣習等宗教上の事項について、いかなる形においても調停し、又は干渉してはならない。

第九章　補則

（報告及び質問）

第七十八条の二　所轄庁は、宗教法人について次の各号の一に該当する疑いがあると認めるときは、この法律を施行するため必要な限度において、当該宗教法人の業務又は事業の管理運営に関する事項に関し、当該宗教法人に対し報告を求め、又は当該職員に当該宗教法人の代表役員、責任役員その他の関係者に対し質問させることができる。この場合において、当該職員が質問するために当該宗教法人の施設に立ち入るときは、当該宗教法人の代表役員、責任役員その他の関係者の同意を得なければならない。

一　当該宗教法人が行う公益事業以外の事業について第六条第二項の規定に違反する事実があること。

・・・

三　当該宗教法人について第八十一条第一項第一号から第四号までの一に該当する事由があること。

２　前項の規定により報告を求め、又は当該職員に質問させようとする場合においては、所轄庁は、当該所轄庁が文部科学大臣であるときはあらかじめ宗教法人審議会に諮問してその意見を聞き、当該所轄庁が都道府県知事であるときはあらかじめ文部科学大臣を通じて宗教法人審議会の意見を聞かなければならない。

３　前項の場合においては、文部科学大臣は、報告を求め、又は当該職員に質問させる事項及び理由を宗教法人審議会に示して、その意見を聞かなければならない。

４　所轄庁は、第一項の規定により報告を求め、又は当該職員に質問させる場合には、宗教法人の宗教上の特性及び慣習を尊重し、信教の自由を妨げることがないように特に留意しなければならない。

５　第一項の規定により質問する当該職員は、その身分を示す証明書を携帯し、宗教法人の代表役員、責任役員その他の関係者に提示しなければならない。

６　第一項の規定による権限は、犯罪捜査のために認められたものと解釈してはならない。

（解散命令）

第八十一条　裁判所は、宗教法人について左の各号の一に該当する事由があると認めたときは、所轄庁、利害関係人若しくは検察官の請求により又は職権で、その解散を命ずることができる。

一　法令に違反して、著しく公共の福祉を害すると明らかに認められる行為をしたこと。

二　第二条に規定する宗教団体の目的を著しく逸脱した行為をしたこと又は一年以上にわたつてその目的のための行為をしないこと。

・・・

３　第一項の規定による裁判には、理由を付さなければならない。

４　裁判所は、第一項の規定による裁判をするときは、あらかじめ当該宗教法人の代表役員若しくはその代務者又は当該宗教法人の代理人及び同項の規定による裁判の請求をした所轄庁、利害関係人又は検察官の陳述を求めなければならない。

５　第一項の規定による裁判に対しては、当該宗教法人又は同項の規定による裁判の請求をした所轄庁、利害関係人若しくは検察官に限り、即時抗告をすることができる。この場合において、当該即時抗告が当該宗教法人の解散を命ずる裁判に対するものであるときは、執行停止の効力を有する。

６　裁判所は、第一項の規定による裁判が確定したときは、その解散した宗教法人の主たる事務所の所在地の登記所に解散の登記の嘱託をしなければならない。

［１］宗教法人解散命令に対する抗告棄却決定に対する特別抗告

　　　（平成8年1月30日，最高裁判所第一小法廷，決定，棄却）

［２］法人の解散命令等に関する規定（消費者庁）

カーリングの原理（５）

カーリングの原理（４）からの続き

村田次郎さんの結論を整理すると次のように理解できる。

図：カーリングストーン（リングモデル）に働く力

広義の摩擦力は，速度方向と逆向きに働く動摩擦力と，衝突/固着点の回りに働く旋回力から成り立つ。その向きは図のように表わされる。衝突/固着の瞬間の撃力の寄与を除く旋回力はストーンの中心から衝突/固着点までの位置ベクトルと同じ向きになり，動摩擦力は進行方向の速度ベクトル$\bm{u}$と回転方向の速度ベクトル$\bm{w}$を合成したものと逆方向になる。

これらの力をリング上で平均すれば対称性によっていずれの場合もその合力は0になる。ところで，y軸の正負の領域間でアイスシートに対するストーンの速度に非対称性がある。上半面の速度$v_{-}$は下半面の速度$v_{+}$より小さくなる。旋回力と動摩擦力はともに速度が遅い方が大きいとする。これによって，旋回力では，ストーンの自転方向にカールする力の成分が残るが，同摩擦力では，打ち消し合いが依然として生ずることになる。

$y$軸方向の力を打ち消しあうペアが，同じ速度領域なのか異なった速度領域なのかがポイントだということになる。

（付）旋回力の図で，中心に示した矢印は衝突/固着時の撃力による重心の運動量変化に対応する力を現したものである。この項も，速度の非対称性の有無にかかわらず，キャンセルする項になっている。

付録（TikZでの変数や反復の方法について少しだけ勉強した）：

\begin{tikzpicture}
\tikzstyle{every node}=[font = \large];
\draw[step=1.0, dotted] (-7.5,-3.5) grid (7.5,3.5);
\draw[->] (-7.5,0) -- (-0.5,0) node[below left]{$x$};
\draw[->] ( 0.5,0) -- ( 7.5,0) node[below left]{$x$};
\draw[->] (-4,-2.5) -- (-4,2.5) node[left]{$y$};
\draw[->] ( 4,-2.5) -- ( 4,2.5) node[left]{$y$};
\draw (-4,0) circle(2) circle(1pt) node[below left]{${\rm O}_1$};
\draw ( 4,0) circle(2) circle(1pt) node[below left]{${\rm O}_2$};
\draw[<- br="" gray="" thick="" ultra=""> \node at (-3,0) [above right]{$u$};
\node at (-4,-3) [red, above]{旋回力の方向};
\draw[->, thick, gray] (-4,-0.5) arc [start angle=-90, end angle=90, radius=0.5] node[left]{$w$};
\draw[<- 4="" 5="" br="" gray="" thick="" ultra=""> \node at (5, 0) [above right]{$u$};
\node at (4,-3) [blue, above]{動摩擦力の方向};
\draw[->, thick, gray] ( 4,-0.5) arc [start angle=-90, end angle=90, radius=0.5] node[left]{$w$};
\draw[dashed, blue] (6,1.4)--(2,1.4);
\draw[dashed, blue] (6,-1.4)--(2,-1.4);
\draw[dashed, red] (-5.4,-2)--(-5.4,2);
\draw[dashed, red] (-2.6,-2)--(-2.6,2);
\node[purple] at (-3.3, 0.5) {2,4};
\node[purple] at (-3.3,-0.5) {1,3};

\foreach \t/\s/\r/\d\q in {1/10/-1/0/{-}, 3/5/0/0/{-}, 5/-5/0/0.4/{+}, 7/-10/-1/-0.4/{+}}
{
\tikzmath{
integer \n, \m;
\n = (\t+1)/2; \m=abs(4-\t)-2;
\x1 = -4+2*cos(\t*45); \y1 = 2*sin(\t*45);
\u1 = -4+2.5*cos(\t*45); \v1 = 2.5*sin(\t*45);
\w1 = -4+1.5*cos(\t*45); \z1 = 1.5*sin(\t*45);
\x2 = 4+2*cos(\t*45); \y2 = 2*sin(\t*45);
\u2 = 4+(2.5+\r)*cos(\t*45+\s); \v2 = (2.5+\r)*sin(\t*45+\s);
\w2 = 4+1.5*cos(\t*45); \z2 = 1.5*sin(\t*45);
}
\filldraw[red] (\x1,\y1) circle(2pt) (\w1,\z1) node{$\n\ v_\q$};
\draw[->, thick, red] (\x1,\y1)--(\u1,\v1);
\draw[->, thick, purple] (-4,0)--({\m*(\u1-\x1)-4},{\m*(\v1-\y1)});
\filldraw[blue] (\x2,\y2) circle(2pt) (\w2,\z2) node{$\n\ v_\q$};
\draw[->, thick, blue] (\x2,\y2)--(\u2,{\v2+\d});
}
\end{tikzpicture}

カーリングの原理（４）

カーリングの原理（３）からの続き 

問題を現象論的な微分方程式で扱えれば簡単なのだが，動摩擦力だけではカールがうまく説明できない。もちろん，進行方向前後で動摩擦係数が異なり，進行方向後部の摩擦が大きいとすればよいのだけれど，物理的な理由づけが難しい。

村田さんの精密実験の結論は「摩擦＝微小な衝突の重ね合わせとして，離散的な摩擦支点中心の旋回が，左右非対称な頻度で起こる」というもので「従来の左右非対称摩擦説を支持するものであった」とまとめられている。

実験データの分析や結論は問題ないが，最後のまとめ方がちょっと納得できない。というのも普通の動摩擦力であれば，物体の速度と逆方向に働くのだけれど，摩擦支点中心の旋回力＝動摩擦力（旋回）の向きは通常の動摩擦力（並進）とは異なる。もちろん広義の摩擦ということではくくることはできるだろうが，現象論的には動摩擦力（並進）および動摩擦力（旋回）として区別して取り扱う方がよいと考えられる。

そこで，$x$軸方向に速度$u$で進む質量$M$，半径$R$の環状ストーンモデルを考える。このリングの進行方向から$\theta$の部分とアイスシートの間に衝突/固着点 P が生じ，これを中心に重心が角速度$\Omega$で旋回運動を始めた場合の旋回力の大きさと方向を求めてみる。

図：衝突/固着点Pの回りの旋回力を導くための概念図

衝突/固着点Pで働く抗力はこの点の回りの角運動量に寄与しないので，点Pの回りの衝突/固着前後の角運動量が保存される。ただしリングの重心の回りの自転角速度は小さいのでその寄与は無視した。保存される角運動量は，$L= M u R \sin \theta = I_{\rm P} \Omega  = 2MR^2 \Omega$である。したがって，$\Omega = \frac{u \sin \theta}{2R}$となる。ここで，$I_{\rm P}= 2MR^2 \ $はP点の回りのリングの慣性能率である。

衝突/固着後に重心は運動量$\ \bm{p}=M\bm{V}=M R\Omega(\sin\theta, -\cos \theta)\ $で動き出す。なお，運動量の大きさは$\ p=MR\Omega\ $である。微小時間$\Delta t$の後，重心は$\ \Omega \Delta t\ $だけ反時計回りに回転し，これにともなって，運動量ベクトルは$\ \Delta {\bm p}=p \Omega \Delta t (\cos\theta, \sin\theta) = M R \Omega ^2  \Delta t (\cos\theta, \sin\theta)\ $だけ変化する。

したがって，衝突/固着時の旋回力は $\bm{F}=\dfrac{\Delta {\bm p}}{\Delta t} = M R \Omega^2 (\cos\theta, \sin \theta) = M \frac{u^2 \sin^2\theta}{4R} (\cos\theta, \sin\theta)\ $である。これについても，幾何学的な対称性によって，このまま一様に角度平均をとればその寄与は0になってしまい，ストーンのカールの原因にはなれない。

そこで，衝突/固着の確率がアイスシートに対する衝突/固着点Pの速度$v(\theta)$に逆比例すると考え，無次元化した$\frac{u}{v(\theta)}$をかける。つまり，旋回力$\ \bm{F}=M \frac{u^3 \sin^2\theta}{4R v(t)} (\cos\theta, \sin\theta)\ $を導入すれば，リングの回転方向にカールすることが説明できることになる。

前回と同様に，Mathematicaで試してみると次のようになる。

mg\[Mu] = 2; mcl = 1; ux = 1; uy = 0; w = 0.1; r = 0.1; q = 0.944;
v[t_] := Sqrt[w^2 - 2*w *(ux*Sin[t] - uy*Cos[t]) + ux^2 + uy^2]
k[t_] := mg\[Mu]/v[t]^1.5 *(1 - q*Cos[t])
F[t_] := -k[t]*{ ux - w * Sin[t], uy + w*Cos[t]}
G[t_] := mcl ux^3/(4 r*v[t])*Sin[t]^3
T[t_] := -r * k[t]*(w + uy*Cos[t] - ux* Sin[t])
Plot[{v[t], F[t], G[t]}, {t, 0, 2 Pi}, PlotRange -> {-3, 3}]
NIntegrate[{F[t], G[t]}, {t, 0, 2 Pi}]
{{-12.5427, 0.592575}, 0.592504}

図：旋回力と前後非対称の動摩擦力の比較

（注）旋回力と同じ効果を与える動摩擦力の前後非対称係数 q を入れた場合

［１］The Physics of Sliding Cylinders and Curling Rocks（A. R. Penner）

カーリングの原理（３）

カーリングの原理（２）からの続き 

経験的にわかっていることは，(1) ストーンは反時計（左）回りで左にカールし，時計（右）回りで右にカールする，(2) カールする度合いは，回転の大小で余り変わらない，(3) ストーンの並進運動が止まるときに回転も止まる，であるらしい。

理論的には，大きく分けて２つの立場がある。(1) 左右の摩擦係数の違いによる：並進速度と回転速度が加わる部分はアイスシートに対する速度が大きく摩擦係数が小さいが，並進速度から回転速度が引かれる部分はアイスシートに対する速度が小さく摩擦係数が大きい。これがカールの原因となる。(2) 前後の摩擦係数の違いによる：(1) では，円環モデルをとるかぎり横方向の摩擦力は対称性から打ち消しあうので，カールするための原動力がでてこない。ところが，進行方向前後の摩擦係数を変えるとき，後方の摩擦係数が大きければ経験的にしられている回転方向とカール方向の関係が再現できる。

円環モデルの進行方向に対して角度$\theta$の場所にある質量要素のアイスシートに対する速度$v$は，$v = \sqrt{ u_x^2 + u_y^2 + w^2 - 2w (u_y \cos \theta - u_x  \sin \theta ) }\ $ である。つまり，$u_y=0$ならば，$\sin \theta$のみの関数になる。一方，円環モデルの進行方向に対して垂直な方向に現れる摩擦力は，$F_y= - \frac{\mu}{v^p} M g w \ \cos \theta$なので，これを$\theta=0 \sim 2\pi$で積分すれば，必ず0になってしまう。一方，この$\mu$に前後非対称があれば，積分値は0ではなくなる。 

この事情をMathematicaで計算してみると次のようになる。

mg\[Mu] = 2; ux = 1; uy = 0.00; w = 0.1; r = 0.1; q=0.0;
v[t_] := Sqrt[w^2 - 2*w*(ux*Sin[t] - uy*Cos[t]) + ux^2 + uy^2]
k[t_] := mg\[Mu]/v[t]^1.5 *(1 - q*Cos[t]) 
F[t_] := -k[t]*{ux - w*Sin[t], uy + w*Cos[t]}
T[t_] := -r * k[t] * (w + uy*Cos[t] - ux*Sin[t])
Plot[{v[t], F[t], T[t]}, {t, 0, 2 Pi}, PlotRange -> {-3, 3}]
NIntegrate[{F[t], T[t]}, {t, 0, 2 Pi}] 

{{-12.7325, 6.27728*10^-7}, -0.221247}

図：角度の関数としての{速度，摩擦力_x，摩擦力_y，トルク}

・速度方向の逆を向く摩擦力の係数 k[t] にある Cos[t] の項は前後非対称を示す。

・進行方向の速度にy成分u_yを持たせると，その逆方向の摩擦力が働く。

　　　　　　　　付録：HTMLにおける表の練習（四字熟語から）

　１．羊頭狗肉　２．鶏鳴狗盗　３．鶏口牛後　４．沈魚落雁　５．窮鼠噛猫
　６．烏兎匆匆　７．籠鳥檻猿　８．鶏群一鶴　９．鱸膾蓴羮　10．鯨飲馬食

カーリングの原理（２）

カーリングの原理（１）からの続き 

富山大学の対馬勝年先生が，2013年に「氷雪のトライポロジー」というまとまったレポートを出していた。 ただ，カーリングがカールする根拠としてあげた左右に錘のついた棒のモデルや角度方向の摩擦係数の議論は理解できなかった。

そこで，カーリングのカールに関するこれまでの議論を少し復習してみる。

カーリングのストーンの質量は，$M=20{\rm \ kg}\ $であり，氷上に接するのはランニングバンドとよばれる狭い円環部分である。その半径は$\ R=0.1 {\rm \ m}\ $だ。そこで，ストーンを円環によってモデル化すると，中心の周りの慣性モーメントは，$I = M R^2 = 0.2 {\rm \ kg m^2}\ $となる。ストーンの初速度は，$u_0 = 2 {\rm \ m/s}$，回転を与えた場合の初角速度は，$\omega_0 = 1 {\rm \ rad/s}\ $とする。つまり回転方向の初速度は，$w_0=R \omega = 0.1 {\rm \ m/s}\ $となる。

摩擦のメカニズムを，動摩擦力$\bm{F}\ $ によって現象論的に表現すると，その力は，ストーンと氷の接点の相対速度ベクトル$\bm{v}\ $とは逆向きで，大きさが垂直抗力に比例するものとなる。その比例定数が動摩擦係数 $\ \mu\ $になり，必要ならばこれに速度依存性を導入する。つまり，$\bm{F} =- \mu(v) \ Mg \ \hat{\bm{v}} = -\dfrac{\mu(v)}{v}\ Mg\  \bm{v} \rightarrow -\dfrac{\mu(v)}{v^p}\ Mg\  \bm{v}$。

なお，動摩擦係数の値を$\mu = 0.01 \ $のオーダーとすれば，動摩擦力の大きさは，$F = \mu \ M g = 2 {\rm \ N \ } $となる。ストーンの初期運動エネルギー$K_0$が，停止するまでに摩擦力がする仕事 $F d$と等しいと置けば，$K_0 = \frac{M}{2}u_0^2 = F d \ $から 停止距離は $ \ d= \frac{K_0}{F} = 20 {\rm \ m\ }$ である。

図：カーリングストーンの円環モデル

図の角度$\ \theta \ $ の位置の円環要素$ \delta M(\theta) $の氷に対する相対速度ベクトルは，$\bm{v} = (v_x, v_y)  = (u_x - w \sin \theta, \  u_y + w \cos \theta) \ $であり，その大きさは，$v = \sqrt{v_x^2+v_y^2} = \sqrt{ u_x^2 + u_y^2 + w^2 - 2w (u_y \cos \theta - u_x  \sin \theta ) }\ $である。したがって，摩擦力は，$\bm{F} = - \frac{\mu}{v} M g \  (u_x - w \sin \theta, \  u_y + w \cos \theta) \ $ となる。

また，この円環要素に働く摩擦力のトルクの大きさは，

$ N =  (\bm{R} \times \bm{F})_z = R_x F_y - R_y F_x = R\  ( \cos \theta F_y - \sin \theta F_x) $

$  \quad = - \frac{\mu}{v} R M g \ (w + u_y \cos \theta - u_x \sin \theta )$

［１］カーリングと氷物性（前野紀一，2010）

［２］カーリングストーンの動力学（前野紀一，2013）

［３］カーリング・ストーンはなぜ曲がるか？新理論:エッジ・モデル（前野紀一，2018）
［４］カーリング・ストーンの曲がりの説明について（対馬勝年，2011）
［５］カーリングストーンと同様にカール運動する物体（森克馬，対馬勝年，2018）

カーリングの原理（１）

立教大学の村田次郎さんの名前は，TRIUMFの時間反転実験の記事で知った。偏極したリチウム8のベータ崩壊からでてくる電子の縦偏極との相関をみるというものだ。σ・(ｊ×ｐ) という項の係数Rを調べたところ，時間反転を破る効果は見えていない。最近では余剰次元探索の実験を手がけている。

その村田さんが，カーリングのストーンが曲がる原理を精密実験によって解き明かしたという論文が出た。軽井沢のアイスパークで自分自身が投げた122回のデータを，コンパクトデジタルカメラと三脚だけを用いて測定し，余剰次元実験で開発した画像処理型変位計の技術でミクロン単位で分析した結果である。

その結果は，カーリング石の下面が氷と歯車の様に嚙み合って旋廻する現象が基本であること。摩擦支点の形成確率が，速度に依存することから，速度依存の動摩擦係数が導かれること。動摩擦係数の速度依存性により，左右非対称に旋廻の中心が形成されることがこの謎の答えであることを解き明かしたというものだった。

これによって，時計回り（反時計回り）のストーンは進行方向に向かって右側（左側）に曲がることが説明される。まあ，野球のボールの回転方向と，曲がる方向の関係はマグヌス効果で説明され，カーリングの原理とは違うものの同じ対応関係になっているので，それほど違和感がないかもしれない。 

図：a curling stone rotating, sliding and curving on the ice sheet with pebble olympic fine detailed photo shot from the sealing of curling hall（memeplexによる）

亀趺

朝の散歩で，Googleマップを眺めていたら前栽のあたりの名所旧跡マークで「亀趺（きふ）」というのがでてきた。さっそくいってみた。

前栽駅の北東徒歩7分にある黄檗宗の雲井寺といっても，無人の小さなお堂とお稲荷さんなどがあるだけだ。その一角に亀趺がある。雲井寺の地蔵堂というのが前栽駅前にあり，こちらの方は前をよく通るのでお馴染だった。

亀趺というのは石碑の台石の一種で，大亀の形をしている。これは実は亀ではなく，龍の九子の一つで龍になれなかった「贔屓（ひいき）」という想像上の霊獣だ。贔屓は「一生懸命努力して力を出すさま」を意味するとされた。やがて「特別に便宜を図ったり，力添えをする」という今の意味に転じた。

亀趺はもともと中国の貴族階級の風習だったが，江戸時代に日本でも取り入れられた。日本の亀趺（平勢隆郎）によれば，大名家の墓や神格の顕彰に用いられるようだが，なぜ，前栽の小さな寺にあったのだろう。

亀趺の上の碑の三面に銘文があったので，解読を試みたが，花崗岩が風化しているので，ちょっと全部はわからなかった（たぶん15%くらいは誤りかもしれない）。とくに第三面はほとんどだめだったので二面分だけ。

建立放蕩納妙連経

毎冩一字一華一香

合唱三拝喜捨一銭

皆成彿道以此功徳

八彿出世寺受授記

多宝證明韻言善哉

修瑜枷地右縁無縁

共雪昔雲財法二也

写真：天理前栽雲井寺の亀趺（2022.9.26撮影）

［１］亀趺を持つ石碑の系譜（藤井正直）
［２］亀趺を持つ石碑の系譜（二）（藤井正直）
［３］亀趺を持つ石碑の系譜（三）（藤井正直）

［４］日本近世の亀趺碑（平勢隆郎）

［５］日本近世の亀趺碑−その続（平勢隆郎）

［６］ひかり拓本（奈良文化財研究所）

山所池

物理学科同窓会からの続き

10月8日に久々に開催された同窓会の報告を関係者に送ったところ，当日欠席のＳ正泰さんからメールが来て「蛍池にあった池は学校になったみたいですね」とあった。正泰さんは蛍池東町に下宿しており，蛍池中町の文化住宅に住んでいた自分から最近距離の同級生だった。

幸荘という名前の文化住宅は阪急宝塚線西の蛍ヶ池公園の北西に隣接していて，大家さんは千林大宮に住んでいた。幸荘の近所にＯさんという父の会社の人の知人が住んでいたので，その伝手で，この文化住宅に住むことになった。家賃は当初1.1万円で（数年後オイルショックのあとに値上げ騒動があった），隣接した一軒家に住んでいた管理人のＳさんのところに通帳を持って毎月払い込みに行くのだ。

蛍ヶ池公園の南に道を挟んで山所池（やまんじょいけ）という戦国時代（1595）から続く由緒正しい大きなため池があり，これが正泰さんの指摘していた池だった。阪急蛍池駅に最も近い池なので，これが蛍池だと勘違いする人もいたが，そうではない。

蛍池のほうは，阪急蛍池駅東400mあたりの刀根山病院（元：刀根山療養所，現：国立病院機構大阪刀根山医療センター）にあるひとまわり小さな池の方らしい。当時は，どこにあるのかまったく気にも留めていなかったのだけど。

蛍ヶ池公園では，夏になると盆踊り大会が開催されていて，河内音頭などを踊った記憶がある。山所池は，大半が埋め立てられて，1986年に豊中市立第十八中学校となった。私たちは1981年には池田市石橋，1984年には池田市緑丘に住んでいたので，当時の状況はよく知らない。

図：蛍ヶ池公園=左上，山所池=左中，蛍池=右下（googlemapから引用）

物理学科同窓会（１）

阪大理物アラS51同窓会2022が10月8日土曜日に新大阪で開かれた。S51は卒業年の昭和51年，アラはアラフォーのアラなので，卒業年度には分散がある。当時の物理学科のクラス定員は40名，入学年度は1972年（昭和47年），今年でもう50年になる。

台風やコロナの影響でしばらく中断していて4年ぶりの対面の同窓会。去年一昨年はzoomを使っての開催だった。まだ，コロナの影響はおさまっていないので，出席者は9名，ビュフッェ形式ではなく，大きなアクリル板に仕切られた和洋折衷膳＋フリードリンク（＝ほとんどビールをサーブしてもらう）。

まだ健康寿命（男性70.4歳，女性73.6歳）に達していないので，皆おおむね元気そうだった。そういえば，健康寿命というものが存在することを教えてくれた大石君が，今は，幸福寿命だといっていたのだけれど，怪しいウェルビーイングキャンペーンが頭にあったので，ちゃんと説明を聞くのを怠ってしまった。そこで，自力で調べてみた。

日本語で幸福寿命を検索しても，腸内細菌とホルモンの宣伝しか出てこない。こりゃだめだ。その中で目を引いたのが昨年の日経ビジネスの記事のキーワード「幸福感」は60代以降に上昇傾向だ。ダートマス大学のブランチフラワー教授が145カ国の年齢別主観的幸福度の研究によれば、先進国でも発展途上国でも，欧米でもアジアでも，主観的幸福度は次の図のように「U字型」になっている。

図：年齢別の主観的幸福度（Blanchflowerの論文から引用）

内閣府の2019年の調査でも，60歳以上で主観的満足度が急上昇するとなっている。むむむ。どうなのだろうか。

［１］LWC指標利活用ハンドブック（スマートシティインスティチュート）

［２］世界幸福度調査World Happiness Report2020の概要（産業精神保健研究機構）

［３］国民の幸せな人生（well-being）を政策目標に（第一生命経済研究所）

箱の数

鈴木寛太郎の番組で紹介されていた奈良県立医大の入試問題がおもしろかった。

図１：すべての四角形の数え上げ問題

問題は簡単で，図１のような階段型の図においてとりうるすべての可能な四角形の数を数えるというものだ。このぐらいならば，順番に数えつくすこともできなくはない。

縦（横）が $n$ マスの場合の四角形の数を$S_n$とすると，$S_1=1,\  S_2=5\ $である。次に，$S_{n}-S_{n-1}$の漸化式を，$n \ge 2$の場合に求めてみよう。$S_{n-1}$ブロックの各列の上に1マスが加わっているので，これを含んで新たに加わる四角形の数は，左から1列目の１マスに対して$1 \cdot n$個，2列目の１マスに対して$2 \cdot (n-1)$個，…，$n$列目が$n \cdot 1$個となる。

これらの和が，$S_{n}-S_{n-1}$になり次式で与えられる。$\displaystyle S_{n}-S_{n-1} = \sum_{k=1}^n k \cdot  (n-k+1) = \sum_{k=1}^n \{ (n+1) k -k^2 \} $

$= \dfrac{n(n+1)^2}{2} - \dfrac{n(n+1)(2n+1)}{6} = \dfrac{n(n+1)(n+2)}{6}$

そこで，$S_k-S_{k-1}= \dfrac{k(k+1)(k+2)}{6}\ $の両辺を$k=2$から$k=n$までそれぞれ加えると，

$\displaystyle S_n-S_1 = \sum_{k=2}^n \dfrac{k(k+1)(k+2)}{6} = \dfrac{1}{6}  \sum_{k=2}^n (k^3+3k^2+2k) $

$ S_n = 1 + \dfrac{1}{6} \Bigl\{ \dfrac{n^2(n+1)^2}{4} -1 +3\bigl( \dfrac{n(n+1)(2n+1)}{6} -1 \bigr)  +2\bigl( \dfrac{n(n+1)}{2}-1 \bigr) \Bigr\}$

$ = \dfrac{n(n+1)}{24}\Bigl\{ n(n+1) + 2(2n+1) + 4 \Bigr\} = \dfrac{n(n+1)(n+2)(n+3)}{4!}=_{n+3}C_4$

答えは求まったがとても簡単な式で意味付けができそうだ。鈴木貫太郎の番組の答えをみると小倉悠司さんの素晴らしい解説が紹介されていてオオオッとなった。

図２：四角形数え上げの答えの解釈

与えられた図形の上に次のような$n+3$個の赤点を置く。上から始めて時計回りに2点を縦の分割線，さらに2点を横の分割線として4点選ぶと，1つの四角形に一意的に対応する。つまり，$n+3$個の点から4点選ぶことが四角形の選び方と対応するというわけだ。すごくない。

未来をあきらめない（１）

いや，ほとんどあきらめてるわけだけれど。だいだい，このリノベーションの時代になんで，東京オリンピック（2020）と大阪万博（2025）とカジノIRとリニア新幹線なんだ，と誰かがボヤいていたけれども，全くその通りで時代錯誤で先につながらない利権プロジェクトばかりが蠢いている。統一教会さえ完全に排除できれば，日韓トンネルの方がまだマシと思わず口走ってしまいそうになる。

ソフトイーサの開発で有名な，日本トップクラスのIT技術者である登大遊（1984-）さんが，デジタル庁 デジタル臨時行政調査会作業部会 テクノロジーベースの規制改革推進委員会（第1回）で，「テクノロジーマップ、技術カタログの在り方について」という資料を提出していた。外観はお役人が大好きなポンチ絵テクノロジー満載のコテコテ画像だったので，どうかなと思いつつ読んでみると，久々に心が洗われた。

いつの間にかIT（情報技術）という標語が，動詞も含めたパッケージとしてのDX（デジタルトランスフォーメ—ション）という無意味なカタカナ掛け声に変わってしまった。そのまったく進化しない日本型組織の強固な殻をどうやって内側から破るかという戦略が書かれている。

そこにあるのは，組織内の，(A) 技術研究的な人材（同僚制・独立責任・専門性重視・試行錯誤主義）と (B) 経営事務的な人材（官僚制・指揮命令・上意下達・計画主義）の特徴を捉え，IT技術イノベーションを，(A)に如何に浸透させて，それを(B)の力で展開拡散するかを具体的に検討したものだ。まさに，細胞に浸入して増殖するウイルスのイメージだ。

技術研究的な人材は，試行錯誤・ 業務革新を担い，自分の責任で頭脳を働かせることができる。一方，経営事務的な人材は，組織的な集団思考と決定に頼って仕事をし，大規模化・組織化・運用を担う。この技術研究的な人材(A)をターゲットとして，適切な情報を注入すれば，業務革新的な情報が経営事務的な人材(B)が担う組織に提供され，フィードバックを受けつつ前進できるという作戦だ。

(A) の技術研究的な人材に届く情報の例として，1990年代のコンピュータ雑誌群があげられていた。自分にとっては，1980-1990年代のそれなのだが，本質的に変わらない。登さんが指摘するように，本当に良質で大量の技術情報があふれていた。そのころは，帰宅時の書店通いで，ほとんどの主要雑誌を毎日のように立ち読みすることができた。さらに，日経バイトと日経パソコンとTheBASICとSoftware DesignとMacの雑誌を数冊，定期購読して浴びるように情報をインプットしていた時代だ。

登さんのすごいところは，アイディアをすぐに実装できるところだ。Git + GitLab + 静的 HTML 生成器という開発実例を示している。MarkDownを前提としたHTM重視なのはとてもリーズナブルなのだけれど，おじさんとしては，pdf化されたパッケージもないと安心できない。歳を取ってしまったからだろう。

これを読んで，日本の大学改革がなぜ失敗だったかがよくわかった。大学組織の(A)技術研究的な人材や環境を破壊して，(B) 経営事務的な組織論理を貫徹させようとしているからだ。一方，初中等教育のGIGAスクールの分析は難しい。(A)と(B)が縮退した組織だから。教師の二面性をどう捉えるか。

ノーベル物理学賞2022

先週発表された今年のノーベル物理学賞は，量子力学の基礎が対象だった。アスペやツァイリンガーはこれまでも予想していた人が多かったので（クラウザーの代わりに日本人や他の科学者を含めるパターンだ）昨年のような意外感がなく，当たりましたという声があちこちから聞こえてきた。

NHKは，量子もつれというキーワードを看板にして，その応用面の量子情報科学という新しい分野の開拓につながる大きな貢献をしたというニュースに仕立て上げていた。あたかも，青色発光ダイオードやリチウムイオンバッテリーの話のような組立で，ちょっと違和感がある。

ノーベル財団のページに書かれた2022年のノーベル物理学賞の受賞理由は，"for experiments with entangled photons, establishing the violation of Bell inequalities and pioneering quantum information science" （量子もつれのある光子の実験によって，ベル不等式の破れを確立し，量子情報科学を開拓したことに対して）である。

重要なのは「ベル不等式の破れを確立」したところだ。もし，ベル不等式を最初に提案した，ジョン・スチュアート・ベル（1928-1990）が存命だったら，確実にノーベル賞をもらっていたはずだ。この実験は，量子力学の基礎における最重要テーマである局所実在をどう捉えるかに関わるもので，相対論でいえば，マイケルソン・モーレーの実験に相当する。

数年前，甲南大学で物理教育学会近畿支部主催の大学入試問題検討会があったとき，駿台予備校の牛尾健一さんが参加していた。帰りに岡本駅前の居酒屋でちょっと一杯ということになり，学生のころ微分形式やクリフォード代数を勉強しておけばおもしろかったはずなのにとか，1970年代当時の量子力学の教科書にはなんでベル不等式がなかったのかとか，言いたい放題の雑談をして楽しかった。

1970年代には，ベル不等式やCHSH不等式はわかっていたし，クラウザーの実験も始まっていたけれど，量子力学の観測問題とひとくくりにされたコンテンツは，教科書の外に追いやられていた。

［１］量子論とベルの不等式の破れ（東大理学部物理学科学生展示）

［２］日経サイエンス2022年ノーベル物理学賞（吉田彩）

［３］「ベルの不等式の破れ」について知るための約 100 分（田崎晴明）

［４］量子の不可解な偶然−訳者解説（木村元・筒井泉）

［５］ノーベル物理学賞2022年の解説（日本物理学会）

［６］2022年ノーベル物理学賞解説（彩恵りり）

［７］ベル不等式の意味（K. Sugiyama）

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

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

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. と思ったが・・・）。

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

Ｊアラート

Ｊアラートとは，全国瞬時警報システムのことである。総務省消防庁が運用しており，弾道ミサイル情報，緊急地震速報，津波警報など，対処に時間的余裕のない事態に関する情報を，携帯電話等に配信される緊急速報メール，市町村防災行政無線等により，国から住民まで瞬時に伝達するシステムだ。

10月4日の朝，NHKの7時のニュースを見ていたら，Ｊアラート発出にともない画面が急に変わって，北朝鮮からのミサイル発射にかかわるニュースが延々と長時間垂れ流された。午前7時27分には，北海道と東京都の諸島部（都道府県名なしの町村名表記だけ）にアラートが出され，午前7時28分には，北海道が消えて，青森県と東京都の諸島部になった。

北朝鮮から一発のミサイルが飛んでくると仮定して，NHKは東北・北海道と東京諸島部に同時に警報が出されているという事態をおかしいと考えないのか。十分に吟味されていない情報がコメントもなしにそのまま右から左へと提供されている。案の定，警戒不要の東京諸島部には誤ってＪアラートが発出されていたとの説明が後日あった。

そもそも，Ｊアラートが発出された時刻は，丁度ミサイルが津軽海峡の約800km（国際宇宙ステーションの軌道高度400kmの2倍）上を数km/sで通過している時間だ。このタイミングでアラートを発出する必然性がほとんどないにもかかわらず（防衛省も内閣府もわかっているはずだろう），安全な場所に避難して下さいとさんざん煽っていた。いったいどういうこと。これでは，内閣支持率が下がって国会で問題が発生したときにいつも都合よくミサイルが飛んでくるとか，軍備増強・憲法改正の環境作りのための過剰宣撫といわれても仕方がない。

発射時刻・場所 7時22分，北朝鮮慈江道舞坪里

通過時刻・場所 7時29分，青森県津軽海峡

落下時刻・場所 7時44分，釜石市の東3200kmの太平洋

飛距離 4600km，到達高度 1000km，速度マッハ17

J-アラート 午前7時27分 北海道＋東京都諸島部

J-アラート 午前7時29分 青森県＋東京都諸島部

水平方向の平均速度vx0は，vx0 = 4600/(22*60) = 3.5 km/sなので，津軽海峡上空までは，6−7分，1360±100 kmとなる。4600-1360=3140 kmだから，残りの距離が日本から落下地点までの距離3200kmとほぼ一致する。

一方，到達高度hが1000km，飛行時間Tが 1320sである。一様重力場における斜方投射モデルを採用して，鉛直方向の初速度をvy0，有効重力加速度をg とする。h = vy0 (T/2) -g/2 (T/2)^2 = g/2 (T/2)^2 が成り立つので，g = 2h/660^2 = 0.0046 km/s^2，vy0 = 2h/660 = 3.0 km/s となる。初速度v0が v0^2=vx0^2+vy0^2 を満たすので，v0 = 4.6 km/s （マッハ13.5）だ。この単純なモデルでは，マッハ17=5.8 km/s はうまく再現できなかったし，有効重力加速度gの値がもっともらしいのかどうかも微妙である。