はじめに

　この記事は前回の、非円形摩擦車の記事：
　p5.jsで非円形摩擦車を作成する
の続きです。非円形摩擦車を外側に作る場合は前回ニュートン法とか交えて詳しく説明したんですが、外側に作る場合だったので、内側にも作れるのを説明したいと思いました。
　まず、この記事では二つのことがポイントになります。今回内側の図形に使うのは周期がTAUより大きい曲線です。具体的には5TAUです（$10\pi$）。そうでないと内側で摩擦する図形を作れないからです。また、それを外側にも作ります。もうひとつは、内側のため、ニュートン法の式がちょっと変わることですね。

コード全文

　INNER star wheel remake

/*
  数値の算出元：https://editor.p5js.org/dark_fox/sketches/X4N0iyFCO
  ニュートン法を内側から外側に向かって使う
  ニュートンで書き直しました。お疲れさま。

  まずTAUを5*TAUにする
  これにより0～0.2で一周する
  ちなみに0～0.4で一周させるには0～1で(5/2)*TAUにすればいい
  なぜなら0～0.4でTAUだから。必然的に0～1の要請は(5/2)*TAUになる。
  このコードでは5*TAUにしているので0～0.2で一周してTAUになる。それで卵型。
  こういう遊びができるのは被積分関数が周期0.2だから。でないとおかしなことになる。

  つまり被積分関数の周期が1より小さく、整数分の1であることで、こういうことができる
  わけですね。
*/

/*
  非円形摩擦車
  参考：
  https://www.youtube.com/watch?v=uqlpZN6eX2o
  https://www.youtube.com/watch?v=xbONXxrTIcc
*/

const master = {};
const WHEEL_NUM = 5;
const CANVAS_SCALE = Math.floor(Math.min(window.innerWidth, window.innerHeight)*0.99);
const CANVAS_WIDTH = CANVAS_SCALE;
const CANVAS_HEIGHT = CANVAS_SCALE;
const SCALE_FACTOR = CANVAS_SCALE*0.1;

// インテグラル0からtまでのfという関数の1000刻みでのt:0～1の値を格納するための関数
function getIntegralArray(f,dt=1000){
  const result=[0];
  let sum=0;
  for(let i=0;i<dt;i++){
    sum+=(1/dt)*f(i/dt);
    result.push(sum);
  }
  // xに対してresult[x/dt]が疑似的にfの積分になる
  // 間の場合は補間する
  // たとえばfloor(x/dt)やfract(x/dt)を取るなどする
  // 0～dtに入るように適切に制限する
  return result;
}

// newton法関連

// radius関数の0～1におけるdetail刻みでの最大値を出すだけ
// なぜ必要なのかは後述
function maxCurveRadius(radius,dt=1000){
  let r=-Infinity;
  for(let i=0;i<dt;i++){
    r=Math.max(radius(i/dt),r);
  }
  return r;
}

// innerの場合はこっちが要る。つまり最小半径。
function minCurveRadius(radius,dt=1000){
  let r=Infinity;
  for(let i=0;i<dt;i++){
    r=Math.min(radius(i/dt),r);
  }
  return r;
}

function getIntegralOuter(radius,diffTheta,a,dt=1000){
  let sum=0;
  for(let i=0; i<dt; i++){
    sum+=(1/dt)*diffTheta(i/dt)*radius(i/dt)/(a-radius(i/dt));
  }
  return sum;
}

// 内側版。符号が逆。
function getIntegralInner(radius,diffTheta,a,dt=1000){
  let sum=0;
  for(let i=0; i<dt; i++){
    sum+=(1/dt)*diffTheta(i/dt)*radius(i/dt)/(radius(i/dt)-a);
  }
  return sum;
}

// Powerは共通で使う。
function getIntegralPower(radius,diffTheta,a,dt=1000){
  let sum=0;
  for(let i=0; i<dt; i++){
    sum+=(1/dt)*diffTheta(i/dt)*radius(i/dt)/Math.pow(a-radius(i/dt),2);
  }
  return sum;
}

// 精度やばいです
function applyNewtonOuter(radius, diffTheta, target = Math.PI*2, dt=1000){
  const r = maxCurveRadius(radius,dt);
  // r+0.1とかr+0.01
  let initialValueDiff = 0.1;
  while(1){
    if(getIntegralOuter(radius, diffTheta, r+initialValueDiff) > target) break;
    initialValueDiff *= 0.1;
  }
  let a0 = r + initialValueDiff;
  for(let k=0; k<20; k++){
    a0 += (getIntegralOuter(radius,diffTheta,a0,dt)-target)/getIntegralPower(radius,diffTheta,a0,dt);
    console.log(a0);
  }
  return a0;
}

