ang2pix_ring.c

/* -----------------------------------------------------------------------------
 *
 *  Copyright (C) 1997-2010 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_ring.c
 *
 */

/* Standard Includes */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "chealpix.h"

/*
 * gives the pixel number ipix (RING)
 * corresponding to angles theta and phi
 *
 */
void
ang2pix_ring(const long nside, double theta, double phi, long *ipix)
{

	int			nl2,
				nl4,
				ncap,
				npix,
				jp,
				jm,
				ipix1;
	double		z,
				za,
				tt,
				tp,
				tmp;
	int			ir,
				ip,
				kshift;

	double		piover2 = 0.5 * M_PI;
	double		twopi = 2.0 * M_PI;
	double		z0 = 2.0 / 3.0;

	z = cos(theta);
	za = fabs(z);
	if (phi >= twopi)
		phi = phi - twopi;
	if (phi < 0.)
		phi = phi + twopi;
	tt = phi / piover2;			/* ! in [0,4) */

	nl2 = 2 * nside;
	nl4 = 4 * nside;
	ncap = nl2 * (nside - 1);	/* ! number of pixels in the north polar cap */
	npix = 12 * nside * nside;

	if (za <= z0)
	{

		jp = (int) floor(nside * (0.5 + tt - z * 0.75));	/* index of ascending
															 * edge line */
		jm = (int) floor(nside * (0.5 + tt + z * 0.75));	/* index of descending
															 * edge line */

		ir = nside + 1 + jp - jm;	/* ! in {1,2n+1} (ring number counted from
									 * z=2/3) */
		kshift = 0;
		if (fmod(ir, 2) == 0.)
			kshift = 1;			/* ! kshift=1 if ir even, 0 otherwise */

		ip = (int) floor((jp + jm - nside + kshift + 1) / 2) + 1;	/* ! in {1,4n} */
		if (ip > nl4)
			ip = ip - nl4;

		ipix1 = ncap + nl4 * (ir - 1) + ip;
	}
	else
	{

		tp = tt - floor(tt);	/* !MOD(tt,1.d0) */
		tmp = sqrt(3. * (1. - za));

		jp = (int) floor(nside * tp * tmp); /* ! increasing edge line index */
		jm = (int) floor(nside * (1. - tp) * tmp);	/* ! decreasing edge line
													 * index */

		ir = jp + jm + 1;		/* ! ring number counted from the closest pole */
		ip = (int) floor(tt * ir) + 1;	/* ! in {1,4*ir} */
		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;			/* ! in {0, npix-1} */

}