diff --git a/libsharp/sharp_geomhelpers.c b/libsharp/sharp_geomhelpers.c index 520ff25..1653534 100644 --- a/libsharp/sharp_geomhelpers.c +++ b/libsharp/sharp_geomhelpers.c @@ -289,6 +289,53 @@ void sharp_make_hw_geom_info (int nrings, int ppring, double phi0, DEALLOC(stride_); } +/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */ +void sharp_make_fejer2_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+1); + int *nph=RALLOC(int,nrings); + double *phi0_=RALLOC(double,nrings); + ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings); + int *stride_=RALLOC(int,nrings); + + int n=nrings+1; + SET_ARRAY(weight,0,n,0.); + weight[0]=2.; + for (int k=1; k<=(n/2-1); ++k) + weight[2*k-1]=2./(1.-4.*k*k); + weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1.; + real_plan plan = make_real_plan(n); + real_plan_backward_fftpack(plan,weight); + kill_real_plan(plan); + for (int m=0; m