preamble = () => `
void forward(in float lambda, in float phi, inout float x, inout float y) {
float fu = 1.4, k = 12.;
if (lambda + phi < -fu) {
float u = (lambda - phi + 1.6) * (lambda + phi + fu) / 8.;
lambda += u;
phi -= 0.8 * u * sin(phi + halfPi);
}
// Hammer
float A = 1.68, B = 2.;
float cy = cos(phi), cxcy = cos(lambda / B) * cy;
float K = sqrt(2. / (1. + cxcy));
x = A * K * cy * sin(lambda / B);
y = K * sin(phi);
float d = (1. - cos(lambda * phi)) / k;
if (y < 0.) x *= (1. + d);
if (y > 0.) y *= (1. + d / 1.5 * x * x);
}
void newton2d(in float x, in float y, inout float a, inout float b, in float eps) {
float err2 = 100., da, db, px, py, pax, pay, pbx, pby, tx, ty, eps2 = 1e-8;
bool done = false;
for (int i = 0; i < 15; i++) if (!done) {
forward(a, b, px, py);
// diffs
tx = px - x,
ty = py - y;
if (abs(tx) + abs(ty) < eps2) {
done = true; // we're there!
continue;
}
// backtrack if we overshot
float h = tx * tx + ty * ty;
if (h > err2) {
a -= da /= 2.;
b -= db /= 2.;
continue;
}
err2 = h;
// partial derivatives
// **TODO: there's probably a better solution with dFdx(px)**
float ea = (a > 0. ? -1. : 1.) * eps,
eb = (b > 0. ? -1. : 1.) * eps;
forward(a + ea, b, pax, pay);
forward(a, b + eb, pbx, pby);
float dxa = (pax - px) / ea,
dya = (pay - py) / ea,
dxb = (pbx - px) / eb,
dyb = (pby - py) / eb,
// determinant
D = dyb * dxa - dya * dxb,
// newton step — or half-step for small D
l = (abs(D) < 0.5 ? 0.5 : 1.) / D;
da = (ty * dxb - tx * dyb) * l;
db = (tx * dya - ty * dxa) * l;
a += da;
b += db;
if (abs(da) + abs(db) < eps2) done = true; // we're crawling
}
}
void newton2d(in float x, in float y, inout float a, inout float b) {
newton2d(x, y, a, b, 1e-6);
}
`