Initial commit

1 files changed, 492 insertions, 0 deletions

diff --git a/pkg/images/imcoords/src/healpix.x b/pkg/images/imcoords/src/healpix.x
new file mode 100644
index 00000000..1156607c
--- /dev/null
+++ b/pkg/images/imcoords/src/healpix.x
@@ -0,0 +1,492 @@
+include	<math.h>
+
+define	MTYPES	"|nest|ring|"
+define	NEST	1
+define	RING	2
+
+define	NS_MAX		8192
+define	TWOTHIRDS	0.66666666667
+
+
+# ANG2PIX -- Compute the HEALPix map row from a spherical coordinate.
+#
+# It is up to the caller to know the coordinate type, map type, and
+# resolution for the map.
+#
+# The returned row is 1 indexed.
+
+procedure ang2row (row, lng, lat, mtype, nside)
+
+int	row			#O Table row
+double	lng			#I Longitude (deg)
+double	lat			#I Latitude (deg)
+int	mtype			#I HEALPix map type
+int	nside			#I Resolution parameter
+
+int	ipix
+double	phi, theta
+errchk	ang2pix_nest, ang2pix_ring
+
+begin
+	# Check parameters and call appropriate procedure.
+
+	if (nside < 1 || nside > NS_MAX)
+	    call error (1, "nside out of range")
+
+	if (lat < -90D0 || lat > 90D0)
+	    call error (2, "latitude out of range")
+
+	phi = DEGTORAD (lng)
+	theta = DEGTORAD (90D0 - lat)
+
+	switch (mtype) {
+	case NEST:
+	    call ang2pix_nest (nside, theta, phi, ipix)
+	case RING:
+	    call ang2pix_ring (nside, theta, phi, ipix)
+	default:
+	    call error (3, "unknown HEALPix map type")
+	}
+
+	row = ipix + 1
+end
+
+
+# PIX2ANG -- Compute spherical coordinate from HEALPix map row.
+#
+# It is up to the caller to know the coordinate type, map type, and
+# resolution for the map.
+
+procedure row2ang (row, lng, lat, mtype, nside)
+
+int	row			#I Table row (1 indexed)
+double	lng			#O Longitude (deg)
+double	lat			#O Latitude (deg)
+int	mtype			#I HEALPix map type
+int	nside			#I Resolution parameter
+
+int	ipix
+double	phi, theta
+errchk	pix2ang_nest, pix2ang_ring
+
+begin
+	# Check input parameters and call appropriate procedure.
+
+	if (nside < 1 || nside > NS_MAX)
+	    call error (1, "nside out of range")
+
+	if (row < 1 || row > 12*nside*nside)
+	    call error (1, "row out of range")
+
+	ipix = row - 1
+
+	switch (mtype) {
+	case NEST:
+	    call pix2ang_nest (nside, ipix, theta, phi)
+	case RING:
+	    call pix2ang_ring (nside, ipix, theta, phi)
+	default:
+	    call error (3, "unknown HEALPix map type")
+	}
+
+	lng = RADTODEG (phi)
+	lat = 90D0 - RADTODEG (theta)
+end
+
+
+# The following routines are SPP translations of the HEALPix software from
+# the authors identified below.  If it matters, the C version was used
+# though the translation is not necessarily exact.  Comments were
+# largely removed.
+#
+# I'm not sure if the arguments to the floor function in the original
+# can be negative.  Assuming not I just do an integer truncation.
+
+# -----------------------------------------------------------------------------
+#
+#  Copyright (C) 1997-2008 Krzysztof M. Gorski, Eric Hivon,
+#                          Benjamin D. Wandelt, Anthony J. Banday,
+#                          Matthias Bartelmann,
+#                          Reza Ansari & Kenneth M. Ganga
+#
+#
+#  This file is part of HEALPix.
+#
+#  HEALPix is free software; you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published
+#  by
+#  the Free Software Foundation; either version 2 of the License, or
+#  (at your option) any later version.
+#
+#  HEALPix is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with HEALPix; if not, write to the Free Software
+#  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
+#  02110-1301  USA
+#
+#  For more information about HEALPix see http://healpix.jpl.nasa.gov
+#
+#----------------------------------------------------------------------------- 
+
+	
+# ANG2PIX_NEST -- Compute HEALPix index for a nested map.
+
+procedure ang2pix_nest (nside, theta, phi, ipix)
+
+int	nside			#I Resolution parameter
+double	theta			#I Latitude (rad from pole)
+double	phi			#I Longitude (rad)
+int	ipix			#O HEALPix index
+
+double	z, za, tt, tp, tmp
+int	face_num, jp, jm
+long	ifp, ifm
+int	ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt
+int	 x2pix[128], y2pix[128]
+int	 setup_done
+
+errchk	mk_xy2pix
+
+data	setup_done/NO/
+
+begin
+	if (setup_done == NO) {
+	    call mk_xy2pix (x2pix, y2pix)
+	    setup_done = YES
+	}
+
+	z  = cos (theta)
+	za = abs (z)
+	if (phi >= TWOPI)
+	    phi = phi -  TWOPI
+	if (phi < 0.)
+	    phi = phi +  TWOPI
+	tt = phi / HALFPI
+	
+	if (za <= TWOTHIRDS) {
+	    jp = int (NS_MAX * (0.5 + tt - z * 0.75))
+	    jm = int (NS_MAX * (0.5 + tt + z * 0.75))
+	    
+	    ifp = jp / NS_MAX
+	    ifm = jm / NS_MAX
+	    
+	    if (ifp==ifm)
+	        face_num = mod (ifp, 4) + 4
+	    else if (ifp<ifm)
+	        face_num = mod (ifp, 4)
+	    else
+	        face_num = mod (ifm, 4) + 8
+	    
+	    ix = mod (jm, NS_MAX)
+	    iy = NS_MAX - mod (jp, NS_MAX) - 1
+	} else {
+	    ntt = int (tt)
+	    if (ntt >= 4)
+	        ntt = 3
+	    tp = tt - ntt
+	    tmp = sqrt (3. * (1. - za))
+	    
+	    jp = int (NS_MAX * tp * tmp)
+	    jm = int (NS_MAX * (1. - tp) * tmp)
+	    jp = min (jp, NS_MAX-1)
+	    jm = min (jm, NS_MAX-1)
+	    
+	    if (z >= 0) {
+		face_num = ntt
+		ix = NS_MAX - jm - 1
+		iy = NS_MAX - jp - 1
+	    } else {
+		face_num = ntt + 8
+		ix =  jp
+		iy =  jm
+	    }
+	}
+	
+	ix_low = mod (ix, 128) + 1
+	ix_hi  = ix / 128 + 1
+	iy_low = mod (iy, 128) + 1
+	iy_hi  = iy / 128 + 1
+
+	ipf = (x2pix[ix_hi] + y2pix[iy_hi]) * (128 * 128) +
+	    (x2pix[ix_low] + y2pix[iy_low])
+	ipf = ipf / (NS_MAX/nside)**2
+	ipix = ipf + face_num * nside**2
+end
+
+	
+# ANG2PIX_RING -- Compute HEALPix index for a ring map.
+
+procedure ang2pix_ring (nside, theta, phi, ipix)
+
+int	nside			#I Resolution parameter
+double	theta			#I Latitude (rad from pole)
+double	phi			#I Longitude (rad)
+int	ipix			#O HEALPix index
+
+int	nl2, nl4, ncap, npix, jp, jm, ipix1
+double	z, za, tt, tp, tmp
+int	ir, ip, kshift
+  
+begin
+	z = cos (theta)
+	za = abs (z)
+	if ( phi >= TWOPI)
+	    phi = phi - TWOPI
+	if (phi < 0.)
+	    phi = phi + TWOPI
+	tt = phi / HALFPI
+
+	nl2 = 2 * nside
+	nl4 = 4 * nside
+	ncap  = nl2 * (nside - 1)
+	npix  = 12 * nside * nside
+
+	if (za <= TWOTHIRDS) {
+
+	    jp = int (nside * (0.5 + tt - z * 0.75))
+	    jm = int (nside * (0.5 + tt + z * 0.75))
+
+	    ir = nside + 1 + jp - jm
+	    kshift = 0
+	    if (mod (ir,2) == 0)
+	        kshift = 1
+
+	    ip = int ((jp + jm - nside + kshift + 1) / 2) + 1
+	    if (ip > nl4)
+	        ip = ip - nl4
+
+	    ipix1 = ncap + nl4 *  (ir - 1) + ip 
+	} else {
+
+	    tp = tt - int (tt)
+	    tmp = sqrt (3. *  (1. - za))
+
+	    jp = int (nside * tp * tmp)
+	    jm = int (nside * (1. - tp) * tmp)
+
+	    ir = jp + jm + 1
+	    ip = int (tt * ir) + 1
+	    if (ip > 4*ir)
+	        ip = ip - 4 * ir
+
+	    ipix1 = 2 * ir * (ir - 1) + ip
+	    if (z<=0.) {
+	      ipix1 = npix - 2 * ir * (ir + 1) + ip
+	    }
+	}
+	ipix = ipix1 - 1
+end
+
+
+# PIX2ANG_NEST -- Translate HEALpix nested row to spherical coordinates.
+
+procedure pix2ang_nest (nside, ipix, theta, phi)
+
+int	nside			#I Resolution parameter
+int	ipix			#I HEALPix index
+double	theta			#O Latitude (rad from pole)
+double	phi			#O Longitude (rad)
+
+int	npface, face_num
+int	ipf, ip_low, ip_trunc, ip_med, ip_hi
+int	ix, iy, jrt, jr, nr, jpt, jp, kshift, nl4
+double	z, fn, fact1, fact2
+
+int	pix2x[1024], pix2y[1024]
+
+int	jrll[12], jpll[12], setup_done
+data	jrll/2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4/
+data	jpll/1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7/
+data	setup_done/NO/
+
+begin
+	if (setup_done == NO) {
+	    call mk_pix2xy (pix2x,pix2y)
+	    setup_done = YES
+	}
+
+	fn = 1. * nside
+	fact1 = 1. / (3. * fn * fn)
+	fact2 = 2. / (3. * fn)
+	nl4   = 4 * nside
+
+	npface = nside * nside
+
+	face_num = ipix / npface + 1
+	ipf = mod (ipix, npface)
+
+	ip_low = mod (ipf, 1024) + 1
+	ip_trunc = ipf / 1024 
+	ip_med = mod (ip_trunc, 1024) + 1
+	ip_hi  = ip_trunc / 1024 + 1
+
+	ix = 1024*pix2x[ip_hi] + 32*pix2x[ip_med] + pix2x[ip_low]
+	iy = 1024*pix2y[ip_hi] + 32*pix2y[ip_med] + pix2y[ip_low]
+
+	jrt = ix + iy
+	jpt = ix - iy
+
+	jr =  jrll[face_num] * nside - jrt - 1
+	nr = nside
+	z  = (2 * nside - jr) * fact2
+	kshift = mod (jr - nside, 2)
+	if( jr < nside) {
+	    nr = jr
+	    z = 1. - nr * nr * fact1
+	    kshift = 0
+	} else if (jr > 3*nside) {
+	    nr = nl4 - jr
+	    z = - 1. + nr * nr * fact1
+	    kshift = 0
+	}
+
+	jp = (jpll[face_num] * nr + jpt + 1 + kshift)/2
+	if (jp > nl4)
+	    jp = jp - nl4
+	if (jp < 1)
+	    jp = jp + nl4
+
+	theta = acos(z)
+	phi = (jp - (kshift+1)*0.5) * (HALFPI / nr)
+end
+
+
+# PIX2ANG_RING -- Convert HEALpix pixel to spherical coordinates.
+
+procedure pix2ang_ring (nside, ipix, theta, phi)
+
+int	nside			#I Resolution parameter
+int	ipix			#I HEALPix index
+double	theta			#O Latitude (rad from pole)
+double	phi			#O Longitude (rad)
+
+int	nl2, nl4, npix, ncap, iring, iphi, ip, ipix1
+double	fact1, fact2, fodd, hip, fihip
+
+begin
+	npix = 12 * nside * nside
+	ipix1 = ipix + 1
+	nl2 = 2 * nside
+	nl4 = 4 * nside
+	ncap = 2 * nside * (nside - 1)
+	fact1 = 1.5 * nside
+	fact2 = 3.0 * nside * nside
+
+	if (ipix1 <= ncap) {
+
+	    hip = ipix1 / 2.
+	    fihip = int (hip)
+	    iring = int (sqrt (hip - sqrt (fihip))) + 1
+	    iphi = ipix1 - 2 * iring * (iring - 1)
+
+	    theta = acos (1. - iring * iring / fact2)
+	    phi = (iphi - 0.5) * PI / (2. * iring)
+
+	} else if (ipix1 <= nl2 * (5 * nside + 1)) {
+
+	    ip    = ipix1 - ncap - 1
+	    iring = (ip / nl4) + nside
+	    iphi  = mod (ip, nl4) + 1
+
+	    fodd  = 0.5 * (1 + mod (iring + nside, 2))
+	    theta = acos ((nl2 - iring) / fact1)
+	    phi   = (iphi - fodd) * PI / (2. * nside)
+
+	} else {
+
+	    ip    = npix - ipix1 + 1
+	    hip   = ip/2.
+	    
+	    fihip = int (hip)
+	    iring = int (sqrt (hip - sqrt (fihip))) + 1
+	    iphi  = 4. * iring + 1 - (ip - 2. * iring * (iring-1))
+
+	    theta = acos (-1. + iring * iring / fact2)
+	    phi   = (iphi - 0.5) * PI / (2. * iring)
+
+	}
+end
+
+
+# MK_XY2PIX
+#
+# Sets the array giving the number of the pixel lying in (x,y)
+# x and y are in {1,128}
+# the pixel number is in {0,128**2-1}
+#
+# if  i-1 = sum_p=0  b_p * 2^p
+# then ix = sum_p=0  b_p * 4^p
+# iy = 2*ix
+# ix + iy in {0, 128**2 -1}
+
+procedure mk_xy2pix (x2pix, y2pix)
+
+int	x2pix[128], y2pix[128]
+
+int	i, j, k, ip, id
+
+begin
+	do i = 1, 128
+	    x2pix[i] = 0
+
+	do i = 1, 128 {
+	    j = i - 1
+	    k = 0
+	    ip = 1
+	    while (j != 0) {
+		id = mod (j, 2)
+		j  = j / 2
+		k  = ip * id + k
+		ip = ip * 4
+	    }
+	    x2pix[i] = k
+	    y2pix[i] = 2 * k
+	}
+end
+
+
+# MK_PIX2XY
+#
+# Constructs the array giving x and y in the face from pixel number
+# for the nested (quad-cube like) ordering of pixels.
+# 
+# The bits corresponding to x and y are interleaved in the pixel number.
+# One breaks up the pixel number by even and odd bits.
+
+procedure mk_pix2xy (pix2x, pix2y)
+
+int	pix2x[1024], pix2y[1024]
+
+int	kpix, jpix, ix, iy, ip, id
+
+begin
+
+	do kpix = 1, 1024
+	    pix2x[kpix] = 0
+
+	do kpix = 1, 1024 {
+	    jpix = kpix - 1
+	    ix = 0
+	    iy = 0
+	    ip = 1 
+	    while (jpix != 0) {
+		id = mod (jpix, 2)
+		jpix = jpix / 2
+		ix = id * ip + ix
+		  
+		id = mod (jpix, 2)
+		jpix = jpix / 2
+		iy = id * ip + iy
+		  
+		ip = 2 * ip
+	    }
+
+	    pix2x[kpix] = ix
+	    pix2y[kpix] = iy
+	}
+
+end