From 94c6a99ad85208e4fb1b013252334a4495b3cb46 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Fri, 23 Nov 2012 13:28:51 +0100
Subject: [PATCH] allow any number of rings in sharp_make_ecp_geom_info()

---
 libsharp/sharp_geomhelpers.c | 49 ++++++++++++------------------------
 libsharp/sharp_geomhelpers.h |  3 +--
 2 files changed, 17 insertions(+), 35 deletions(-)

diff --git a/libsharp/sharp_geomhelpers.c b/libsharp/sharp_geomhelpers.c
index 7ecdd85..b1d0142 100644
--- a/libsharp/sharp_geomhelpers.c
+++ b/libsharp/sharp_geomhelpers.c
@@ -158,21 +158,6 @@ static void gauss_legendre_tbl(int n, double* x, double* w)
     }
   }
 
-static void makeweights (int bw, double *weights)
-  {
-  const double pi = 3.141592653589793238462643383279502884197;
-  const double fudge = pi/(4*bw);
-  for (int j=0; j<2*bw; ++j)
-    {
-    double tmpsum = 0;
-    for (int k=0; k<bw; ++k)
-      tmpsum += 1./(2*k+1) * sin((2*j+1)*(2*k+1)*fudge);
-    tmpsum *= sin((2*j+1)*fudge);
-    tmpsum *= 2./bw;
-    weights[j] = tmpsum;
-    }
-  }
-
 void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
   int stride_lon, int stride_lat, sharp_geom_info **geom_info)
   {
@@ -207,7 +192,8 @@ void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
   DEALLOC(stride_);
   }
 
-void sharp_make_ecp_geom_info (int nrings, int nphi, double phi0,
+/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
+void sharp_make_ecp_geom_info (int nrings, int ppring, double phi0,
   int stride_lon, int stride_lat, sharp_geom_info **geom_info)
   {
   const double pi=3.141592653589793238462643383279502884197;
@@ -219,17 +205,21 @@ void sharp_make_ecp_geom_info (int nrings, int nphi, double phi0,
   ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
   int *stride_=RALLOC(int,nrings);
 
-  UTIL_ASSERT((nrings&1)==0,
-    "Even number of rings needed for equidistant grid!");
-  makeweights(nrings/2,weight);
-  for (int m=0; m<nrings; ++m)
+  for (int m=0; m<(nrings+1)/2; ++m)
     {
-    theta[m] = (m+0.5)*pi/nrings;
-    nph[m]=nphi;
-    phi0_[m]=phi0;
+    theta[m]=pi*(m+0.5)/nrings;
+    theta[nrings-1-m]=pi-theta[m];
+    nph[m]=nph[nrings-1-m]=ppring;
+    phi0_[m]=phi0_[nrings-1-m]=phi0;
     ofs[m]=(ptrdiff_t)m*stride_lat;
-    stride_[m]=stride_lon;
-    weight[m]*=2*pi/nphi;
+    ofs[nrings-1-m]=(ptrdiff_t)((nrings-1-m)*stride_lat);
+    stride_[m]=stride_[nrings-1-m]=stride_lon;
+    double wgt=0;
+    int n=nrings;
+    for (int j=1; j<=n/2; ++j)
+      wgt += cos(((j*(2*m+1))%(2*n))*(pi/nrings))/(4.*j*j-1.);
+    wgt = 2./n*(1.-2*wgt)*2*pi/nph[m];
+    weight[m]=weight[nrings-1-m]=wgt;
     }
 
   sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL, weight,
@@ -243,13 +233,6 @@ void sharp_make_ecp_geom_info (int nrings, int nphi, double phi0,
   DEALLOC(stride_);
   }
 
-static inline 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 Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
 void sharp_make_hw_geom_info (int nrings, int ppring, double phi0,
   int stride_lon, int stride_lat, sharp_geom_info **geom_info)
@@ -263,7 +246,7 @@ void sharp_make_hw_geom_info (int nrings, int ppring, double phi0,
   ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
   int *stride_=RALLOC(int,nrings);
 
-  for (int m=0; m<=nrings/2; ++m)
+  for (int m=0; m<(nrings+1)/2; ++m)
     {
     theta[m]=pi*m/(nrings-1.);
     if (theta[m]<1e-15) theta[m]=1e-15;
diff --git a/libsharp/sharp_geomhelpers.h b/libsharp/sharp_geomhelpers.h
index 80a7ac2..991665d 100644
--- a/libsharp/sharp_geomhelpers.h
+++ b/libsharp/sharp_geomhelpers.h
@@ -25,7 +25,7 @@
 /*! \file sharp_geomhelpers.h
  *  SHARP helper function for the creation of grid geometries
  *
- *  Copyright (C) 2006-2011 Max-Planck-Society
+ *  Copyright (C) 2006-2012 Max-Planck-Society
  *  \author Martin Reinecke
  */
 
@@ -69,7 +69,6 @@ void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
     is \a stride_lat.
     \note The spacing of pixel centers is equidistant in colatitude and
       longitude.
-    \note \a nrings must be an even number.
     \note The sphere is pixelized in a way that the colatitude of the first ring
       is \a 0.5*(pi/nrings) and the colatitude of the last ring is
       \a pi-0.5*(pi/nrings). There are no pixel centers at the poles.