// innerの場合、符号が逆になる。色々逆になる。
function applyNewtonInner(radius, diffTheta, target = Math.PI*2, dt=1000){
  const r = minCurveRadius(radius,dt);
  // r+0.1とかr+0.01
  let initialValueDiff = 0.1;
  while(1){
    if(getIntegralInner(radius, diffTheta, r-initialValueDiff) > target) break;
    initialValueDiff *= 0.1;
  }
  let a0 = r - initialValueDiff;
  for(let k=0; k<20; k++){
    a0 -= (getIntegralInner(radius,diffTheta,a0,dt)-target)/getIntegralPower(radius,diffTheta,a0,dt);
    console.log(a0);
  }
  return a0;
}

function setup() {
  createCanvas(CANVAS_WIDTH, CANVAS_HEIGHT);

  const r1 = (t) => 2+0.4*cos(5*TAU*t);
  // ここをいじることで不等速回転にできる？
  const theta1 = (t) => TAU*t;
  const diffTheta1 = (t) => TAU;
  master.r1 = r1;
  master.theta1 = theta1;

  const detail = 800;

  const a2 = applyNewtonOuter(r1, diffTheta1, 5*TAU, detail);
  const f2 = (t) => diffTheta1(t)*r1(t)/(a2-r1(t));
  const intF2 = getIntegralArray(f2,detail);
  //console.log(TAU-intF[800]); // ほぼTAU

  const r2 = (t) => a2-r1(t);
  // 要するにtheta2はintFなんですが値が離散的なので（配列で与えられるので）
  // 前後の値を算出したうえで補間で出しています
  const theta2 = (t) => {
    const x = fract(t);
    const nx = floor(x*detail);
    const fx = fract(x*detail);
    const leftValue = intF2[nx];
    const rightValue = intF2[nx+1];
    return leftValue*(1-fx) + rightValue*fx;
  }
  master.a2 = a2;
  master.r2 = r2;
  master.theta2 = theta2;

  const a3 = applyNewtonInner(r1, diffTheta1, 5*TAU);
  const f3 = (t) => diffTheta1(t)*r1(t)/(r1(t)-a3);
  const intF3 = getIntegralArray(f3,detail);

  const r3 = (t) => r1(t)-a3;
  const theta3 = (t) => {
    const x = fract(t);
    const nx = floor(x*detail);
    const fx = fract(x*detail);
    const leftValue = intF3[nx];
    const rightValue = intF3[nx+1];
    return leftValue*(1-fx) + rightValue*fx;
  }
  master.a3 = a3;
  master.r3 = r3;
  master.theta3 = theta3;

  const pathDetail = 400;

  const path1 = new Path2D();
  for(let h=0;h<=pathDetail;h++){
    const t = h/pathDetail;
    const x = SCALE_FACTOR*r1(t)*cos(theta1(t));
    const y = SCALE_FACTOR*r1(t)*sin(theta1(t));
    if(h===0){
      path1.moveTo(x,y);
    }else{
      path1.lineTo(x,y);
    }
  }
  master.path1 = path1;

  const path2 = new Path2D();
  for(let h=0;h<=pathDetail;h++){
    const t = h/pathDetail;
    const x = SCALE_FACTOR*r2(t)*cos(theta2(t));
    const y = SCALE_FACTOR*r2(t)*sin(theta2(t));
    if(h===0){
      path2.moveTo(x,y);
    }else{
      path2.lineTo(x,y);
    }
  }
  master.path2 = path2;

  const path3 = new Path2D();
  for(let h=0;h<=pathDetail;h++){
    const t = h/pathDetail;
    const x = SCALE_FACTOR*r3(t)*cos(theta3(t));
    const y = SCALE_FACTOR*r3(t)*sin(theta3(t));
    if(h===0){
      path3.moveTo(x,y);
    }else{
      path3.lineTo(x,y);
    }
  }
  master.path3 = path3;
}

function draw(){
  background(0);
  const {
    a2,r1,r2,theta1,theta2,path1,path2,
    a3,theta3,path3
  } = master;
  
  noStroke();
  const elapsedTime = millis()/3000;

  fill(255,128,64);
  push();
  translate(CANVAS_WIDTH*0.5, CANVAS_HEIGHT*0.5);
  // 全体を反時計回りに回しているのでtimeのあとの右側の位置は-timeの時の位置が来る。
  // マイナスになるのはそういうこと。
  // 座標軸をマイナスに回してそれに基づいてpathで描画しているところをイメージすればわかる。
  rotate(-theta1(elapsedTime));
  drawingContext.fill(path1);
  pop();

  fill(255,128,64,128);
  for(let k=0;k<WHEEL_NUM;k++){
    const phi = k*TAU/WHEEL_NUM;
    push();
    // 中心位置に移る計算
    translate(CANVAS_WIDTH*0.5 + a2 * SCALE_FACTOR*cos(phi), CANVAS_HEIGHT*0.5 + a2*SCALE_FACTOR*sin(phi));
    // 鏡映変換で、軸に垂直な直線による折り返しを実行する。
    applyMatrix(-cos(2*phi),-sin(2*phi),-sin(2*phi),cos(2*phi),0,0);
    rotate(-theta2(elapsedTime+k/WHEEL_NUM) - theta1(elapsedTime) + theta1(elapsedTime+k/WHEEL_NUM));
    drawingContext.fill(path2);
    pop();
  }

  fill(255);
  
  for(let k=0;k<WHEEL_NUM;k++){
    const phi = k*TAU/WHEEL_NUM;
    push();
    // 中心位置に移る計算
    translate(CANVAS_WIDTH*0.5 + a3 * SCALE_FACTOR*cos(phi), CANVAS_HEIGHT*0.5 + a3*SCALE_FACTOR*sin(phi));
    // 内側の場合はそのまま
    rotate(-theta3(elapsedTime+k/WHEEL_NUM) - theta1(elapsedTime) + theta1(elapsedTime+k/WHEEL_NUM));
    drawingContext.fill(path3);
    pop();
  }
}

