更新履歴

2025/03/28 おまけとしてJupyter Notebookと周辺ファイルを作成．
2025/03/27 公開

前略

まあ，まずはこれを見てみてくれ；

こいつを作った．この記事では，これを作った経緯と導出過程を説明する．

Qiitaでの初めての投稿です！内容が内容なのでどこにも公開する場所ないし，なんかいいところないかなーって思ってQiitaにたどり着きました．
たいよろ！

前例

Googleでおっぱい　関数と検索してみると，トップにDesmosのグラフが出てきた；

Qiitaでも同じようなことをやっている記事があった；

どれもとてもいい形状に見えるし，ぱっと見て「うおおおおおっぱい！」と言いたくなるような形だ．魅力的だし，時点$t$をいじれば揺れるやつなんか，アイデアに富んでてすごい，
一方で，これらの数式は非常にアドホックな方法で構築されていて，論理的ではない．その数式のその項はいったいどんな意味を持って構築されているのか質問したくなる．

今度はmath breast functionと検索してみたら，Quoraにてカテナリー曲線を用いたモデルが出てきた；

そうそう！こういうやつ！物理的な根拠に基づく，こういう導出過程を踏んで構築されたもののほうが，私としては意味のある数式に見えるし，説明力もある．
でも，この形は全然魅力的じゃない．おっぱいとは，もっと魅惑的で，むっちりとしているはずだ．分析モデルとしては有用だと思うけど，こんなのは全然おっぱいじゃないよね．

3Dのものも見たいと思って，3D breast functionとかで調べたら，こんな式が出てきた；

$$z=\exp(-((x-4)^2+(y-4)^2)^2/1000) + \exp(-((x+4)^2+(y+4)^2)^2/1000) + 0.1\exp(-((x+4)^2+(y+4)^2)^2)+0.1\exp(-((x-4)^2+(y- 4)^2)^2)$$

……なんだこれは．たしかに3Dだ．2つの膨らみもある．でも，これはただ単にガウス関数を足し合わせたものであって，おっぱいとは似て非なるものだ．

あとおもしろいのもあった，それは

$$x=1$$

これは2次元でも3次元でも共通する数式であり，非常に簡潔だ．小さければ頂上の高さは$0$，勾配も$0$．数式的にも美しいし，いわゆる「絶壁」を明確に表現している．

しかし，私は大きいおっぱいのほうが好きだ．巨乳派なのである．
(あと，たっぱもおしりもでかいほうが私は好きなのだが，本稿では弁明を差し控える)

モチベーション

前項で挙げた条件は満たしたい．したがって論理的な導出過程と魅力的な形状を両立し，さらに3次元で大きいおっぱいを表現したいと考えている．

したがって，この条件を換言すれば以下のような数式を構築することになる．

物理学に基づく導出
現実的な力学仮定
パラメトリックな方程式
魅力的な形状
3次元での表現
大きいおっぱい！←←←←←←←

一方で，魅力的なおっぱいは複雑な数式になることは，前項の例を見ても明らかだ．なので，解析的にではなく数値的に形状を導出することも考えた(例えば有限要素法による数値シミュレーション)．
しかしながら，私は数式でおっぱいを表したいのだ．だからこそ私はおっぱい　関数などとGoogleで検索したわけだ．

ただ……大きな問題がある．私は物理学者ではないし，物理を専門的に学んだこともない．したがって，基本のところから着実に攻めて行こうと思う．

前提条件

乳首について

多くの数式をザッピングしていくと，魅力的な数式によるおっぱいは乳首が付いていて，魅力的ではない数式のおっぱいは乳首が付いていない．また，乳首が付いている関数は乳首項というおまけ的な項が加算されていて，おおよそ以下の形で表現されている；

$$y \propto t_n\exp{\left( -\frac{(x-x_n)^2}{2c^2} \right)}$$

この項は通称ガウスの乳首関数というようだ．実際，$a=\frac{1}{\sqrt{2\pi \sigma^2}}, c^2=\sigma^2$とすれば，ガウス関数の定義そのままだ．

一方で，そのままのガウス関数はちょっと面倒な点がある．
もしも乳首を細くしたいとして$c$を小さくする．すると$\frac{1}{\sqrt{2\pi\sigma^2}}$が大きくなる．一方で$\exp{\left( -\frac{(x-x_n)^2}{2\sigma^2} \right)}$の項は$x \rightarrow\infty,-\infty$で最小$0$, $x=x_n$で最大$1$だから，結果的に「乳首を細くしたいだけなのに，乳首の高さも変わってしまう」問題が発生する．これはパラメトライゼーションとしてはあまり適切ではない．したがって$t_n$の項は定数のほうが好ましい．
また，仮に少し角を立たせたいとしても，この関数形ではコントロールできない．

したがって，以下の基本乳首関数を導入する．

$$\tag{1} n(y,z,t,s,k)=t\exp{\left(- \left(\frac{y^2+z^2}{s^2}\right)^k \right)}$$

すなわち，元の関数に指数$k$を追加したものになる．また，この関数は明らかに$y=z=0$で最大値$t \cdot 1=t$を取る．
もし$y$軸方向に$y_{\text{shift}}$だけ平行移動したいのであれば，$x=n(y-y_{\text{shift}},z,t,s,k)$を評価してやれば良い．

脂肪体について

もし乳首と本体，すなわち脂肪体を分けて評価するなら，脂肪体の関数形を定める必要がある．

脂肪体の関数形についてなんかいいアイデアないかなーと探していたら，ワコール人間科学研究所が著者の『～重力からバストを守る～「重力によるバストの動きと皮膚研究」』を見つけた；

このPDFの第5セクションにて，無重力状態でのバストの3D計測が実施されている．このPDFに明示的に書かれているように

無重力下では～(略)～バスト形状は、半球状
の丸いきれいな形をしている。

と書かれている．すなわち，どんな大きなおっぱいさえも，無重力なら半球状に近似でき，逆説的に言えばあの魅惑的な形状は重力によって形作られていることになる．

すなわち

脂肪体は無重力なら近似的に半球として扱える
脂肪体は重力によって変形してあの形状になる

半球の数式は簡単だ．半径$1$の球の数式は
$$x^2+y^2+z^2=1$$
なのだから，その半分を表せばいい．都合により$x$について解くと，
$$x=\pm\sqrt{1-(y^2+z^2)}$$
となる．これもまた都合により$x \ge 0$の半球を考えると，基本脂肪体関数は

$$\tag{2} f(y,z)=\sqrt{1-(y^2+z^2)}$$

と記述できる．

次は重力について考えたい．しかしながら，私は前述の通り物理学しょしんしゃだ．過去に習ったはずの運動方程式もすっかり忘れてしまったし，運動エネルギーってなんなん？レベルだ．そもそも専門が違うため，どうすればいいかよくわからなかった．

しかし，ふと脳裏に電撃が走った．私は重力について考えたいのではなく，重力による変位を考えたいのだ．ある点$P=(x,y,z)$について，それだけ$z$方向に動くかどうか，それをモデル化すれば良い．

いろいろ考えを巡らせていくと，このモデルは「$x=0$の点で固定された長さ$L$の棒が，どんな$z$軸上で変位をもたらすか？」を説明するモデルだと気づいた．私はこの旨をcopilotくんに入力して，本稿のキモとなる断面曲線方程式に出会った．

1. 棒モデル(Stick Model)

1.1. 断面曲線方程式の導出

棒の長さを$L$,重量を$m$とする．この時，質量密度$\rho$は
$$\rho=\frac{m}{L}$$
で表される．こいつを使うと，重力分布$q(x)$が
$$q(x)=\rho g = \frac{mg}{L}$$
で表される．$mg$は物体に作用する力であるから，重力分布とは，単位長さあたりの物体に作用する力とも解釈できる．

ヤング率$E$, 断面二次モーメント$I$を定義する……え？なにそれ？

ヤング率ってなんぞや

ヤング率とは「物体が受ける力と変形量の比例関係」らしい．
調べてみたところ，材料ごとに決まっているようだ．すなわちこれは「素材の強さ」を表しているに違いない．

断面二次モーメントってなんぞや

断面二次モーメントとは，「断面形状によって決まる部材の曲げにくさ」らしい．とりあえず，今考えているのはただの棒だから，定数だろう．

……ともあれ，これは物性値のようなものだろう．一旦は定数として扱えるはずなので，そのままにしておこう．

すると，たわみ量を表すたわみ関数$v(x)$は以下の断面曲線方程式と呼ばれる微分方程式の解である；

$$\frac{d}{dx} V(x)=\frac{d^2}{dx^2}M(x)=\frac{d^2}{dx^2}\left( EI \frac{d^2}{dx^2} v(x) \right) = q(x)$$

