From 4f73d2e800e9d0bca48fcc4fead20ef7fea511f0 Mon Sep 17 00:00:00 2001 From: Grzegorz Kucmierz Date: Fri, 6 Mar 2026 01:39:02 +0000 Subject: [PATCH] Initial commit: BigInt port of healpixjs --- .gitea/workflows/npm-publish.yaml | 24 + .gitignore | 4 + README.md | 45 ++ index.js | 8 + package-lock.json | 13 + package.json | 33 ++ src/CircleFinder.js | 68 +++ src/Constants.js | 9 + src/Fxyf.js | 57 +++ src/Healpix.js | 781 ++++++++++++++++++++++++++++++ src/Hploc.js | 197 ++++++++ src/Pointing.js | 32 ++ src/RangeSet.js | 66 +++ src/Vec3.js | 117 +++++ src/Xyf.js | 11 + src/Zphi.js | 9 + src/index.js | 12 + src/pstack.js | 45 ++ 18 files changed, 1531 insertions(+) create mode 100644 .gitea/workflows/npm-publish.yaml create mode 100644 .gitignore create mode 100644 README.md create mode 100644 index.js create mode 100644 package-lock.json create mode 100644 package.json create mode 100644 src/CircleFinder.js create mode 100644 src/Constants.js create mode 100644 src/Fxyf.js create mode 100644 src/Healpix.js create mode 100644 src/Hploc.js create mode 100644 src/Pointing.js create mode 100644 src/RangeSet.js create mode 100644 src/Vec3.js create mode 100644 src/Xyf.js create mode 100644 src/Zphi.js create mode 100644 src/index.js create mode 100644 src/pstack.js diff --git a/.gitea/workflows/npm-publish.yaml b/.gitea/workflows/npm-publish.yaml new file mode 100644 index 0000000..0734af9 --- /dev/null +++ b/.gitea/workflows/npm-publish.yaml @@ -0,0 +1,24 @@ +name: npm-publish + +on: + push: + tags: + - 'v*' + workflow_dispatch: + +jobs: + publish: + runs-on: self-hosted + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-node@v4 + with: + node-version: '20.x' + registry-url: 'https://registry.npmjs.org' + - run: npm ci + - name: Publish to npm + run: | + echo "//registry.npmjs.org/:_authToken=${NODE_AUTH_TOKEN}" > .npmrc + npm publish --access public + env: + NODE_AUTH_TOKEN: ${{ secrets.NPM_AUTH_TOKEN }} diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0d545c9 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +node_modules/ +.npmrc +*.log +.DS_Store diff --git a/README.md b/README.md new file mode 100644 index 0000000..ee81903 --- /dev/null +++ b/README.md @@ -0,0 +1,45 @@ +# @gkucmierz/healpixjs-bigint + +![npm (scoped)](https://img.shields.io/npm/v/@gkucmierz/healpixjs-bigint) +![NPM License](https://img.shields.io/npm/l/%40gkucmierz%2Fhealpixjs-bigint) + +A true mathematically sound **BigInt** port of the popular HEALPix spatial mapping algorithm. + +## Motivation: Breaking the 32-bit ceiling 🚀 +The original `healpixjs` engine is a fantastic library. However, under the hood it computes Z-order curves and bitwise coordinate spreads (like `xyf2pix`) using standard JavaScript bitwise operators (`<<`, `>>`, `|`, `&`). + +By specification, JavaScript forces bitwise operators onto **32-bit signed integers**. This means any `Nside` resolution that generates spatial face indices greater than `$2,147,483,647$` mathematically overflows causing the map grids to fracture or the CPU to loop infinitely. + +**HEALPix Order <= 13** is the absolute limit for standard JavaScript math. + +`@gkucmierz/healpixjs-bigint` rewrites the core geometry algorithms using native `BigInt` (e.g., `1n << 20n`). This allows you to construct and query mapping grids at massive depths (virtually tested up to Order 29) without ever triggering `NaN` rendering freezes or out-of-bounds array overflows. + +## Usage + +```bash +npm install @gkucmierz/healpixjs-bigint +``` + +```javascript +// Native ESM architecture +import { Healpix, Pointing } from '@gkucmierz/healpixjs-bigint'; + +// Constructing an massive grid (e.g. Order 20 -> Nside 1048576) +const hp = new Healpix(1n << 20n); + +const ptg = new Pointing(1.57, 0.5); // (theta, phi) + +// Radius is in radians. Fact is the oversampling multiplier (e.g. 4) +const ranges = hp.queryDiscInclusive(ptg, 0.05, 4); + +// ranges.r contains pure BigInt boundaries spanning massive precision +console.log(ranges.r); // BigInt64Array [158671400n, 158671408n, ...] +``` + +## Live Implementation Example + +This library powers extreme sub-meter zoom depths on the [geo-words](https://gitea.7u.pl/gkucmierz/geo-words) interactive 3D map. + +## Credits + +All credit for the structural JavaScript foundation algorithms goes to the original authors of `healpixjs`. This fork modifies it exclusively for precision integer performance in WebGL/Canvas pipelines. diff --git a/index.js b/index.js new file mode 100644 index 0000000..70f98d5 --- /dev/null +++ b/index.js @@ -0,0 +1,8 @@ +// Export the main BigInt modules +export * from './src/Healpix.js'; +export * from './src/Pointing.js'; +export * from './src/RangeSet.js'; +export * from './src/Vec3.js'; +export * from './src/Hploc.js'; +export * from './src/Constants.js'; +export * from './src/Fxyf.js'; diff --git a/package-lock.json b/package-lock.json new file mode 100644 index 0000000..0825074 --- /dev/null +++ b/package-lock.json @@ -0,0 +1,13 @@ +{ + "name": "@gkucmierz/healpixjs-bigint", + "version": "1.0.0", + "lockfileVersion": 3, + "requires": true, + "packages": { + "": { + "name": "@gkucmierz/healpixjs-bigint", + "version": "1.0.0", + "license": "MIT" + } + } +} diff --git a/package.json b/package.json new file mode 100644 index 0000000..a8d73b3 --- /dev/null +++ b/package.json @@ -0,0 +1,33 @@ +{ + "name": "@gkucmierz/healpixjs-bigint", + "version": "1.0.0", + "description": "True BigInt port of healpixjs supporting orders beyond the 32-bit (<=14) limit", + "keywords": [ + "healpix", + "astronomy", + "math", + "bigint", + "64-bit", + "3d", + "sphere", + "projection", + "geography", + "mapping", + "precision" + ], + "files": [ + "index.js", + "src/" + ], + "main": "index.js", + "type": "module", + "scripts": { + "test": "echo \"Error: no test specified\" && exit 1" + }, + "repository": { + "type": "git", + "url": "git+https://gitea.7u.pl/gkucmierz/healpixjs-bigint.git" + }, + "author": "Grzegorz Kućmierz", + "license": "MIT" +} \ No newline at end of file diff --git a/src/CircleFinder.js b/src/CircleFinder.js new file mode 100644 index 0000000..35639e6 --- /dev/null +++ b/src/CircleFinder.js @@ -0,0 +1,68 @@ +import { Vec3 } from './Vec3.js'; +export class CircleFinder { + /** + * @param point: Vec3 + */ + constructor(point) { + let np = point.length; + //HealpixUtils.check(np>=2,"too few points"); + if (!(np >= 2)) { + console.log("too few points"); + return; + } + this.center = point[0].add(point[1]); + this.center.normalize(); + this.cosrad = point[0].dot(this.center); + for (let i = 2; i < np; ++i) { + if (point[i].dot(this.center) < this.cosrad) { // point outside the current circle + this.getCircle(point, i); + } + } + } + ; + /** + * @parm point: Vec3 + * @param q: int + */ + getCircle(point, q) { + this.center = point[0].add(point[q]); + this.center.normalize(); + this.cosrad = point[0].dot(this.center); + for (let i = 1; i < q; ++i) { + if (point[i].dot(this.center) < this.cosrad) { // point outside the current circle + this.getCircle2(point, i, q); + } + } + } + ; + /** + * @parm point: Vec3 + * @param q1: int + * @param q2: int + */ + getCircle2(point, q1, q2) { + this.center = point[q1].add(point[q2]); + this.center.normalize(); + this.cosrad = point[q1].dot(this.center); + for (let i = 0; i < q1; ++i) { + if (point[i].dot(this.center) < this.cosrad) { // point outside the current circle + this.center = (point[q1].sub(point[i])).cross(point[q2].sub(point[i])); + this.center.normalize(); + this.cosrad = point[i].dot(this.center); + if (this.cosrad < 0) { + this.center.flip(); + this.cosrad = -this.cosrad; + } + } + } + } + ; + getCenter() { + return new Vec3(this.center.x, this.center.y, this.center.z); + } + getCosrad() { + return this.cosrad; + } + ; +} +//# sourceMappingURL=CircleFinder.js.map \ No newline at end of file diff --git a/src/Constants.js b/src/Constants.js new file mode 100644 index 0000000..79f4a73 --- /dev/null +++ b/src/Constants.js @@ -0,0 +1,9 @@ +export class Constants { +} +// static halfpi = Math.PI/2.; +Constants.halfpi = 1.5707963267948966; +Constants.inv_halfpi = 2. / Math.PI; +/** The Constant twopi. */ +Constants.twopi = 2 * Math.PI; +Constants.inv_twopi = 1. / (2 * Math.PI); +//# sourceMappingURL=Constants.js.map \ No newline at end of file diff --git a/src/Fxyf.js b/src/Fxyf.js new file mode 100644 index 0000000..3af6725 --- /dev/null +++ b/src/Fxyf.js @@ -0,0 +1,57 @@ +/** + * Partial porting to Javascript of Fxyf.java from Healpix3.30 + */ +import { Hploc } from './Hploc.js'; +export class Fxyf { + constructor(x, y, f) { + this.fx = x; + this.fy = y; + this.face = f; + // coordinate of the lowest corner of each face + this.jrll = new Uint8Array([2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4]); + this.jpll = new Uint8Array([1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7]); + this.halfpi = Math.PI / 2.; + } + toHploc() { + let loc = new Hploc(); + let jr = this.jrll[this.face] - this.fx - this.fy; + let nr; + if (jr < 1) { + nr = jr; + let tmp = nr * nr / 3.; + loc.z = 1 - tmp; + if (loc.z > 0.99) { + loc.sth = Math.sqrt(tmp * (2.0 - tmp)); + loc.have_sth = true; + } + } + else if (jr > 3) { + nr = 4 - jr; + let tmp = nr * nr / 3.; + loc.z = tmp - 1; + if (loc.z < -0.99) { + loc.sth = Math.sqrt(tmp * (2.0 - tmp)); + loc.have_sth = true; + } + } + else { + nr = 1; + loc.z = (2 - jr) * 2.0 / 3.; + } + let tmp = this.jpll[this.face] * nr + this.fx - this.fy; + if (tmp < 0) { + tmp += 8; + } + if (tmp >= 8) { + tmp -= 8; + } + loc.phi = (nr < 1e-15) ? 0 : (0.5 * this.halfpi * tmp) / nr; + return loc; + } + ; + toVec3() { + return this.toHploc().toVec3(); + } + ; +} +//# sourceMappingURL=Fxyf.js.map \ No newline at end of file diff --git a/src/Healpix.js b/src/Healpix.js new file mode 100644 index 0000000..157b592 --- /dev/null +++ b/src/Healpix.js @@ -0,0 +1,781 @@ +"use strict"; +import { CircleFinder } from "./CircleFinder.js"; +import { Constants } from "./Constants.js"; +import { Fxyf } from "./Fxyf.js"; +import { Hploc } from "./Hploc.js"; +import { Pointing } from "./Pointing.js"; +import { pstack } from "./pstack.js"; +import { RangeSet } from "./RangeSet.js"; +import { Vec3 } from "./Vec3.js"; +import { Xyf } from "./Xyf.js"; +import { Zphi } from "./Zphi.js"; +/** + * Partial porting to Javascript of HealpixBase.java from Healpix3.30 + */ +// import Fxyf from './Fxyf.js'; +// import Hploc from './Hploc.js'; +// import Xyf from './Xyf.js'; +// import Vec3 from './Vec3.js'; +// import Pointing from './Pointing.js'; +// import CircleFinder from './CircleFinder.js'; +// import Zphi from './Zphi.js'; +// import pstack from './pstack.js'; +// import Constants from './Constants.js'; +// import RangeSet from './RangeSet.js'; +export class Healpix { + constructor(nside_in) { + this.order_max = 29; + this.inv_halfpi = 2.0 / Math.PI; + this.twothird = 2.0 / 3.; + // console.log("twothird "+this.twothird); + // this.ns_max=1L< 0) { + this.nside = nside_in; + this.npface = BigInt(this.nside) * BigInt(this.nside); + this.npix = 12n * this.npface; + this.order = this.nside2order(this.nside); + this.nl2 = 2 * this.nside; + this.nl3 = 3 * this.nside; + this.nl4 = 4 * this.nside; + this.fact2 = 4.0 / Number(this.npix); + this.fact1 = (this.nside << 1) * this.fact2; + this.ncap = 2 * this.nside * (this.nside - 1); // pixels in each polar cap + // console.log("order: "+this.order); + // console.log("nside: "+this.nside); + } + this.bn = []; + this.mpr = []; + this.cmpr = []; + this.smpr = []; + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // TODO INFINITE LOOP!!!!!! FIX ITTTTTTTTTT + // Uncaught RangeError: Maximum call stack size exceeded + // MOVED TO computeBn() + // for (let i=0; i <= this.order_max; ++i) { + // this.bn[i]=new Healpix(1< 0) && (ix < nsm1) && (iy > 0) && (iy < nsm1)) { + let fpix = Math.floor(face_num << (2 * this.order)); + let px0 = this.spread_bits(ix); + let py0 = this.spread_bits(iy) << 1; + let pxp = this.spread_bits(ix + 1); + let pyp = this.spread_bits(iy + 1) << 1; + let pxm = this.spread_bits(ix - 1); + let pym = this.spread_bits(iy - 1) << 1; + result[0] = fpix + pxm + py0; + result[1] = fpix + pxm + pyp; + result[2] = fpix + px0 + pyp; + result[3] = fpix + pxp + pyp; + result[4] = fpix + pxp + py0; + result[5] = fpix + pxp + pym; + result[6] = fpix + px0 + pym; + result[7] = fpix + pxm + pym; + } + else { + for (let i = 0; i < 8; ++i) { + let x = ix + this.xoffset[i]; + let y = iy + this.yoffset[i]; + let nbnum = 4; + if (x < 0) { + x += this.nside; + nbnum -= 1; + } + else if (x >= this.nside) { + x -= this.nside; + nbnum += 1; + } + if (y < 0) { + y += this.nside; + nbnum -= 3; + } + else if (y >= this.nside) { + y -= this.nside; + nbnum += 3; + } + let f = this.facearray[nbnum][face_num]; + if (f >= 0) { + let bits = this.swaparray[nbnum][face_num >>> 2]; + if ((bits & 1) > 0) { + x = Math.floor(this.nside - x - 1); + } + if ((bits & 2) > 0) { + y = Math.floor(this.nside - y - 1); + } + if ((bits & 4) > 0) { + let tint = x; + x = y; + y = tint; + } + result[i] = this.xyf2nest(x, y, f); + } + else { + result[i] = -1; + } + } + } + return result; + } + ; + nside2order(nside) { + return ((nside & (nside - 1)) != 0) ? -1 : Math.log2(nside); + } + ; + nest2xyf(ipix) { + ipix = BigInt(ipix); + let pix = ipix & (this.npface - 1n); + let xyf = new Xyf(Number(this.compress_bits(pix)), Number(this.compress_bits(pix >> 1n)), Number(ipix >> (2n * BigInt(this.order)))); + return xyf; + } + ; + xyf2nest(ix, iy, face_num) { + return (BigInt(face_num) << (2n * BigInt(this.order))) + + this.spread_bits(BigInt(ix)) + (this.spread_bits(BigInt(iy)) << 1n); + } + ; + loc2pix(hploc) { + let z = hploc.z; + let phi = hploc.phi; + let za = Math.abs(z); + let tt = this.fmodulo((phi * this.inv_halfpi), 4.0); // in [0,4) + let pixNo; + if (za <= this.twothird) { // Equatorial region + let temp1 = this.nside * (0.5 + tt); + let temp2 = this.nside * (z * 0.75); + let jp = Math.floor(temp1 - temp2); // index of ascending edge line + let jm = Math.floor(temp1 + temp2); // index of descending edge line + let ifp = Math.floor(jp >>> this.order); // in {0,4} + let ifm = Math.floor(jm >>> this.order); + let face_num = Math.floor((ifp == ifm) ? (ifp | 4) : ((ifp < ifm) ? ifp : (ifm + 8))); + let ix = Math.floor(jm & (this.nside - 1)); + let iy = Math.floor(this.nside - (jp & (this.nside - 1)) - 1); + pixNo = this.xyf2nest(ix, iy, face_num); + } + else { // polar region, za > 2/3 + let ntt = Math.min(3, Math.floor(tt)); + let tp = tt - ntt; + let tmp = ((za < 0.99) || (!hploc.have_sth)) ? + this.nside * Math.sqrt(3 * (1 - za)) : + this.nside * hploc.sth / Math.sqrt((1.0 + za) / 3.); + let jp = Math.floor(tp * tmp); // increasing edge line index + let jm = Math.floor((1.0 - tp) * tmp); // decreasing edge line index + if (jp >= this.nside) { + jp = this.nside - 1; // for points too close to the boundary + } + if (jm >= this.nside) { + jm = this.nside - 1; + } + if (z >= 0) { + pixNo = this.xyf2nest(Math.floor(this.nside - jm - 1), Math.floor(this.nside - jp - 1), ntt); + } + else { + pixNo = this.xyf2nest(Math.floor(jp), Math.floor(jm), ntt + 8); + } + } + return pixNo; + } + ; + /** Returns the normalized 3-vector corresponding to the center of the + supplied pixel. + @param pix long the requested pixel number. + @return the pixel's center coordinates. */ + pix2vec(pix) { + return this.pix2loc(pix).toVec3(); + } + ; + /** Returns the Zphi corresponding to the center of the supplied pixel. + @param pix the requested pixel number. + @return the pixel's center coordinates. */ + pix2zphi(pix) { + return this.pix2loc(pix).toZphi(); + } + pix2ang(pix, mirror) { + return this.pix2loc(pix).toPointing(mirror); + } + /** + * @param pix long + * @return Hploc + */ + pix2loc(pix) { + let loc = new Hploc(undefined); + let xyf = this.nest2xyf(pix); + let jr = ((this.jrll[xyf.face]) << this.order) - xyf.ix - xyf.iy - 1; + let nr; + if (jr < this.nside) { + nr = jr; + let tmp = (nr * nr) * this.fact2; + loc.z = 1 - tmp; + if (loc.z > 0.99) { + loc.sth = Math.sqrt(tmp * (2. - tmp)); + loc.have_sth = true; + } + } + else if (jr > this.nl3) { + nr = this.nl4 - jr; + let tmp = (nr * nr) * this.fact2; + loc.z = tmp - 1; + if (loc.z < -0.99) { + loc.sth = Math.sqrt(tmp * (2. - tmp)); + loc.have_sth = true; + } + } + else { + nr = this.nside; + loc.z = (this.nl2 - jr) * this.fact1; + } + let tmp = (this.jpll[xyf.face]) * nr + xyf.ix - xyf.iy; + // assert(tmp<8*nr); // must not happen + if (tmp < 0) { + tmp += 8 * nr; + } + loc.phi = (nr == this.nside) ? 0.75 * Constants.halfpi * tmp * this.fact1 : (0.5 * Constants.halfpi * tmp) / nr; + // loc.setPhi((nr == this.nside) ? 0.75 * Constants.halfpi * tmp * this.fact1 : (0.5 * Constants.halfpi * tmp)/nr); + return loc; + } + ; + za2vec(z, a) { + const sin_theta = Math.sqrt(1 - z * z); + const X = sin_theta * Math.cos(a); + const Y = sin_theta * Math.sin(a); + return new Vec3(X, Y, z); + } + ang2vec(theta, phi) { + const z = Math.cos(theta); + return this.za2vec(z, phi); + } + vec2ang(v) { + const { z, a } = this.vec2za(v.getX(), v.getY(), v.getZ()); + return { theta: Math.acos(z), phi: a }; + } + vec2za(X, Y, z) { + const r2 = X * X + Y * Y; + if (r2 == 0) + return { z: z < 0 ? -1 : 1, a: 0 }; + else { + const PI2 = Math.PI / 2; + const a = (Math.atan2(Y, X) + PI2) % PI2; + z /= Math.sqrt(z * z + r2); + return { z, a }; + } + } + ang2pix(ptg, mirror) { + return this.loc2pix(new Hploc(ptg)); + } + ; + fmodulo(v1, v2) { + if (v1 >= 0) { + return (v1 < v2) ? v1 : v1 % v2; + } + var tmp = v1 % v2 + v2; + return (tmp === v2) ? 0.0 : tmp; + } + ; + compress_bits(v) { + v = BigInt(v); + let raw = (v & 0x5555n) | ((v & 0x55550000n) >> 15n); + let compressed = BigInt(this.ctab[Number(raw & 0xffn)]) | (BigInt(this.ctab[Number(raw >> 8n)]) << 4n); + return compressed; + } + ; + spread_bits(v) { + v = BigInt(v); + return BigInt(this.utab[Number(v & 0xffn)]) + | (BigInt(this.utab[Number((v >> 8n) & 0xffn)]) << 16n) + | (BigInt(this.utab[Number((v >> 16n) & 0xffn)]) << 32n) + | (BigInt(this.utab[Number((v >> 24n) & 0xffn)]) << 48n); + } + ; + /** + * Returns a range set of pixels that overlap with the convex polygon + * defined by the {@code vertex} array. + *

