gkucmierz/healpixjs-bigint
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Initial commit: BigInt port of healpixjs

Browse Source
This commit is contained in:
Grzegorz Kućmierz
2026-03-06 01:39:02 +00:00
commit 4f73d2e800
18 changed files with 1531 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
.gitea/workflows/npm-publish.yaml Normal file
View File

@@ -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 }}

4
.gitignore vendored Normal file
View File

@@ -0,0 +1,4 @@
node_modules/
.npmrc
*.log
.DS_Store

45
README.md Normal file
View File

@@ -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.

8
index.js Normal file
View File

@@ -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';

13
package-lock.json generated Normal file
View File

@@ -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"
}
}
}

33
package.json Normal file
View File

@@ -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"
}

68
src/CircleFinder.js Normal file
View File

@@ -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

9
src/Constants.js Normal file
View File

@@ -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

57
src/Fxyf.js Normal file
View File

@@ -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

781
src/Healpix.js Normal file
View File

@@ -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<<order_max;
this.ns_max = Math.pow(2, this.order_max);
this.ctab = new Uint16Array([
0, 1, 256, 257, 2, 3, 258, 259, 512, 513, 768, 769, 514, 515, 770, 771, 4, 5, 260, 261, 6, 7, 262,
263, 516, 517, 772, 773, 518, 519, 774, 775, 1024, 1025, 1280, 1281, 1026, 1027, 1282, 1283,
1536, 1537, 1792, 1793, 1538, 1539, 1794, 1795, 1028, 1029, 1284, 1285, 1030, 1031, 1286,
1287, 1540, 1541, 1796, 1797, 1542, 1543, 1798, 1799, 8, 9, 264, 265, 10, 11, 266, 267, 520,
521, 776, 777, 522, 523, 778, 779, 12, 13, 268, 269, 14, 15, 270, 271, 524, 525, 780, 781, 526,
527, 782, 783, 1032, 1033, 1288, 1289, 1034, 1035, 1290, 1291, 1544, 1545, 1800, 1801, 1546,
1547, 1802, 1803, 1036, 1037, 1292, 1293, 1038, 1039, 1294, 1295, 1548, 1549, 1804, 1805,
1550, 1551, 1806, 1807, 2048, 2049, 2304, 2305, 2050, 2051, 2306, 2307, 2560, 2561, 2816,
2817, 2562, 2563, 2818, 2819, 2052, 2053, 2308, 2309, 2054, 2055, 2310, 2311, 2564, 2565,
2820, 2821, 2566, 2567, 2822, 2823, 3072, 3073, 3328, 3329, 3074, 3075, 3330, 3331, 3584,
3585, 3840, 3841, 3586, 3587, 3842, 3843, 3076, 3077, 3332, 3333, 3078, 3079, 3334, 3335,
3588, 3589, 3844, 3845, 3590, 3591, 3846, 3847, 2056, 2057, 2312, 2313, 2058, 2059, 2314,
2315, 2568, 2569, 2824, 2825, 2570, 2571, 2826, 2827, 2060, 2061, 2316, 2317, 2062, 2063,
2318, 2319, 2572, 2573, 2828, 2829, 2574, 2575, 2830, 2831, 3080, 3081, 3336, 3337, 3082,
3083, 3338, 3339, 3592, 3593, 3848, 3849, 3594, 3595, 3850, 3851, 3084, 3085, 3340, 3341,
3086, 3087, 3342, 3343, 3596, 3597, 3852, 3853, 3598, 3599, 3854, 3855
]);
this.utab = new Uint16Array([0, 1, 4, 5, 16, 17, 20, 21, 64, 65, 68, 69, 80, 81, 84, 85, 256, 257, 260, 261, 272, 273, 276, 277,
320, 321, 324, 325, 336, 337, 340, 341, 1024, 1025, 1028, 1029, 1040, 1041, 1044, 1045, 1088,
1089, 1092, 1093, 1104, 1105, 1108, 1109, 1280, 1281, 1284, 1285, 1296, 1297, 1300, 1301,
1344, 1345, 1348, 1349, 1360, 1361, 1364, 1365, 4096, 4097, 4100, 4101, 4112, 4113, 4116,
4117, 4160, 4161, 4164, 4165, 4176, 4177, 4180, 4181, 4352, 4353, 4356, 4357, 4368, 4369,
4372, 4373, 4416, 4417, 4420, 4421, 4432, 4433, 4436, 4437, 5120, 5121, 5124, 5125, 5136,
5137, 5140, 5141, 5184, 5185, 5188, 5189, 5200, 5201, 5204, 5205, 5376, 5377, 5380, 5381,
5392, 5393, 5396, 5397, 5440, 5441, 5444, 5445, 5456, 5457, 5460, 5461, 16384, 16385, 16388,
16389, 16400, 16401, 16404, 16405, 16448, 16449, 16452, 16453, 16464, 16465, 16468, 16469,
16640, 16641, 16644, 16645, 16656, 16657, 16660, 16661, 16704, 16705, 16708, 16709, 16720,
16721, 16724, 16725, 17408, 17409, 17412, 17413, 17424, 17425, 17428, 17429, 17472, 17473,
17476, 17477, 17488, 17489, 17492, 17493, 17664, 17665, 17668, 17669, 17680, 17681, 17684,
17685, 17728, 17729, 17732, 17733, 17744, 17745, 17748, 17749, 20480, 20481, 20484, 20485,
20496, 20497, 20500, 20501, 20544, 20545, 20548, 20549, 20560, 20561, 20564, 20565, 20736,
20737, 20740, 20741, 20752, 20753, 20756, 20757, 20800, 20801, 20804, 20805, 20816, 20817,
20820, 20821, 21504, 21505, 21508, 21509, 21520, 21521, 21524, 21525, 21568, 21569, 21572,
21573, 21584, 21585, 21588, 21589, 21760, 21761, 21764, 21765, 21776, 21777, 21780, 21781,
21824, 21825, 21828, 21829, 21840, 21841, 21844, 21845]);
this.jrll = new Int16Array([2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4]);
this.jpll = new Int16Array([1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7]);
this.xoffset = new Int16Array([-1, -1, 0, 1, 1, 1, 0, -1]);
this.yoffset = new Int16Array([0, 1, 1, 1, 0, -1, -1, -1]);
this.facearray = [
new Int16Array([8, 9, 10, 11, -1, -1, -1, -1, 10, 11, 8, 9]),
new Int16Array([5, 6, 7, 4, 8, 9, 10, 11, 9, 10, 11, 8]),
new Int16Array([-1, -1, -1, -1, 5, 6, 7, 4, -1, -1, -1, -1]),
new Int16Array([4, 5, 6, 7, 11, 8, 9, 10, 11, 8, 9, 10]),
new Int16Array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]),
new Int16Array([1, 2, 3, 0, 0, 1, 2, 3, 5, 6, 7, 4]),
new Int16Array([-1, -1, -1, -1, 7, 4, 5, 6, -1, -1, -1, -1]),
new Int16Array([3, 0, 1, 2, 3, 0, 1, 2, 4, 5, 6, 7]),
new Int16Array([2, 3, 0, 1, -1, -1, -1, -1, 0, 1, 2, 3]) // N
];
// questo forse deve essere un UInt8Array. Viene usato da neighbours
this.swaparray = [
new Int16Array([0, 0, 3]),
new Int16Array([0, 0, 6]),
new Int16Array([0, 0, 0]),
new Int16Array([0, 0, 5]),
new Int16Array([0, 0, 0]),
new Int16Array([5, 0, 0]),
new Int16Array([0, 0, 0]),
new Int16Array([6, 0, 0]),
new Int16Array([3, 0, 0]) // N
];
if (nside_in <= this.ns_max && nside_in > 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<<i);
// this.mpr[i]=bn[i].maxPixrad();
// this.cmpr[i]=Math.cos(mpr[i]);
// this.smpr[i]=Math.sin(mpr[i]);
// }
}
computeBn() {
for (let i = 0; i <= this.order_max; ++i) {
this.bn[i] = new Healpix(1 << i);
this.mpr[i] = this.bn[i].maxPixrad();
this.cmpr[i] = Hploc.cos(this.mpr[i]);
this.smpr[i] = Hploc.sin(this.mpr[i]);
}
}
getNPix() {
return this.npix;
}
;
getBoundaries(pix) {
let points = new Array();
let xyf = this.nest2xyf(pix);
let dc = 0.5 / this.nside;
let xc = (xyf.ix + 0.5) / this.nside;
let yc = (xyf.iy + 0.5) / this.nside;
points[0] = new Fxyf(xc + dc, yc + dc, xyf.face).toVec3();
points[1] = new Fxyf(xc - dc, yc + dc, xyf.face).toVec3();
points[2] = new Fxyf(xc - dc, yc - dc, xyf.face).toVec3();
points[3] = new Fxyf(xc + dc, yc - dc, xyf.face).toVec3();
return points;
}
;
/** Returns a set of points along the boundary of the given pixel.
* Step 1 gives 4 points on the corners. The first point corresponds
* to the northernmost corner, the subsequent points follow the pixel
* boundary through west, south and east corners.
*
* @param pix pixel index number
* @param step the number of returned points is 4*step
* @return {@link Vec3} for each point
*/
getBoundariesWithStep(pix, step) {
// var points = new Array();
let points = new Array();
let xyf = this.nest2xyf(pix);
let dc = 0.5 / this.nside;
let xc = (xyf.ix + 0.5) / this.nside;
let yc = (xyf.iy + 0.5) / this.nside;
let d = 1.0 / (this.nside * step);
for (let i = 0; i < step; i++) {
points[i] = new Fxyf(xc + dc - i * d, yc + dc, xyf.face).toVec3();
points[i + step] = new Fxyf(xc - dc, yc + dc - i * d, xyf.face).toVec3();
points[i + 2 * step] = new Fxyf(xc - dc + i * d, yc - dc, xyf.face).toVec3();
points[i + 3 * step] = new Fxyf(xc + dc, yc - dc + i * d, xyf.face).toVec3();
}
return points;
}
;
getPointsForXyfNoStep(x, y, face) {
// let nside = Math.pow(2, this.order);
let points = new Array();
let xyf = new Xyf(x, y, face);
let dc = 0.5 / this.nside;
let xc = (xyf.ix + 0.5) / this.nside;
let yc = (xyf.iy + 0.5) / this.nside;
points[0] = new Fxyf(xc + dc, yc + dc, xyf.face).toVec3();
points[1] = new Fxyf(xc - dc, yc + dc, xyf.face).toVec3();
points[2] = new Fxyf(xc - dc, yc - dc, xyf.face).toVec3();
points[3] = new Fxyf(xc + dc, yc - dc, xyf.face).toVec3();
return points;
}
getPointsForXyf(x, y, step, face) {
let nside = step * Math.pow(2, this.order);
let points = new Array();
let xyf = new Xyf(x, y, face);
let dc = 0.5 / nside;
let xc = (xyf.ix + 0.5) / nside;
let yc = (xyf.iy + 0.5) / nside;
points[0] = new Fxyf(xc + dc, yc + dc, xyf.face).toVec3();
points[1] = new Fxyf(xc - dc, yc + dc, xyf.face).toVec3();
points[2] = new Fxyf(xc - dc, yc - dc, xyf.face).toVec3();
points[3] = new Fxyf(xc + dc, yc - dc, xyf.face).toVec3();
return points;
}
/** Returns the neighboring pixels of ipix.
This method works in both RING and NEST schemes, but is
considerably faster in the NEST scheme.
@param ipix the requested pixel number.
@return array with indices of the neighboring pixels.
The returned array contains (in this order)
the pixel numbers of the SW, W, NW, N, NE, E, SE and S neighbor
of ipix. If a neighbor does not exist (this can only happen
for the W, N, E and S neighbors), its entry is set to -1. */
neighbours(ipix) {
let result = new Int32Array(8);
let xyf = this.nest2xyf(ipix);
let ix = xyf.ix;
let iy = xyf.iy;
let face_num = xyf.face;
var nsm1 = this.nside - 1;
if ((ix > 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.
* <p>
* This method is more efficient in the RING scheme.
* <p>
* 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. <p>
This method is more efficient in the RING scheme. <p>
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<<order_max)>=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

197
src/Hploc.js Normal file
View File

@@ -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

32
src/Pointing.js Normal file
View File

@@ -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

66
src/RangeSet.js Normal file
View File

@@ -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

117
src/Vec3.js Normal file
View File

@@ -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

11
src/Xyf.js Normal file
View File

@@ -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

9
src/Zphi.js Normal file
View File

@@ -0,0 +1,9 @@
export class Zphi {
/** Creation from individual components */
constructor(z_, phi_) {
this.z = z_;
this.phi = phi_;
}
;
}
//# sourceMappingURL=Zphi.js.map

12
src/index.js Normal file
View File

@@ -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

45
src/pstack.js Normal file
View File

@@ -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