ここで，$V(x)$はせん断力，$M(x)$は曲げモーメントを表す……えっなにそれ？

せん断力ってなんぞや

調べてみたところ，一言で言うなら，「棒の中でずれようとする力」のことらしい．おそらく，実際はずれは発生しないけど，内部的には力が発生していて，それがせん断力ということなのだろう．

曲げモーメントってなんぞや

これも調べてみたところ，「棒を曲げようとする力の大きさ」のことらしい．これは非常に直感的だ．なぜなら，棒を曲げようとしたら内部的には形状を歪めようとする力が発生しているはずだからだ．

さて，とりあえず微分方程式を解いていこう．本稿では解き方については説明しないが，やっていることといえば愚直に積分を繰り返すだけだ．解は以下の通りである；
$$v(x) = \frac{mg}{24LEI}x^4+\frac{C_1}{6EI}x^3+\frac{C_2}{2EI}x^2+C_3x+C_4$$

そう，積分しただけだから，積分定数がくっついている．こいつらをなんとかしなければ，実際には使えない．

境界条件というものが微分方程式に関する問題には定められていることが多い．固定された棒の場合の境界条件は以下の通り

$v(0)=0$
$\frac{d}{dx}v(x) \big|_{x=0}=0$
$\frac{d^2}{dx^2}v(x) \big|_{x=L}=0$
$\frac{d^3}{dx^3}v(x) \big|_{x=L}=0$

こいつらについて意味を確認しておこう．

1.について

$v$が$x$の変化による$z$の変化量であることを思い出そう．これについて，$v(0)=0$というのは，$x=0$では$z$の変位は$0$ということだ．つまり，「根本では棒は固定されている」という条件になる．

2.について

$v$が$x$の変化による$z$の変位なら，$\frac{d}{dx}v(x)$は$z$の変位の変化量を表している．すなわち，$x=0$の点では$v(x)$の傾きは$0$．つまり「根本では棒は水平である」という条件になる．

3.について

断面曲線方程式を再びにらめっこする．すると以下の関係が導かれるだろう；
$$\frac{d^2}{dx^2}v(x)=\frac{M(x)}{EI} \propto M(x)$$
つまり，たわみ関数の二階微分は曲げモーメントに比例する．$x=L$でこれが$0$ということは，分子が$0$になることに他ならない．つまり，この条件は「端っこでは棒を曲げようとする力はない」という条件になる．

4.について

もう１回断面曲線方程式を確認すると，以下の関係もあるだろう；
$$\frac{d^3}{dx^3}v(x)=\frac{d}{dx}\frac{M(x)}{EI}=\frac{V(x)}{EI} \propto V(x)$$
つまり，たわみ関数の三階微分はせん断力に比例する．$x=L$でこれが$0$ということは，この条件は「端っこでは棒の中でずれようとする力はない」という条件になる．

特に，3,4,の条件は自由端条件というらしい．覚えておこう．

さて，未知変数は4つ，条件も4つということで，積分定数$C_1,C_2,C_3,C_4$を具体的に計算できる．これも計算していくと，求めるたわみ関数
$$\tag{3} v(x)=\frac{mg}{24LEI}x^4-\frac{mg}{6EI}x^3+\frac{mgL}{4EI}x^2$$
が定まった．

1.2. 数式の組み立て

とりあえずパーツは集まった．

基本脂肪体関数(再掲)

$$\tag{2} f(y,z)=\sqrt{1-(y^2+z^2)}$$

たわみ関数(再掲)

$$\tag{3} v(x)=\frac{mg}{24LEI}x^4-\frac{mg}{6EI}x^3+\frac{mgL}{4EI}x^2$$

基本乳首関数(再掲)

$$\tag{1} n(y,z,t,s,k)=t\exp{\left(- \left(\frac{y^2+z^2}{s^2}\right)^k \right)}$$

この関数を組み合わせて，おっぱいを形作っていこう．

まずは3次元上に脂肪体を描画しよう．つまり，
$$x=f(y,z)$$
ここまではただ単に半球を3次元上に描画しただけ．

今度はこれに対してたわみを付加しよう．たわみ関数$v(x)$は$z$の変位量を$x$で説明するものだ．ここで$y$軸方向に$y_{\text{shift}}$だけ移動させたいなら，関数の引数の$y$の代わりに$y-y_{\text{shift}}$を与えてやればよかった．今は$z$軸方向に$v(x)$だけ移動させたいのだから，引数の$z$の代わりに$z-v(x)$を与えてやればよい．つまり，
$$x=f(y,z-v(x))$$
となる．また，想定長さは$L$なのだから，$x$軸方向に$L$倍しよう．すなわち，
$$x=Lf(y,z-v(x))$$
となる．最後に乳首関数を使って，乳頭と乳輪を描こう．第二項は乳頭で，第三項は乳輪だ．乳首を描画する$y,z$座標を$y_{\text{n}},z_{\text{n}}$としよう．先程の平行移動と組み合わせれば，
$$x=Lf(y,z-v(x))+n(y-y_{\text{n}},z-v(x)-z_{\text{n}},t_1,s_1,k_1)+n(y-y_{\text{n}},z-v(x)-z_{\text{n}},t_2,s_2,k_2)$$

……完成だ．私の第一号のおっぱい関数！

1.3. 結果

……出来はともかくとして，とりあえず結果を見て欲しい．
使用したパラメータは以下の通りである；

\left( \begin{array}\\
    L\\
    m\\
    g\\
    E\\
    I\\
    t_1\\
    s_1\\
    k_1\\
    t_2\\
    s_2\\
    k_2\\
    y_n\\
    z_n
    \end{array} \right) = \left( \begin{array}\\
    1\\
    0.006\\
    -9.8\\
    0.1\\
    0.1\\
    0.07\\
    0.1\\
    2.5\\
    0.05\\
    0.25\\
    3.3\\
    0\\
    0.24
    \end{array} \right)

……なんだこの饅頭は．

いろいろ突っ込みどころがある．しかし確実に言えることがある．
これはおっぱいではない．

球体としての形状を失っている
そもそも変形のしかたがおかしい

特に，根本の下端の形状が気になる．

こいつを見てみると，$x=0$の時点で，重力によって脂肪体が傾いていない．図の赤線がそれを示している．

私が思う大きなおっぱいは，紫色のラインのような変位を持つはずである．すなわち

$x=0$では明らかに負の傾きを持つ
ある$x$において変位の変化量，すなわち$\frac{d}{dx}v(x)$は正に転ずる(可能性がある)．

原因1. 境界条件の設定ミス

私は境界条件に$\frac{d}{dx}v(x) \big|_{x=0}=0$を仮定したのだから．仮定どおりの動きをしていると思えば，至極当然である．

1 .について解決するのは簡単だ．境界条件として，何らかの定数$d$を用いて

$$\frac{d}{dx}v(x) \big|_{x=0}=d \neq 0$$

を仮定してやれば良い．

原因2. 曲げモーメント条件を満たさない

2 . の条件は物理学的にはどんな意味があるのだろうか？

1階微分が負から正に変化する関数は，2階微分は正だ．ここでたわみ関数の2階微分が何者だったかを再び考えると，曲げモーメント$M(x)$を$EI$でスケーリングしたものだった．$E,I$は両方とも正だから，$M(x)$が正であることを暗黙的に示している．

ちなみに，今のたわみ関数における曲げモーメントは
$$M(x)=EI\frac{d^2}{dx^2} v(x)=\frac{mg}{2L}x^2-mgx+\frac{mg}{2}L$$
である．$g$が負であることを考えれば，この関数は$x$について上に凸の2次関数だ．つまり，$\frac{d}{dx}M(x)=0$の点が最大値を取る．具体的に解いてみると，
$$\frac{mg}{L}x=mg \Leftrightarrow x=L$$
つまり，端っこで曲げモーメントが最大．この時の曲げモーメントは
$$M(L)=\frac{mg}{2}L-mgL+\frac{mg}{2}L=0$$
したがって，最大が$0$の上に凸の関数は，その地点から離れた瞬間に負になる．つまり．このモデルは私の求める要件を満たしていない．

では，曲げモーメント$M(x)$が常に正というのは具体的にはどういう条件なのだろうか？曲げモーメントはせん断力の積分として定義されていて，せん断力は重力分布の積分として定義されていた；
$$\frac{d}{dx}M(x)=V(x) \Leftrightarrow M(x)=\int V(x) dx =\int \int q(x) dx dx$$
今のモデルについて考えると，$q(x)=\frac{mg}{L} \lt 0$なのだから，いくら積分をしたところで，負の値の積分は負．そしてその値の積分は負だ．したがって，なんとかして重力分布を正にするような形式に書き換える必要がある．しかし，正の重力分布というのはいわば反重力や外力に対応するものだ．これはどういった理由付けが可能かと考えたところ，以下の論文を見つけた．

