add support for Huffenberger & Wandelt grid

libsharp/sharp_geomhelpers.c

@@ -25,7 +25,7 @@
 /*! \file sharp_geomhelpers.c
  *  Spherical transform library
  *
- *  Copyright (C) 2006-2011 Max-Planck-Society
+ *  Copyright (C) 2006-2012 Max-Planck-Society
  *  Copyright (C) 2007-2008 Pavel Holoborodko (for gauss_legendre_tbl)
  *  \author Martin Reinecke \author Pavel Holoborodko
  */
@@ -241,3 +241,55 @@ void sharp_make_ecp_geom_info (int nrings, int nphi, double phi0,
   DEALLOC(ofs);
   DEALLOC(stride_);
   }
+
+static double eps(int j, int J)
+  {
+  if ((j==0)||(j==J)) return 0.5;
+  if ((j>0)&&(j<J)) return 1.;
+  return 0.;
+  }
+
+/* Weights from Keiner & Potts: "Fast evaluation of quadrature formulae
+   on the sphere", 2000 */
+/* nrings MUST be odd */
+void sharp_make_hw_geom_info (int nrings, int ppring, double phi0,
+  int stride_lon, int stride_lat, sharp_geom_info **geom_info)
+  {
+  const double pi=3.141592653589793238462643383279502884197;
+
+  double *theta=RALLOC(double,nrings);
+  double *weight=RALLOC(double,nrings);
+  int *nph=RALLOC(int,nrings);
+  double *phi0_=RALLOC(double,nrings);
+  ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
+  int *stride_=RALLOC(int,nrings);
+
+  UTIL_ASSERT((nrings&1)==1,"nrings must be an odd number");
+  int lmax=(nrings-1)/2;
+
+  for (int m=0; m<nrings; ++m)
+    {
+    theta[m]=pi*m/(nrings-1.);
+    if (theta[m]<1e-15) theta[m]=1e-15;
+    if (theta[m]>pi-1e-15) theta[m]=pi-1e-15;
+    nph[m]=ppring;
+    phi0_[m]=phi0;
+    ofs[m]=(ptrdiff_t)m*stride_lat;
+    stride_[m]=stride_lon;
+    double wgt=0;
+    double prefac=4*pi*eps(m,2*lmax)/lmax;
+    for (int l=0; l<=lmax; ++l)
+      wgt+=eps(l,lmax)/(1-4*l*l)*cos((pi*m*l)/lmax);
+    weight[m]=prefac*wgt/nph[m];
+    }
+
+  sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL, weight,
+    geom_info);
+
+  DEALLOC(theta);
+  DEALLOC(weight);
+  DEALLOC(nph);
+  DEALLOC(phi0_);
+  DEALLOC(ofs);
+  DEALLOC(stride_);
+  }

libsharp/sharp_geomhelpers.h

@@ -75,6 +75,22 @@ void sharp_make_gauss_geom_info (int nrings, int nphi, int stride_lon,
 void sharp_make_ecp_geom_info (int nrings, int nphi, double phi0,
   int stride_lon, int stride_lat, sharp_geom_info **geom_info);
 
+/*! Creates a geometry information describing an ECP map with \a nrings
+    iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
+    pixel in each ring is \a phi0 (in radians). The index difference between
+    two adjacent pixels in an iso-latitude ring is \a stride_lon, the index
+    difference between the two start pixels in consecutive iso-latitude rings
+    is \a stride_lat.
+    \note The spacing of pixel centers is equidistant in colatitude and
+      longitude.
+    \note \a nrings must be an odd number.
+    \note The sphere is pixelized in a way that the colatitude of the first ring
+      is \a 0 and that of the last ring is \a pi.
+    \note This is the grid used by Huffenberger & Wandelt 2010
+    \ingroup geominfogroup */
+void sharp_make_hw_geom_info (int nrings, int ppring, double phi0,
+  int stride_lon, int stride_lat, sharp_geom_info **geom_info);
+
 #ifdef __cplusplus
 }
 #endif

libsharp/sharp_test.c

@@ -193,7 +193,7 @@ int main(int argc, char **argv)
   MPI_Init(NULL,NULL);
 #endif
   sharp_module_startup("sharp_test",argc,7,
-    "<healpix|ecp|gauss> <lmax> <nside|nphi> <niter> <spin> <ntrans>",1);
+    "<healpix|ecp|gauss|hw> <lmax> <nside|nphi> <niter> <spin> <ntrans>",1);
 
   int lmax=atoi(argv[2]);
   int niter=atoi(argv[4]);
@@ -237,6 +237,17 @@ int main(int argc, char **argv)
     check_accuracy(tinfo,lmax,lmax,npix,spin,ntrans,niter);
     sharp_destroy_geom_info(tinfo);
     }
+  else if (strcmp(argv[1],"hw")==0)
+    {
+    int nrings=2*lmax+1;
+    int ppring=atoi(argv[3]);
+    ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
+    printf("\nTesting HW grid (%d rings, %d pixels/ring, %ld pixels)\n",
+          nrings,ppring,(long)npix);
+    sharp_make_hw_geom_info (nrings, ppring, 0., 1, ppring, &tinfo);
+    check_accuracy(tinfo,lmax,lmax,npix,spin,ntrans,niter);
+    sharp_destroy_geom_info(tinfo);
+    }
   else
     UTIL_FAIL("unknown grid geometry");