コピペが...多いです。リファクタリングの練習にうってつけですね（こら！）。
実行結果：

基本となる摩擦車の関数について

　今回用いるのは次の関数です。

  const r1 = (t) => 2+0.4*cos(5*TAU*t);
  // ここをいじることで不等速回転にできる？
  const theta1 = (t) => TAU*t;
  const diffTheta1 = (t) => TAU;

極座標で星形を作っています。特に工夫はありません。半径の関数の周期は今回0.2となっています。前回は1でした。関数のグラフとしての周期は1ですが、半径に関しては0.2です。ここがポイントです。

向きが逆になる

　要請の式がちょっと異なります。前回は

r_1(t)+r_2(t)=a,~~~~~r_1(t)\dot{\theta}_1(t) = r_2(t)\dot{\theta}_2(t)

でしたが、内側の場合こうなります：

r_1(t)-r_3(t)=a,~~~~~r_1(t)\dot{\theta}_1(t) = r_3(t)\dot{\theta}_3(t).

ゆえに、計算において符号が若干異なるんですが、似たような計算で似たようなことが成り立ちます。まあ要するにこれを満たすr_3を作ればいいです。符号が逆なので、微分は一致します。

\dot{r}_1(t) = \dot{r}_3(t).

積分の式

　積分の式も逆になります。前回同様、

\theta_2(t) = \int_0^t \frac{\dot{\theta}_1(\tau)r_1(\tau)}{a-r_1(\tau)}d\tau

です。なお被積分関数は違う文字の方がいいと思ったので$\tau$にしてあります。これに対して、目的の内側の関数については、

\theta_3(t) = \int_0^t \frac{\dot{\theta}_1(\tau)r_1(\tau)}{r_1(\tau)-a}d\tau

となります。ここで$\theta_1$は微分、つまり定数の形なので、被積分関数は周期0.2となっています。ゆえに、前回は

\theta_2(1)=2\pi

を要請しましたが、今回は0.2の倍数であればOKです。そこで今回は大胆に、

\theta_2(0.2)=2\pi,~~~~\theta_3(0.2)=2\pi

を要請としています（0.4や0.6でもOKなので試してみてください）。つまり、

\theta_2(1)=10\pi,~~~~\theta_3(1)=10\pi

というわけですね。

ニュートン法の適用について

　$\theta_2$の方に関しては、ニュートン法の適用の仕方は完全に一緒です。なので、$\theta_3$の方だけ説明します。

　今回は最小半径を取っています。あそこが最小半径よりちょっと小さいと積分値がちゃんと有限になること、最小半径に小さい方から近づいていくと正の無限大に発散することがわかるでしょうか。この辺は数学的感覚が要るんですが、まあ何となくわかると思います。

　で、0のとき負の有限値になります。微分するとその結果は前回と同じ感じで、今回は正になるので単調増加です。そういうわけでこんなグラフになり、a軸とただ1点で交わるので、これをニュートン法で求めることになります。その辺りの計算を、符号が違うだけで前回とほぼおんなじロジックで計算しています。説明は以上です。

　また、$10\pi$なのでこうなりますが、もし前回と同じ$2\pi$の場合、0で原点を通ってしまい、内側で小さくなって接することができません。ゆえにこの手法では周期的であることがマストです。それでこのようなコードになっています。

描画部分について

　描画部分もほぼ同じです。唯一違うのは、鏡映変換をしないことだけです。

  fill(255);
  
  for(let k=0;k<WHEEL_NUM;k++){
    const phi = k*TAU/WHEEL_NUM;
    push();
    // 中心位置に移る計算
    translate(CANVAS_WIDTH*0.5 + a3 * SCALE_FACTOR*cos(phi), CANVAS_HEIGHT*0.5 + a3*SCALE_FACTOR*sin(phi));
    // 内側の場合はそのまま
    rotate(-theta3(elapsedTime+k/WHEEL_NUM) - theta1(elapsedTime) + theta1(elapsedTime+k/WHEEL_NUM));
    drawingContext.fill(path3);
    pop();
  }

なぜなら回る方向が一緒だからです。それで内側になります。これで説明は終わりです。

おわりに

　これを使って非円形歯車が作れるらしいです...作れる人すごいです。
　ここまでお読みいただいてありがとうございました。