この論文によれば，乳房には重力のほかに，皮膚による張力が存在し，重力だけでは垂れ下がってしまうはずの脂肪体を，張力によって引っ張り上げているというのだ．すなわち，おっぱいには常に反重力が発生している．(は？)

……もとより，とりあえず重力分布$q(x)$については更改の余地がある．

あと，

正面から見たら丸すぎる

問題も考えられる．

原因3. 脂肪体関数の定式化の誤り

これは単純に，WacoalのPDFに従って，脂肪体を完全な半球としてしまったからだと考える．つまり，脂肪体関数についても修正を行う必要があるようだ．

2. 変動重力分布モデル(Variable Gravity Distribution Model)

2.1. 脂肪体関数の調整

棒モデル(Stick Model)では，脂肪体関数は完全な半球だった．しかし，魅力的なおっぱいとしては完全な半球じゃないほうが良い．
(あと，個人的にはなんかもうちょっと細長くて，斜めに歪んだ形のほうが魅力的だ)．

したがって，脂肪体関数を以下の形状に変えよう；

$$\tag{4} f(y,z)=\sqrt{1-((y_sy)^2+z^2+ay_syz)}$$

ここで，$a$は歪みパラメータ，$y_s$は$y$軸拡縮パラメータである．

一応，断面がどんな変化をするか確認してみよう．もし$y_s=1, a=0$なら，$x=0$において，断面は半径$1$の真円になり，その面積は$\pi$だ．おっぱいだけにね．

まずは$a \neq 0$としてみよう．この時，$x=0$において断面の式は

\begin{align*}
0 &= \sqrt{1-(y^2+z^2+ayz)}\\
\Leftrightarrow y^2+z^2+ayz &= 1
\end{align*}

こいつをよーく確認してみると，$a$を動かすと，$z=y$方向に引き伸ばされたり，$z=-y$方向に引き伸ばされたりするようだ．そして，楕円の長軸もしくは短軸は$z=y$ないし$z=-y$の交点と原点との距離に等しい．

これも真面目に計算すると，$y \ge 0$だけに解を制限すれば，
$$y=\sqrt{\frac{1}{2+a}},\sqrt{\frac{1}{2-a}}$$
という解が得られるはずだ．これはただの交点．実際は原点との距離なので，直角2等辺三角形の辺の比を使えば，長軸，短軸はそれぞれ

$$r_{\text{long}}=\sqrt{\frac{2}{2-a}},r_{\text{short}}=\sqrt{\frac{2}{2+a}}$$
となり，楕円の面積の公式から，
$$r_{\text{long}}r_{\text{short}}\pi=\frac{2\pi}{\sqrt{4-a^2}}$$
となる．

こんどは$y_s \neq 1$としよう．でも，これって単純に$y$軸方向に$\frac{1}{y_s}$倍するだけだ．つまり，$x=0$における断面の面積は
$$\frac{\pi}{y_s}$$
となる．

また，脂肪体関数を変更するのに合わせて，乳首関数も変えておこう．具体的には，
$$\tag{5} n(y,z,t,s,k)=t\exp{\left(- \left(\frac{(y_sy)^2+z^2+ay_syz}{s^2}\right)^k \right)}$$
となる．

2.2. 重力分布の調整

2.2.1. 変動重力分布

棒モデルにおいては，重力分布は$\frac{mg}{L}$と，$x$の地点によらず定数だった．
しかし，実際は棒を変形させているわけではなく，球体を変形させている．
換言するならば，$yz$平面と平行に輪切りにしていった時，断面の単位厚さあたりの重量は$x$によって異なるべきだ．
また，輪切りの厚さが微小なら，その重量は断面の面積に比例するはずだ．

つまり，断面の面積を重力分布に組み込むことで，$x$の地点ごとの重さの違いを表現しようという試みである．

簡単化のために半径$1$の球を，$x$軸に並行な方向に$L$倍だけ長くした楕円体を考えよう．つまり，こんな式で表される楕円体を考える．

$$\left(\frac{x}{L}\right)^2+y^2+z^2=1$$
これを輪切りにしていった時，$x=0$の地点では面積は半径$1$の円の面積，すなわち$\pi$になる．逆に言えば，半径は$1$だ．

楕円体を$y=0$を所与のものとして，$z$について解こう．半径について考えたいから，再び非負の範囲の$z$だけ考えればいいから，
$$z=\sqrt{1-\left(\frac{x}{L}\right)^2}$$
となる．半径が求まれば，あとは$2$乗して$\pi$をかければ良い．ただ，正直$\pi$は邪魔だ．おっぱい関数を求めている私がいうのもアレだが，$\pi$は定数倍しているだけで，本質的な成分ではない．したがって，$x$が任意の地点の断面積$S(x)$は以下の関係を満たしている，ということを表すと；

$$S(x) \propto 1- \left(\frac{x}{L} \right)^2$$

となる．

断面あたりの重力分布と比例しているという仮定のもとでは，新たな重力分布は
$$q(x)=\frac{mg}{L} \left( 1- \left(\frac{x}{L} \right)^2 \right)$$
となる．ただ，調整度合いをパラメータとして扱いたいから，パラメータ$\delta$を導入して，
$$q(x)=\frac{mg}{L} \left( 1- \delta \left(\frac{x}{L} \right)^2 \right)$$
としよう．

2.2.2 反重力項

重力分布が$x$によって可変になったのだが，まだ曲げモーメント条件を満たさない．曲げモーメント条件を満たすには，重力分布関数に反重力項を加える必要がある．

脂肪体の先端の部分に$1$の力がかかっているとする．力の伝達は指数関数で近似できると仮定すると，例えば以下のような反重力項が考えられる；

$$\exp \left( \frac{x}{L}-1 \right)$$
これを重力分布関数に突っ込みたいが，元の$q(x)$とスケールを合わせたい．したがって私はアドホックに以下の形にすることにした(正直，ここに関しては再考の余地がある．なんかいいアイデアがあったら教えて下さい)；

$$\frac{\gamma}{L^2}\exp \left( \frac{x}{L}-1 \right)$$
ここで$\gamma$はパラメータである．

2.3 最終的な重力分布関数とたわみ関数

上の調整を組み合わせて，新たな重力分布関数
$$q(x) = \frac{mg}{L} \left( 1-\delta \left(\frac{x}{L}\right)^2 \right) + \frac{\gamma}{L^2} \exp \left( \frac{x}{L}-1 \right)$$
を定義する．こいつを使って，断面曲線方程式を解いていこう．解き方については先程言った通りなので省略する．ただし，境界条件について，

$v(0)=0$
$\frac{d}{dx}v(x) \big|_{x=0}=d$
$\frac{d^2}{dx^2}v(x) \big|_{x=L}=0$
$\frac{d^3}{dx^3}v(x) \big|_{x=L}=0$

を用いる．棒モデルのときは，2つ目の境界条件の右辺が$0$だったが，根本の角度を調整したいのでパラメータ$d$を導入してコントロール可能にしている．

結局のところ，解は
$$\tag{6} v(x)=\frac{L^{2} \gamma e^{-1 + \frac{x}{L}}}{E I} - \frac{L^{2} \gamma}{e E I} + \frac{x^{2} \left(- L \delta m g + 2 L m g\right)}{8 E I} + \frac{x \left(e E I d - L \gamma\right)}{e E I} + \frac{m g x^{4}}{24 E I L} + \frac{x^{3} \left(L \delta m g - 3 L m g - 3 \gamma\right)}{18 E I L} - \frac{\delta m g x^{6}}{360 E I L^{3}}$$
となる．

2.4. 数式の組み立て

組み立ての方法は先の棒モデルの時と変わらない．

脂肪体関数(再掲)

$$\tag{4} f(y,z)=\sqrt{1-((y_sy)^2+z^2+ay_syz)}$$

変動重力分布たわみ関数(再掲)

$$\tag{6} v(x)=\frac{L^{2} \gamma e^{-1 + \frac{x}{L}}}{E I} - \frac{L^{2} \gamma}{e E I} + \frac{x^{2} \left(- L \delta m g + 2 L m g\right)}{8 E I} + \frac{x \left(e E I d - L \gamma\right)}{e E I} + \frac{m g x^{4}}{24 E I L} + \frac{x^{3} \left(L \delta m g - 3 L m g - 3 \gamma\right)}{18 E I L} - \frac{\delta m g x^{6}}{360 E I L^{3}}$$

乳首関数(再掲)

$$\tag{5} n(y,z,t,s,k)=t\exp{\left(- \left(\frac{(y_sy)^2+z^2+ay_syz}{s^2}\right)^k \right)}$$

