From 7d23ebcc39411219224c29bc098178e762c638ea Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Thu, 22 Nov 2012 23:33:20 +0100
Subject: [PATCH] add MW geometry

---
 libsharp/sharp_geomhelpers.c | 31 +++++++++++++++++++++++++++++++
 libsharp/sharp_geomhelpers.h | 18 +++++++++++++++++-
 2 files changed, 48 insertions(+), 1 deletion(-)

diff --git a/libsharp/sharp_geomhelpers.c b/libsharp/sharp_geomhelpers.c
index 14deb00..7ecdd85 100644
--- a/libsharp/sharp_geomhelpers.c
+++ b/libsharp/sharp_geomhelpers.c
@@ -291,3 +291,34 @@ void sharp_make_hw_geom_info (int nrings, int ppring, double phi0,
   DEALLOC(ofs);
   DEALLOC(stride_);
   }
+
+void sharp_make_mw_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);
+  int *nph=RALLOC(int,nrings);
+  double *phi0_=RALLOC(double,nrings);
+  ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
+  int *stride_=RALLOC(int,nrings);
+
+  for (int m=0; m<nrings; ++m)
+    {
+    theta[m]=pi*(2.*m+1.)/(2.*nrings-1.);
+    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;
+    }
+
+  sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL, NULL,
+    geom_info);
+
+  DEALLOC(theta);
+  DEALLOC(nph);
+  DEALLOC(phi0_);
+  DEALLOC(ofs);
+  DEALLOC(stride_);
+  }
diff --git a/libsharp/sharp_geomhelpers.h b/libsharp/sharp_geomhelpers.h
index ba534d0..80a7ac2 100644
--- a/libsharp/sharp_geomhelpers.h
+++ b/libsharp/sharp_geomhelpers.h
@@ -87,11 +87,27 @@ void sharp_make_ecp_geom_info (int nrings, int nphi, double phi0,
       longitude.
     \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
+    \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);
 
+/*! Creates a geometry information describing a 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 The sphere is pixelized in a way that the colatitude of the first ring
+      is \a pi/(2*nrings-1) and that of the last ring is \a pi.
+    \note This is the grid introduced by McEwen & Wiaux 2011.
+    \note This function does \e not define any quadrature weights.
+    \ingroup geominfogroup */
+void sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
+  int stride_lon, int stride_lat, sharp_geom_info **geom_info);
+
 #ifdef __cplusplus
 }
 #endif