斜方投射(空気抵抗あり)の飛距離を求めたい

斜方投射はきっと多くの方が通る道でしょう。スポーツをしている人々も大学生も斜方投射したことがない人の方が少ないのではないでしょうか。
ですが、空気抵抗ありの斜方投射で飛距離を求めるのはかなりの難易度です。その方法をここに記します。

まずは式から

このページより、$x$を飛距離、$y$を初期条件の高さ、$k$を空気抵抗係数、$m$を質量、$g$を重力加速度として

X=\frac{kx}{mv_0\cos θ}

Y=-\frac{k^2y}{m^2g}

Z=1+\frac{kv_0\sin θ}{mg}

Y=XZ+\log(1−X)

ということがわかります。指数をとると

e^Y=(1−X)e^{XZ}

正負を反転させて

-e^Y=(X-1)e^{XZ}

$e^{-Z}$倍して

-e^{Y-Z}=(X-1)e^{Z(X-1)}

$Z$倍して

-Ze^{Y-Z}=Z(X-1)e^{Z(X-1)}

ここでランバート(ランベルト)の$W$関数という$xe^x$の逆関数(初等関数ではない)を使い

W(-Ze^{Y-Z})=Z(X-1)

$X$についてとくと

X=1+\frac{W(-Ze^{Y-Z})}{Z}

$X=\frac{kx}{mv_0\cos θ}$より

x=\frac{mv_0\cos θ(W(-Ze^{Y-Z})+Z)}{kZ}

Y=-\frac{k^2y}{m^2g}

Z=1+\frac{kv_0\sin θ}{mg}

$Y$と$Z$に代入するのは式が長くなるのでやめます

いざプログラム

言語はJavaScriptにします

function Lambert_W_0(x) { //ランバート(ランベルト)のW関数(0の方)
	var ret = 0;
	for (var i = 0;i < 100;i++) {
		var power = Math.pow(Math.E,ret);
		var little_error = ret*power-x;
		ret = ret-little_error/(power*(ret+1)-(ret+2)*little_error/2/(ret+1));
		if (ret < -1) {
			ret = -1;
		}
	}
	return ret;
}

function compute_x(rad,m,v,g,k,y) {
	if (k == 0) { //空気抵抗なしの時
		var A = g/(v*v*Math.cos(rad)*Math.cos(rad));
		var B = Math.tan(rad);
		return (B+Math.sqrt(B*B+2*A*y))/A;
	}
	if (m == 0) { //NaNが出ないように
		return 0;
	}
	var Y = -k*k*y/(m*m*g);
	var Z = 1+k*v*Math.sin(rad)/(m*g);
	var T = Lambert_W_0(-Z*Math.pow(Math.E,Y-Z))+Z;
	return m*v*Math.cos(rad)*T/(Z*k);
}

シミュレーションでおそらく正しく動作することを確認済み