+ * This method is more efficient in the RING scheme. + *

+ * This method may return some pixels which don't overlap with the polygon + * at all. The higher {@code fact} is chosen, the fewer false positives are + * returned, at the cost of increased run time. + * + * @param vertex + * an array containing the vertices of the requested convex + * polygon. + * @param fact + * The overlapping test will be done at the resolution + * {@code fact*nside}. For NESTED ordering, {@code fact} must be + * a power of 2, else it can be any positive integer. A typical + * choice would be 4. + * @return the requested set of pixel number ranges + */ + queryPolygonInclusive(vertex, fact) { + let inclusive = (fact != 0); + let nv = vertex.length; + // let ncirc = inclusive ? nv+1 : nv; + if (!(nv >= 3)) { + console.log("not enough vertices in polygon"); + return; + } + let vv = new Array(); + for (let i = 0; i < nv; ++i) { + vv[i] = Vec3.pointing2Vec3(vertex[i]); + } + let normal = new Array(); + let flip = 0; + let index = 0; + let back = false; + while (index < vv.length) { + let first = vv[index]; + let medium = null; + let last = null; + if (index == vv.length - 1) { + last = vv[1]; + medium = vv[0]; + } + else if (index == vv.length - 2) { + last = vv[0]; + medium = vv[index + 1]; + } + else { + medium = vv[index + 1]; + last = vv[index + 2]; + } + normal[index] = first.cross(medium).norm(); + let hnd = normal[index].dot(last); + if (index == 0) { + flip = (hnd < 0.) ? -1 : 1; + let tmp = new Pointing(first); // TODO not used + back = false; + } + else { + let flipThnd = flip * hnd; + if (flipThnd < 0) { + let tmp = new Pointing(medium); + vv.splice(index + 1, 1); + normal.splice(index, 1); + back = true; + index -= 1; + continue; + } + else { + let tmp = new Pointing(first); + back = false; + } + } + normal[index].scale(flip); + index += 1; + } + nv = vv.length; + let ncirc = inclusive ? nv + 1 : nv; + let rad = new Array(ncirc); + rad = rad.fill(Constants.halfpi); + // rad = rad.fill(1.5707963267948966); + // let p = "1.5707963267948966"; + // rad = rad.fill(parseFloat(p)); + if (inclusive) { + let cf = new CircleFinder(vv); + normal[nv] = cf.getCenter(); + rad[nv] = Hploc.acos(cf.getCosrad()); + } + return this.queryMultiDisc(normal, rad, fact); + } + ; + /** + * For NEST schema only + * + * @param normal: + * Vec3[] + * @param rad: + * Float32Array + * @param fact: + * The overlapping test will be done at the resolution + * {@code fact*nside}. For NESTED ordering, {@code fact} must be + * a power of 2, else it can be any positive integer. A typical + * choice would be 4. + * @return RangeSet the requested set of pixel number ranges + */ + queryMultiDisc(norm, rad, fact) { + this.computeBn(); + let inclusive = (fact != 0); + let nv = norm.length; + // HealpixUtils.check(nv==rad.lengt0,"inconsistent input arrays"); + if (!(nv == rad.length)) { + console.error("inconsistent input arrays"); + return; + } + let res = new RangeSet(4 << 1); + // Removed code for Scheme.RING + let oplus = 0; + if (inclusive) { + if (!(Math.pow(2, this.order_max - this.order) >= fact)) { + console.error("invalid oversampling factor"); + } + if (!((fact & (fact - 1)) == 0)) { + console.error("oversampling factor must be a power of 2"); + } + oplus = this.ilog2(fact); + } + let omax = this.order + oplus; // the order up to which we test + // TODO: ignore all disks with radius>=pi + // let crlimit = new Float32Array[omax+1][nv][3]; + let crlimit = new Array(omax + 1); + let o; + let i; + for (o = 0; o <= omax; ++o) { // prepare data at the required orders + crlimit[o] = new Array(nv); + let dr = this.bn[o].maxPixrad(); // safety distance + for (i = 0; i < nv; ++i) { + crlimit[o][i] = new Float64Array(3); + crlimit[o][i][0] = (rad[i] + dr > Math.PI) ? -1 : Hploc.cos(rad[i] + dr); + crlimit[o][i][1] = (o == 0) ? Hploc.cos(rad[i]) : crlimit[0][i][1]; + crlimit[o][i][2] = (rad[i] - dr < 0.) ? 1. : Hploc.cos(rad[i] - dr); + } + } + let stk = new pstack(12 + 3 * omax); + for (let i = 0n; i < 12n; i++) { // insert the 12 base pixels in reverse + // order + stk.push(11n - i, 0); + } + while (stk.size() > 0) { // as long as there are pixels on the stack + // pop current pixel number and order from the stack + let pix = stk.ptop(); + let o = stk.otop(); + stk.pop(); + let pv = this.bn[o].pix2vec(pix); + let zone = 3; + for (let i = 0; (i < nv) && (zone > 0); ++i) { + let crad = pv.dot(norm[i]); + for (let iz = 0; iz < zone; ++iz) { + if (crad < crlimit[o][i][iz]) { + zone = iz; + } + } + } + if (zone > 0) { + this.check_pixel(o, omax, zone, res, pix, stk, inclusive); + } + } + return res; + } + ; + /** Integer base 2 logarithm. + @param arg + @return the largest integer {@code n} that fulfills {@code 2^n<=arg}. + For negative arguments and zero, 0 is returned. */ + ilog2(arg) { + let max = Math.max(arg, 1); + return 31 - Math.clz32(max); + } + ; + /** Computes the cosine of the angular distance between two z, phi positions + on the unit sphere. */ + cosdist_zphi(z1, phi1, z2, phi2) { + return z1 * z2 + Hploc.cos(phi1 - phi2) * Math.sqrt((1.0 - z1 * z1) * (1.0 - z2 * z2)); + } + /** + * @param int o + * @param int omax + * @param int zone + * @param RangeSet pixset + * @param long pix + * @param pstack stk + * @param boolean inclusive + */ + check_pixel(o, omax, zone, pixset, pix, stk, inclusive) { + if (zone == 0) + return; + if (o < this.order) { + if (zone >= 3) { // output all subpixels + let sdist = 2n * BigInt(this.order - o); // the "bit-shift distance" between map orders + pixset.append1(pix << sdist, ((pix + 1n) << sdist)); + } + else { // (zone>=1) + for (let i = 0n; i < 4n; ++i) { + stk.push(4n * pix + 3n - i, o + 1); // add children + } + } + } + else if (o > this.order) { // this implies that inclusive==true + if (zone >= 2) { // pixel center in shape + pixset.append(pix >> (2n * BigInt(o - this.order))); // output the parent pixel at order + stk.popToMark(); // unwind the stack + } + else { // (zone>=1): pixel center in safety range + if (o < omax) { // check sublevels + for (let i = 0n; i < 4n; ++i) { // add children in reverse order + stk.push(4n * pix + 3n - i, o + 1); // add children + } + } + else { // at resolution limit + pixset.append(pix >> (2n * BigInt(o - this.order))); // output the parent pixel at order + stk.popToMark(); // unwind the stack + } + } + } + else { // o==order + if (zone >= 2) { + pixset.append(pix); + } + else if (inclusive) { // and (zone>=1) + if (this.order < omax) { // check sublevels + stk.mark(); // remember current stack position + for (let i = 0n; i < 4n; ++i) { // add children in reverse order + stk.push(4n * pix + 3n - i, o + 1); // add children + } + } + else { // at resolution limit + pixset.append(pix); // output the pixel + } + } + } + } + /** Returns the maximum angular distance between a pixel center and its + corners. + @return maximum angular distance between a pixel center and its + corners. */ + maxPixrad() { + let zphia = new Zphi(2. / 3., Math.PI / this.nl4); + let xyz1 = this.convertZphi2xyz(zphia); + let va = new Vec3(xyz1[0], xyz1[1], xyz1[2]); + let t1 = 1. - 1. / this.nside; + t1 *= t1; + let zphib = new Zphi(1 - t1 / 3, 0); + let xyz2 = this.convertZphi2xyz(zphib); + let vb = new Vec3(xyz2[0], xyz2[1], xyz2[2]); + return va.angle(vb); + } + ; + /** + * this is a workaround replacing the Vec3(Zphi) constructor. + */ + convertZphi2xyz(zphi) { + let sth = Math.sqrt((1.0 - zphi.z) * (1.0 + zphi.z)); + let x = sth * Hploc.cos(zphi.phi); + let y = sth * Hploc.sin(zphi.phi); + let z = zphi.z; + return [x, y, z]; + } + ; + /** Returns a range set of pixels which overlap with a given disk.

