Aquila-Consortium/libsharp2
0
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

add support for a Fejer 2 type grid

Browse source
This commit is contained in:
Martin Reinecke 2012-11-27 16:32:34 +01:00
parent ccf358f9d4
commit a9bc412c7f
2 changed files with 62 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

47
libsharp/sharp_geomhelpers.c
View file

@ -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<nrings; ++m)
weight[m]=weight[m+1];
for (int m=0; m<(nrings+1)/2; ++m)
{
theta[m]=pi*(m+1)/(nrings+1.);
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;
ofs[nrings-1-m]=(ptrdiff_t)((nrings-1-m)*stride_lat);
stride_[m]=stride_[nrings-1-m]=stride_lon;
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*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_);
}
void sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
int stride_lon, int stride_lat, sharp_geom_info **geom_info)
{

15
libsharp/sharp_geomhelpers.h
View file

@ -91,6 +91,21 @@ void sharp_make_ecp_geom_info (int nrings, int nphi, double phi0,
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 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 The sphere is pixelized in a way that the colatitude of the first ring
is \a pi/(nrings+1) and that of the last ring is \a pi-pi/(nrings+1).
\note This is the grid used by Huffenberger & Wandelt 2010.
\ingroup geominfogroup */
void sharp_make_fejer2_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