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

sharp_make_subset_healpix_geom_info

Browse source
This commit is contained in:
Dag Sverre Seljebotn 2015-03-30 12:17:30 +02:00
parent cf7de6ba70
commit c3b6277545
2 changed files with 33 additions and 7 deletions
Download patch file Download diff file Expand all files Collapse all files

28
libsharp/sharp_geomhelpers.c
View file

@ -34,14 +34,14 @@
#include "sharp_geomhelpers.h"
#include "c_utils.h"
#include "ls_fft.h"
#include <stdio.h>
void sharp_make_weighted_healpix_geom_info (int nside, int stride,
const double *weight, sharp_geom_info **geom_info)
void sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
const int *rings, const double *weight, sharp_geom_info **geom_info)
{
const double pi=3.141592653589793238462643383279502884197;
ptrdiff_t npix=(ptrdiff_t)nside*nside*12;
ptrdiff_t ncap=2*(ptrdiff_t)nside*(nside-1);
int nrings=4*nside-1;
double *theta=RALLOC(double,nrings);
double *weight_=RALLOC(double,nrings);
@ -49,9 +49,10 @@ void sharp_make_weighted_healpix_geom_info (int nside, int stride,
double *phi0=RALLOC(double,nrings);
ptrdiff_t *ofs=RALLOC(ptrdiff_t,nrings);
int *stride_=RALLOC(int,nrings);
ptrdiff_t curofs=0, checkofs; /* checkofs used for assertion introduced when adding rings arg */
for (int m=0; m<nrings; ++m)
{
int ring=m+1;
int ring = (rings==NULL)? (m+1) : rings[m];
ptrdiff_t northring = (ring>2*nside) ? 4*nside-ring : ring;
stride_[m] = stride;
if (northring < nside)
@ -59,7 +60,7 @@ void sharp_make_weighted_healpix_geom_info (int nside, int stride,
theta[m] = 2*asin(northring/(sqrt(6.)*nside));
nph[m] = 4*northring;
phi0[m] = pi/nph[m];
ofs[m] = 2*northring*(northring-1)*stride;
checkofs = 2*northring*(northring-1)*stride;
}
else
{
@ -71,14 +72,21 @@ void sharp_make_weighted_healpix_geom_info (int nside, int stride,
phi0[m] = 0;
else
phi0[m] = pi/nph[m];
ofs[m] = (ncap + (northring-nside)*nph[m])*stride;
checkofs = (ncap + (northring-nside)*nph[m])*stride;
ofs[m] = curofs;
}
if (northring != ring) /* southern hemisphere */
{
theta[m] = pi-theta[m];
ofs[m] = (npix - nph[m])*stride - ofs[m];
checkofs = (npix - nph[m])*stride - checkofs;
ofs[m] = curofs;
}
weight_[m]=4.*pi/npix*((weight==NULL) ? 1. : weight[northring-1]);
if (rings==NULL) {
UTIL_ASSERT(curofs==checkofs, "Bug in computing ofs[m]");
}
ofs[m] = curofs;
curofs+=nph[m];
}
sharp_make_geom_info (nrings, nph, ofs, stride_, phi0, theta, weight_,
@ -92,6 +100,12 @@ void sharp_make_weighted_healpix_geom_info (int nside, int stride,
DEALLOC(stride_);
}
void sharp_make_weighted_healpix_geom_info (int nside, int stride,
const double *weight, sharp_geom_info **geom_info)
{
sharp_make_subset_healpix_geom_info(nside, stride, 4 * nside - 1, NULL, weight, geom_info);
}
static inline double one_minus_x2 (double x)
{ return (fabs(x)>0.1) ? (1.+x)*(1.-x) : 1.-x*x; }

12
libsharp/sharp_geomhelpers.h
View file

@ -38,6 +38,18 @@
extern "C" {
#endif
/*! Creates a geometry information describing a HEALPix map with an
Nside parameter \a nside. \a weight contains the relative ring
weights and must have \a 2*nside entries. The rings array contains
the indices of the rings, with 1 being the first ring at the north
pole; if NULL then we take them to be sequential. Pass 4 * nside - 1
as nrings and NULL to rings to get the full HEALPix grid.
\note if \a weight is a null pointer, all weights are assumed to be 1.
\note if \a rings is a null pointer, take all rings
\ingroup geominfogroup */
void sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
const int *rings, const double *weight, sharp_geom_info **geom_info);
/*! Creates a geometry information describing a HEALPix map with an
Nside parameter \a nside. \a weight contains the relative ring
weights and must have \a 2*nside entries.