これらを用いて，1本の数式を組み立てよう．乳首の座標の決定方法についても棒モデルのときと同じだ．すなわち，
$$x=Lf(y,z-v(x))+n(y-y_{\text{n}},z-v(x)-z_{\text{n}},t_1,s_1,k_1)+n(y-y_{\text{n}},z-v(x)-z_{\text{n}},t_2,s_2,k_2)$$

……完成だ．私の第2️号のおっぱい関数！

2.5. 結果

使用したパラメータは以下の通り；

\left( \begin{array}\\
    L\\
    m\\
    g\\
    E\\
    I\\
    d\\
    \delta\\
    \gamma\\
    a\\
    y_s\\
    t_1\\
    s_1\\
    k_1\\
    t_2\\
    s_2\\
    k_2\\
    y_n\\
    z_n
    \end{array} \right) = \left( \begin{array}\\
    1\\
    0.006\\
    -9.8\\
    0.1\\
    0.1\\
    -0.8\\
    0.28\\
    0.1\\
    0.6\\
    1.1\\
    0.07\\
    0.1\\
    2.5\\
    0.05\\
    0.25\\
    3.3\\
    0\\
    0.24
    \end{array} \right)

でもって，これが結果である．

ぱっと見，非常にいい感じだ．特に……

ここが良いよね．横からの形状，正面からの形状はすごく良い．

.

.

.

.

.

.

.

.

察しの良い方なら気づいているだろう．私は上からの図は今まで見せていなかった．なので，上からの図を提示しよう．

……問題点が伝わらないだろうか？なら，$L$を$1$から$1.35$に変更して，適切な$\gamma$を設定した図と比べてみよう；

きっと伝わったはずだ．伝わらないはずがない．だってこんな記事を読んでいる時点で，あなたもおっぱいに興味があるはずだからだ．

2.6. 問題点

わかりやすく言語化しよう．$L$を拡大するという操作は，脂肪体の仮想半径を大きくするという操作に対応する．すなわち，厳密には$x$軸方向に$L$倍するという操作はおかしい．

$x$軸方向に$L$倍スケールされた半球は，楕円体の$x$軸方向の大きさが変わるだけである．つまり，あくまでも半球は半球のままなのだ(図の赤線)；

しかし，$L$が大きくなったらこんなふうに脂肪体が必ず膨れ上がるのだろうか？否，そんな事はありえない．ボリューム感はないが高さはある人もいるだろうし，ボリューム感はあっても高さがない人もいるだろう．

したがって，ボリューム感を扱うために，新しいパラメータを導入することにしよう．ボリューム感を適切にモデル化できれば，図の紫色の線のような脂肪体の形状にもできるはずだ．脂肪体関数を適切な形に変形して，それに基づいて重力分布関数を改め，再び微分方程式を解き，新しいたわみ関数を導出することにしよう．

3. 膨張脂肪体関数+変動重力分布モデル(Inflatable Fat Function+Variable Gravity Distribution model, IFF+VGD model)

名前がめちゃくちゃ長くなったので，今度から変動重力分布モデルのことをVGDとして，今から定義する膨張脂肪体関数をIFFと呼称することにしよう．

3.1. 膨張脂肪体関数の導出

もともとの脂肪体関数は以下のようなものだった；

脂肪体関数(再掲)

$$\tag{4} f(y,z)=x=\sqrt{1-((y_sy)^2+z^2+ay_syz)}$$
これに関して，$0 \le x$の範囲で膨張可能(Inflatable)な脂肪体関数を導出する．

パラメータ$b$を考えておく．このパラメータは脂肪体関数の膨張具合を示すパラメータだ．これを用いて膨張脂肪体関数を導出していこう．

3.1.1. 想定半径の変化

脂肪体関数の想定半径は$1+\frac{b}{2}$とする．すなわち，$1=1^2$を$\left(1+ \frac{b}{2} \right)^2$で置き換えて，
$$x=\sqrt{\left(1+ \frac{b}{2} \right)^2-((y_sy)^2+z^2+ay_syz)}$$

3.1.2. 極大値・極小値を取る座標の変化

$y,z$が極大値・極小値を取る座標は$x=\frac{b}{2}$であるとする．もともとの球体は$x=0$がこの座標だったため，この操作は換言するならば$x$軸方向に$\frac{b}{2}$だけ平行移動することにほかならない．しかるに，
$$x-\frac{b}{2}=\sqrt{\left(1+ \frac{b}{2} \right)^2-((y_sy)^2+z^2+ay_syz)}$$

3.1.3. 根本での半径の制限

どんな$b$の値だとしても，$x=0$の断面の半径は$1$であるとする．すなわち，両辺2乗して
$$\left( x-\frac{b}{2} \right)^2=\left(1+ \frac{b}{2} \right)^2-((y_sy)^2+z^2+ay_syz)$$
これを$x=0$で評価するから，
$$\begin{align*}\left( 0-\frac{b}{2} \right)^2 &=\left(1+ \frac{b}{2} \right)^2-((y_sy)^2+z^2+ay_syz)\
\Leftrightarrow 1+b&=((y_sy)^2+z^2+ay_syz)\end{align*}$$
したがって，今の$x=0$での断面の半径は$\sqrt{1+b}$である．これを調整するために，元の式の$((y_sy)^2+z^2+ay_syz)$の項に$1+b$をかける操作を行う．
すなわち，変形後の式は
$$\left( x-\frac{b}{2} \right)^2=\left(1+ \frac{b}{2} \right)^2-(1+b)((y_sy)^2+z^2+ay_syz)$$

3.1.4.(≈π) x=1が上限

どんな$b$の値だとしても，$x=1$が脂肪体関数の最大長さである．
現在の想定半径は$1+\frac{b}{2}$であり，$x$は$x=\frac{b}{2}$だけ平行移動していたのだから，現在の$x$の上限は
$$1+\frac{b}{2}+\frac{b}{2}=1+b$$
なら，$x$軸方向に$\frac{1}{1+b}$倍したいなら，$x$を$1+b$倍すればいいので，
$$\left( (1+b)x-\frac{b}{2} \right)^2=\left(1+ \frac{b}{2} \right)^2-(1+b)((y_sy)^2+z^2+ay_syz)$$
整理して，
$$\tag{7} f(x,y,z)=\frac{1}{\sqrt{1+b}}\sqrt{1+bx-((y_sy)^2+z^2+ay_syz)}$$
これが膨張脂肪体関数である．

3.2. 重力分布関数の変化，それによるたわみ関数の変化

重力分布関数関数はもともと以下の形だった；
$$q(x) = \frac{mg}{L} \left( 1-\delta \left(\frac{x}{L}\right)^2 \right) + \frac{\gamma}{L^2} \exp \left( \frac{x}{L}-1 \right)$$
ここの第一項のカッコの中身は，断面の面積の変化によるものだった．では，脂肪体関数が膨張可能なものになった場合，断面積はどんなふうになるだろうか？

再び$y=0$を仮定して解いていこう．すなわち
$$x=\frac{1}{\sqrt{1+b}}\sqrt{1+bx-z^2}$$
を$z$について解く．解いていくと，再び$z$が正負両方出てくるが，興味があるのはそれの2乗なので，正直どっちでもいい．また，今は$x=1$が上限だと考えているが，実際は$x=L$が上限なので，$x$を$\frac{x}{L}$で置き換えれば，
$$z^2=1+b\left(\frac{x}{L}\right)-(1+b)\left(\frac{x}{L}\right)^2$$
これが断面積に比例し，断面あたりの重力分布に比例する．したがって，重力分布関数は
$$q(x) = \frac{mg}{L} \left( 1+\delta \left(b\left(\frac{x}{L}\right)-(1+b)\left(\frac{x}{L}\right)^2\right) \right) + \frac{\gamma}{L^2} \exp \left( \frac{x}{L}-1 \right)$$
となり，これを使って断面曲線方程式を解く．

したがって，IFF-VGDたわみ関数は
$$\tag{8} v(x)=\frac{L^{2} \gamma e^{-1 + \frac{x}{L}}}{E I} - \frac{L^{2} \gamma}{e E I} + \frac{x^{2} \left(L b \delta m g - 3 L \delta m g + 6 L m g\right)}{24 E I} + \frac{x \left(e E I d - L \gamma\right)}{e E I} + \frac{m g x^{4}}{24 E I L} + \frac{x^{3} \left(- L b \delta m g + 2 L \delta m g - 6 L m g - 6 \gamma\right)}{36 E I L} + \frac{b \delta m g x^{5}}{120 E I L^{2}} + \frac{\delta m g x^{6} \left(- b - 1\right)}{360 E I L^{3}}$$
となる．

