LevelS SHT library  3.50
sharp_geomhelpers.c
Go to the documentation of this file.
1 /*
2  * This file is part of libsharp.
3  *
4  * libsharp is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * libsharp is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with libsharp; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18 
19 /*
20  * libsharp is being developed at the Max-Planck-Institut fuer Astrophysik
21  * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
22  * (DLR).
23  */
24 
25 /*! \file sharp_geomhelpers.c
26  * Spherical transform library
27  *
28  * Copyright (C) 2006-2018 Max-Planck-Society
29  * \author Martin Reinecke
30  */
31 
32 #include <math.h>
33 #include "sharp_geomhelpers.h"
34 #include "sharp_legendre_roots.h"
35 #include "c_utils.h"
36 #include "pocketfft/pocketfft.h"
37 
38 void sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
39  const int *rings, const double *weight, sharp_geom_info **geom_info)
40  {
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);
44 
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; /* checkofs used for assertion introduced when adding rings arg */
52  for (int m=0; m<nrings; ++m)
53  {
54  int ring = (rings==NULL)? (m+1) : rings[m];
55  ptrdiff_t northring = (ring>2*nside) ? 4*nside-ring : ring;
56  stride_[m] = stride;
57  if (northring < nside)
58  {
59  theta[m] = 2*asin(northring/(sqrt(6.)*nside));
60  nph[m] = 4*northring;
61  phi0[m] = pi/nph[m];
62  checkofs = 2*northring*(northring-1)*stride;
63  }
64  else
65  {
66  double fact1 = (8.*nside)/npix;
67  double costheta = (2*nside-northring)*fact1;
68  theta[m] = acos(costheta);
69  nph[m] = 4*nside;
70  if ((northring-nside) & 1)
71  phi0[m] = 0;
72  else
73  phi0[m] = pi/nph[m];
74  checkofs = (ncap + (northring-nside)*nph[m])*stride;
75  ofs[m] = curofs;
76  }
77  if (northring != ring) /* southern hemisphere */
78  {
79  theta[m] = pi-theta[m];
80  checkofs = (npix - nph[m])*stride - checkofs;
81  ofs[m] = curofs;
82  }
83  weight_[m]=4.*pi/npix*((weight==NULL) ? 1. : weight[northring-1]);
84  if (rings==NULL) {
85  UTIL_ASSERT(curofs==checkofs, "Bug in computing ofs[m]");
86  }
87  ofs[m] = curofs;
88  curofs+=nph[m];
89  }
90 
91  sharp_make_geom_info (nrings, nph, ofs, stride_, phi0, theta, weight_,
92  geom_info);
93 
94  DEALLOC(theta);
95  DEALLOC(weight_);
96  DEALLOC(nph);
97  DEALLOC(phi0);
98  DEALLOC(ofs);
99  DEALLOC(stride_);
100  }
101 
102 void sharp_make_weighted_healpix_geom_info (int nside, int stride,
103  const double *weight, sharp_geom_info **geom_info)
104  {
105  sharp_make_subset_healpix_geom_info(nside, stride, 4 * nside - 1, NULL, weight, geom_info);
106  }
107 
108 void sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
109  int stride_lon, int stride_lat, sharp_geom_info **geom_info)
110  {
111  const double pi=3.141592653589793238462643383279502884197;
112 
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);
119 
120  sharp_legendre_roots(nrings,theta,weight);
121  for (int m=0; m<nrings; ++m)
122  {
123  theta[m] = acos(-theta[m]);
124  nph[m]=nphi;
125  phi0_[m]=phi0;
126  ofs[m]=(ptrdiff_t)m*stride_lat;
127  stride_[m]=stride_lon;
128  weight[m]*=2*pi/nphi;
129  }
130 
131  sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight,
132  geom_info);
133 
134  DEALLOC(theta);
135  DEALLOC(weight);
136  DEALLOC(nph);
137  DEALLOC(phi0_);
138  DEALLOC(ofs);
139  DEALLOC(stride_);
140  }
141 
142 /* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
143 void sharp_make_fejer1_geom_info (int nrings, int ppring, double phi0,
144  int stride_lon, int stride_lat, sharp_geom_info **geom_info)
145  {
146  const double pi=3.141592653589793238462643383279502884197;
147 
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);
154 
155  weight[0]=2.;
156  for (int k=1; k<=(nrings-1)/2; ++k)
157  {
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);
160  }
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);
165 
166  for (int m=0; m<(nrings+1)/2; ++m)
167  {
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]);
176  }
177 
178  sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight,
179  geom_info);
180 
181  DEALLOC(theta);
182  DEALLOC(weight);
183  DEALLOC(nph);
184  DEALLOC(phi0_);
185  DEALLOC(ofs);
186  DEALLOC(stride_);
187  }
188 
189 /* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
190 void sharp_make_cc_geom_info (int nrings, int ppring, double phi0,
191  int stride_lon, int stride_lat, sharp_geom_info **geom_info)
192  {
193  const double pi=3.141592653589793238462643383279502884197;
194 
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);
201 
202  int n=nrings-1;
203  SET_ARRAY(weight,0,nrings,0.);
204  double dw=-1./(n*n-1.+(n&1));
205  weight[0]=2.+dw;
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);
212  weight[n]=weight[0];
213 
214  for (int m=0; m<(nrings+1)/2; ++m)
215  {
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]);
225  }
226 
227  sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight,
228  geom_info);
229 
230  DEALLOC(theta);
231  DEALLOC(weight);
232  DEALLOC(nph);
233  DEALLOC(phi0_);
234  DEALLOC(ofs);
235  DEALLOC(stride_);
236  }
237 
238 /* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
239 void sharp_make_fejer2_geom_info (int nrings, int ppring, double phi0,
240  int stride_lon, int stride_lat, sharp_geom_info **geom_info)
241  {
242  const double pi=3.141592653589793238462643383279502884197;
243 
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);
250 
251  int n=nrings+1;
252  SET_ARRAY(weight,0,n,0.);
253  weight[0]=2.;
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];
262 
263  for (int m=0; m<(nrings+1)/2; ++m)
264  {
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]);
273  }
274 
275  sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, weight,
276  geom_info);
277 
278  DEALLOC(theta);
279  DEALLOC(weight);
280  DEALLOC(nph);
281  DEALLOC(phi0_);
282  DEALLOC(ofs);
283  DEALLOC(stride_);
284  }
285 
286 void sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
287  int stride_lon, int stride_lat, sharp_geom_info **geom_info)
288  {
289  const double pi=3.141592653589793238462643383279502884197;
290 
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);
296 
297  for (int m=0; m<nrings; ++m)
298  {
299  theta[m]=pi*(2.*m+1.)/(2.*nrings-1.);
300  if (theta[m]>pi-1e-15) theta[m]=pi-1e-15;
301  nph[m]=ppring;
302  phi0_[m]=phi0;
303  ofs[m]=(ptrdiff_t)m*stride_lat;
304  stride_[m]=stride_lon;
305  }
306 
307  sharp_make_geom_info (nrings, nph, ofs, stride_, phi0_, theta, NULL,
308  geom_info);
309 
310  DEALLOC(theta);
311  DEALLOC(nph);
312  DEALLOC(phi0_);
313  DEALLOC(ofs);
314  DEALLOC(stride_);
315  }
#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)
Definition: sharp.c:198
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)
#define DEALLOC(ptr)
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)

Generated on Mon Dec 10 2018 10:24:20 for LevelS SHT library