うごくものをうごかしたい

なんとなく書きたくなったので書きます。10割が算数の話だしまぁ「認識されないディティール」の話なので意味はないです。自分用のメモみたいなもん。

ひも

f:id:phi16_ind:20210612162407p:plain

CRTの後ろのこの子は、Capが元々垂れてたケーブルにしていたものの、CRTを動かすと動いてしまうのでやばいよね、ということで動的にいい感じにする必要がありました。

f:id:phi16_ind:20210612162312p:plain

そこで元のメッシュはsubdiv付き線分にしてもらって、いい感じにするということでいい感じにしました。

まず紐の垂れる形状は懸垂線というやつで、 $y=\cosh(x)$ の一部になることが知られています。

ではどの一部なのかが問題です。

今回の制約条件は「端点の位置」「紐の長さ」です。紐が平面にあると仮定すると、紐の形状を決定するのは「(一様) 拡大」「平行移動」の3つのパラメータ。

端点の位置を $(0,0)$ と $(1,v)$ として正規化して、紐の形状を $f(x)=\beta+\gamma\cosh((x-\alpha)/\gamma)$ で表すことにします。正規化された状態での (固定の) 紐の長さを $L$ とします。

与えられた制約から:

始点 $f(0)=0$
終点 $f(1)=v$
弧長 $\gamma\sinh((1-\alpha)/\gamma) - \gamma\sinh((0-\alpha)/\gamma) = L$

なので、解きます。解けません。

諦めて二分法にします。 $\gamma$ を固定した際の弧長を計算してぱたぱた動かしていきます。

$\gamma$ から $\alpha$ が計算できそうなことがわかったのでそこから倒します。がちゃがちゃ叩いたら次のようになりました。

$A=\alpha/\gamma, B = 1/\gamma, C = v/\gamma$ とすると
$\cosh(B-A) - \cosh(0-A) = C$ より
$e^A = (e^B\sqrt{e^{-3B}(e^B(C^2-2) + e^{2B} + 1)}-C)\times e^B/(e^B-1)$

よって $\gamma$ を固定すると $\alpha$ が出るのでそこから弧長 $\gamma(\sinh(B-A) - \sinh(0-A) )$ が計算できます。

後はぺちぺちやるだけです。ただこの時点で $\gamma$ よりも $1/\gamma$ の方が単位がまともそうなことに気づいたので、 $1/\gamma$ に対して二分法をやります。

float ic = 3, dc = ic/2;
float c = 1/ic, ac = 0;
for(int i=0;i<12;i++) {
    float E = exp(1/c);
    float C = v/c;
    float eA = (E*sqrt(exp(-3/c)*(E*(C*C-2) + E*E + 1)) - C) * E/(E-1);
    ac = log(eA);
    float l = c*(sinh(1/c-ac) - sinh(0/c-ac));
    if(l < L) ic += dc;
    else ic -= dc;
    c = 1/ic;
    dc /= 2.0;
}

初期値は実験的に決めて $1/\gamma = 3$ になりました。これで c つまり $\gamma$ がわかったので、後は素直に答えが出ます。

さて曲線の形状がわかったので実際に曲げます。単なる平行移動では角度が変わらないので太さ一定になりません。よってちゃんと接空間を取って回転させる必要があります。

今回の曲線は (正規化された状態で) $(t, \beta+\gamma\cosh((t-\alpha)/\gamma) )$ と表されます。

よって接ベクトルは $(1, \sinh((t-\alpha)/\gamma) )$ (の方向) です。あとbinormalはY軸と直行するように設定しました。後はやるだけです。

全体としてこんな感じになってます。

float dist = distance(crtPoint.xz, fixPoint.xz);
float scale = (crtPoint.y - fixPoint.y) / dist;
float3 su = computeCurve(scale, 1.2/dist); // これがさっきのめんどくさいやつ

float3 target = lerp(fixPoint, crtPoint, float3(t, (su.y+su.z*cosh((t-su.x)/su.z))/scale, t));
// calculate tangent
float3 d = (crtPoint - fixPoint) * float3(1, sinh((t-su.x)/su.z)/scale, 1);
d = normalize(d);
float3 T = normalize(d);
float3 up = float3(0,1,0);
float3 B = normalize(cross(up,T));
float3 N = cross(T,B);
perp = mul(transpose(float3x3(B,N,T)), perp);

v.vertex.xyz = target + perp;

ちなみに元々ついてる斜めのケーブルを真っ直ぐに (tangent = +Z) にするために最初適当な逆回転を入れています。そしてケーブルの中心軸が target、そこからのズレを perp として、回転するときは perp だけを回す感じですね。

そういえば弧長パラメータは取ってないです。めんどくさくて。多分わからんし。

あとCRTが動いたときの揺れはUdonから渡ってきた適当な値を元にして target 側に足し合わせてるだけです (重みは $1-(2t-1)^2$ らしい)。

正しくはないけどさすがにわからないと思うのでOKだと思います。

前後移動

あとは大した話ではないのだけど。

f:id:phi16_ind:20210612171613p:plain

この曲線は「 $\tanh$ の一部分を切り出して貼っつけて繰り返して $\cos$ に突っ込んだもの」です。ええと。

結論から言うとこれなんだけど。

f:id:phi16_ind:20210612203439p:plain

まず $\tanh$ を適当に伸ばしたりずらしたものがあって。

f:id:phi16_ind:20210612203758p:plain

ひっくり返して貼っ付けます。で、そうすると端点の傾きも一致するのでループできます。

f:id:phi16_ind:20210612203923p:plain

出力のY座標が「整数付近が広い」のと、「半整数付近も広い」のが重要。

あとは $x\mapsto\cos(2\pi x)$ に突っ込むとこう。

f:id:phi16_ind:20210612204117p:plain

つまり整数らへんは1に、半整数らへんは-1になるわけです。

この曲線のパラメータは「傾きのきつさ (f1, 0より大きい値) 」「デューティ比っぽいの (f2, 0から1の間)」で制御できてて、こう、いい感じというわけです。

そういえばコントローラの画面にグラフを映すためにコレをやっています。これをやるために微分もついでに計算する必要がありました。

float2 crtF(float2 p, float x) {
    float s = p.x;
    float b = p.y;
    float v0 = tanh((0-b)*s), v1 = tanh((1-b)*s);
    float u = frac(x/2)*2-1;
    float ug = 1;
    float v = sign(u) * (tanh((abs(u)-b)*s)-v0)/(v1-v0);
    float vg = ug * s * pow(1 / cosh((abs(u)-b)*s), 2) / (v1-v0);
    float w = floor(x) + v*0.5+0.5;
    float wg = vg*0.5;
    float2 o = float2(cos(UNITY_PI*2*w), -UNITY_PI*2*wg*sin(UNITY_PI*2*w));
    o.x += tanh(o.y) * p.x / 15.0f * 0.1; // grad bump
    return o;
}