3.3. 数式の組み立て

ここまでやってると，なんかプラモデルを組み立ててる気分になる．パーツとパーツを組み合わせて新しいパーツを作り，でもってそれを組み合わせてロボットとかをつくる．
しかし，私が組み立てているのはおっぱいの数式だ．文字式だらけで面食らうかもしれないが，私は至って真面目に，少年心をもっておっぱいの数式を組み立てている．

膨張脂肪体関数(再掲)

$$\tag{7} f(x,y,z)=\frac{1}{\sqrt{1+b}}\sqrt{1+bx-((y_sy)^2+z^2+ay_syz)}$$

IFF-VGDたわみ関数(再掲)

$$\tag{8} v(x)=\frac{L^{2} \gamma e^{-1 + \frac{x}{L}}}{E I} - \frac{L^{2} \gamma}{e E I} + \frac{x^{2} \left(L b \delta m g - 3 L \delta m g + 6 L m g\right)}{24 E I} + \frac{x \left(e E I d - L \gamma\right)}{e E I} + \frac{m g x^{4}}{24 E I L} + \frac{x^{3} \left(- L b \delta m g + 2 L \delta m g - 6 L m g - 6 \gamma\right)}{36 E I L} + \frac{b \delta m g x^{5}}{120 E I L^{2}} + \frac{\delta m g x^{6} \left(- b - 1\right)}{360 E I L^{3}}$$

乳首関数(再掲)

$$\tag{5} n(y,z,t,s,k)=t\exp{\left(- \left(\frac{(y_sy)^2+z^2+ay_syz}{s^2}\right)^k \right)}$$

これらを用いて，1本の数式を組み立てよう．ここで，脂肪体関数の引数に$x$が増えたことに注意すると，脂肪体関数は拡大しないといけないから，$x$を$\frac{x}{L}$で置き換えて，
$$x=Lf\left(\frac{x}{L},y,z-v(x)\right)+n(y-y_{\text{n}},z-v(x)-z_{\text{n}},t_1,s_1,k_1)+n(y-y_{\text{n}},z-v(x)-z_{\text{n}},t_2,s_2,k_2)$$

……完成だ．私の第3号のおっぱい関数！

3.4. 結果

使用したパラメータは以下の通り；

\left( \begin{array}\\
    L\\
    m\\
    g\\
    E\\
    I\\
    d\\
    \delta\\
    \gamma\\
    a\\
    y_s\\
    b\\
    t_1\\
    s_1\\
    k_1\\
    t_2\\
    s_2\\
    k_2\\
    y_n\\
    z_n
    \end{array} \right) = \left( \begin{array}\\
    1\\
    0.006\\
    -9.8\\
    0.1\\
    0.1\\
    -0.8\\
    0.28\\
    0.1\\
    0.6\\
    1.1\\
    0\\
    0.07\\
    0.1\\
    2.5\\
    0.05\\
    0.25\\
    3.3\\
    0\\
    0.24
    \end{array} \right)

ここで$b=0$としたのは，一旦もとのVGDモデルとの整合性を比較するためである．$b=0$なら，VGDモデルと一致する．

見たところ，大丈夫そうだ．
ここまで式が複雑になってくると，だんだん自分の計算に自信がなくなってくる．

では，$b \gt 0$にした場合は一体どうなるか？当然，仮定する重力分布の値も大きくなるから，垂れ下がるのは目に見えている．したがって適切に$\gamma$を動かして，脂肪体を支えることにしよう．

これは$b=1$にしたときの横から見た図である．実は，$\gamma$の値は変えていない．なぜかはわからないけど，$\gamma$を変えなくてもいい感じの形状を保っている．

上から見た図．

非常に良い感じの曲線を描いている．

．

．

．

．

．

……うん，非常にいい感じだ．ただ一つの点を除いては．

3.5. 問題点

私は今まで

重力分布$q(x)$の改変によるたわみ関数の変更
脂肪体関数$f$の変更

を行ってきた．逆に言えば，乳首関数については何も手を付けていなかった．

この図を見て欲しい．何かとてつもなく不自然な点があるはずだ．
特に，乳頭部を良く確認してほしい．

．

．

．

．

．

．

そう，なぜか垂れているのだ(図の赤線)．それもそのはず，乳首関数は以下の形で乳首を描画していたことを思い出そう；

$$x=Lf\left(\frac{x}{L},y,z-v(x)\right)+n(y-y_{\text{n}},z-v(x)-z_{\text{n}},t_1,s_1,k_1)+n(y-y_{\text{n}},z-v(x)-z_{\text{n}},t_2,s_2,k_2)$$
この式によれば，乳首の$z$座標は，たわみ関数とともに動くことがわかる．じゃあ，これがいけないのかと思ってこれを
$$x=Lf\left(\frac{x}{L},y,z-v(x)\right)+n(y-y_{\text{n}},z-z_{\text{n}},t_1,s_1,k_1)+n(y-y_{\text{n}},z-z_{\text{n}},t_2,s_2,k_2)$$
とやってみても，乳首の$z$座標が変わるだけで，乳頭の角度は変わらない．

理想的な形は紫色の線のような乳頭部だ．数学的に表現するなら，乳首は脂肪体上の点において，法線ベクトルと並行であるはずだ．しかし，今の乳首関数はただのガウス関数を拡縮したもの．すなわち，現在の乳首関数のピークの法線ベクトルは明らかに$(1,0,0)$であり，つまりは$x$軸と平行である．だから，例えば$z_n$を$z_n=-1.25$などとすると

こうなる．こんなものは乳首ではない．

4. IFF-VGD+最適乳首モデル(IFF-VGD-optN model)

計算において混乱しないために$X_{\text{pre}}$は$X$の変形適用前のもので，$X_{\text{post}}$は$X$の変形適用後のものとする．また，点$P$の$x,y,z$成分をそれぞれ$P.x,P.y,P.z$と表記することにする．

4.1. たわみ前脂肪体関数における乳首座標の計算

乳首の法線ベクトルを脂肪体の法線ベクトルに合わせるためには，そもそも乳首が位置する場所における脂肪体の法線ベクトルが必要であり，それを計算するには乳首が脂肪体上のどの点にあるのかを求めなければならない．

ただ，正直なところ，パラメータに対してあるべき乳首の位置というのは複雑な関係があるから，パラメータを使って表現することはちょっと難しいだろう．そもそも，任意に選んだ点が脂肪体上にあるとは限らないため，いろいろ矛盾が発生してしまう．

したがって，まずはたわみ関数$v(x)$による変形前の脂肪体関数$x=f(x,y,z)$における乳首座標$P_{\text{pre}}$を計算することにしよう．

とりあえず，乳首の位置はコントロールしたい．でも，2次元的にではなく，3次元的に決定したい．そして，その位置というのは確実に脂肪体上にある必要がある．

こういった条件から，私は変形前の脂肪体の外側から原点への直線と，脂肪体の交点を求める方法を思いついた．(私はこれを思いついた時，乳首光線法と名前をつけようと思ったのだが，あまりにもありきたりすぎるのでやめた．)

任意に3次元空間上の点$P_{\text{ray}}=(x_{\text{ray}},y_{\text{ray}},z_{\text{ray}})$を定める．この時，原点$O=(0,0,0)$とその点を結ぶ直線は，以下の平面の交線であるはずだ；

$$y=\frac{y_{\text{ray}}}{x_{\text{ray}}}x$$
$$z=\frac{z_{\text{ray}}}{x_{\text{ray}}}x$$
この2つの式に，変形前の脂肪体関数の式
$$x=f_{\text{pre}}(x,y,z)$$
を加えると，なんと3本の式が手に入る．未知変数は$x,y,z$の3つ，等式の数も3つ．したがって，この問題は解くことができる．

