36 #include "pocketfft/pocketfft.h" 39 const int *rings,
const double *weight, sharp_geom_info **geom_info)
41 const double pi=3.141592653589793238462643383279502884197;
42 ptrdiff_t npix=(ptrdiff_t)nside*nside*12;
43 ptrdiff_t ncap=2*(ptrdiff_t)nside*(nside-1);
45 double *theta=
RALLOC(
double,nrings);
46 double *weight_=
RALLOC(
double,nrings);
47 int *nph=
RALLOC(
int,nrings);
48 double *phi0=
RALLOC(
double,nrings);
49 ptrdiff_t *ofs=
RALLOC(ptrdiff_t,nrings);
50 int *stride_=
RALLOC(
int,nrings);
51 ptrdiff_t curofs=0, checkofs;
52 for (
int m=0; m<nrings; ++m)
54 int ring = (rings==NULL)? (m+1) : rings[m];
55 ptrdiff_t northring = (ring>2*nside) ? 4*nside-ring : ring;
57 if (northring < nside)
59 theta[m] = 2*asin(northring/(sqrt(6.)*nside));
62 checkofs = 2*northring*(northring-1)*stride;
66 double fact1 = (8.*nside)/npix;
67 double costheta = (2*nside-northring)*fact1;
68 theta[m] = acos(costheta);
70 if ((northring-nside) & 1)
74 checkofs = (ncap + (northring-nside)*nph[m])*stride;
77 if (northring != ring)
79 theta[m] = pi-theta[m];
80 checkofs = (npix - nph[m])*stride - checkofs;
83 weight_[m]=4.*pi/npix*((weight==NULL) ? 1. : weight[northring-1]);
85 UTIL_ASSERT(curofs==checkofs,
"Bug in computing ofs[m]");
103 const double *weight, sharp_geom_info **geom_info)
109 int stride_lon,
int stride_lat, sharp_geom_info **geom_info)
111 const double pi=3.141592653589793238462643383279502884197;
113 double *theta=
RALLOC(
double,nrings);
114 double *weight=
RALLOC(
double,nrings);
115 int *nph=
RALLOC(
int,nrings);
116 double *phi0_=
RALLOC(
double,nrings);
117 ptrdiff_t *ofs=
RALLOC(ptrdiff_t,nrings);
118 int *stride_=
RALLOC(
int,nrings);
121 for (
int m=0; m<nrings; ++m)
123 theta[m] = acos(-theta[m]);
126 ofs[m]=(ptrdiff_t)m*stride_lat;
127 stride_[m]=stride_lon;
128 weight[m]*=2*pi/nphi;
144 int stride_lon,
int stride_lat, sharp_geom_info **geom_info)
146 const double pi=3.141592653589793238462643383279502884197;
148 double *theta=
RALLOC(
double,nrings);
149 double *weight=
RALLOC(
double,nrings);
150 int *nph=
RALLOC(
int,nrings);
151 double *phi0_=
RALLOC(
double,nrings);
152 ptrdiff_t *ofs=
RALLOC(ptrdiff_t,nrings);
153 int *stride_=
RALLOC(
int,nrings);
156 for (
int k=1; k<=(nrings-1)/2; ++k)
158 weight[2*k-1]=2./(1.-4.*k*k)*cos((k*pi)/nrings);
159 weight[2*k ]=2./(1.-4.*k*k)*sin((k*pi)/nrings);
161 if ((nrings&1)==0) weight[nrings-1]=0.;
162 rfft_plan plan = make_rfft_plan(nrings);
163 rfft_backward(plan,weight,1.);
164 destroy_rfft_plan(plan);
166 for (
int m=0; m<(nrings+1)/2; ++m)
168 theta[m]=pi*(m+0.5)/nrings;
169 theta[nrings-1-m]=pi-theta[m];
170 nph[m]=nph[nrings-1-m]=ppring;
171 phi0_[m]=phi0_[nrings-1-m]=phi0;
172 ofs[m]=(ptrdiff_t)m*stride_lat;
173 ofs[nrings-1-m]=(ptrdiff_t)((nrings-1-m)*stride_lat);
174 stride_[m]=stride_[nrings-1-m]=stride_lon;
175 weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(nrings*nph[m]);
191 int stride_lon,
int stride_lat, sharp_geom_info **geom_info)
193 const double pi=3.141592653589793238462643383279502884197;
195 double *theta=
RALLOC(
double,nrings);
196 double *weight=
RALLOC(
double,nrings);
197 int *nph=
RALLOC(
int,nrings);
198 double *phi0_=
RALLOC(
double,nrings);
199 ptrdiff_t *ofs=
RALLOC(ptrdiff_t,nrings);
200 int *stride_=
RALLOC(
int,nrings);
204 double dw=-1./(n*n-1.+(n&1));
206 for (
int k=1; k<=(n/2-1); ++k)
207 weight[2*k-1]=2./(1.-4.*k*k) + dw;
208 weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1. -dw*((2-(n&1))*n-1);
209 rfft_plan plan = make_rfft_plan(n);
210 rfft_backward(plan,weight,1.);
211 destroy_rfft_plan(plan);
214 for (
int m=0; m<(nrings+1)/2; ++m)
216 theta[m]=pi*m/(nrings-1.);
217 if (theta[m]<1e-15) theta[m]=1e-15;
218 theta[nrings-1-m]=pi-theta[m];
219 nph[m]=nph[nrings-1-m]=ppring;
220 phi0_[m]=phi0_[nrings-1-m]=phi0;
221 ofs[m]=(ptrdiff_t)m*stride_lat;
222 ofs[nrings-1-m]=(ptrdiff_t)((nrings-1-m)*stride_lat);
223 stride_[m]=stride_[nrings-1-m]=stride_lon;
224 weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
240 int stride_lon,
int stride_lat, sharp_geom_info **geom_info)
242 const double pi=3.141592653589793238462643383279502884197;
244 double *theta=
RALLOC(
double,nrings);
245 double *weight=
RALLOC(
double,nrings+1);
246 int *nph=
RALLOC(
int,nrings);
247 double *phi0_=
RALLOC(
double,nrings);
248 ptrdiff_t *ofs=
RALLOC(ptrdiff_t,nrings);
249 int *stride_=
RALLOC(
int,nrings);
254 for (
int k=1; k<=(n/2-1); ++k)
255 weight[2*k-1]=2./(1.-4.*k*k);
256 weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1.;
257 rfft_plan plan = make_rfft_plan(n);
258 rfft_backward(plan,weight,1.);
259 destroy_rfft_plan(plan);
260 for (
int m=0; m<nrings; ++m)
261 weight[m]=weight[m+1];
263 for (
int m=0; m<(nrings+1)/2; ++m)
265 theta[m]=pi*(m+1)/(nrings+1.);
266 theta[nrings-1-m]=pi-theta[m];
267 nph[m]=nph[nrings-1-m]=ppring;
268 phi0_[m]=phi0_[nrings-1-m]=phi0;
269 ofs[m]=(ptrdiff_t)m*stride_lat;
270 ofs[nrings-1-m]=(ptrdiff_t)((nrings-1-m)*stride_lat);
271 stride_[m]=stride_[nrings-1-m]=stride_lon;
272 weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
287 int stride_lon,
int stride_lat, sharp_geom_info **geom_info)
289 const double pi=3.141592653589793238462643383279502884197;
291 double *theta=
RALLOC(
double,nrings);
292 int *nph=
RALLOC(
int,nrings);
293 double *phi0_=
RALLOC(
double,nrings);
294 ptrdiff_t *ofs=
RALLOC(ptrdiff_t,nrings);
295 int *stride_=
RALLOC(
int,nrings);
297 for (
int m=0; m<nrings; ++m)
299 theta[m]=pi*(2.*m+1.)/(2.*nrings-1.);
300 if (theta[m]>pi-1e-15) theta[m]=pi-1e-15;
303 ofs[m]=(ptrdiff_t)m*stride_lat;
304 stride_[m]=stride_lon;
#define RALLOC(type, num)
void sharp_make_fejer2_geom_info(int nrings, int ppring, double phi0, int stride_lon, int stride_lat, sharp_geom_info **geom_info)
#define SET_ARRAY(ptr, i1, i2, val)
void sharp_make_subset_healpix_geom_info(int nside, int stride, int nrings, const int *rings, const double *weight, sharp_geom_info **geom_info)
void sharp_make_geom_info(int nrings, const int *nph, const ptrdiff_t *ofs, const int *stride, const double *phi0, const double *theta, const double *wgt, sharp_geom_info **geom_info)
void sharp_make_cc_geom_info(int nrings, int ppring, double phi0, int stride_lon, int stride_lat, sharp_geom_info **geom_info)
void sharp_make_gauss_geom_info(int nrings, int nphi, double phi0, int stride_lon, int stride_lat, sharp_geom_info **geom_info)
void sharp_make_weighted_healpix_geom_info(int nside, int stride, const double *weight, sharp_geom_info **geom_info)
#define UTIL_ASSERT(cond, msg)
void sharp_make_mw_geom_info(int nrings, int ppring, double phi0, int stride_lon, int stride_lat, sharp_geom_info **geom_info)
void sharp_legendre_roots(int n, double *x, double *w)
void sharp_make_fejer1_geom_info(int nrings, int ppring, double phi0, int stride_lon, int stride_lat, sharp_geom_info **geom_info)