+ This method is more efficient in the RING scheme.

+ This method may return some pixels which don't overlap with + the polygon at all. The higher {@code fact} is chosen, the fewer false + positives are returned, at the cost of increased run time. + @param ptg the angular coordinates of the disk center + @param radius the radius (in radians) of the disk + @param fact The overlapping test will be done at the resolution + {@code fact*nside}. For NESTED ordering, {@code fact} must be a power + of 2, else it can be any positive integer. A typical choice would be 4. + @return the requested set of pixel number ranges */ + queryDiscInclusive(ptg, radius, fact) { + this.computeBn(); + let inclusive = (fact != 0); + let pixset = new RangeSet(); + if (radius >= Math.PI) { // disk covers the whole sphere + pixset.append1(0n, this.npix); + return pixset; + } + let oplus = 0; + if (inclusive) { + // HealpixUtils.check ((1L<=fact,"invalid oversampling factor"); + if (!((fact & (fact - 1)) == 0)) { + console.error("oversampling factor must be a power of 2"); + } + oplus = this.ilog2(fact); + } + let omax = Math.min(this.order_max, this.order + oplus); // the order up to which we test + let vptg = Vec3.pointing2Vec3(ptg); + let crpdr = new Array(omax + 1); + let crmdr = new Array(omax + 1); + let cosrad = Hploc.cos(radius); + let sinrad = Hploc.sin(radius); + for (let o = 0; o <= omax; o++) { // prepare data at the required orders + let dr = this.mpr[o]; // safety distance + let cdr = this.cmpr[o]; + let sdr = this.smpr[o]; + crpdr[o] = (radius + dr > Math.PI) ? -1. : cosrad * cdr - sinrad * sdr; + crmdr[o] = (radius - dr < 0.) ? 1. : cosrad * cdr + sinrad * sdr; + } + let stk = new pstack(12 + 3 * omax); + for (let i = 0n; i < 12n; i++) { // insert the 12 base pixels in reverse order + stk.push(11n - i, 0); + } + while (stk.size() > 0) { // as long as there are pixels on the stack + // pop current pixel number and order from the stack + let pix = stk.ptop(); + let curro = stk.otop(); + stk.pop(); + let pos = this.bn[curro].pix2zphi(pix); + // cosine of angular distance between pixel center and disk center + let cangdist = this.cosdist_zphi(vptg.z, ptg.phi, pos.z, pos.phi); + if (cangdist > crpdr[curro]) { + let zone = (cangdist < cosrad) ? 1 : ((cangdist <= crmdr[curro]) ? 2 : 3); + this.check_pixel(curro, omax, zone, pixset, pix, stk, inclusive); + } + } + return pixset; + } +} +//# sourceMappingURL=Healpix.js.map \ No newline at end of file diff --git a/src/Hploc.js b/src/Hploc.js new file mode 100644 index 0000000..16c75b8 --- /dev/null +++ b/src/Hploc.js @@ -0,0 +1,197 @@ +import { Pointing } from './Pointing.js'; +import { Vec3 } from './Vec3.js'; +import { Zphi } from './Zphi.js'; +export class Hploc { + constructor(ptg) { + Hploc.PI4_A = 0.7853981554508209228515625; + Hploc.PI4_B = 0.794662735614792836713604629039764404296875e-8; + Hploc.PI4_C = 0.306161699786838294306516483068750264552437361480769e-16; + Hploc.M_1_PI = 0.3183098861837906715377675267450287; + if (ptg) { + this.sth = 0.0; + this.have_sth = false; + this.z = Hploc.cos(ptg.theta); + this._phi = ptg.phi; + if (Math.abs(this.z) > 0.99) { + this.sth = Hploc.sin(ptg.theta); + this.have_sth = true; + } + } + } + setZ(z) { + this.z = z; + } + ; + get phi() { + return this._phi; + } + ; + set phi(phi) { + this._phi = phi; + } + ; + setSth(sth) { + this.sth = sth; + } + ; + toPointing(mirror) { + const st = this.have_sth ? this.sth : Math.sqrt((1.0 - this.z) * (1.0 + this.z)); + return new Pointing(null, false, Hploc.atan2(st, this.z), this._phi); + } + toVec3() { + var st = this.have_sth ? this.sth : Math.sqrt((1.0 - this.z) * (1.0 + this.z)); + var vector = new Vec3(st * Hploc.cos(this.phi), st * Hploc.sin(this.phi), this.z); + // var vector = new Vec3(st*Math.cos(this.phi),st*Math.sin(this.phi),this.z); + return vector; + } + ; + toZphi() { + return new Zphi(this.z, this.phi); + } + static sin(d) { + let u = d * Hploc.M_1_PI; + let q = Math.floor(u < 0 ? u - 0.5 : u + 0.5); + let x = 4.0 * q; + d -= x * Hploc.PI4_A; + d -= x * Hploc.PI4_B; + d -= x * Hploc.PI4_C; + if ((q & 1) != 0) { + d = -d; + } + return this.sincoshelper(d); + } + ; + static cos(d) { + // let u = d * Hploc.M_1_PI - 0.5; + let u = d * Hploc.M_1_PI - 0.5; + // u -= 0.5; + let q = 1 + 2 * Math.floor(u < 0 ? u - 0.5 : u + 0.5); + let x = 2.0 * q; + let t = x * Hploc.PI4_A; + d = d - t; + d -= x * Hploc.PI4_B; + d -= x * Hploc.PI4_C; + if ((q & 2) == 0) { + d = -d; + } + return Hploc.sincoshelper(d); + } + ; + static sincoshelper(d) { + let s = d * d; + let u = -7.97255955009037868891952e-18; + u = u * s + 2.81009972710863200091251e-15; + u = u * s - 7.64712219118158833288484e-13; + u = u * s + 1.60590430605664501629054e-10; + u = u * s - 2.50521083763502045810755e-08; + u = u * s + 2.75573192239198747630416e-06; + u = u * s - 0.000198412698412696162806809; + u = u * s + 0.00833333333333332974823815; + u = u * s - 0.166666666666666657414808; + return s * u * d + d; + } + ; + /** This method calculates the arc sine of x in radians. The return + value is in the range [-pi/2, pi/2]. The results may have + maximum error of 3 ulps. */ + static asin(d) { + return Hploc.mulsign(Hploc.atan2k(Math.abs(d), Math.sqrt((1 + d) * (1 - d))), d); + } + ; + /** This method calculates the arc cosine of x in radians. The + return value is in the range [0, pi]. The results may have + maximum error of 3 ulps. */ + static acos(d) { + return Hploc.mulsign(Hploc.atan2k(Math.sqrt((1 + d) * (1 - d)), Math.abs(d)), d) + (d < 0 ? Math.PI : 0); + } + ; + static mulsign(x, y) { + let sign = Hploc.copySign(1, y); + return sign * x; + } + ; + static copySign(magnitude, sign) { + return sign < 0 ? -Math.abs(magnitude) : Math.abs(magnitude); + // let finalsign = 1; + // if (Object.is(finalsign , -0)){ + // sign = -1; + // }else if (Object.is(finalsign , 0)){ + // sign = 1; + // }else { + // sign = Math.sign(finalsign); + // } + // return finalsign * magnitude; + } + static atanhelper(s) { + let t = s * s; + let u = -1.88796008463073496563746e-05; + u = u * t + (0.000209850076645816976906797); + u = u * t + (-0.00110611831486672482563471); + u = u * t + (0.00370026744188713119232403); + u = u * t + (-0.00889896195887655491740809); + u = u * t + (0.016599329773529201970117); + u = u * t + (-0.0254517624932312641616861); + u = u * t + (0.0337852580001353069993897); + u = u * t + (-0.0407629191276836500001934); + u = u * t + (0.0466667150077840625632675); + u = u * t + (-0.0523674852303482457616113); + u = u * t + (0.0587666392926673580854313); + u = u * t + (-0.0666573579361080525984562); + u = u * t + (0.0769219538311769618355029); + u = u * t + (-0.090908995008245008229153); + u = u * t + (0.111111105648261418443745); + u = u * t + (-0.14285714266771329383765); + u = u * t + (0.199999999996591265594148); + u = u * t + (-0.333333333333311110369124); + return u * t * s + s; + } + ; + static atan2k(y, x) { + let q = 0.; + if (x < 0) { + x = -x; + q = -2.; + } + if (y > x) { + let t = x; + x = y; + y = -t; + q += 1.; + } + return Hploc.atanhelper(y / x) + q * (Math.PI / 2); + } + ; + /** This method calculates the arc tangent of y/x in radians, using + the signs of the two arguments to determine the quadrant of the + result. The results may have maximum error of 2 ulps. */ + static atan2(y, x) { + let r = Hploc.atan2k(Math.abs(y), x); + r = Hploc.mulsign(r, x); + if (Hploc.isinf(x) || x == 0) { + r = Math.PI / 2 - (Hploc.isinf(x) ? (Hploc.copySign(1, x) * (Math.PI / 2)) : 0); + } + if (Hploc.isinf(y)) { + r = Math.PI / 2 - (Hploc.isinf(x) ? (Hploc.copySign(1, x) * (Math.PI * 1 / 4)) : 0); + } + if (y == 0) { + r = (Hploc.copySign(1, x) == -1 ? Math.PI : 0); + } + return Hploc.isnan(x) || Hploc.isnan(y) ? NaN : Hploc.mulsign(r, y); + } + ; + /** Checks if the argument is a NaN or not. */ + static isnan(d) { + return d != d; + } + ; + /** Checks if the argument is either positive or negative infinity. */ + static isinf(d) { + return Math.abs(d) === +Infinity; + } + ; +} +Hploc.PI4_A = 0.7853981554508209228515625; +Hploc.PI4_B = 0.794662735614792836713604629039764404296875e-8; +Hploc.PI4_C = 0.306161699786838294306516483068750264552437361480769e-16; +Hploc.M_1_PI = 0.3183098861837906715377675267450287; +//# sourceMappingURL=Hploc.js.map \ No newline at end of file diff --git a/src/Pointing.js b/src/Pointing.js new file mode 100644 index 0000000..d65b2c3 --- /dev/null +++ b/src/Pointing.js @@ -0,0 +1,32 @@ +import { Hploc } from './Hploc.js'; +export class Pointing { + /** + * + * @param {*} vec3 Vec3.js + * @param {*} mirror + * @param {*} in_theta radians + * @param {*} in_phi radians + */ + constructor(vec3, mirror, in_theta, in_phi) { + if (vec3 != null) { + this.theta = Hploc.atan2(Math.sqrt(vec3.x * vec3.x + vec3.y * vec3.y), vec3.z); + if (mirror) { + this.phi = -Hploc.atan2(vec3.y, vec3.x); + } + else { + this.phi = Hploc.atan2(vec3.y, vec3.x); + } + if (this.phi < 0.0) { + this.phi = this.phi + 2 * Math.PI; + } + if (this.phi >= 2 * Math.PI) { + this.phi = this.phi - 2 * Math.PI; + } + } + else { + this.theta = in_theta; + this.phi = in_phi; + } + } +} +//# sourceMappingURL=Pointing.js.map \ No newline at end of file diff --git a/src/RangeSet.js b/src/RangeSet.js new file mode 100644 index 0000000..23e1160 --- /dev/null +++ b/src/RangeSet.js @@ -0,0 +1,66 @@ +export class RangeSet { + /** + * @param int cap: initial capacity + */ + constructor(cap) { + if (cap < 0) + console.error("capacity must be positive"); + this.r = new BigInt64Array(cap << 1); + this.sz = 0; + } + ; + /** Append a single-value range to the object. + @param val value to append */ + append(val) { + this.append1(val, val + 1n); + } + ; + /** Append a range to the object. + @param a first long in range + @param b one-after-last long in range */ + append1(a, b) { + if (a >= b) + return; + if ((this.sz > 0) && (a <= this.r[this.sz - 1])) { + if (a < this.r[this.sz - 2]) + console.error("bad append operation"); + if (b > this.r[this.sz - 1]) + this.r[this.sz - 1] = b; + return; + } + // this.ensureCapacity(this.sz+2); + let cap = this.sz + 2; + if (this.r.length < cap) { + let newsize = Math.max(2 * this.r.length, cap); + let rnew = new BigInt64Array(newsize); + rnew.set(this.r); + this.r = rnew; + } + this.r[this.sz] = a; + this.r[this.sz + 1] = b; + this.sz += 2; + } + ; + /** Make sure the object can hold at least the given number of entries. + * @param cap int + * */ + ensureCapacity(cap) { + if (this.r.length < cap) + this.resize(Math.max(2 * this.r.length, cap)); + } + ; + /** + * @param newsize int + */ + resize(newsize) { + if (newsize < this.sz) + console.error("requested array size too small"); + if (newsize == this.r.length) + return; + let rnew = new BigInt64Array(newsize); + rnew.set(this.r.slice(0, this.sz)); + this.r = rnew; + } + ; +} +//# sourceMappingURL=RangeSet.js.map \ No newline at end of file diff --git a/src/Vec3.js b/src/Vec3.js new file mode 100644 index 0000000..ba4dca6 --- /dev/null +++ b/src/Vec3.js @@ -0,0 +1,117 @@ +/** + * Partial porting to Javascript of Vec3.java from Healpix3.30 + */ +import { Hploc } from './Hploc.js'; +import { Pointing } from './Pointing.js'; +export class Vec3 { + constructor(in_x, in_y, in_z) { + if (in_x instanceof Pointing) { + let ptg = in_x; + let sth = Hploc.sin(ptg.theta); + this.x = sth * Hploc.cos(ptg.phi); + this.y = sth * Hploc.sin(ptg.phi); + this.z = Hploc.cos(ptg.theta); + } + else { + this.x = in_x; + this.y = in_y; + this.z = in_z; + } + } + getX() { + return this.x; + } + ; + getY() { + return this.y; + } + ; + getZ() { + return this.z; + } + ; + /** Scale the vector by a given factor + @param n the scale factor */ + scale(n) { + this.x *= n; + this.y *= n; + this.z *= n; + } + ; + /** Vector cross product. + @param v another vector + @return the vector cross product between this vector and {@code v} */ + cross(v) { + return new Vec3(this.y * v.z - v.y * this.z, this.z * v.x - v.z * this.x, this.x * v.y - v.x * this.y); + } + ; + /** Vector addition + * @param v the vector to be added + * @return addition result */ + add(v) { + return new Vec3(this.x + v.x, this.y + v.y, this.z + v.z); + } + ; + /** Normalize the vector */ + normalize() { + let d = 1. / this.length(); + this.x *= d; + this.y *= d; + this.z *= d; + } + ; + /** Return normalized vector */ + norm() { + let d = 1. / this.length(); + return new Vec3(this.x * d, this.y * d, this.z * d); + } + ; + /** Vector length + @return the length of the vector. */ + length() { + return Math.sqrt(this.lengthSquared()); + } + ; + /** Squared vector length + @return the squared length of the vector. */ + lengthSquared() { + return this.x * this.x + this.y * this.y + this.z * this.z; + } + ; + /** Computes the dot product of the this vector and {@code v1}. + * @param v1 another vector + * @return dot product */ + dot(v1) { + return this.x * v1.x + this.y * v1.y + this.z * v1.z; + } + ; + /** Vector subtraction + * @param v the vector to be subtracted + * @return subtraction result */ + sub(v) { + return new Vec3(this.x - v.x, this.y - v.y, this.z - v.z); + } + ; + /** Angle between two vectors. + @param v1 another vector + @return the angle in radians between this vector and {@code v1}; + constrained to the range [0,PI]. */ + angle(v1) { + return Hploc.atan2(this.cross(v1).length(), this.dot(v1)); + } + /** Invert the signs of all components */ + flip() { + this.x *= -1.0; + this.y *= -1.0; + this.z *= -1.0; + } + static pointing2Vec3(pointing) { + let sth = Hploc.sin(pointing.theta); + let x = sth * Hploc.cos(pointing.phi); + let y = sth * Hploc.sin(pointing.phi); + let z = Hploc.cos(pointing.theta); + return new Vec3(x, y, z); + } + ; +} +//# sourceMappingURL=Vec3.js.map \ No newline at end of file diff --git a/src/Xyf.js b/src/Xyf.js new file mode 100644 index 0000000..ec2f03a --- /dev/null +++ b/src/Xyf.js @@ -0,0 +1,11 @@ +/** + * Partial porting to Javascript of Xyf.java from Healpix3.30 + */ +export class Xyf { + constructor(x, y, f) { + this.ix = x; + this.iy = y; + this.face = f; + } +} +//# sourceMappingURL=Xyf.js.map \ No newline at end of file diff --git a/src/Zphi.js b/src/Zphi.js new file mode 100644 index 0000000..7137315 --- /dev/null +++ b/src/Zphi.js @@ -0,0 +1,9 @@ +export class Zphi { + /** Creation from individual components */ + constructor(z_, phi_) { + this.z = z_; + this.phi = phi_; + } + ; +} +//# sourceMappingURL=Zphi.js.map \ No newline at end of file diff --git a/src/index.js b/src/index.js new file mode 100644 index 0000000..adc6e57 --- /dev/null +++ b/src/index.js @@ -0,0 +1,12 @@ +export { Constants } from "./Constants.js"; +export { pstack } from "./pstack.js"; +export { CircleFinder } from './CircleFinder.js'; +export { Fxyf } from './Fxyf.js'; +export { Healpix } from './Healpix.js'; +export { Pointing } from './Pointing.js'; +export { RangeSet } from './RangeSet.js'; +export { Vec3 } from './Vec3.js'; +export { Xyf } from './Xyf.js'; +export { Zphi } from './Zphi.js'; +export { Hploc } from './Hploc.js'; +//# sourceMappingURL=index.js.map \ No newline at end of file diff --git a/src/pstack.js b/src/pstack.js new file mode 100644 index 0000000..55b5c2b --- /dev/null +++ b/src/pstack.js @@ -0,0 +1,45 @@ +export class pstack { + /** Creation from individual components */ + constructor(sz) { + this.p = new BigInt64Array(sz); + this.o = new Int32Array(sz); + this.s = 0; + this.m = 0; + } + ; + /** + * @param p long + * @param o int + */ + push(p_, o_) { + this.p[this.s] = p_; + this.o[this.s] = o_; + ++this.s; + } + ; + pop() { + --this.s; + } + ; + popToMark() { + this.s = this.m; + } + ; + size() { + return this.s; + } + ; + mark() { + this.m = this.s; + } + ; + otop() { + return this.o[this.s - 1]; + } + ; + ptop() { + return this.p[this.s - 1]; + } + ; +} +//# sourceMappingURL=pstack.js.map \ No newline at end of file