解は
$$x_{\text{pre}}=\frac{x_{\text{ray}} \left(b x_{\text{ray}} + \sqrt{4 a y_{s} y_{\text{ray}} z_{\text{ray}} + b^{2} x_{\text{ray}}^{2} + 4 b x_{\text{ray}}^{2} + 4 x_{\text{ray}}^{2} + 4 y_{\text{ray}}^{2} y_{s}^{2} + 4 z_{\text{ray}}^{2}}\right)}{2 \left(a y_{s} y_{\text{ray}} z_{\text{ray}} + b x_{\text{ray}}^{2} + x_{\text{ray}}^{2} + y_{\text{ray}}^{2} y_{s}^{2} + z_{\text{ray}}^{2}\right)}$$
$$y_{\text{pre}}=\frac{y_{\text{ray}} \left(b x_{\text{ray}} + \sqrt{4 a y_{s} y_{\text{ray}} z_{\text{ray}} + b^{2} x_{\text{ray}}^{2} + 4 b x_{\text{ray}}^{2} + 4 x_{\text{ray}}^{2} + 4 y_{\text{ray}}^{2} y_{s}^{2} + 4 z_{\text{ray}}^{2}}\right)}{2 \left(a y_{s} y_{\text{ray}} z_{\text{ray}} + b x_{\text{ray}}^{2} + x_{\text{ray}}^{2} + y_{\text{ray}}^{2} y_{s}^{2} + z_{\text{ray}}^{2}\right)}$$
$$z_{\text{pre}}=\frac{z_{\text{ray}} \left(b x_{\text{ray}} + \sqrt{4 a y_{s} y_{\text{ray}} z_{\text{ray}} + b^{2} x_{\text{ray}}^{2} + 4 b x_{\text{ray}}^{2} + 4 x_{\text{ray}}^{2} + 4 y_{\text{ray}}^{2} y_{s}^{2} + 4 z_{\text{ray}}^{2}}\right)}{2 \left(a y_{s} y_{\text{ray}} z_{\text{ray}} + b x_{\text{ray}}^{2} + x_{\text{ray}}^{2} + y_{\text{ray}}^{2} y_{s}^{2} + z_{\text{ray}}^{2}\right)}$$
となる．(ちょっと気持ち悪い式だ)

この値を使って，たわみ関数適用前の乳首座標$P_{\text{pre}}=(x_{\text{pre}},y_{\text{pre}},z_{\text{pre}})$を得る．

4.2. たわみ関数の適用と変形後脂肪体における乳首座標の計算

脂肪体関数に対してたわみ関数を適用して，変形後脂肪体関数を
$$x=f_{\text{post}}(x,y,z)=Lf_{\text{pre}}\left(\frac{x}{L},y,z-v(x)\right)$$
としよう．でもって，たわみ関数適用後の乳首座標$P_{\text{post}}$は以下のように求められる；
$$P_{\text{post}}=(LP_{\text{pre}}.x,P_{\text{pre}}.y,P_{\text{pre}}.z+v(LP_{\text{pre}}.x))$$

4.3. 変形後脂肪体関数における乳首座標における法線計算

さて，ここから変形後脂肪体関数の勾配を計算しよう．

変形後脂肪体関数は以下で定義される関数だったことを思い出そう；
$$x=f_{\text{post}}(x,y,z)$$
両辺について確認してみると，両方とも$x$が含まれている．しかも，右辺には関数の引数として$x$が含まれている．したがってこの関数は陰関数だ．

陰関数$F(x,y,z)=0$における勾配は，
$$\nabla F(x,y,z)=\left( \frac{\partial}{\partial x} F(x,y,z),\frac{\partial}{\partial y} F(x,y,z),\frac{\partial}{\partial z} F(x,y,z) \right)$$

であった．勾配ベクトルは陰関数の接平面に垂直なのだから，勾配ベクトルに座標を代入すればそれが法線ベクトルになる．

ではこれを実際にやってみよう．$F$を
$$F(x,y,z)=x-f_{\text{post}}(x,y,z)$$
とする．これで勾配ベクトルを計算する．

……と言いたいのだが，勾配ベクトルはめちゃくちゃ複雑な式になるはずだ．したがって，数式処理プログラムをつかって解いてもらおう(私はPythonのsympyを用いた)．

例として$\frac{\partial}{\partial x} F(x,y,z)$を示すと

$$\frac{\partial}{\partial x} F(x,y,z) = - \frac{L \left(- \frac{a y y_{s} \left(- \frac{L \gamma e^{-1 + \frac{x}{L}}}{E I} - \frac{x \left(L b \delta g m - 3 L \delta g m + 6 L g m\right)}{12 E I} - \frac{e E I d - L \gamma}{e E I} - \frac{g m x^{3}}{6 E I L} - \frac{x^{2} \left(- L b \delta g m + 2 L \delta g m - 6 L g m - 6 \gamma\right)}{12 E I L} - \frac{b \delta g m x^{4}}{24 E I L^{2}} - \frac{\delta g m x^{5} \left(- b - 1\right)}{60 E I L^{3}}\right)}{2} - \frac{\left(- \frac{2 L \gamma e^{-1 + \frac{x}{L}}}{E I} - \frac{x \left(L b \delta g m - 3 L \delta g m + 6 L g m\right)}{6 E I} - \frac{2 \left(e E I d - L \gamma\right)}{e E I} - \frac{g m x^{3}}{3 E I L} - \frac{x^{2} \left(- L b \delta g m + 2 L \delta g m - 6 L g m - 6 \gamma\right)}{6 E I L} - \frac{b \delta g m x^{4}}{12 E I L^{2}} - \frac{\delta g m x^{5} \left(- b - 1\right)}{30 E I L^{3}}\right) \left(z - \frac{L^{2} \gamma e^{-1 + \frac{x}{L}}}{E I} + \frac{L^{2} \gamma}{e E I} - \frac{x^{2} \left(L b \delta g m - 3 L \delta g m + 6 L g m\right)}{24 E I} - \frac{x \left(e E I d - L \gamma\right)}{e E I} - \frac{g m x^{4}}{24 E I L} - \frac{x^{3} \left(- L b \delta g m + 2 L \delta g m - 6 L g m - 6 \gamma\right)}{36 E I L} - \frac{b \delta g m x^{5}}{120 E I L^{2}} - \frac{\delta g m x^{6} \left(- b - 1\right)}{360 E I L^{3}}\right)}{2} + \frac{b}{2 L}\right)}{\sqrt{b + 1} \sqrt{- a y y_{s} \left(z - \frac{L^{2} \gamma e^{-1 + \frac{x}{L}}}{E I} + \frac{L^{2} \gamma}{e E I} - \frac{x^{2} \left(L b \delta g m - 3 L \delta g m + 6 L g m\right)}{24 E I} - \frac{x \left(e E I d - L \gamma\right)}{e E I} - \frac{g m x^{4}}{24 E I L} - \frac{x^{3} \left(- L b \delta g m + 2 L \delta g m - 6 L g m - 6 \gamma\right)}{36 E I L} - \frac{b \delta g m x^{5}}{120 E I L^{2}} - \frac{\delta g m x^{6} \left(- b - 1\right)}{360 E I L^{3}}\right) - y^{2} y_{s}^{2} - \left(z - \frac{L^{2} \gamma e^{-1 + \frac{x}{L}}}{E I} + \frac{L^{2} \gamma}{e E I} - \frac{x^{2} \left(L b \delta g m - 3 L \delta g m + 6 L g m\right)}{24 E I} - \frac{x \left(e E I d - L \gamma\right)}{e E I} - \frac{g m x^{4}}{24 E I L} - \frac{x^{3} \left(- L b \delta g m + 2 L \delta g m - 6 L g m - 6 \gamma\right)}{36 E I L} - \frac{b \delta g m x^{5}}{120 E I L^{2}} - \frac{\delta g m x^{6} \left(- b - 1\right)}{360 E I L^{3}}\right)^{2} + 1 + \frac{b x}{L}}} + 1$$
……気持ち悪い！
いくら解析的に解けているとはいえ，いくらなんでもひどすぎる！

……まあいい．とりあえず，計算はできた．しかし陰関数の微分の結果を羅列することは本質的ではない．したがって，うまいことやって勾配ベクトルを
$$\nabla F(x,y,z)=\left( \frac{\partial}{\partial x} F(x,y,z),\frac{\partial}{\partial y} F(x,y,z),\frac{\partial}{\partial z} F(x,y,z) \right)$$
として計算出来たとしよう．でもって，この勾配ベクトルに$(x,y,z)=P_{\text{post}}$を代入してやれば，たわみ後脂肪体関数上の，乳首座標における法線ベクトル
$$R=\nabla F(P_{\text{post}}.x,P_{\text{post}}.y,P_{\text{post}}.z)$$
が計算できた．

4.4. 回転軸と回転量の計算

法線ベクトル$R$が計算出来たのだから，今の乳首関数のピークの法線ベクトル$(1,0,0)$を$R$の方向に向けられれば，回転が定義できる．

回転については回転行列というものが用いられるが，その一部にロドリゲスの回転公式というものがある．

