"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