昔語りの老害などと仰らずに、温故知新の心でお付き合いください。
１．はじめに
アナログ プラニメータは、図１に示すように非常にシンプルな器械ですが、Ｏ部を机上に固定し、Ｔ部の照準マークで図形の外周をなぞるだけで、その図形の面積を測ることができます。学生時代の実験で一度使ったことがありますが、その原理については分からないまま、ウン十年も過ぎてしまいました。最近、部屋の整理中に、以前にネットオークションで入手していた死蔵品が見つかったため、これを機会に動作原理の探求をしてみることにしました。
図１ ポーラー プラニメータ（プラニメータの代表機種）
（部位の命名は自己流です。以下同様）
同じアナログ機器でも、「図を拡大縮小するパンタグラフ」や「掛算や割算のための計算尺」であれば、幾何や対数の知識だけで原理は容易に理解できます。しかし、プラニメータはなかなかの難物です。1850年代、我ら電気屋が神と崇める J.C.Maxwell も面積計の研究に手を染め、図２のような複雑な装置を論文発表しています。しかし、その直後に、J.Amsler が先行して発明していた図１の Amsler ポーラー プラニメータの存在を知り、そのシンプルさや理論の明快さに感服して兜を脱いだ（意訳）という話が伝わっています。偉人さえをも唸らせたこのツール、凡人が簡単に発明したり理解できるものではなさそうなことは窺い知ることができます。
図２ Maxwell のプラニメータ（Platometer）
出典: Pre-Amsler Planimeters
https://persweb.wabash.edu/facstaff/footer/planimeter/PreAmsler.htm
それでも、ネットを検索すれば、原理を説明した多くの文書を拾い出すことができます。しかし、省略が多いためか、何となく分かったような気分にはなれるものの、どこか引っかかるところが残ります。自信をもって他人に説明できるほどのレベルには達しません。曖昧なところを突かれれば、直ぐにでも行き詰まって立ち往生してしまいそうです。そこで、既存の文書の後追いは止め、手持ちのプラニメータをモデルにして、自分が納得できるまで、曖昧さが残らないように、ゼロから証明を試みようと思い立ちました。
探求のための補助ツールとして、MATLAB (R2019a) の Symbolic Math Toolbox を利用しています。ただし、結果を共有するだけであれば、この環境が必須というわけではありません。
２．プラニメータの各部構造
図３に各部の拡大写真をまとめています。
図３ プラニメータの各部拡大写真
図１では全部が一体になっているように見えますが、アーム１と２は簡単に分離でき、個別に扱うことができます。アーム１のＯ端の下部には短い針が出ており、これを机の面に刺すことで固定できます。しかし、片持ち梁としてアーム１の全体を支えるだけの力はなく、アーム２と組み合わされるまでは、Ｊ端は机上に倒れ落ちた状態になっています。また、アーム１はＯ端を中心にして回転させて使うのですが、ここには軸受けのような可動部はありません。アーム全体が金属の一体ものです。机の面との摩擦は気にせず、Ｏ端の土台ごと回転させる割り切った構造です。
アーム２は、目盛車Ｗ・補助車・支柱の３点で支えられて自立し、机の面を確実にホールドしています。Ｊ端の上部には竪穴があり、ここに、アーム１のＪ端下部にある「先端が球状になった連結棒」を落とし込みます。この連結部は単純な構造ですが、ガタがないよう精度良く仕上げられ、自在継手として働きます。アーム１の荷重の半分はここで支えられますが、「アーム２を保持している３つの接地点」を浮かせるような捻じりの力は伝わりません。
アーム２のＴ端には、下面に照準マークを彫り込んだ分厚い平凸レンズが取り付けられています。このレンズと机の面との間には若干の隙間があり、両者が擦れ合うことはありません。測定時には、この照準マークを図形の輪郭線に合わせて、時計回りになぞっていきます。簡単な作業のように思えますが、実際にやってみると、なかなかの技術が要ります。その原因は、Ｔ点を移動させるときの反作用が、思い掛けない方向に働くためです。鉛筆で線をなぞるのとは勝手が違います。反作用の外乱に負けないように、レンズをしっかりと保持していなければなりません。力み過ぎによる無理なトルクが他の機構部に伝わらないよう、レンズ外周にある把持部分は可動リング構造になっています。
アーム２のＪ端寄りのＷ点には、プラニメータの最重要部品である目盛車が取り付けられています。この車軸はアームと平行に取り付けられているため、アームと直角な方向に移動したときには、机の面との相互作用により車は自由に回転します。しかし、平行に移動したときには横滑りするだけです。斜め方向に移動したときには、直角方向の成分に相当する量だけ回転します。Ｊ端を左、Ｔ端を右に見るように置いたとき、アームを手前に移動すると目盛が増える方向に回転します。
目盛車には全周を100等分した目盛が刻まれており、さらにバーニア目盛も付いているので、1/2000 回転の分解能で読み取れます。なお、大きな図形を測定するときには、1回転以内の変化だけでは収まり切らなくなるため、ウォームギアの目盛によって回転回数も分かるようになっています。目盛車はピボット軸受けで支持されているため摩擦トルクは非常に小さく、勢い良く回すと、小さな慣性モーメントにも拘わらず、長いこと回り続けます。そのため、回転方向への移動に対して車がスリップすることはなく、測定精度が落ちる心配はなさそうです。
補助車の車軸は、目盛車のそれとは直交するように取り付けられています。しかし、少々ガタもあり、精度にはあまり気を遣っていないようです。目盛があるわけでもなく、あえて車にする必要もなさそうに思えます。それでも車にしているのは、たぶん、反作用によるトレース乱れを減らそうとする配慮なのでしょう。残るもう一か所の接地点である支柱の簡素な作りとは対照的です。
ここで使っているプラニメータは、比較的簡単な構造のため以上が全てです。もっと高級な製品では、図形の縮尺に合わせてアーム２の長さを調節できたり、目盛をワンタッチでリセットできるボタンが付いていたりと、利便性が向上しています。しかし、ここでは原理の探求が目的ですから、簡素な製品の方が分かりやすくて好都合です。
３．プラニメータの動作
図４は、面積を測定する手順を示したアニメーションです。手持ちのプラニメータの実測パラメータをそのまま使ってシミュレーションしています。
図４ 測定の手順
基本的には、計測系全体を机の上方から見下ろした平面図として描いています。しかし、目盛車だけは特別扱いです。この平面図上に表示されるべき本来の目盛車は、Ｊ点付近にある緑色の短い線だけです。しかし、これだけでは回転の様子が把握できないので、Ｔ点方向から観察した目盛車の円形の立面図も重ね描きしています。
目盛車の外周に接した黒色の細い線は、立面図として描いた机の面のつもりです。これにより、机面から受ける作用によって車が回転している様子が観察できます。なお、測定の開始時と終了時だけは、目盛車の拡大図を表示して指示値を読み取れるようにしています。
図の右側にある青色の複雑な曲線が被測定図形です。まずは、この図形の外側の適当な位置にＯ点を設置します（図形の内側に設置するのは厳禁です）。適当とはいっても、Ｔ点が図形全体をアクセスできるような位置でなければなりません。
次に、被測定図形の輪郭線上の任意の位置に測定の始点を定めます。ここにＴ点を置いた状態で、目盛車の目盛を0位置まで回してリセットします。その後、Ｔ点を図形の輪郭線に沿わせて時計回りになぞって行きます。移動につれて、目盛車は行きつ戻りつしながら回転します。その傍らに表示される $n$ の値は、ウォームギアが示す回転回数のつもりです。被測定図形を一周し、始点の位置まで戻ったところで移動を止めます。そのときの目盛を読んで、次の式で計算すれば面積 $\Delta S$ [$\mathrm{cm}^2$] が求まります。
$$\Delta S=100\times(n+\frac{m}{10})\tag{1}$$
ただし、$n$: ウォームギアの目盛で読んだ回転回数、$m$: 一周を10とみなした目盛車の目盛の読み
ここで100という係数が唐突に現れますが、このような値になるように設計されているからです。詳細は後回しにして、実測パラメータを使った結論だけを先に示すと次のようになります、
$$
\begin{align}
K&=2\,\pi\,r_\mathrm{w}\,r_2=2{\times}3.1416{\times}(19.05/2){\times}167\\
&=9994.5\,[\mathrm{mm}^2]\simeq100\,[\mathrm{cm}^2]\tag{2}
\end{align}
$$
ただし、$r_\mathrm{w}$: 目盛車の半径 $[\mathrm{mm}]$、$r_2$: アーム２の長さ $[\mathrm{mm}]$
プラニメータによるこの測定結果は、別の方法で数値計算した面積と良く一致しています。この面積は、多角形（$N$角形）で表現された被測定図形の整列した座標データを使い、次の式から求めています。
$$
\Delta S=\frac{1}{2}\left|\sum_{n=1}^{N}(\boldsymbol{V}_n{\times}\boldsymbol{V}_{n+1})\right|\tag{3}
$$
ただし、$\boldsymbol{V}_n$: ベクトル $\overrightarrow{\mathsf{OT}_n}$、閉曲線なので $\boldsymbol{V}_{N{+}1}=\boldsymbol{V}_1$
また、実機を使って図４と同一の実寸試料を測定した結果は 232.0 [mm] でした。大した練習もせず、１回だけ試した結果にしてはマアマアの出来です。（いまさら測定名人になったところで稼ぎにはなりません。何回もやってみる気は起きません。）
４．原理の探求
４.１ 試行１（独立変数 θ1, θ2）
ここで使う各変数の意味は図５に示すとおりです。殆どの変数は $\rm x,y$ 座標上で扱いますが、目盛車の移動量だけは $\rm u,v$ 座標上でも表現します。$\rm v$ 軸はアーム２と平行、$\rm u$ 軸はそれと直交した方向で、目盛車の目盛値が大きくなる向きを正とします。図５の縮尺は NTS ですが、手持ちの実機での値は、$r_1{=}164$、$r_2{=}167$、$a{=}16$、$b{=}15.5$、$r_w{=}19.05/2$ （いずれも [mm]）です。ちなみに、質量は 230 [g] です。
図５ 各変数の定義
図中の全ての変数をそのまま取り込んで機械的に立式/変形できれば、文章による冗長な説明を省けて、書き手も読み手も楽になります。しかし、その反面、式の変形は複雑になって難航しそうです。そこで、少々寄り道にはなりますが、変数 $b$ だけは事前に除外しておこうと思います。次段で示すように、$b$ の値によって目盛車の横滑り方向（$\rm v$ 方向）への移動量は変化しますが、回転方向（$\rm u$ 方向）の移動量は影響を受けません。それゆえ、目盛の値を決める $\rm u$ 方向の変化量だけに着目している本ケースでは、以後は $b{=}0$ と考えて検討を進めることにします。これにより、目盛車の軸はアーム２の中心軸と重なり、目盛車の中心点Ｗもこの上に移動します。既に、図４ではこの条件を取り入れています。
図６ 変数 $b$ の影響
図６左のように、アーム２がその角度を保ちながら併進運動する場合には、目盛車の取り付け位置 $b$ によらず、その回転量が同一であることは明らかです。一方、回転運動する図６右のような場合には、次のように考えることができます。この場合、回転の中心は素直にＪ点に置きたいところですが、念のため、より一般的な任意のＣ点を考えることにします。図において、目盛車も含めたアーム２の全体が、Ｃ点を中心にして微小角度 $\delta\theta$ だけ回転したとき、目盛車の中心点の移動距離 $\delta l$ は $f\delta\theta$ となります。また、図中の $\alpha$ の値は $\sin^{{-}1}(g/f)$ で表せます。ここで、$f$ は $g,h,b$ の値で決まる量なので、 $\delta l$ の値はまだ $b$ の影響下にあります。しかし、これを $\rm u$, $\rm v$ 方向に分解すると $\delta u$, $\delta v$ は次のようになります。
$$
\begin{align}
\delta u&=\delta l\!\cdot\!\sin\alpha=f\!\cdot\!\delta\theta\!\cdot\!\sin(\sin^{{-}1}\frac{g}{f})\\
&=f\!\cdot\!\delta\theta\!\cdot\!\frac{g}{f}=g\,\delta\theta\tag{4}
\end{align}
$$
$$
\begin{align}
\delta v&=\delta l\!\cdot\!\cos\alpha=f\!\cdot\!\delta\theta\!\cdot\!\cos(\sin^{{-}1}\frac{g}{f})\\
&=f\!\cdot\!\delta\theta\!\cdot\!\frac{\sqrt{f^2-g^2}}{f}
=\sqrt{f^2-g^2}\ \delta\theta=(b+h)\delta\theta\tag{5}
\end{align}
$$
これより、横滑り量 $\delta v$ は $b$ の影響を受けますが、車の回転方向の移動量 $\delta u$ は $b$ とは無関係になることが分かります。これにより $b{=}0$ と置くことが可能になります。
さて、他の殆どの関連文書では、さらに結果を先取りし、$a$ の値も $a{=}r_2$ などと設定して、式を簡略化しようとしています。しかし、少々飛躍ぎみで、疑義を感じさせる要因にもなります。前述の $b$ を変えても、逐一の計測量は変わりません。しかし、$a$ を変えた場合には最終結果が同じになるだけで、逐一の動作は異なってきます。現時点で $a$ まで任意の値にできると考えるのはフライングぎみです。この記事では、$a$ は与えられた値のままにして検討を進めます。
プラニメータは二次元的な構造なので、交流回路計算に慣れている電気屋としては、まずは、ベクトルよりも複素数を使って考える方が都合が良さそうに感じます。ということで、ジョイント部Ｊ、トレース端Ｔ、目盛車Ｗの各位置の座標をそれぞれ $z_\mathrm{j}{=}x_\mathrm{j}{+}iy_\mathrm{j}$ 、 $z_\mathrm{t}{=}x_\mathrm{t}{+}iy_\mathrm{t}$ 、 $z_\mathrm{w}{=}x_\mathrm{w}{+}iy_\mathrm{w}$ と置くと、パラメータ $\theta_1$, $\theta_2$, $r_1$, $r_2$, $a$ で表現した各位置は次のようになります。ここでは、$\rm x$ 軸を実軸、$\rm y$ 軸を虚軸としています。
$$z_\mathrm{j}=r_1 e^{i\theta_1}\tag{6}$$
$$z_\mathrm{t}=z_\mathrm{j}+r_2 e^{i\theta_2}=r_1 e^{i\theta_1}+r_2 e^{i\theta_2}\tag{7}$$
$$
\begin{align}
z_\mathrm{w}&=z_\mathrm{j}+a e^{i\theta_2}=r_1 e^{i\theta_1}+a e^{i\theta_2}\\
&=z_\mathrm{t}-(r_2-a)e^{i\theta_2}\tag{8}
\end{align}
$$
最終的に必要になるのは、$z_\mathrm{w}$ の $\rm u$ 軸成分です。 $z_\mathrm{w}$ の値を、$\rm u$ 軸が実軸、$\rm v$ 軸が虚軸となるように変換して $z_\mathrm{w'}$ と置けば、この値は次のように表せます。$\rm u, v$ 軸は $\rm x, y$ 軸よりも $\pi/2{-}\theta_2$ だけ遅れているので、$\rm u, v$ 軸上から観察する $z_\mathrm{w}$ の値は、$\rm x, y$ 軸上で $z_\mathrm{w}$ を $\pi/2{-}\theta_2$ だけ進めた値と同じになっています。（ここでいう位相の「進み/遅れ」は、時計の針による表現とは逆です）
$$
\begin{align}
z_\mathrm{w'}&=z_\mathrm{w}e^{i(\frac{\pi}{2}-\theta_2)}={z_\mathrm{t}-(r_2-a)e^{i\theta_2}}e^{i(\frac{\pi}{2}-\theta_2)}\\
&=z_\mathrm{t}e^{i(\frac{\pi}{2}-\theta_2)}-(r_2-a)e^{i\frac{\pi}{2}}=i(z_\mathrm{t}e^{-i\theta_2}-r_2+a)\tag{9}
\end{align}
$$
なにやら $z_\mathrm{w'}$ と $z_\mathrm{t}$ の関係式らしきものは得られましたが、この中の $\theta_2$ は未知数のままです。色々と解決策を検討してみましたが、柔軟性が欠如した頭では容易に先へは進めそうもありません。一旦保留して、別の視点から攻めてみることにします。（ここで得られた多くの式はあとで復活します。無駄にはなっていません。）
４.２ 試行２（独立変数 xt, yt）
先の「試行１」では、$\theta_1$ と $\theta_2$ を独立変数として検討を始めましたが、途中で行き詰ってしまいました。やはり、ブラニメータへの入力値であるＴ点の座標 $x_{\rm t}$, $y_{\rm t}$ を独立変数として扱う方が自然なのでしょう。この場合、図７のように、２つの円の交点を求める計算が必要になります。交点の座標が求められれば、$\theta_1$ と $\theta_2$ も自動的に決まるのであとは簡単に処理できそうです。
図７ Ｔ点の座標からＪ点の座標を求める
しかし、交点は２つあるので、どちらを採用するかを予め考えておく必要があります。この選び方によって、Ｏ点からＴ点を見たとき、2つのアームの作る形が「＜」の字になるか「＞」の字になるかの違いが現れます。これに関連して、まずは、Ｔ点に許容される移動範囲について調べてみます。プラニメータの構造上、機械的な可動範囲は半径が $|r_1{-}r_2|$ と $r_1{+}r_2$ の２つの同心円に囲まれた範囲です。ただし、同心円の円周上は除外します。なぜなら、次の２つの条件下では、Ｔ点は身動きができないロック状態に陥ってしまうからです。
- ２つのアームが完全に畳み込まれた $|r_1{-}r_2|$ の円周上から、「Ｏ点を中心とした放射方向」にＴ点を直進させようとしたとき
- ２つのアームが完全に伸び切った $r_1{+}r_2$ の円周上から、Ｏ点に向かってＴ点を直進させようとしたとき
さらに、これらの状態では、次の理由によって、円の中心線方向への移動だけではなく、円周方向への移動もできなくなります。
Ｔ点が円周方向に移動するためには、２本のアームが１本の直線状態を保ちながら、Ｊ点も半径 $r_1$ の円周上を移動しなければなりません。Ｊ点にこの移動のための力を与えるには、Ｔ点周りにトルクをかけてやる必要があります。しかし、既述のように、Ｔ点の把持部分は可動リング構造になっており、ここからトルクをかけることはできません。
結局、四方八方どちらの方向にも身動きができない状態になります。しかし、いずれにしても、Ｔ点が同心円上に乗らないように注意しておけば、特に意識することなく、計測の開始時点で選んだ「＜」あるいは「＞」の形状が計測終了までそのままに保たれます。
それでは、Ｏ点を原点 [0,0] として、2つの円の式を次のように置き、交点の座標を計算していきます。
$$
\begin{align}
&x^2+y^2=r_1^2\tag{10}\\
&(x-x_{\rm t})^2+(y-y_{\rm t})^2=r_2^2\tag{11}
\end{align}
$$
(11)式を次のように展開します。
$$x^2-2x_{\rm t}x+x_{\rm t}^2+y^2-2y_{\rm t}y+y_{\rm t}^2=r_2^2\tag{12}$$
(10)式から(12)式を引いたものを変形し、交点の $y$ の仮値を求めます。
$$2x_{\rm t}x-x_{\rm t}^2+2y_{\rm t}y-y_{\rm t}^2=r_1^2-r_2^2\tag{13}$$
$$y=-\frac{x_{\rm t}}{y_{\rm t}}x+\frac{x_{\rm t}^2 +y_{\rm t}^2+r_1^2-r_2^2}{2y_{\rm t}}\tag{14}$$
まずは、$y_\rm t{\neq}0$ として、(14)式を(10)式に代入して整理します。
$$x^2+\left(-\frac{x_{\rm t}}{y_{\rm t}}x+\frac{x_{\rm t}^2 +y_{\rm t}^2+r_1^2-r_2^2}{2y_{\rm t}}\right)^2=r_1^2$$
$$
\begin{align}
\left\{1+\left(\frac{x_{\rm t}}{y_{\rm t}}\right)^2\right\}x^2&-\frac{x_{\rm t}}{y_{\rm t}^2}\left(x_{\rm t}^2 +y_{\rm t}^2+r_1^2-r_2^2\right)x\\
&+\left(\frac{x_{\rm t}^2 +y_{\rm t}^2+r_1^2-r_2^2}{2y_{\rm t}}\right)^2-r_1^2=0\tag{15}
\end{align}
$$
少々複雑になったので、$x$ の2次式の各項をまとめて次のように簡略化します。
$$Ax^2+Bx+C=0$$
これより、2つの交点の $x$ の値はそれぞれ次のようになります。
$$x_\oplus=\frac{-B+\sqrt{B^2-4AC}}{2A}\tag{16}$$
$$x_\ominus=\frac{-B-\sqrt{B^2-4AC}}{2A}\tag{17}$$
(14)式を使って、これらに対応する $y$ の値はそれぞれ次のようになります。
$$y_\oplus=-\frac{x_{\rm t}}{y_{\rm t}}x_\oplus+\frac{x_{\rm t}^2 +y_{\rm t}^2+r_1^2-r_2^2}{2y_{\rm t}}\tag{18}$$
$$y_\ominus=-\frac{x_{\rm t}}{y_{\rm t}}x_\ominus+\frac{x_{\rm t}^2 +y_{\rm t}^2+r_1^2-r_2^2}{2y_{\rm t}}\tag{19}$$
最終的には、これらの交点 $(x_\oplus, y_\oplus)$ あるいは $(x_\ominus, y_\ominus)$ がＪ点 $(x_{\rm j}, y_{\rm j})$ になるわけですが、アームの形（「＜」,「＞」）のどちらを採用するかによって $\oplus$, $\ominus$ の選び方が変わってきます。
図８ 解の選択（$y_\rm t{\neq}0$ のとき）
図８は、各種のＴ点（$x_{\rm t}$, $y_{\rm t}$）における $\oplus$ 解, $\ominus$ 解の位置を示した図です。Ｔ点が存在する象限によってアームを色分けしています。(16), (17)式からも分かるように、$\rm x$ 成分が大きい点が $\oplus$ 解、小さい点が $\ominus$ 解となります。また、アームが作る形が「＜」のものは実線、「＞」のものは破線で描いています。これらを対比しながら観察すると次のようになります。
「＜」を採用する場合：
$$
\begin{align}
(x_{\rm j}, y_{\rm j})=
\begin{cases}
(x_\ominus, y_\ominus)&(y_{\rm t}>0 のとき)\\
(x_\oplus, y_\oplus)&(y_{\rm t}<0 のとき)
\end{cases}\tag{20}
\end{align}
$$
「＞」を採用する場合：
$$
\begin{align}
(x_{\rm j}, y_{\rm j})=
\begin{cases}
(x_\oplus, y_\oplus)&(y_{\rm t}>0 のとき)\\
(x_\ominus, y_\ominus)&(y_{\rm t}<0 のとき)
\end{cases}\tag{21}
\end{align}
$$
以上は $y_\rm t{\neq}0$ の場合でしたが、$y_\rm t{=}0$ の場合は次のようになります。
(11)式の $y_\rm t$ を 0 と置き、(10)式から得られる $y^2=r_1^2-x^2$ を代入して変形すると、
$$(x-x_{\rm t})^2+r_1^2-x^2=r_2^2$$
$$-2x_{\rm t}x=r_2^2-r_1^2-x_{\rm t}^2$$
$$x=\frac{x_{\rm t}^2+r_1^2-r_2^2}{2x_{\rm t}}\tag{22}$$
この場合は2つの交点とも $x$ の値は同じです。しかし、 これを(10)式に代入すると、2つの $y$ の値が得られます。
$$y_\oplus=\sqrt{r_1^2-x^2}=\sqrt{r_1^2-\left(\frac{x_{\rm t}^2+r_1^2-r_2^2}{2x_{\rm t}}\right)^2}\tag{23}$$
$$y_\ominus=-\sqrt{r_1^2-\left(\frac{x_{\rm t}^2+r_1^2-r_2^2}{2x_{\rm t}}\right)^2}\tag{24}$$
これを $y_\rm t{\neq}0$ の場合と同様に図示すると図９のようになります。このケースでは、Ｔ点は $\rm x$ 座標上だけに存在します。アームの色はＴ点の $\rm x$ 座標値により塗り分けています。(23),(24)式より、Ｊ点の $\rm y$ 座標値が正の場合に $\oplus$ となります。
図９ 解の選択（$y_\rm t{=}0$ のとき）
この図から、
「＜」を採用する場合：
$$
\begin{align}
(x_{\rm j}, y_{\rm j})=
\begin{cases}
(x, y_\oplus)&(x_{\rm t}>0 のとき)\\
(x, y_\ominus)&(x_{\rm t}<0 のとき)
\end{cases}\tag{25}
\end{align}
$$
「＞」を採用する場合：
$$
\begin{align}
(x_{\rm j}, y_{\rm j})=
\begin{cases}
(x, y_\ominus)&(x_{\rm t}>0 のとき)\\
(x, y_\oplus)&(x_{\rm t}<0 のとき)
\end{cases}\tag{26}
\end{align}
$$
図４のアニメーションの作成では、ここで得られた式を便利に使って役立てています。しかし、複数の「場合分け」もあって取り扱いが非常に煩雑です。このまま原理の証明にまで発展させて行くには苦労がありそうです。
また、アニメーションでの目盛車の最終指示値は、上式から直接求められているわけではありません。被測定図形の輪郭線を短い線分に分割したうえ、各区画での目盛車の微小変化を逐一積算していくという手間をかけています。これから考えると、原理説明のための式の展開は、微小変位量を基本にして行うのが本来の姿かもしれません。そうであれば、上式よりも、「試行１」で使った式の方がシンプルで扱いやすそうな気がしてきます。
４.３ 試行３（微小変位に注目）
ここでは、「試行１」で一時的に放棄した(7),(8),(9)式を蘇らせ、利用しやすい $\rm x,y$ 座標形式に変形してから、微小変位量どうしを関係づける式を導いていきます。
(7)式から必要部分のみを取り出して再掲します。
$$z_\mathrm{t}=r_1 e^{i\theta_1}+r_2 e^{i\theta_2}\tag{7'}$$
指数表現のこの式を、一般の複素座標で書き直すと、
$$
\begin{align}
x_{\rm t}+i y_{\rm t}&=r_1(\cos\theta_1+i\sin\theta_1)+r_2(\cos\theta_2+i\sin\theta_2)\\
&=r_1\cos\theta_1+r_2\cos\theta_2+i(r_1\sin\theta_1+r_2\sin\theta_2)\tag{27}
\end{align}
$$
これを実数部、虚数部に分けて対応づけると、
$$
\begin{align}
\begin{cases}
x_{\rm t}=r_1\cos\theta_1+r_2\cos\theta_2\\
y_{\rm t}=r_1\sin\theta_1+r_2\sin\theta_2
\end{cases}\tag{28}
\end{align}
$$
微小変位量として関係づけるために、これらを全微分すると、
$$
\begin{align}
{\rm d}x_{\rm t}&=\frac{\partial x_{\rm t}}{\partial\theta_1}{\rm d}\theta_1+\frac{\partial x_{\rm t}}{\partial\theta_2}{\rm d}\theta_2\\[1.5ex]
&=-r_1\sin\theta_1{\rm d}\theta_1-r_2\sin\theta_2{\rm d}\theta_2\tag{29}
\end{align}
$$
$$
\begin{align}
{\rm d}y_{\rm t}&=\frac{\partial y_{\rm t}}{\partial\theta_1}{\rm d}\theta_1+\frac{\partial y_{\rm t}}{\partial\theta_2}{\rm d}\theta_2\\[1.5ex]
&=r_1\cos\theta_1{\rm d}\theta_1+r_2\cos\theta_2{\rm d}\theta_2\tag{30}
\end{align}
$$
行列表現で書き換えると、
$$
\begin{align}
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]
&=
\left[
\begin{array}{c}
\partial x_{\rm t}/\partial\theta_1 & \partial x_{\rm t}/\partial\theta_2\\
\partial y_{\rm t}/\partial\theta_1 & \partial y_{\rm t}/\partial\theta_2
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}\theta_1\\
{\rm d}\theta_2
\end{array}
\right]\\[2ex]
&=
\left[
\begin{array}{c}
-r_1\sin\theta_1 & -r_2\sin\theta_2\\
r_1\cos\theta_1 & r_2\cos\theta_2
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}\theta_1\\
{\rm d}\theta_2
\end{array}
\right]
=\boldsymbol{A}
\left[
\begin{array}{c}
{\rm d}\theta_1\\
{\rm d}\theta_2
\end{array}
\right]\tag{31}
\end{align}
$$
先の「試行１」では、$(x_{\rm t}, x_{\rm t})$ から $(\theta_1, \theta_2)$ を求めるところで行き詰まりましたが、微小変位量の関係式として線形化した現状では、逆行列を使って次のように簡単に求めることができます。
$$
\begin{align}
\left[
\begin{array}{c}
{\rm d}\theta_1\\
{\rm d}\theta_2
\end{array}
\right]
&=\boldsymbol{A}^{-1}
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]
=\boldsymbol{B}
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]\\[2ex]
&=
\left[
\begin{array}{c}
\partial\theta_1/\partial x_{\rm t} & \partial\theta_1/\partial y_{\rm t}\\
\partial\theta_2/\partial x_{\rm t} & \partial\theta_2/\partial y_{\rm t}
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]\tag{32}
\end{align}
$$
ここで、$\boldsymbol{B}$ の値は、
$$
\begin{align}
\boldsymbol{A}=
\left[
\begin{array}{c}
a_{11} & a_{12}\\
a_{21} & a_{22}
\end{array}
\right]
\end{align}
$$
として、次のようになります。
$$
\begin{align}
\boldsymbol{B}&=\boldsymbol{A}^{-1}=
\frac{
\mathrm{adj}(\boldsymbol{A})
}
{|\boldsymbol{A}|}=
\frac{1}{a_{11}a_{22}-a_{12}a_{21}}
\left[
\begin{array}{c}
a_{22} & -a_{12}\\
-a_{21} & a_{11}
\end{array}
\right]\\[2ex]
&=
\frac{1}{-r_1r_2\sin(\theta_1-\theta_2)}
\left[
\begin{array}{c}
r_2\cos\theta_2 & r_2\sin\theta_2\\
-r_1\cos\theta_1 & -r_1\sin\theta_1
\end{array}
\right]\tag{33}
\end{align}
$$
2本のアームが完全に伸び切った状態（$\theta_1{=}\theta_2$）や完全に畳み込まれた状態（$|\theta_1{-}\theta_2|{=}\pi$）では、$|\boldsymbol{A}|{=0}$ となって上式は使えなくなります。しかし、既に述べているように、これらの状態は避けて計測することにしているため心配は要りません。
それにしても、${\rm d}\theta_1$ や ${\rm d}\theta_2$ を求めるはずの上記行列 $\boldsymbol{B}$ の中には、未知数の $\theta_1$ や $\theta_2$ が複雑に入り込んでいます。このままで先に進めるのかは疑問です。しかし、諦めが早過ぎるのも問題なので、気を取り直して続けることにします。
次に、(8)式からは次の式を取り出します。
$$z_\mathrm{w}=r_1 e^{i\theta_1}+a e^{i\theta_2}\tag{8'}$$
(8')式は、(7')式の $r_2$ が $a$ に変わっているだけですから、(7')式の場合と同様にして、
$$
\begin{align}
\left[
\begin{array}{c}
{\rm d}x_{\rm w}\\
{\rm d}y_{\rm w}
\end{array}
\right]
&=
\left[
\begin{array}{c}
\partial x_{\rm w}/\partial\theta_1 & \partial x_{\rm w}/\partial\theta_2\\
\partial y_{\rm w}/\partial\theta_1 & \partial y_{\rm w}/\partial\theta_2
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}\theta_1\\
{\rm d}\theta_2
\end{array}
\right]\\[2ex]
&=
\left[
\begin{array}{c}
-r_1\sin\theta_1 & -a\sin\theta_2\\
r_1\cos\theta_1 & a\cos\theta_2
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}\theta_1\\
{\rm d}\theta_2
\end{array}
\right]
=\boldsymbol{C}
\left[
\begin{array}{c}
{\rm d}\theta_1\\
{\rm d}\theta_2
\end{array}
\right]\tag{34}
\end{align}
$$
次に、(9)式から抜き出した式を一部変形して下式が得られます。
$$z_\mathrm{w'}=z_\mathrm{w}e^{i(\frac{\pi}{2}-\theta_2)}=i z_\mathrm{w}e^{-i\theta_2}\tag{9'}$$
ここで、$z_\mathrm{w'}$ は ${\rm u,v}$ 座標上で表した $z_\mathrm{w}$ の値だったことを思い返し、$z_\mathrm{w'}{=}u+iv$ と置き換えて、(9')式を一般複素数形式に書き換えると、
$$
\begin{align}
u+iv&=i(x_{\rm w}+iy_{\rm w})(\cos\theta_2-i\sin\theta_2)\\
&=i(x_{\rm w}\cos\theta_2+y_{\rm w}\sin\theta_2-ix_{\rm w}\sin\theta_2+iy_{\rm w}\cos\theta_2)\\
&=x_{\rm w}\sin\theta_2-y_{\rm w}\cos\theta_2+i(x_{\rm w}\cos\theta_2+y_{\rm w}\sin\theta_2)\tag{35}
\end{align}
$$
これを実数部、虚数部に分けて対応づけると、
$$
\begin{align}
\begin{cases}
u=x_{\rm w}\sin\theta_2-y_{\rm w}\cos\theta_2\\
v=x_{\rm w}\cos\theta_2+y_{\rm w}\sin\theta_2
\end{cases}\tag{36}
\end{align}
$$
行列表現にすると、
$$
\begin{align}
\left[
\begin{array}{c}
u\\
v
\end{array}
\right]=
\left[
\begin{array}{c}
\sin\theta_2 & -\cos\theta_2\\
\cos\theta_2 & \sin\theta_2
\end{array}
\right]
\left[
\begin{array}{c}
x_{\rm w}\\
y_{\rm w}
\end{array}
\right]=
\boldsymbol{D}
\left[
\begin{array}{c}
x_{\rm w}\\
y_{\rm w}
\end{array}
\right]\tag{37}
\end{align}
$$
これは座標変換のための式なので、微小変位量に対してだけではなく常に成り立っています。しかし、ここでは各変数の微小変位量どうしを関係づけようとしているので、上式も次のように書き換えます。
$$
\begin{align}
\left[
\begin{array}{c}
{\rm d}u\\
{\rm d}v
\end{array}
\right]=
\boldsymbol{D}
\left[
\begin{array}{c}
{\rm d}x_{\rm w}\\
{\rm d}y_{\rm w}
\end{array}
\right]\tag{38}
\end{align}
$$
以上、(32),(34),(38)式をまとめると、Ｔ点の $\rm x,y$ 座標上の微小変位に対する、Ｗ点の $\rm u,v$ 座標上の微小変位の関係式が得られます。
$$
\begin{align}
\left[
\begin{array}{c}
{\rm d}u\\
{\rm d}v
\end{array}
\right]=
\boldsymbol{D}
\left[
\begin{array}{c}
{\rm d}x_{\rm w}\\
{\rm d}y_{\rm w}
\end{array}
\right]=
\boldsymbol{D}
\boldsymbol{C}
\left[
\begin{array}{c}
{\rm d}\theta_1\\
{\rm d}\theta_2
\end{array}
\right]=
\boldsymbol{D}
\boldsymbol{C}
\boldsymbol{A}^{-1}
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]\tag{39}
\end{align}
$$
どうも、$\boldsymbol{D}\boldsymbol{C}\boldsymbol{A}^{-1}$ の中に、プラニメータの特性の全てが隠されているようです。順を追って、この値を計算していきます。(37),(34),(33)式より、
$$
\begin{align}
\boldsymbol{D}\boldsymbol{C}&=
\left[
\begin{array}{c}
\sin\theta_2 & -\cos\theta_2\\
\cos\theta_2 & \sin\theta_2
\end{array}
\right]
\left[
\begin{array}{c}
-r_1\sin\theta_1 & -a\sin\theta_2\\
r_1\cos\theta_1 & a\cos\theta_2
\end{array}
\right]\\[2ex]
&=
\left[
\begin{array}{c}
-r_1\sin\theta_1\sin\theta_2{-}r_1\cos\theta_1\cos\theta_2 & -a\sin^2\theta_2{-}a\cos^2\theta_2\\
-r_1\sin\theta_1\cos\theta_2{+}r_1\cos\theta_1\sin\theta_2 & -a\sin\theta_2\cos\theta_2{+}a\sin\theta_2\cos\theta_2
\end{array}
\right]\\[2ex]
&=
\left[
\begin{array}{c}
-r_1\cos(\theta_1{-}\theta_2) & -a\\
-r_1\sin(\theta_1{-}\theta_2) & 0
\end{array}
\right]\tag{40}
\end{align}
$$
$$
\begin{align}
\boldsymbol{D}\boldsymbol{C}\boldsymbol{A}^{-1}&=
\frac{1}{-r_1r_2\sin(\theta_1-\theta_2)}
\left[
\begin{array}{c}
-r_1\cos(\theta_1{-}\theta_2) & -a\\
-r_1\sin(\theta_1{-}\theta_2) & 0
\end{array}
\right]
\left[
\begin{array}{c}
r_2\cos\theta_2 & r_2\sin\theta_2\\
-r_1\cos\theta_1 & -r_1\sin\theta_1
\end{array}
\right]\\[2ex]
&=
\frac{1}{-r_1r_2\sin(\theta_1-\theta_2)}*\\[1.5ex]
&\hspace{2em}*
\left[
\begin{array}{c}
-r_1r_2\cos(\theta_1{-}\theta_2)\cos\theta_2{+}r_1a\cos\theta_1 & -r_1r_2\cos(\theta_1{-}\theta_2)\sin\theta_2{+}r_1a\sin\theta_1\\
-r_1r_2\sin(\theta_1{-}\theta_2)\cos\theta_2 & -r_1r_2\sin(\theta_1{-}\theta_2)\sin\theta_2
\end{array}
\right]\\[2ex]
&=
\frac{1}{\sin(\theta_1-\theta_2)}
\left[
\begin{array}{c}
\cos(\theta_1{-}\theta_2)\cos\theta_2{-}(a/r_2)\cos\theta_1 & \cos(\theta_1{-}\theta_2)\sin\theta_2{-}(a/r_2)\sin\theta_1\\
\sin(\theta_1{-}\theta_2)\cos\theta_2 & \sin(\theta_1{-}\theta_2)\sin\theta_2
\end{array}
\right]\tag{41}
\end{align}
$$
なお、目盛車の回転方向の移動量に関係するのは ${\rm d}u$ だけなので、他方の ${\rm d}v$ だけに寄与する「上式最終行の行列の2行目の要素」は、今後の計算では使用しません。もし使用したとしても、先に強制的に $b{=}0$ と置き換えているので、${\rm d}v$ の値は実機の本来の値とは異なったものになります。無視しておくのが安全です。
ところで、「Ｔ点の微小移動量」と「目盛車の回転方向への微小移動量」を関連づけようとしている(41)式の行列の中には、(33)式よりもさらに複雑に未知数 $\theta_1$, $\theta_2$ が入り込んでいます。(9)式を上回って状況が悪化しているように見えます。しかし、それも厭わず、果敢に検討を進めていこうと思います。
そこで、複雑になった式を簡略化するため、行列 $\boldsymbol{D}\boldsymbol{C}\boldsymbol{A}^{-1}$ の1行1列と1行2列の要素を $H_{\rm x}$, $H_{\rm y}$ と置き換えます。
$$
H_{\rm x}= \frac{\cos(\theta_1{-}\theta_2)\cos\theta_2{-}(a/r_2)\cos\theta_1}{\sin(\theta_1-\theta_2)}\tag{42}
$$
$$
H_{\rm y}= \frac{\cos(\theta_1{-}\theta_2)\sin\theta_2{-}(a/r_2)\sin\theta_1}{\sin(\theta_1-\theta_2)}\tag{43}
$$
すると、(39)式より次の関係が得られます。
$$
{\rm d}u=
(\boldsymbol{D}\boldsymbol{C}\boldsymbol{A}^{-1})_{(1行目)}
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]=
\left[
\begin{array}{c}
H_{\rm x} & \hspace{-0.7em}H_{\rm y}
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]=
\boldsymbol{H}{\cdot}{\rm d}\boldsymbol{s}
\tag{44}
$$
ここで、$\boldsymbol{H}{=}[H_{\rm x}\ H_{\rm y}]$ は、面積の測定が可能な全領域に広がるプラニメータ起源の正体不明のベクトル場、${\rm d}\boldsymbol{s}{=}[{\rm d}x_{\rm t}\ {\rm d}y_{\rm t}]$ は、Ｔ点が「被測定図形の輪郭線 $C$」上を移動するときの微小変位ベクトルです。したがって、$\boldsymbol{H}{\cdot}{\rm d}\boldsymbol{s}$ は、 $\boldsymbol{H}$ の「$C$ の接線方向」成分で重みづけされた微小移動量 $|{\rm d}\boldsymbol{s}|$ になっています。被測定図形を一巡しながら ${\rm d}u$ の値を積算したものが、目盛車の回転方向の総移動量 $\Delta u$ になるはずなので、次の式が成り立ちます。
$$
\Delta u=\int{\rm d}u=\oint_C\boldsymbol{H}{\cdot}{\rm d}\boldsymbol{s}\tag{45}
$$
電磁気学に気合を入れて取り組んだことがある電気屋であれば、この式を見てすぐに、次のストークスの定理が思い浮かびます。
$$
\oint_C\boldsymbol{H}{\cdot}{\rm d}\boldsymbol{s}=\iint_S({\rm rot}\boldsymbol{H}){\cdot}\boldsymbol{n},{\rm d}S\tag{46}
$$
ここで、積分範囲として示している $S$ は、輪郭線 $C$ に囲まれた開曲面、${\rm d}S$ はその面上の微小面積、$\boldsymbol{n}$ は、$S$ 上の各点における法線ベクトルで、その向きは、$C$ の周回積分の方向に対して右ネジが進む向きとします。しかし、自由度があまり大きいのは混乱の原因になるので、この記事では $S$ を平面に限定して次のように決めておきます。
- 周回積分の方向は反時計回り（プラニメータの仕様とは逆ですが、数学的に自然な方向）
- その結果として、法線ベクトルは $\rm z$ の正の向きの単位ベクトル ${\bf k}$
(45),(46)式より、もし $({\rm rot}\boldsymbol{H}){\cdot}\boldsymbol{n}$ が一定値 $R$ であれば次の関係が成立します。
$$
\Delta u=\iint_S({\rm rot}\boldsymbol{H}){\cdot}\boldsymbol{n},{\rm d}S=R\iint_S{\rm d}S=R\Delta S\tag{47}
$$
ただし、$\Delta S$：輪郭線 $C$ に囲まれた面積（被測定図形の面積）
すると、面積 $\Delta S$ は次式により求められることになります。
$$\Delta S=\Delta u/R\tag{48}$$
したがって、プラニメータで面積が測定できることの証明は、 $({\rm rot}\boldsymbol{H}){\cdot}\boldsymbol{n}$ が一定であることを証明することに帰着します。$\boldsymbol{H}$ は二次元場で、$H_{\rm z}{=}0$, $\ \partial\!\star\!\!/\partial z{=}0$, $\ \boldsymbol{n}{=}{\bf k}$ であることを考慮すると、
$$
\begin{align}
({\rm rot}\boldsymbol{H}){\cdot}\boldsymbol{n}&=(\nabla{\times}\boldsymbol{H}){\cdot}\boldsymbol{n}=
\left[
\begin{array}{c}
{\bf i} & {\bf j} & {\bf k}\\
\partial/\partial x & \partial/\partial y & \partial/\partial z\\
H_{\rm x} & H_{\rm y} & H_{\rm z}
\end{array}
\right]{\cdot}\boldsymbol{n}\\[2ex]
&=
\left[
\begin{array}{c}
{\bf i} & {\bf j} & {\bf k}\\
\partial/\partial x & \partial/\partial y & 0\\
H_{\rm x} & H_{\rm y} & 0
\end{array}
\right]{\cdot}{\bf k}
=\left(\frac{\partial H_{\rm y}}{\partial x}-\frac{\partial H_{\rm x}}{\partial y}\right){\bf k}{\cdot}{\bf k}\\[2ex]
&=\frac{\partial H_{\rm y}}{\partial x}-\frac{\partial H_{\rm x}}{\partial y}\tag{49}
\end{align}
$$
しかし、(42), (43)式に示すように $H_{\rm x}$ や $H_{\rm y}$ は、直接 $x, y$ の関数としては表されていないので、次のように $\theta_1$ と $\theta_2$ を仲介させて計算します。
$$
\begin{align}
\frac{\partial H_{\rm y}}{\partial x}-\frac{\partial H_{\rm x}}{\partial y}&=
\left(\frac{\partial H_{\rm y}}{\partial\theta_1}\frac{\partial\theta_1}{\partial x}
{+}\frac{\partial H_{\rm y}}{\partial\theta_2}\frac{\partial\theta_2}{\partial x}\right)\\[1ex]
&\hspace{2em}-
\left(\frac{\partial H_{\rm x}}{\partial\theta_1}\frac{\partial\theta_1}{\partial y}
{+}\frac{\partial H_{\rm x}}{\partial\theta_2}\frac{\partial\theta_2}{\partial y}\right)\tag{50}
\end{align}
$$
ここで使われている $x$, $y$ は、$x_{\rm t}$, $y_{\rm t}$ がとり得る全ての値を表しているので、(50)式の中の $\partial\theta_1/\partial x$, $\ \partial\theta_2/\partial x$, $\ \partial\theta_1/\partial y$, $\ \partial\theta_2/\partial y$ は、(32)式の $\boldsymbol{B}$ の各要素 $\partial\theta_1/\partial x_{\rm t}$, $\ \partial\theta_2/\partial x_{\rm t}$, $\ \partial\theta_1/\partial y_{\rm t}$, $\ \partial\theta_2/\partial y_{\rm t}$ と同じものになります。そこで、
$$
\begin{align}
\boldsymbol{B}=
\left[
\begin{array}{c}
b_{11} & b_{12}\\
b_{21} & b_{22}
\end{array}
\right]
\end{align}
$$
として、(50)式を書き換えると、次のようになります。
$$
\begin{align}
\frac{\partial H_{\rm y}}{\partial x}-\frac{\partial H_{\rm x}}{\partial y}&=
\left(\frac{\partial H_{\rm y}}{\partial\theta_1}b_{11}
{+}\frac{\partial H_{\rm y}}{\partial\theta_2}b_{21}\right)\\[1ex]
&\hspace{2em}-
\left(\frac{\partial H_{\rm x}}{\partial\theta_1}b_{12}
{+}\frac{\partial H_{\rm x}}{\partial\theta_2}b_{22}\right)\tag{51}
\end{align}
$$
これを、MATLAB の Symbolic Math Toolbox を使い下記のプログラムで計算させると、次の結果が得られます。これで、$({\rm rot}\boldsymbol{H}){\cdot}\boldsymbol{n}$ の値は、Ｔ端の可動範囲全域で一定値 $R$ であることが確認できました。
$$
({\rm rot}\boldsymbol{H}){\cdot}\boldsymbol{n}=\frac{\partial H_{\rm y}}{\partial x}-\frac{\partial H_{\rm x}}{\partial y}=-\frac{1}{r_2}=R\tag{52}
$$
% planimeter20.m
% ポーラープラニメータのrotationの数式計算
clear
close all
syms r1 r2 th1 th2 a
A = [-r1*sin(th1) -r2*sin(th2) ; r1*cos(th1) r2*cos(th2)];
B = inv(A);
C = [-r1*sin(th1) -a*sin(th2) ; r1*cos(th1) a*cos(th2)];
D = [sin(th2) -cos(th2) ; cos(th2) sin(th2)];
H = D*C*B;
BB = formula(B); % 行列内の各要素の読み出しを可能にする。
HH = formula(H); % 例えば、B(1,1)はエラーになるが、BB(1,1)ならＯＫ。
Hx = HH(1,1);
Hy = HH(1,2);
Rot_raw = diff(Hy,th1)*BB(1,1)+diff(Hy,th2)*BB(2,1) ...
-diff(Hx,th1)*BB(1,2)-diff(Hx,th2)*BB(2,2);
Rot_sim = simplify(Rot_raw)
>> planimeter20
Rot_sim =
-1/r2
MATLAB は一気に正解を出してくれるので助かりますが、導出過程の確認もせずに鵜呑みにするわけにはいきません。人力による確認結果を添付しておきます（ $\rm\TeX$ 化の気力が湧きません。手書きで勘弁願います。）。AI時代、人間に残される仕事はこのようなものしかなくなってしまうのでしょうか？ まるで、AIの奴隷です。
(52)式で、一定値 $R$ が負の値になっているのは、(47)式の周回積分の向きを、プラニメータの仕様とは逆向きの反時計回りとして計算しているためです。このままでは不自然なので、時計回りの実機の仕様に即した新しい変数 $\Delta u_\circlearrowright$, $R_\circlearrowright$ を導入すると、$\Delta u_\circlearrowright{=}{-}\Delta u$, $R_\circlearrowright{=}{-}R{=}1/r_2$ となります。また、目盛車の回転量 $\Delta N_\circlearrowright$ は $\Delta u_\circlearrowright/(2\pi\,r_{\rm w})$ ですから、面積 $\Delta S$ を表す(48)式は、最終的に次のように変形することができます。
$$\Delta S=\Delta u_\circlearrowright/R_\circlearrowright=r_2\,\Delta u_\circlearrowright=2\,\pi\, r_{\rm w}\,r_2\, \Delta N_\circlearrowright=K \Delta N_\circlearrowright\tag{53}$$
ここで、$K$ は(2)式に相当する値、$\Delta N_\circlearrowright$ は(1)式の $(n+m/10)$ に相当する値になります。（めでたく証明完了）
$a$ と $r_1$ はプラニメータの構造上の重要なパラメータですが、最終結果の(53)式からは消し去られています。意外なことに、計測機能には関係しないパラメータのようです。
多くの文書において、$a$ の値を一方的に決めつけて説明している根拠はここにあります。$a{=}r_2$ とすれば、(42),(43)式は $H_{\rm x}{=}\sin\theta_2$, $H_{\rm y}{=}{-}\cos\theta_2$ と大幅に簡略化されるため、$R$ の値は非常に少ない手数で導くことができます。原理を手短かに伝えようとするきには捨て難い論法です。しかし、納得できる十分な説明もなく、天下り的に押し付けられれば、読者はついて行けなくなります。
また、$r_1$ の値が測定結果に関係しないということは、$r_1{=}\infty$ でも良い、つまり、Ｊ点の軌跡は円弧だけとは限らず、直線でも良いことを示しています。これを利用したものが、類似品のリニア プラニメータで、アーム２のＪ点が直線上を移動するような構造になっています。
多くの文書では、まずは、比較的シンプルなリニアプラニメータで原理を説明し、次いで、ポーラープラニメータに移って、軽く流して説明を終えるという論法をとっています。しかし、当記事では、先にポーラープラニメータを詳細に説明し終えているので、リニアプラニメータはその一種に過ぎないと解釈すれば、重ねての原理説明は不要となります。
この記事では、ストークスの定理を使って説明していますが、一般にはグリーンの定理で説明するのが普通です。結局はどちらも同じことを表現しているのですが、個人的にはグリーンの定理にはそれほど馴染みがなかったため、このような結果になりました。グリーンの定理は、ベクトル解析の知識がなくても使える反面、やや抽象的なため、実際の応用にまでは考えが及び難いという気がします。同じ電磁気学の教科書でも、理学系では両定理が併記されているものが多いのに対し、私が学んだ工学系では、ほぼストークスの定理一本で済ませています。
５．おわりに
以上検討の結果、今まで感じていたモヤモヤがかなりスッキリしてきました。しかし、他の方にとっては、数式が多すぎ、気合を入れないとなかなか読み進めない内容になってしまったかもしれません。
世の中には、今回のグリーンの定理系を使った数式偏重ぎみの証明のほかに、数式が比較的少ない幾何学的な証明方法もあります。しかし、これも図や説明を端折られることが多いため、かえって疑問が生じやすく、完全理解はさらに遠くなります。今後、機会があれば、こちらの証明方法でも省略なく詳細に論じてみたいと思います。
また、この記事では、Ｊ点の軌跡は円弧と直線の2種類に限定して説明しましたが、少し捻って考えてみると、もっと自由な曲線でも良さそうです。これについては、鋭意作成中の別稿で紹介するつもりです。
この記事のタグとして「プラニメータ」を登録したところ、Qiita には今まで実績がなかったようで、まさかの一番乗りです。どうも、会員のみなさんには興味が薄いテーマを持ち込んでしまったようです。なにとぞご容赦ください。
６．プログラム
図５～９は下記の記事を利用して MATLAB で作図しています。実舞台での作品としては初披露ですが、本記事には直接の関係はないため、プログラムは割愛します。
MATLABを 三角定規やコンパスのように使って 図を描く
ここに添付しているのは、図４のアニメーション作成用の MATLAB プログラムだけです（Toolbox は使用していません）。試行２で導いた式を利用しています。
■■ 折り畳んであります ■■
% planimeter19.m
% ポーラー プラニメータ のアニメーション
%
%
% 構造
% ＝＝
%
% Ｙ
% ↑ ●Ｔ（トレースの中心点）
% ｜ ↑
% ｜ ｜ アーム２（ＪＴ）長さr2
% ｜ r2 θ2
% ｜ ｜
% ｜ Ｊ●－－－－－－
% ｜ ↑
% ｜ r1 θ1 アーム１（ＯＪ）長さr1
% ｜ ｜
% ●－－－－－－－－－－－－－→Ｘ
% Ｏ(0,0)
%
%
% アーム２の詳細図
% ＝＝＝＝＝＝＝＝
%
% ・－－－－－－→Ｖ トレースの
% ｜ アーム２（全長: r2） 中心
% ｜
% ｜ Ｊ点 Ｔ点
% ｜ ●＊＊＊＊＊＊＊＊●＊＊＊＊＊＊＊／ 略 ／＊＊＊＊＊●
% ｜ a ↑
% ｜ ＪＷ間の距離 Ｗ点
% ｜ （目盛車の接地点）
% ↓ 目盛車はアーム２と直交する方向Ｕの
% Ｕ 移動量で回る。車の半径: r。
% １回転（2πr移動）で100目盛。
%
% 各部の寸法
% r1=164;
% r2=167;
% a=16;
% r=19.05/2; % 目盛車の半径（ノギスでの計測）
clear
close all
% gifアニメーションファイルの作成のための準備
path=mfilename('fullpath'); % このプログラムのファイル名を取得。
% 拡張子なしのフルパス表現。
filename=[path '_01.gif']; % アニメーションgifのファイル名。
delay=1/24; % 一コマの秒数（24コマ/秒）。
% 最初と最後のコマの秒数は別途指定。
% 実機の各部の寸法
r1=164;
r2=167;
a=16;
r=19.05/2; % 目盛車の半径
% アニメーション用の数値
lead_lag=1; % VjがVtより、進みのとき1、遅れのとき-1。
% （Vt: ベクトルＯ→Ｔ、Vj: ベクトルＯ→Ｊ）
n_start=120; % 計測開始点のfree_curve上の要素番号
rt=7; % 照準レンズの半径
% 被計測図形の概形（スプライン曲線の骨格のx,y座標値）
curve_path=[
0 0
1 0
0.8 0.6
1 1
0 1
0.1 0.5
0 0
];
% 被計測図形の大きさや位置の調整
kx=150; ky=150; dx=120; dy=-70;
curve_path=[curve_path(:,1)*kx+dx curve_path(:,2)*ky+dy];
% 被計測図形のスプライン曲線化データの取得
[X,Y]=free_curve(curve_path,-10,2); % ローカル関数の呼び出し。
% 指定された要素番号の点が始点になるように要素の並びを回転シフト。
X(end)=[]; % 閉曲線では始点と終点の座標が重複するので、
Y(end)=[]; % 終点を一旦削除。
X=[X([n_start:end]) X([1:n_start-1]) X(n_start)]; % 並べ替えの
Y=[Y([n_start:end]) Y([1:n_start-1]) Y(n_start)]; % シフト処理
% プラニメータの仕様に合わせ、要素の並びが時計回りになるよう反転。
X=fliplr(X);
Y=fliplr(Y);
% 被測定図形を、一定長の短線分の折線として再構成する。
% 各折れ点を、Ｔ点の移動ステップとする。
[xt,yt]=divide_line_into_fixed_short_segments(X,Y,2);
% ローカル関数の呼び出し。
hf=figure(1); % 作図用の画面を開く
plot(xt,yt,'Color','#0072BD','LineWidth',1); % 被計測図形を描く。
hold on
ths=atan2d(yt(2)-yt(end-1),xt(2)-xt(end-1))+90; % 被計測図形上に
plot(xt(1)+[-15 15]*cosd(ths), ... % 計測開始の位置を
yt(1)+[-15 15]*sind(ths), 'Color','#0072BD'); % 表示する。
axis equal
axis([-50 300 -100 200]);
grid on
xlabel('[mm]');
ylabel('[mm]');
% figureの上下左右の余白の調整。ローカル関数の呼び出し。
figure_trimmer(gcf,[-10 0 -25 -40]);
% 被計測図形を多角形とみなして、面積を数値計算（比較の基準とする）
area=0;
for m=1:length(xt)-1
area=area+cross([xt(m) yt(m) 0],[xt(m+1) yt(m+1) 0]);
end
area=norm(area)/2; % 数値計算された面積
% アーム結合点Ｊの拘束線を描く。
% 原点を中心とした半径r1の円弧
th=linspace(0,360,201); % 大きい円の描画用の全周の分割角度
xr1=r1*cosd(th);
yr1=r1*sind(th);
plot(xr1,yr1,'Color','#D95319','LineWidth',1);
thc=linspace(0,360,51); % 小さい円の描画用の全周の分割角度
angw=0; % 目盛車の回転角度（初期値）
N=length(xt);
for n=1:N % 始点から終点まで各点(xt,yt)について処理。
% (xt,yt)を中心にした半径r2の円と、
% (0,0)を中心にした半径r1の円の交点の座標を計算。
xsq=xt(n)^2;
ysq=yt(n)^2;
if abs(yt(n))>=0.01 % yt≠0 のとき
% yt≠0に対応する、２つの円の交点(xj1,yj1),(xj2,yj2)の式。
% ２つの交点のどちらを採用するかはlead_lagをもとに決定。
RR=r1^2-r2^2+xsq+ysq;
A=1+xsq/ysq;
B=-xt(n)*RR/ysq;
C=(RR/(2*yt(n)))^2-r1^2;
xj1=(-B+sqrt(B^2-4*A*C))/(2*A);
xj2=(-B-sqrt(B^2-4*A*C))/(2*A);
yj1=RR/(2*yt(n))-xt(n)*xj1/yt(n);
yj2=RR/(2*yt(n))-xt(n)*xj2/yt(n);
if sign(lead_lag)*yt(n)<0
xj=xj1; yj=yj1; % アーム結合点Ｊの位置座標
else
xj=xj2; yj=yj2; % アーム結合点Ｊの位置座標
end
else % yt=0 のとき
% yt=0に対応する、２つの円の交点(xj1,yj1),(xj2,yj2)の式。
% この場合、どちらの交点もx座標値は同じ。
% ２つの交点のどちらを採用するかはlead_lagをもとに決定。
xj1=(xsq+r1^2-r2^2)/(2*xt(n));
xj2=xj1;
yj1=+sqrt(r1^2-xj1^2);
yj2=-sqrt(r1^2-xj2^2);
if sign(lead_lag)*xt(n)>0
xj=xj1; yj=yj1; % アーム結合点Ｊの位置座標
else
xj=xj2; yj=yj2; % アーム結合点Ｊの位置座標
end
end
th1=atan2d(yj,xj); % r1アームの角度
th2=atan2d(yt(n)-yj,xt(n)-xj); % r2アームの角度
xw=xj+a*(xt(n)-xj)/r2; % 目盛車の取り付け位置の座標
yw=yj+a*(yt(n)-yj)/r2;
xwm=xj+(a+r)*(xt(n)-xj)/r2; % 目盛車の立面図における
ywm=yj+(a+r)*(yt(n)-yj)/r2; % 紙面への接地点の座標。
xtl=xj+(r2-rt)*(xt(n)-xj)/r2; % 照準レンズの取り付け位置の座標
ytl=yj+(r2-rt)*(yt(n)-yj)/r2;
% r1,r2アームの描画
hp1=plot([0 xj xtl],[0 yj ytl],'Color','#0072BD','LineWidth',3);
% 原点Ｏ、連結点Ｊの描画
hp2=plot([0 xj],[0 yj],'o', ...
'MarkerSize',7,'Color','none','LineWidth',2, ...
'MarkerFaceColor','#0072BD');
text(-5,-13,'O','FontSize',15,'HorizontalAlign','center');
hp2a=text(xj-5,yj+13,'J','FontSize',15,'HorizontalAlign','center');
% 照準レンズの描画
hp3=plot(xt(n)+rt*cosd(thc),yt(n)+rt*sind(thc), ...
'Color','#0072BD','LineWidth',2);
% 照準点の描画
hp4=plot(xt(n),yt(n),'o','Color','#D95319','MarkerSize',2, ...
'LineWidth',2);
hp4a=text(xt(n)+11,yt(n)+12,'T','FontSize',15, ...
'HorizontalAlign','center');
% 目盛車の平面図（紙面への投影図。短い直線）
hp5=plot([xw+r*cosd(th2+90) xw+r*cosd(th2-90)], ...
[yw+r*sind(th2+90) yw+r*sind(th2-90)],'Color','#0b0', ...
'LineWidth',2);
% 目盛車の立面図（円）
hp6=fill(xw+r*cosd(thc),yw+r*sind(thc),'y', ...
'EdgeColor','#0072BD','LineWidth',1);
xtn=xt(n); ytn=yt(n); % 現在(now)のＴ点の位置の座標
xjn=xj; yjn=yj; % 現在のＪ点の位置の座標
xwn=xw; ywn=yw; % 現在の目盛車の位置の座標
if n>1 % nの繰り返しループの初回以外のループ。
dVw=[xwn-xwo ywn-ywo]; % 目盛車の位置の変化ベクトル
thw=atan2d((ytn+yto)/2-(yjn+yjo)/2,(xtn+xto)/2-(xjn+xjo)/2);
% r2アームの角度の区間平均値(度)
nw=[cosd(thw-90) sind(thw-90)];
% 目盛車の自由移動方向の単位ベクトル
dw=dot(dVw,nw); % 目盛車の移動距離の自由移動方向の成分
angw=angw+360*dw/(2*pi*r); % 目盛車の積算回転角度(度)
else % nの繰り返しループの初回ループ。
thw=th2;
end
% 目盛車の立面図での、紙面との接地線
hp7=plot(xwm+((2*pi*r/4)*[-12 1]+ ...
2*pi*r*angw/360)*cosd(th2+90), ...
ywm+((2*pi*r/4)*[-12 1]+ ...
2*pi*r*angw/360)*sind(th2+90), ...
'Color','k','LineWidth',0.5);
% 目盛車の紙面との接地線上の位置目盛
hp8=plot(xwm+((2*pi*r/4)*[-11:0]+ ...
2*pi*r*angw/360)*cosd(th2+90), ...
ywm+((2*pi*r/4)*[-11:0]+ ...
2*pi*r*angw/360)*sind(th2+90), ...
'ok','MarkerSize',2,'MarkerFaceColor','k');
rot=floor(angw/360); % 目盛車の回転回数（整数。ウォームギア表示器
% の表示値のつもり。現在の回転回数が
% rot～rot+1[回転]の範囲にあることを示す）
thwi=thw+180; % r2アームのＴ点からＪ点向きの角度
% 後ろに隠れてしまった目盛車の平面図を前面に引き出す。
uistack(hp5,'top');
% 目盛車の立面図上の四角形の模様の描画
wwx=xw+9*(cosd(thwi+angw+[0 1 2 3 0]*90));
wwy=yw+9*(sind(thwi+angw+[0 1 2 3 0]*90));
hp9=plot(wwx,wwy,'Color','#0072BD','LineWidth',1);
% 目盛車の立面図上の0点目盛の赤線の描画
hp10=plot(xw+[0 10*cosd(thwi+angw)],yw+[0 10*sind(thwi+angw)], ...
'r','LineWidth',2);
% 目盛車の回転回数の数値の表示
hp11=text(xw+27*cosd(thw+100),yw+27*sind(thw+100), ...
['n=' num2str(rot)],'FontSize',15, ...
'HorizontalAlign','center','Rotation',thwi-90);
% アニメの最初と最後には、
% 目盛車の拡大立面図を重ね描きしたコマを追加する。
x0=80;
y0=-45;
if n==1 | n==N
if n==1 % アニメの最初の画面 n==1
thb=0; % 目盛車の回転角度
Pale=[ 90 290 290 90 90
-25 -25 40 40 -25]; % 背景をマスクする領域
else % アニメの最後の画面 n==N
thb=angw;
Pale=[ 90 290 290 90 90
-90 -90 40 40 -90]; % 背景をマスクする領域
end
% 説明文の背景が目立たないようにマスクする。
h19=fill(Pale(1,:),Pale(2,:),[1 1 1],'FaceAlpha',0.6, ...
'EdgeColor','none');
% 目盛車の拡大平面図
h20=plot(x0+[-1 1]*35*cosd(th2+90), ...
y0+[-1 1]*35*sind(th2+90), ...
'Color','#0072BD','LineWidth',4,'Color','#0a0');
% 目盛車の拡大立面図の円全面を黄色にする。
h21=fill(x0+35*cosd(thc),y0+35*sind(thc),'y', ...
'EdgeColor','#0072BD','LineWidth',2);
% 目盛車の拡大立面図の内接四角形の赤味を増す。
h22=fill(x0+35*cosd(thb+th2+[180 90 0 -90 -180 ]), ...
y0+35*sind(thb+th2+[180 90 0 -90 -180 ]), ...
[1 0.85 0],'EdgeColor','none');
% 後ろに隠れてしまった目盛車の平面図を前面に引き出す。
uistack(h20,'top');
% 目盛車の拡大立面図の周辺のr2アームの一部
h23=plot([x0-[35 65 35]*cosd(th2) NaN x0+[35 50 35]*cosd(th2)], ...
[y0-[35 65 35]*sind(th2) NaN y0+[35 50 35]*sind(th2)], ...
'Color','#0072BD','LineWidth',4);
% 目盛読み取りの指標とする三角形
h24=fill(x0-43*cosd(th2)+8*cosd(th2+[0 120 240 0]), ...
y0-43*sind(th2)+8*sind(th2+[0 120 240 0]), ...
[0.8500 0.3250 0.0980],'EdgeColor','none');
% 目盛の読み取り値を表示
h25=text(x0-55*cosd(th2-8),y0-55*sind(th2-8), ...
['m=' num2str(10*mod(angw,360)/360,'%4.2f')], ...
'FontSize',15,'HorizontalAlign','left', ...
'Rotation',thwi-90);
% 説明文
if n==1
h26(1)=text(100,20,'目盛をリセットして再計測','FontSize',15, ...
'HorizontalAlign','left');
h26(2)=text(180,-5,'スタート !!','FontSize',20, ...
'HorizontalAlign','left','Color','#D95319');
elseif n==N
h26(1)=text(100,20,'計測完了 !!','FontSize',20, ...
'HorizontalAlign','left','Color','#D95319');
h26(2)=text(130,-10,'回転回数=(n+m/10)','FontSize',15, ...
'HorizontalAlign','left');
h26(3)=text(200,-30, ...
['=' num2str(rot+mod(angw,360)/360,'%6.3f')], ...
'FontSize',15,'HorizontalAlign','left');
h26(4)=text(145+3,-50,['面積=100×' ...
num2str(rot+mod(angw,360)/360,'%6.3f')], ...
'FontSize',15,'HorizontalAlign','left');
h26(5)=text(180+3,-70,['=' ...
num2str(100*(rot+mod(angw,360)/360),'%6.1f') ...
'[cm^2]'],'FontSize',15,'HorizontalAlign','left');
end
% 目盛車の拡大立面図の詳細目盛
for nn=1:10
if nn==1
% 目盛の0位置を示す赤線
h27=plot(x0+[32 0 10]*cosd(-(nn-1)*36+th2+180+thb), ...
y0+[32 0 10]*sind(-(nn-1)*36+th2+180+thb), ...
'r','LineWidth',4);
end
% 目盛の数値
h28(nn)=text(x0+20*cosd(-(nn-1)*36+th2+180+thb), ...
y0+20*sind(-(nn-1)*36+th2+180+thb), ...
num2str(nn-1),'FontSize',14, ...
'HorizontalAlign','center', ...
'Rotation',-(nn-1)*36+th2+90+thb);
% 目盛の点
h29(nn)=plot(x0+32*cosd(-(nn-1)*36+th2+180+thb), ...
y0+32*sind(-(nn-1)*36+th2+180+thb),'ko', ...
'MarkerSize',3,'MarkerFaceColor','k');
end
% 最初と最後のコマだけに表示するオブジェクトのハンドルをまとめる。
h_add=[h19 h20 h21 h22 h23 h24 h25 h26 h27 h28 h29];
if n == 1 % アニメの最初のコマの画像を書き出す。
drawnow
frames=getframe(hf); % figure画像を取得。
[AA, map] = rgb2ind(frame2im(frames), 256); % 画像形式変換
imwrite(AA, map, filename, 'gif', 'DelayTime', 3, ...
'LoopCount',Inf);
% 表示は3秒間。１コマ目で繰返し再生回数の設定。
delete(h_add);
end
if n==N % アニメの最後の画像として描画されている目盛車の
% 拡大図はいったん見えなくしておく。その前に、
% 拡大図なしのコマを書き出す必要があるため。
for nn=1:length(h_add)
h_add(nn).Visible='off';
end
end
end
% 用済みの現在値を、前回値として保存。
xto=xtn; yto=ytn; % 前回(old)のＴ点の位置の座標
xjo=xjn; yjo=yjn; % 前回のＪ点の位置の座標
xwo=xwn; ywo=ywn; % 前回の目盛車の位置の座標
drawnow
% gifアニメーションファイルとして出力
frames=getframe(hf); % figure画像を取得。
[AA, map] = rgb2ind(frame2im(frames), 256); % 画像形式変換
imwrite(AA, map, filename, 'gif', 'DelayTime', delay, ...
'WriteMode', 'append');
% ２コマ目以降には、"追記"の指示が必要。
if n==N % 目盛車の詳細図を重ね描きして最後の一コマにする。
% いったん隠してあった目盛車の拡大図を見えるようにする。
for nn=1:length(h_add)
h_add(nn).Visible='on';
end
drawnow;
% 拡大図つきの最後の１コマを書き出す。表示時間は7秒間。
frames=getframe(hf); % figure画像を取得。
[AA, map] = rgb2ind(frame2im(frames), 256); % 画像形式変換
imwrite(AA, map, filename, 'gif', 'DelayTime', 7, ...
'WriteMode', 'append');
end
if n<N % コマごとに変化する画像は、そのつど削除。
delete([hp1 hp2 hp2a hp3 hp4 hp4a hp5 hp6 hp7 hp8 hp9 hp10 hp11]);
end
end
% 事前の数値計算の面積と、プラニメータによる面積を比較
disp(' ');
disp(['多角形としての数値計算面積：' num2str(area) '[平方mm]'])
disp(['プラニメータによる計測面積：' ...
num2str((angw/360)*2*pi*r*r2) '[平方mm]'])
disp(' ');
% ====================================
% スプライン曲線で、歪んだ円や開曲線を描くためのローカル関数。
% ====================================
function [X,Y]=free_curve(curve_path,angs,kang)
% スプライン曲線で描いた歪んだ円や開曲線の座標を取得する。
% 始点,終点に同一点を指定したときは閉曲線になる。
% 【入力】
% curve_path: 線を通過させる点の[x,y]座標の対を縦に並べた配列。
% [x1,y1 ; ... ; xn,yn]。
% angs: 始点/終点での線の傾き [degree]。NaNのときは自由端。
% スカラーのときは、始点/終点とも同一の傾き。
% 2要素の配列のときは、それぞれ始点,終点での傾き。
% kang: angsの効果が及ぶ範囲の広さ（適量の目安は1.0）
% 【出力】
% [X,Y]: 結果として得られた折線のx座標列とy座標列
if length(angs)==1
ang1=angs(1);
ang2=ang1;
else
ang1=angs(1);
ang2=angs(2);
end
if isnan(ang1)
dx1=0;
dy1=0;
else
dx1=cosd(ang1);
dy1=sind(ang1);
end
if isnan(ang2)
dx2=0;
dy2=0;
else
dx2=cosd(ang2);
dy2=sind(ang2);
end
N=size(curve_path,1); % 線の経路を指定する点の数
% スプライン曲線の長さの概算見積もり
xx=[];
yy=[];
for n=1:N
xx=[xx curve_path(n,1)];
yy=[yy curve_path(n,2)];
end
L=cumsum(hypot(diff(xx),diff(yy)));
len=L(end);
% 媒介変数
t=linspace(0,1,N)*len*kang;
xo=[dx1 xx dx2];
yo=[dy1 yy dy2];
% スプライン曲線として平滑化
Lsp = spline(t,[xo;yo]);
XY = ppval(Lsp, linspace(0,max(t),201));
X=XY(1,:); % 平滑化された線の
Y=XY(2,:); % 折点のx,y座標値
end
% =============================================-
% 折線を「等長の短線分」の連なりとして再構成するためのローカル関数。
% =============================================-
function [Xd,Yd]=divide_line_into_fixed_short_segments(X,Y,ds)
% 提示された折線（旧折線）を「等長の短線分」の連なりとして構成された
% 新折線に変換する。ただし、旧折線の終端に残った線分だけは、長さを
% 他と揃えることはできず、短い長さになる。
% 【入力】
% X,Y: 旧折線の折点群のx,yデータの配列。
% ds: 入力された折線は、この長さの線分で等分割された折線に
% 変わる。ただし、旧折線の終端に残った線分だけは長さを
% 揃えることはできず、ds未満の値になる。
% 始端と終端だけは旧折線と新折線で共用されるが、その他の
% 折点については、両折線間で関連性はない。
% 新折点は、旧折線の線上には乗るが、旧折点が新折線の
% 線上に乗るとは限らない。乗らないことの方が多い。
% 【出力】
% Xd,Yd: dsの長さに等分割された新折線の折点群のx,yデータの配列。
% 処理済みの旧折線データを逐次消去更新できるよう、別の変数にコピー
xc=X;
yc=Y;
Xd=X(1); % 新折線の折点群のx,y座標を並べた配列の初期値。
Yd=Y(1); % 折線の始点は、旧折線と新折線で同じ。
while true
x1=xc(1); % 旧折線の、処理中の線分の始点の座標
y1=yc(1);
x2=xc(2); % 〃 終点 〃
y2=yc(2);
dx=x2-x1; % 旧折線内の処理中の線分の
dy=y2-y1; % 始点・終点間の座標値の隔たり。
L21=hypot(dx,dy); % 旧折線内の処理中の線分の（残り）長さ。
if L21>ds % dsが旧折線の処理中の線分の残り長さより短い
% とき（次の新折点は残り線分の上にある）
% P1=(x1,y1), P2=(x2,y2), Px=(x3,y3)
%
% ds= |Px-P1| 新折線を構成する一本の短線分の固定長。
% L21=|P2-P1| 旧折線の処理中の線分の長さ
%
% P1 Px P2
% ● ds ● ●
% ↑ 新折線の ↑
% ｜ 次の折点 旧折線の
% ｜ 処理中の線分の終点
% 旧折線の処理中の
% 線分の始点。
% （P2,P1間の長さが、旧折線の本来の線分の長さよりも短い
% 残り線分のときは、新折線の決定済の最新折点P0と等しい）
x3=x1+dx*ds/L21; % 新折線の次の折点のx座標
y3=y1+dy*ds/L21; % 〃 y座標
Xd=[Xd x3]; % 新折点の座標を新折線の配列に追加する。
Yd=[Yd y3];
xc(1)=x3; % 旧折線の処理中の線分の始点を新折点に変更。
yc(1)=y3;
else % dsが旧折線の処理中の線分の残り長さより長い
% とき（次の新折点は旧折線の残り線分上には
% 存在しない）
% 旧折線の、今まで処理対象にしていた線分上には、探している新折点
% は存在しないことが分かったので、旧折線の処理対象線分を次の
% 線分に移動する。
% この場合、新しい折点を求める計算が少々複雑になる。この様子を
% 下図に示す。
% ケース１の場合は、そのまま、余弦定理を用いて新折点を決定でき
% るが、ケース２の場合には、旧折線の処理対象線分をさらに先に
% 移動して、その始点を新たなP1として置き直し、ケース１の状態
% になるまで繰り返す。
%
% P0=(x0,y0), P1=(x1,y1), P2=(x2,y2), Px=(x3,y3)
%
% ds= |Px-P0| 新折線を構成する一本の短線分の固定長。
% L1= |P0-P1| P0-P1の長さ
% Lx= |Px-P1| Px-P1の長さ
% L21=|P2-P1| 旧折線の処理中の線分の長さ
% θ= ∠P0.P1.P2
% = ∠P0.P1.Px
%
% 余弦定理 ds^2=L1^2+Lx^2-2*L1*Lx*cosθ より、Lxを求める。
%
% ==== ケース１（新折線の次の折点はPxで決定） ====
%
% P1: 旧折線の次の処理対象線分の始点
% ●
% L1 Lx Px: P2-P1線分上にあって
% ● ds ● P0点からdsだけ
% P0: 新折線の 離れた点
% 決定済の最新折点 ● P2
% P1点に対応する終点
%
% ==== ケース２（新折線の次の折点は不明。P2よりも先にある） ====
%
% P1
% ●
% L1 ● P2
% P0 ● ds ● Px
%
%
% 旧折線の処理中の線分の終点が、旧折線の終端に達していた場合。
% この先に固定長dsを満たす点はないので、新折線の終点を旧折線の
% 終点と一致させて処理を終わる。
if x2==X(end) & y2==Y(end)
Xd=[Xd X(end)];
Yd=[Yd Y(end)];
return
end
while true
xc(1)=[]; % 旧折線の用済みの折点（そこを始点とする
yc(1)=[]; % 線分上に新折点は存在しない）を削除。
% その結果、旧折線の未処理の全折点xc,ycの
% 番号は一つだけ繰り上がる。
x1=xc(1); % 旧折線のこれから処理する線分の始点の座標
y1=yc(1);
x2=xc(2); % 〃 終点 〃
y2=yc(2);
dx=x2-x1; % 旧折線のこれから処理する線分の始点・終点
dy=y2-y1; % 間の座標値の隔たり。
L21=hypot(dx,dy); % 旧折線のこれから処理する線分の（残り）長さ。
x0=Xd(end); % 新折線の発見済の最新の折点の座標
y0=Yd(end);
% 上図のL1の値
L1=hypot(x1-x0,y1-y0);
% 上図のθ（ここではth）
th=atan2d(dy,dx)-atan2d(y0-y1,x0-x1);
th=abs(mod(th+180,360)-180); % 0～+180°の範囲に変換
% 上図のLx
Lx=L1*cosd(th)+sqrt(L1^2*(cosd(th)^2)+ds^2-L1^2);
if Lx<L21 % 上図の「ケース１」の場合
x3=x1+dx*Lx/L21; % 新折線の、新たに決定した折点。
y3=y1+dy*Lx/L21;
Xd=[Xd x3]; % 新たな折点を、新折線の配列に追加する。
Yd=[Yd y3];
xc(1)=x3; % 旧折線の処理中の線分の始点を、新折線の
yc(1)=y3; % 新たな折点に変更。
break;
else % 上図の「ケース２」の場合
% 旧折線の処理中の線分の終点が、旧折線の終端に達していた場合。
% この先に固定長dsを満たす点はないので、新折線の終点を旧折線の
% 終点と一致させて処理を終わる。
if x2==X(end) & y2==Y(end)
Xd=[Xd X(end)];
Yd=[Yd Y(end)];
return
else % まだ旧折線の終端には達していないとき。
continue; % 旧折線の処理対象線分を一つ先に進める。
end
end % end of if Lx<L21
end % end of while true
end % end of if L21>ds
end % end of while true
end % end of function
% ==========================================
% figureの上下左右の余白調整用のローカル関数
% ==========================================
function figure_trimmer(hf,KK)
% 【入力】
% hf: 余白調整したいfigureのハンドル（呼び出し側では常にgcf固定）
% KK: 余白の変更量のベクトル。 単位は[pixel]。
% （[上余白 下余白 左余白 右余白]。増量は+、減量は-）
%
% 【出力】
% 戻り値は無し
posf=hf.Position; % 余白変更前のfigureのposition
% [左端 下端 幅 高さ]、単位は[pixel]。
naxe=size(hf.Children,1); % figure内のChildren（axesやcolorbar）
% などの総数。
% 全Childrenのposition行列を作成
% 各行がそれぞれのChildrenに対応。
% 各列は[左端 下端 幅 高さ]、単位は[normalized]。
Posa=[]; % position行列の初期値
for n=1:naxe
axes(hf.Children(n))
hn(n)=gca;
Posa=[Posa ; hn(n).Position];
end
% 各Childrenの、余白変更前のfigure内でのpixel単位でのPositionを計算。
Posp(:,[1 3])=Posa(:,[1 3])*posf(3);
Posp(:,[2 4])=Posa(:,[2 4])*posf(4);
% 余白変更後のfigureのpixel単位でのPositionを計算。
posfx=posf+[-KK(3) -KK(2) KK(3)+KK(4) KK(1)+KK(2)];
% 余白変更後のfigure内での各Childrenのpixel単位でのPositionを計算する。
Pospx(:,1)=Posp(:,1)+KK(3);
Pospx(:,2)=Posp(:,2)+KK(2);
Pospx(:,3)=Posp(:,3);
Pospx(:,4)=Posp(:,4);
% 余白変更後の各Childrenのnormalized単位でのPositionを計算する。
Posax(:,[1 3])=Pospx(:,[1 3])/posfx(3);
Posax(:,[2 4])=Pospx(:,[2 4])/posfx(4);
% figureとChildrenのpositionを更新
hf.Position=posfx;
for nn=1:naxe
hn(nn).Position=Posax(nn,:);
end
end