import { fix, rsqrt } from "./Math.js";
const eps = 1e-6;
function assert(condition, message) {
if (!condition) {
throw new Error(message ?? "Assertion Error");
}
}
function wrap(ang) {
if (ang < 0) {
ang += 2 * Math.PI;
}
return ang;
}
function XY2ang(XY, nside, dx, dy) {
const xy = XY2xy(XY, nside);
const xyz = xy2xyz(xy, nside, dx, dy);
return xyz2ang(xyz);
}
function ang2XY(ang, nside) {
const xyz = ang2xyz(ang);
const xy = xyz2xy(xyz, nside);
return xy2XY(xy, nside);
}
function XY2xy(XY, nside) {
assert(nside > 0);
assert(XY >= 0);
assert(XY < 12 * nside * nside);
const hp = fix(XY / (nside * nside)),
hp_frac = XY % (nside * nside);
assert(hp >= 0);
assert(hp < 12);
const x = fix(hp_frac / nside);
assert(x >= 0);
assert(x < nside);
const y = hp_frac % nside;
assert(y >= 0);
assert(y < nside);
return { hp: hp, x: x, y: y };
}
function xy2XY(xy, nside) {
assert(nside > 0);
assert(xy.hp >= 0);
assert(xy.hp < 12);
assert(xy.x >= 0);
assert(xy.x < nside);
assert(xy.y >= 0);
assert(xy.y < nside);
return (xy.hp * nside * nside) + (xy.x * nside) + xy.y;
}
function xy2xyz(xy, nside, dx, dy) {
let z, phi;
let equatorial = true,
zfactor = 1;
let hp = xy.hp,
x = xy.x + dx,
y = xy.y + dy;
if (hp <= 3) {
if (x + y > nside) {
equatorial = false;
zfactor = 1;
}
}
if (hp >= 8) {
if (x + y < nside) {
equatorial = false;
zfactor = -1;
}
}
if (equatorial) {
let zoff = 0,
phioff = 0;
x /= nside;
y /= nside;
if (hp <= 3) {
phioff = 1;
} else if (hp <= 7) {
zoff = -1;
hp -= 4;
} else if (hp <= 11) {
phioff = 1;
zoff = -2;
hp -= 8;
} else {
assert(false);
}
z = 2 / 3 * (x + y + zoff);
phi = Math.PI / 4 * (x - y + phioff + 2 * hp);
} else {
let phi_t;
if (zfactor == -1) {
const tmp = x;
x = y;
y = tmp;
x = nside - x;
y = nside - y;
}
if (y == nside && x == nside) {
phi_t = 0;
} else {
phi_t = Math.PI * (nside - y) / (2 * ((nside - x) + (nside - y)));
}
if (phi_t < Math.PI / 4) {
z = 1 - Math.pow(
Math.PI * (nside - x) / ((2 * phi_t - Math.PI) * nside), 2) / 3;
} else {
z = 1 - Math.pow(
Math.PI * (nside - y) / (2 * phi_t * nside), 2) / 3;
}
z *= zfactor;
assert(0 <= Math.abs(z) && Math.abs(z) <= 1);
if (hp >= 8) {
phi = Math.PI / 2 * (hp - 8) + phi_t;
} else {
phi = Math.PI / 2 * hp + phi_t;
}
}
const r = Math.sqrt(1 - z * z);
phi = wrap(phi);
return { x: r * Math.cos(phi), y: r * Math.sin(phi), z: z };
}
function xyz2xy(xyz, nside) {
assert(nside > 0);
let hp, x, y, dx, dy;
let _x, _y;
const phi = wrap(Math.atan2(xyz.y, xyz.x)),
phi_t = phi % (0.5 * Math.PI);
assert(phi_t >= 0);
if ((xyz.z >= 2 / 3) || (xyz.z <= -2 / 3)) {
let north, zfactor;
if (xyz.z >= 2 / 3) {
north = true;
zfactor = 1;
} else {
north = false;
zfactor = -1;
}
const kx = rsqrt(
(1 - xyz.z * zfactor) * 3
* Math.pow(nside * (2 * phi_t - Math.PI) / Math.PI, 2));
const ky = rsqrt(
(1 - xyz.z * zfactor) * 3
* Math.pow(nside * 2 * phi_t / Math.PI, 2));
if (north) {
_x = nside - kx;
_y = nside - ky;
} else {
_x = ky;
_y = kx;
}
x = Math.min(nside - 1, Math.floor(_x));
y = Math.min(nside - 1, Math.floor(_y));
assert(x >= 0);
assert(x < nside);
assert(y >= 0);
assert(y < nside);
dx = _x - x;
dy = _y - y;
let sector = (phi - phi_t) / (0.5 * Math.PI),
offset = Math.round(sector);
assert(Math.abs(sector - offset) < eps);
offset = ((offset % 4) + 4) % 4;
assert(offset >= 0);
assert(offset <= 3);
if (north) {
hp = offset;
} else {
hp = 8 + offset;
}
} else {
const zunits = (xyz.z + 2 / 3) / (4 / 3),
phiunits = phi_t / (0.5 * Math.PI);
const u1 = zunits + phiunits,
u2 = zunits - phiunits + 1;
assert(u1 >= 0);
assert(u1 <= 2);
assert(u2 >= 0);
assert(u2 <= 2);
_x = u1 * nside;
_y = u2 * nside;
let sector = (phi - phi_t) / (0.5 * Math.PI),
offset = Math.round(sector);
assert(Math.abs(sector - offset) < eps);
offset = ((offset % 4) + 4) % 4;
assert(offset >= 0);
assert(offset <= 3);
if (_x >= nside) {
_x -= nside;
if (_y >= nside) {
_y -= nside;
hp = offset;
} else {
hp = ((offset + 1) % 4) + 4;
}
} else {
if (_y >= nside) {
_y -= nside;
hp = offset + 4;
} else {
hp = 8 + offset;
}
}
assert(_x >= -eps);
assert(_x < nside + eps);
assert(_y >= -eps);
assert(_y < nside + eps);
x = Math.max(0, Math.min(nside - 1, Math.floor(_x)));
y = Math.max(0, Math.min(nside - 1, Math.floor(_y)));
assert(x >= 0);
assert(x < nside);
assert(y >= 0);
assert(y < nside);
dx = _x - x;
dy = _y - y;
}
return { hp: hp, x: x, y: y, dx: dx, dy: dy };
}
function xyz2lon(xyz) {
return wrap(Math.atan2(xyz.y, xyz.x));
}
function xyz2lat(xyz) {
return Math.asin(xyz.z);
}
function xyz2ang(xyz) {
return { lon: xyz2lon(xyz), lat: xyz2lat(xyz) };
}
function ang2x(ang) {
return Math.cos(ang.lat) * Math.cos(ang.lon);
}
function ang2y(ang) {
return Math.cos(ang.lat) * Math.sin(ang.lon);
}
function ang2z(ang) {
return Math.sin(ang.lat);
}
function ang2xyz(ang) {
return { x: ang2x(ang), y: ang2y(ang), z: ang2z(ang) };
}
function boundary(XY, nside, n) {
const edge = new Array(n + 1)
.fill()
.map((_, i) => i / n);
const outline = new Array()
.concat(edge.map(dy => XY2ang(XY, nside, 0, dy)))
.concat(edge.map(dx => XY2ang(XY, nside, dx, 1)))
.concat(edge.reverse().map(dy => XY2ang(XY, nside, 1, dy)))
.concat(edge.map(dx => XY2ang(XY, nside, dx, 0)))
.map(d => {
if (d.lon == 0) {
d.lon += eps;
}
return d;
});
return outline;
}
export const HealPix = {
XY2ang: XY2ang,
ang2XY: ang2XY,
boundary: boundary,
};