ロドリゲスの回転公式は，回転軸を表す単位ベクトル$C$を用いて，角度$q$だけ回転させる場合，変換前の座標を$(x,y,z)$，変換後の座標を$(x',y',z')$とすると，以下の関係があることを述べている；

\left( \begin{array}\\
    x'\\
    y'\\
    z'
    \end{array} \right) = \left( \begin{array}\\
    (C.x)^2(1-\cos(q))+\cos(q) & C.xC.y(1-\cos(q))-C.z\sin(q) & C.xC.z(1-\cos(q))+C.y\sin(q)\\
    C.xC.y(1-\cos(q))+C.z\sin(q) & (C.y)^2(1-\cos(q))+\cos(q) & C.yC.z(1-\cos(q))-C.x\sin(q)\\
    C.xC.z(1-\cos(q))-C.y\sin(q) & C.yC.z(1-\cos(q))+C.x\sin(q) & (C.z)^2(1-\cos(q))+\cos(q)
    \end{array} \right)
    \left( \begin{array}\\
    x\\
    y\\
    z
    \end{array} \right)

したがって

回転軸$C$
回転量$q$
が分かれば，回転を定義できるということになる．

$C$について考えよう．$C$はベクトル$(1,0,0)$を$R$に向けるための回転軸だ．すなわち，$C$は$(1,0,0)$と$R$の両方に直交している単位ベクトルである必要がある．幸い，2つのベクトルに直交するベクトルは，そのベクトルの外積として表せる．また，$C$は単位ベクトルなのだから，そのノルムは$1$だ．したがって，$C$は

$$C_{\text{raw}}=R \times (1,0,0)$$
$$C=\frac{C_{\text{raw}}}{\sqrt{C_{\text{raw}} \cdot C_{\text{raw}}}}$$
のように計算できる．ここで，$\times$は外積，$\cdot$は内積である．

次に回転量を求めよう．でも，その前に$C$の中身を確認してみよう．$C_{\text{raw}}$を具体的に計算してみると，

\begin{align*}
C_{\text{raw}} &= R \times (1,0,0)\\
&=(R.y \cdot 0 - R.z \cdot 0, R.z \cdot 1-R.x \cdot 0, R.x \cdot 0 - R.y \cdot 1)\\
&=(0,R.z,-R.y)
\end{align*}

おや？$C_{\text{raw}}.x$は$0$だ．すなわち，単位ベクトルとしてスケーリングされた$C$についても$x$成分は$0$．すなわち$C.x=0$がどんな法線ベクトル$R$を取ってきても満たされている．

単位ベクトルをいくら回転しても，ノルムが違うベクトルには一致しないから，$R$を単位ベクトルにしておこう．つまり，
$$
R_{\text{normalized}} = \frac{R}{\sqrt{R \cdot R}}
$$
としておく．

これを基に，ロドリゲスの回転公式とにらめっこしよう．具体的に$(x',y',z')=R_{\text{normalized}},C.x=0,(x,y,z)=(1,0,0)$をロドリゲスの回転公式に代入すれば，

\begin{align*}
\left( \begin{array}\\
    R_{\text{normalized}}.x\\
    R_{\text{normalized}}.y\\
    R_{\text{normalized}}.z
    \end{array} \right)
    &= \left( \begin{array}\\
    \cos(q) & -C.z\sin(q) & C.y\sin(q)\\
    C.z\sin(q) & (C.y)^2(1-\cos(q))+\cos(q) & C.yC.z(1-\cos(q))\\
    -C.y\sin(q) & C.yC.z(1-\cos(q)) & (C.z)^2(1-\cos(q))+\cos(q)
    \end{array} \right)
    \left( \begin{array}\\
    1\\
    0\\
    0
    \end{array} \right)\\
    &=\left( \begin{array}\\
    \cos(q)\\
    C.z\sin(q)\\
    -C.y\sin(q)
    \end{array} \right)
    
\end{align*}

となる．1行目に着目すれば，回転量$q$は以下のように計算できる；
$$
R_{\text{normalized}}.x=cos(q) \Leftrightarrow q=\cos^{-1}(R_{\text{normalized}}.x)
$$
以上より，回転軸$C$と回転量$q$が求められた．これを具体的に回転公式に代入して，回転後の座標系，すなわち$x'(x,y,z),y'(x,y,z),z'(x,y,z)$が計算できる．

4.5. 乳首関数の回転と座標系の調整

これに関しては非常につまらない話だ．忘れ去られていた乳首関数を再び書けば，
$$\tag{5} n_{\text{pre}}(y,z,t,s,k)=x=t\exp{\left(- \left(\frac{(y_sy)^2+z^2+ay_syz}{s^2}\right)^k \right)}$$
であるから，$x,y,z$を$x'(x,y,z),y'(x,y,z),z'(x,y,z)$に置き換えることで，回転後の乳首関数が手に入る．左辺の$x$も忘れずに置換しよう．すなわち，
$$\tag{9} x'(x,y,z)=t\exp{\left(- \left(\frac{(y_s y'(x,y,z))^2+z'(x,y,z)^2+ay_sy'(x,y,z)z'(x,y,z)}{s^2}\right)^k \right)}$$
とすれば，回転後の乳首関数を定義できる，しかしながら，脂肪体関数の左辺は$x$であり，回転後の乳首関数の左辺は$x'(x,y,z) \neq x$だ．したがって，右辺同士を足し合わせて$x=$の形にしても，うまいこと行かないはずだ．

$x'(x,y,z)$を具体的に計算してみよう．再びロドリゲスの回転公式に立ち返り，$C.x=0$を忘れずに計算すれば，
$$x'(x,y,z)=\cos(q)x-C.z \sin(q)y+C.y\sin(q)z$$
である．すなわち，これを$x$について解けば，回転後の乳首関数を脂肪体関数の世界に持ってこれる．

両辺から$(-C.z \sin(q)y+C.y\sin(q)z)$を引き，その後$\cos(q)$で割ればいい．ただし，$t$については乳首関数のスケーラーであるから，カッコの外だ．また，実際は乳首座標は$P_{\text{post}}$なのだから，回転後の乳首関数は
$$\tag{10} n_{\text{post}}(x,y,z,t,s,k)=x=\frac{t}{\cos(q)} \left(n_{\text{pre}}(y'(x-P_{\text{post}}.x,y-P_{\text{post}}.y,z-P_{\text{post}}.z),z'(x-P_{\text{post}}.x,y-P_{\text{post}}.y,z-P_{\text{post}}.z),1,s,k)- (-C.z \sin(q)(y-P_{\text{post}}.y)+C.y\sin(q)(z-P_{\text{post}}.z))\right)$$
となる．

4.6. 乳首領域関数

さて，パーツは揃った．変形後の脂肪体関数，回転後の乳首関数．ではこれを足し合わせれば完成するか？というと，実はそうもいかなかったのだ(誤算)．

回転後の乳首関数をプロットしてみよう．

このプロットは，$z_{\text{ray}} \gt 0$とした時の回転後の乳首関数$n_{\text{post}}(x,y,z,1,0.1,2)$をプロットしたものだ．想定通り，乳首関数は上方向に回転している．

しかし問題は乳首領域ではなく，非乳首領域についてである．最終的な数式は，脂肪体関数と乳首関数を$x$軸方向に足し合わせることで完成する．しかし，よく見てみて欲しい．

非乳首領域が$x \approx 0$を満たしていないのだ(図の赤線)．つまり，このまま足し合わせると，脂肪体の形状も歪んでしまう．
私は脂肪体の法線ベクトルが乳首の法線ベクトルであることにしたかったのに，これでは前提が崩れてしまう．これでは意味がない．

理想的な非乳首領域は紫色の線で示したとおりだ．つまり，非乳首領域を$0$に持っていくための何らかの追加的な項が必要だ．

少し考えると，乳首領域ならば$1$,非乳首領域なら$0$になるような関数を定めてやれば良い．なにかないか何かないか……と悩んで試行錯誤していると，私はすでにそれを導いている事に気づいた．

$(9)$式を再び見て欲しい；
$$\tag{9} x'(x,y,z)=t\exp{\left(- \left(\frac{(y_s y'(x,y,z))^2+z'(x,y,z)^2+ay_sy'(x,y,z)z'(x,y,z)}{s^2}\right)^k \right)}$$
この式は，非乳首領域が$x'(x,y,z) \approx 0$で表されている．逆に言えば，乳首領域では$x'(x,y,z)=t$になる．裏を返せば，$t=1$ならば，右辺は乳首領域なら$1$に近い値を取り，非乳首領域なら$0$に近い値を取る．つまり，この右辺を$t=1$と置いたものは目的を満たせる可能性がある．

具体的にプロットして確認してみよう．

うん．たしかに非乳首領域では$x \approx 0$のようだ．でもこれは乳首関数の代わりにはならない．なぜならば

$z_{\text{ray}}$を変化させると，乳首の長さが変わってしまうからだ．すなわち，真の乳首項は
$$(乳首項)=(回転後の乳首関数)\cdot (乳首領域関数)$$
という積の形で表される．実際の乳首座標との整合性を考えれば，乳首領域関数を以下のようになる；

