[Yt-svn] yt: Adding in ang2pix and pix2ang; this enables us to construct ...
hg at spacepope.org
hg at spacepope.org
Thu Feb 10 07:10:11 PST 2011
hg Repository: yt
details: yt/rev/a763d0bcb2df
changeset: 3729:a763d0bcb2df
user: Matthew Turk <matthewturk at gmail.com>
date:
Thu Feb 10 10:09:48 2011 -0500
description:
Adding in ang2pix and pix2ang; this enables us to construct our own Lat/Lon
grids, which means we can plot the results of healpix-style integration using
matplotlib's projection transformations.
diffstat:
yt/utilities/_amr_utils/VolumeIntegrator.pyx | 39 ++++++-
yt/utilities/_amr_utils/healpix_ang2pix_nest.c | 137 +++++++++++++++++++++
yt/utilities/_amr_utils/healpix_interface.pxd | 2 +
yt/utilities/_amr_utils/healpix_pix2ang_nest.c | 157 +++++++++++++++++++++++++
4 files changed, 334 insertions(+), 1 deletions(-)
diffs (truncated from 367 to 300 lines):
diff -r 489e024f45fc -r a763d0bcb2df yt/utilities/_amr_utils/VolumeIntegrator.pyx
--- a/yt/utilities/_amr_utils/VolumeIntegrator.pyx Wed Feb 09 17:54:21 2011 -0500
+++ b/yt/utilities/_amr_utils/VolumeIntegrator.pyx Thu Feb 10 10:09:48 2011 -0500
@@ -106,7 +106,6 @@
healpix_interface.vec2pix_nest(nside, v, &ipix)
return ipix
-
def arr_vec2pix_nest(long nside,
np.ndarray[np.float64_t, ndim=1] x,
np.ndarray[np.float64_t, ndim=1] y,
@@ -124,6 +123,44 @@
tr[i] = ipix
return tr
+def hp_pix2ang_nest(long nside, long ipnest):
+ cdef double theta, phi
+ healpix_interface.pix2ang_nest(nside, ipnest, &theta, &phi)
+ return (theta, phi)
+
+def arr_pix2ang_nest(long nside, np.ndarray[np.int64_t, ndim=1] aipnest):
+ cdef int n = aipnest.shape[0]
+ cdef int i
+ cdef long ipnest
+ cdef np.ndarray[np.float64_t, ndim=2] tr = np.zeros((n, 2), dtype='float64')
+ cdef double theta, phi
+ for i in range(n):
+ ipnest = aipnest[i]
+ healpix_interface.pix2ang_nest(nside, ipnest, &theta, &phi)
+ tr[i,0] = theta
+ tr[i,1] = phi
+ return tr
+
+def hp_ang2pix_nest(long nside, double theta, double phi):
+ cdef long ipix
+ healpix_interface.ang2pix_nest(nside, theta, phi, &ipix)
+ return ipix
+
+def arr_ang2pix_nest(long nside,
+ np.ndarray[np.float64_t, ndim=1] atheta,
+ np.ndarray[np.float64_t, ndim=1] aphi):
+ cdef int n = atheta.shape[0]
+ cdef int i
+ cdef long ipnest
+ cdef np.ndarray[np.int64_t, ndim=1] tr = np.zeros(n, dtype='int64')
+ cdef double theta, phi
+ for i in range(n):
+ theta = atheta[i]
+ phi = aphi[i]
+ healpix_interface.ang2pix_nest(nside, theta, phi, &ipnest)
+ tr[i] = ipnest
+ return tr
+
cdef class star_kdtree_container:
cdef kdtree_utils.kdtree *tree
cdef public np.float64_t sigma
diff -r 489e024f45fc -r a763d0bcb2df yt/utilities/_amr_utils/healpix_ang2pix_nest.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/_amr_utils/healpix_ang2pix_nest.c Thu Feb 10 10:09:48 2011 -0500
@@ -0,0 +1,137 @@
+/* -----------------------------------------------------------------------------
+ *
+ * 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_nest.c */
+
+/* Standard Includes */
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/* Local Includes */
+#include "healpix_vectors.h"
+
+void ang2pix_nest( const long nside, double theta, double phi, long *ipix) {
+
+ /* =======================================================================
+ * subroutine ang2pix_nest(nside, theta, phi, ipix)
+ * =======================================================================
+ * gives the pixel number ipix (NESTED) corresponding to angles theta and phi
+ *
+ * the computation is made to the highest resolution available (nside=8192)
+ * and then degraded to that required (by integer division)
+ * this doesn't cost more, and it makes sure that the treatement of round-off
+ * will be consistent for every resolution
+ * =======================================================================
+ */
+
+ double z, za, z0, tt, tp, tmp;
+ int face_num,jp,jm;
+ long ifp, ifm;
+ int ix, iy, ix_low, ix_hi, iy_low, iy_hi, ipf, ntt;
+ double piover2 = 0.5*M_PI, pi = M_PI, twopi = 2.0*M_PI;
+ int ns_max = 8192;
+ static int x2pix[128], y2pix[128];
+ static char setup_done = 0;
+
+ if( nside<1 || nside>ns_max ) {
+ fprintf(stderr, "%s (%d): nside out of range: %ld\n", __FILE__, __LINE__, nside);
+ exit(0);
+ }
+ if( theta<0. || theta>pi ) {
+ fprintf(stderr, "%s (%d): theta out of range: %f\n", __FILE__, __LINE__, theta);
+ exit(0);
+ }
+ if( !setup_done ) {
+ mk_xy2pix(x2pix,y2pix);
+ setup_done = 1;
+ }
+
+ z = cos(theta);
+ za = fabs(z);
+ z0 = 2./3.;
+ if( phi>=twopi ) phi = phi - twopi;
+ if( phi<0. ) phi = phi + twopi;
+ tt = phi / piover2; /* in [0,4[ */
+
+ if( za<=z0 ) { /* equatorial region */
+
+ /* (the index of edge lines increase when the longitude=phi goes up) */
+ jp = (int)floor(ns_max*(0.5 + tt - z*0.75)); /* ascending edge line index */
+ jm = (int)floor(ns_max*(0.5 + tt + z*0.75)); /* descending edge line index */
+
+ /* finds the face */
+ ifp = jp / ns_max; /* in {0,4} */
+ ifm = jm / ns_max;
+
+ if( ifp==ifm ) face_num = (int)fmod(ifp,4) + 4; /* faces 4 to 7 */
+ else if( ifp<ifm ) face_num = (int)fmod(ifp,4); /* (half-)faces 0 to 3 */
+ else face_num = (int)fmod(ifm,4) + 8; /* (half-)faces 8 to 11 */
+
+ ix = (int)fmod(jm, ns_max);
+ iy = ns_max - (int)fmod(jp, ns_max) - 1;
+ }
+ else { /* polar region, za > 2/3 */
+
+ ntt = (int)floor(tt);
+ if( ntt>=4 ) ntt = 3;
+ tp = tt - ntt;
+ tmp = sqrt( 3.*(1. - za) ); /* in ]0,1] */
+
+ /* (the index of edge lines increase when distance from the closest pole
+ * goes up)
+ */
+ /* line going toward the pole as phi increases */
+ jp = (int)floor( ns_max * tp * tmp );
+
+ /* that one goes away of the closest pole */
+ jm = (int)floor( ns_max * (1. - tp) * tmp );
+ jp = (int)(jp < ns_max-1 ? jp : ns_max-1);
+ jm = (int)(jm < ns_max-1 ? jm : ns_max-1);
+
+ /* finds the face and pixel's (x,y) */
+ if( z>=0 ) {
+ face_num = ntt; /* in {0,3} */
+ ix = ns_max - jm - 1;
+ iy = ns_max - jp - 1;
+ }
+ else {
+ face_num = ntt + 8; /* in {8,11} */
+ ix = jp;
+ iy = jm;
+ }
+ }
+
+ ix_low = (int)fmod(ix,128);
+ ix_hi = ix/128;
+ iy_low = (int)fmod(iy,128);
+ iy_hi = iy/128;
+
+ ipf = (x2pix[ix_hi]+y2pix[iy_hi]) * (128 * 128)+ (x2pix[ix_low]+y2pix[iy_low]);
+ ipf = (long)(ipf / pow(ns_max/nside,2)); /* in {0, nside**2 - 1} */
+ *ipix =(long)( ipf + face_num*pow(nside,2)); /* in {0, 12*nside**2 - 1} */
+}
diff -r 489e024f45fc -r a763d0bcb2df yt/utilities/_amr_utils/healpix_interface.pxd
--- a/yt/utilities/_amr_utils/healpix_interface.pxd Wed Feb 09 17:54:21 2011 -0500
+++ b/yt/utilities/_amr_utils/healpix_interface.pxd Thu Feb 10 10:09:48 2011 -0500
@@ -32,3 +32,5 @@
cdef extern from "healpix_vectors.h":
int pix2vec_nest(long nside, long ipix, double *v)
void vec2pix_nest(long nside, double *vec, long *ipix)
+ void pix2ang_nest(long nside, long ipix, double *theta, double *phi)
+ void ang2pix_nest(long nside, double theta, double phi, long *ipix)
diff -r 489e024f45fc -r a763d0bcb2df yt/utilities/_amr_utils/healpix_pix2ang_nest.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/utilities/_amr_utils/healpix_pix2ang_nest.c Thu Feb 10 10:09:48 2011 -0500
@@ -0,0 +1,157 @@
+/* -----------------------------------------------------------------------------
+ *
+ * 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
+ *
+ *----------------------------------------------------------------------------- */
+/* Standard Includes */
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/* Local Includes */
+#include "healpix_vectors.h"
+
+void pix2ang_nest( long nside, long ipix, double *theta, double *phi) {
+
+ /*
+ c=======================================================================
+ subroutine pix2ang_nest(nside, ipix, theta, phi)
+ c=======================================================================
+ c gives theta and phi corresponding to pixel ipix (NESTED)
+ c for a parameter nside
+ c=======================================================================
+ */
+
+ int npix, 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;
+ double piover2=0.5*M_PI;
+ int ns_max=8192;
+
+ static int pix2x[1024], pix2y[1024];
+ // common /pix2xy/ pix2x, pix2y
+
+ int jrll[12], jpll[12];// ! coordinate of the lowest corner of each face
+ // data jrll/2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4/ ! in unit of nside
+ // data jpll/1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7/ ! in unit of nside/2
+ jrll[0]=2;
+ jrll[1]=2;
+ jrll[2]=2;
+ jrll[3]=2;
+ jrll[4]=3;
+ jrll[5]=3;
+ jrll[6]=3;
+ jrll[7]=3;
+ jrll[8]=4;
+ jrll[9]=4;
+ jrll[10]=4;
+ jrll[11]=4;
+ jpll[0]=1;
+ jpll[1]=3;
+ jpll[2]=5;
+ jpll[3]=7;
+ jpll[4]=0;
+ jpll[5]=2;
+ jpll[6]=4;
+ jpll[7]=6;
+ jpll[8]=1;
+ jpll[9]=3;
+ jpll[10]=5;
+ jpll[11]=7;
+
+
+ if( nside<1 || nside>ns_max ) {
+ fprintf(stderr, "%s (%d): nside out of range: %ld\n", __FILE__, __LINE__, nside);
+ exit(0);
+ }
+ npix = 12 * nside*nside;
More information about the yt-svn
mailing list