$$\tag{11} V(x,y,z,s,k)=n_{\text{pre}}(y'(x-P_{\text{post}}.x,y-P_{\text{post}}.y,z-P_{\text{post}}.z),z'(x-P_{\text{post}}.x,y-P_{\text{post}}.y,z-P_{\text{post}}.z),1,s,k)$$

4.7. 数式の組み立て

さて，最後のステップだ．とはいえ，今までと大きく変わらない．違うのは，乳首項が回転後乳首関数と乳首領域関数の積になっただけだ．すなわち，
$$x=f_{\text{post}}(x,y,z)+ n_{\text{post}}(x,y,z,t_1,s_1,k_1) V(x,y,z,s_1,k_1)+ n_{\text{post}}(x,y,z,t_2,s_2,k_2) V(x,y,z,s_2,k_2)$$
……完成だ．私の第4号のおっぱい関数！

4.8. 結果

使用したパラメータは以下の通り；

\left( \begin{array}\\
    L\\
    m\\
    g\\
    E\\
    I\\
    d\\
    \delta\\
    \gamma\\
    a\\
    y_s\\
    b\\
    t_1\\
    s_1\\
    k_1\\
    t_2\\
    s_2\\
    k_2\\
    x_{\text{ray}}\\
    y_{\text{ray}}\\
    z_{\text{ray}}
    \end{array} \right) = \left( \begin{array}\\
    1\\
    0.006\\
    -9.8\\
    0.1\\
    0.1\\
    -0.8\\
    0.28\\
    0.1\\
    0.6\\
    1.1\\
    1\\
    0.07\\
    0.1\\
    2.5\\
    0.05\\
    0.25\\
    3.3\\
    2\\
    0\\
    -2.5
    \end{array} \right)

脂肪体については特段の変更はないから，乳首部分だけを比較しよう．

IFF+VGD model

IFF+VGD+optN model

いい感じに乳頭の法線ベクトルが脂肪体の法線ベクトルと合っている．

でも個人的にはもっと乳首は上についていて欲しい．でもって，もっとでかい方が良いし，なんならボリューム感も合ったほうがいい．要するに，この形は私好みじゃない．

したがって私は自分にとって最適なパラメータを探すことにした．いろいろいじっていくと，以下のパラメータが自分にとっていいな～と思うおっぱいになった；

\left( \begin{array}\\
    L\\
    m\\
    g\\
    E\\
    I\\
    d\\
    \delta\\
    \gamma\\
    a\\
    y_s\\
    b\\
    t_1\\
    s_1\\
    k_1\\
    t_2\\
    s_2\\
    k_2\\
    x_{\text{ray}}\\
    y_{\text{ray}}\\
    z_{\text{ray}}
    \end{array} \right) = \left( \begin{array}\\
    1.45\\
    0.006\\
    -9.8\\
    0.1\\
    0.1\\
    -0.8\\
    0.28\\
    0.132\\
    0.6\\
    1.1\\
    1\\
    0.07\\
    0.1\\
    2.5\\
    0.05\\
    0.3\\
    4.22\\
    2\\
    0.55\\
    0.3
    \end{array} \right)

すなわち

$L$が大きい: 大きいおっぱい
$\gamma$が大きい: 反重力や張力が大きい
$y_{\text{ray}},z_{\text{ray}} \gt 0$: 乳首は外側にあって，上側にある
$s_2$が大きい: 乳輪が大きい
$k_2$が大きい: 乳輪がくっきりしている
という結果になった(何の分析なんだ？)．

モデルの応用

1. 乳揺れ

前例において，時点$t$を動かすとおっぱいが揺れるというものを再現したものがあった．これを見て私は感動したわけだが，IFF+VGD+optNモデルの枠組みなら，これを扱える．

1.1. 方法

おっぱいが揺れるというのは，すなわちたわみ関数$v(x)$が変化することに他ならない．断面曲線方程式のパラメータは$L,m,g,E,I,d,\delta,\gamma$であったが，この内$L,m,E,I,\delta$は脂肪体の構造的を表すパラメータであるから，あまり動かしたくない．つまり，$g,d,\gamma$を動かすことになる．

各パラメータの意味を再び確認する

$g$は重力加速度で，脂肪体の全体に均一に掛かる加速度だ．
$d$は根本の傾きで，脂肪体の根本の角度だ．
$\gamma$は反重力パラメータで，脂肪体の先端に強くかかり，根本に行くほど弱くなる力のパラメータだ．

つまり，どんな揺らし方がいいかによって，動かすパラメータも変わる．

私は$\gamma$を動かすのがいい気がしたので，$\gamma=\bar{\gamma}+0.01\sin(\pi t)$として，$t \in [-1,1]$で動かしてみた
(マウスカーソルが映っちゃったけど，これは単純に録画ミスです．ごめんね)

揺れた！嬉しい！！

1.2. 問題点

たしかに，たわみ関数のパラメータを動かすことで乳揺れは再現できる．しかし，そもそもたわみ関数は静的なたわみ量を表現しているものであるはずなので，これのパラメータを三角関数で周期的に変化させたとしても，本当の揺れを再現していることにはならないはずだ．

おそらく，時間を通じた揺れを再現したいなら，

$t=0$に力かなにかがかかる
$t$が増加するにつれて揺れが収まっていく
という物理法則を記述しなければいけない．

おそらく，バネか何かを仮定して，その揺れが逓減していくようなモデルを考えればいいだろう．

でも，今のところは面倒なのでパスする．

2. 小さいおっぱい

これも前例についてなのだが，私が簡潔で(数式的に)美しいとしたおっぱい関数
$$x=1$$
で形容される小さいおっぱいについて考えよう．

2.1. 方法

このおっぱい関数$x=1$について，乳首項がない．また，厚さも0だ．つまり，IFF+VGD+optNモデルで記述するには，
$$L=0,t_1=0,t_2=0$$
とすればいい．

でも，よく考えて欲しい．小さなおっぱいは本当に絶壁なのか？　ということだ．
小さいおっぱいはたしかに小さい．しかしながら，確かに脂肪体は存在していて，小さいながらもわずかに膨らんでいるはずだ．

正直私は巨乳派なので，貧乳派の理想のおっぱいはよくわからない．しかし，IFF-VGD+optNモデルは小さいおっぱいも表現できそうなので，再びパラメータをいじって再現してみた．すると
$$L=0.3,b=-0.3$$
とすることで，確かな膨らみが存在するが，小さくてボリューム感がないおっぱいが再現できた．

貧乳派にとっての理想的な貧乳についてよく知らないが，おそらくこんな感じだろうか？
(誓って言うが，私は巨乳派である一方で，小さい胸を否定していたり軽視しているわけではなく，ただ単によく知らないだけなのだ．したがって，小さい胸の魅力や形状の特徴について教えてくれる人が居たら，喜んで話を聞かせてもらうつもりだ．)

2.2. 問題点

これは膨張可能な脂肪体関数の定式化についての問題だ．仮に脂肪体が球状ではなく，円錐状やそれに似た尖った形状だとする．これを再現するために$b$を$0$より小さくしていこう．すると，$b=-1$になると脂肪体が消滅する．

これはそもそもの導出において，想定半径が$1+\frac{b}{2}$で，$y,z$が極大値・極小値を取る，すなわち仮想球の中心は$x=\frac{b}{2}$という仮定を置いたからだ．もし$b=-1$になったら，想定半径が$\frac{1}{2}$なのにもかかわらず，$x=-\frac{1}{2}$が中心．つまり$x=0$では何も存在しない．つまり，IFFモデルの枠組みでは，完全に丸みのないおっぱいは再現不可能である．脂肪体関数の関数形を考えることも必要だと思うが，これも面倒なので保留する．

まとめ

本稿では，既存のおっぱい関数の発展型として，物理法則に基づき，魅力的で，3Dなおっぱい関数を導出した．導出はすごく面白かったし，想定通りのおっぱいになったときはとても嬉しいものだった．

一方で，既存のモデルと同じように，このモデルもまた変な仮定やロジックが含まれているはずだ．なぜなら私は物理学しょしんしゃだからである．したがって，このモデルについて，「ここが変！」とか，「ここを修正したい！」と思う場合，修正は大いに歓迎する．

おまけ

ローカルでもおっぱい関数が楽しめるように，Jupyter Notebookならびにその周辺ファイルを作成した．

このレポジトリでは，モデルの解析解が.tex形式で省略無しで書かれている(空白なしで1,781,222文字になる超巨大な数式だ)．また，ローカルで導出・プロットもできるJupyter Notebookもつけさせていただいた．

でも，それだけじゃない．おっぱいを3Dモデルとして保存できるのだ．つまり，時間をかけて計算したおっぱいをずっと眺められる．

またなにか思いついたら，レポジトリを更新していこうと思う．