sharp_geomhelpers.cc 12.1 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
/*
 *  This file is part of libsharp2.
 *
 *  libsharp2 is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  libsharp2 is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with libsharp2; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/* libsharp2 is being developed at the Max-Planck-Institut fuer Astrophysik */

/*! \file sharp_geomhelpers.c
 *  Spherical transform library
 *
24
 *  Copyright (C) 2006-2020 Max-Planck-Society
Martin Reinecke's avatar
Martin Reinecke committed
25
26
27
 *  \author Martin Reinecke
 */

Martin Reinecke's avatar
Martin Reinecke committed
28
29
#include <cmath>
#include <vector>
Martin Reinecke's avatar
Martin Reinecke committed
30
31
32
33
#include "libsharp2/sharp_geomhelpers.h"
#include "mr_util/gl_integrator.h"
#include "mr_util/fft.h"
#include "mr_util/error_handling.h"
34
#include "mr_util/math_utils.h"
Martin Reinecke's avatar
Martin Reinecke committed
35

Martin Reinecke's avatar
Martin Reinecke committed
36
using namespace std;
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
using namespace mr;

sharp_standard_geom_info::sharp_standard_geom_info(size_t nrings, const size_t *nph, const ptrdiff_t *ofs,
  ptrdiff_t stride_, const double *phi0, const double *theta, const double *wgt)
  : ring(nrings), stride(stride_)
  {
  size_t pos=0;

  nphmax_=0;

  for (size_t m=0; m<nrings; ++m)
    {
    ring[m].theta = theta[m];
    ring[m].cth = cos(theta[m]);
    ring[m].sth = sin(theta[m]);
    ring[m].weight = (wgt != nullptr) ? wgt[m] : 1.;
    ring[m].phi0 = phi0[m];
    ring[m].ofs = ofs[m];
    ring[m].nph = nph[m];
    if (nphmax_<nph[m]) nphmax_=nph[m];
    }
  sort(ring.begin(), ring.end(),[](const Tring &a, const Tring &b)
    { return (a.sth<b.sth); });
  while (pos<nrings)
    {
    pair_.push_back(Tpair());
    pair_.back().r1=pos;
    if ((pos<nrings-1) && approx(ring[pos].cth,-ring[pos+1].cth,1e-12))
      {
      if (ring[pos].cth>0)  // make sure northern ring is in r1
        pair_.back().r2=pos+1;
      else
        {
        pair_.back().r1=pos+1;
        pair_.back().r2=pos;
        }
      ++pos;
      }
    else
      pair_.back().r2=size_t(~0);
    ++pos;
    }

  sort(pair_.begin(), pair_.end(), [this] (const Tpair &a, const Tpair &b)
    {
    if (ring[a.r1].nph==ring[b.r1].nph)
    return (ring[a.r1].phi0 < ring[b.r1].phi0) ? true :
      ((ring[a.r1].phi0 > ring[b.r1].phi0) ? false :
        (ring[a.r1].cth>ring[b.r1].cth));
    return ring[a.r1].nph<ring[b.r1].nph;
    });
  }
Martin Reinecke's avatar
Martin Reinecke committed
89

Martin Reinecke's avatar
Martin Reinecke committed
90
template<typename T> void sharp_standard_geom_info::tclear(T *map) const
91
92
93
94
  {
  for (const auto &r: ring)
    {
    if (stride==1)
Martin Reinecke's avatar
Martin Reinecke committed
95
      memset(&map[r.ofs],0,r.nph*sizeof(T));
96
97
    else
      for (size_t i=0;i<r.nph;++i)
Martin Reinecke's avatar
Martin Reinecke committed
98
        map[r.ofs+ptrdiff_t(i)*stride]=T(0);
99
100
101
    }
  }

Martin Reinecke's avatar
Martin Reinecke committed
102
void sharp_standard_geom_info::clear_map (const any &map) const
103
  {
Martin Reinecke's avatar
Martin Reinecke committed
104
105
  if (map.type()==typeid(double *)) tclear(any_cast<double *>(map));
  else if (map.type()==typeid(float *)) tclear(any_cast<float *>(map));
Martin Reinecke's avatar
Martin Reinecke committed
106
  else MR_fail("bad map data type");
107
108
  }

Martin Reinecke's avatar
Martin Reinecke committed
109
template<typename T> void sharp_standard_geom_info::tadd(bool weighted, size_t iring, const double *ringtmp, T *map) const
110
  {
Martin Reinecke's avatar
Martin Reinecke committed
111
  T *MRUTIL_RESTRICT p1=&map[ring[iring].ofs];
112
113
  double wgt = weighted ? ring[iring].weight : 1.;
  for (size_t m=0; m<ring[iring].nph; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
114
    p1[ptrdiff_t(m)*stride] += T(ringtmp[m]*wgt);
115
116
  }
//virtual
Martin Reinecke's avatar
Martin Reinecke committed
117
void sharp_standard_geom_info::add_ring(bool weighted, size_t iring, const double *ringtmp, const any &map) const
118
  {
Martin Reinecke's avatar
Martin Reinecke committed
119
120
  if (map.type()==typeid(double *)) tadd(weighted, iring, ringtmp, any_cast<double *>(map));
  else if (map.type()==typeid(float *)) tadd(weighted, iring, ringtmp, any_cast<float *>(map));
Martin Reinecke's avatar
Martin Reinecke committed
121
  else MR_fail("bad map data type");
122
  }
Martin Reinecke's avatar
Martin Reinecke committed
123
template<typename T> void sharp_standard_geom_info::tget(bool weighted, size_t iring, const T *map, double *ringtmp) const
124
  {
Martin Reinecke's avatar
Martin Reinecke committed
125
  const T *MRUTIL_RESTRICT p1=&map[ring[iring].ofs];
126
127
  double wgt = weighted ? ring[iring].weight : 1.;
  for (size_t m=0; m<ring[iring].nph; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
128
    ringtmp[m] = p1[ptrdiff_t(m)*stride]*wgt;
129
130
  }
//virtual
Martin Reinecke's avatar
Martin Reinecke committed
131
void sharp_standard_geom_info::get_ring(bool weighted, size_t iring, const any &map, double *ringtmp) const
132
  {
Martin Reinecke's avatar
Martin Reinecke committed
133
134
135
136
  if (map.type()==typeid(const double *)) tget(weighted, iring, any_cast<const double *>(map), ringtmp);
  else if (map.type()==typeid(double *)) tget(weighted, iring, any_cast<double *>(map), ringtmp);
  else if (map.type()==typeid(const float *)) tget(weighted, iring, any_cast<const float *>(map), ringtmp);
  else if (map.type()==typeid(float *)) tget(weighted, iring, any_cast<float *>(map), ringtmp);
Martin Reinecke's avatar
Martin Reinecke committed
137
  else MR_assert(false,"bad map data type",map.type().name());
138
  }
Martin Reinecke's avatar
Martin Reinecke committed
139

Martin Reinecke's avatar
Martin Reinecke committed
140
141
unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (size_t nside, ptrdiff_t stride, size_t nrings,
  const size_t *rings, const double *weight)
Martin Reinecke's avatar
Martin Reinecke committed
142
143
  {
  const double pi=3.141592653589793238462643383279502884197;
Martin Reinecke's avatar
Martin Reinecke committed
144
145
  size_t npix=nside*nside*12;
  size_t ncap=2*nside*(nside-1);
Martin Reinecke's avatar
Martin Reinecke committed
146

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
147
  vector<double> theta(nrings), weight_(nrings), phi0(nrings);
148
  vector<size_t> nph(nrings);
149
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
150
  ptrdiff_t curofs=0, checkofs; /* checkofs used for assertion introduced when adding rings arg */
Martin Reinecke's avatar
Martin Reinecke committed
151
  for (size_t m=0; m<nrings; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
152
    {
Martin Reinecke's avatar
Martin Reinecke committed
153
154
    auto ring = (rings==nullptr)? (m+1) : rings[m];
    size_t northring = (ring>2*nside) ? 4*nside-ring : ring;
Martin Reinecke's avatar
Martin Reinecke committed
155
156
157
158
159
    if (northring < nside)
      {
      theta[m] = 2*asin(northring/(sqrt(6.)*nside));
      nph[m] = 4*northring;
      phi0[m] = pi/nph[m];
Martin Reinecke's avatar
Martin Reinecke committed
160
      checkofs = ptrdiff_t(2*northring*(northring-1))*stride;
Martin Reinecke's avatar
Martin Reinecke committed
161
162
163
164
      }
    else
      {
      double fact1 = (8.*nside)/npix;
Martin Reinecke's avatar
Martin Reinecke committed
165
      double costheta = 2*nside-northring*fact1;
Martin Reinecke's avatar
Martin Reinecke committed
166
167
168
169
170
171
      theta[m] = acos(costheta);
      nph[m] = 4*nside;
      if ((northring-nside) & 1)
        phi0[m] = 0;
      else
        phi0[m] = pi/nph[m];
Martin Reinecke's avatar
Martin Reinecke committed
172
      checkofs = ptrdiff_t(ncap + (northring-nside)*nph[m])*stride;
Martin Reinecke's avatar
Martin Reinecke committed
173
174
175
176
177
      ofs[m] = curofs;
      }
    if (northring != ring) /* southern hemisphere */
      {
      theta[m] = pi-theta[m];
Martin Reinecke's avatar
Martin Reinecke committed
178
      checkofs = ptrdiff_t(npix - nph[m])*stride - checkofs;
Martin Reinecke's avatar
Martin Reinecke committed
179
180
      ofs[m] = curofs;
      }
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
181
182
    weight_[m]=4.*pi/npix*((weight==nullptr) ? 1. : weight[northring-1]);
    if (rings==nullptr) {
Martin Reinecke's avatar
Martin Reinecke committed
183
184
185
      MR_assert(curofs==checkofs, "Bug in computing ofs[m]");
    }
    ofs[m] = curofs;
Martin Reinecke's avatar
Martin Reinecke committed
186
    curofs+=ptrdiff_t(nph[m]);
Martin Reinecke's avatar
Martin Reinecke committed
187
188
    }

Martin Reinecke's avatar
Martin Reinecke committed
189
  return make_unique<sharp_standard_geom_info>(nrings, nph.data(), ofs.data(), stride, phi0.data(), theta.data(), weight_.data());
Martin Reinecke's avatar
Martin Reinecke committed
190
191
  }

Martin Reinecke's avatar
Martin Reinecke committed
192
unique_ptr<sharp_geom_info> sharp_make_weighted_healpix_geom_info (size_t nside, ptrdiff_t stride,
Martin Reinecke's avatar
Martin Reinecke committed
193
  const double *weight)
Martin Reinecke's avatar
Martin Reinecke committed
194
  {
Martin Reinecke's avatar
Martin Reinecke committed
195
  return sharp_make_subset_healpix_geom_info(nside, stride, 4*nside-1, nullptr, weight);
Martin Reinecke's avatar
Martin Reinecke committed
196
197
  }

Martin Reinecke's avatar
Martin Reinecke committed
198
199
unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (size_t nrings, size_t nphi, double phi0,
  ptrdiff_t stride_lon, ptrdiff_t stride_lat)
Martin Reinecke's avatar
Martin Reinecke committed
200
201
202
  {
  const double pi=3.141592653589793238462643383279502884197;

203
204
  vector<size_t> nph(nrings, nphi);
  vector<double> phi0_(nrings, phi0);
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
205
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
206
207
208
209

  mr::GL_Integrator integ(nrings);
  auto theta = integ.coords();
  auto weight = integ.weights();
Martin Reinecke's avatar
Martin Reinecke committed
210
  for (size_t m=0; m<nrings; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
211
212
    {
    theta[m] = acos(-theta[m]);
Martin Reinecke's avatar
Martin Reinecke committed
213
    ofs[m]=ptrdiff_t(m)*stride_lat;
Martin Reinecke's avatar
Martin Reinecke committed
214
215
216
    weight[m]*=2*pi/nphi;
    }

Martin Reinecke's avatar
Martin Reinecke committed
217
  return make_unique<sharp_standard_geom_info>(nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), weight.data());
Martin Reinecke's avatar
Martin Reinecke committed
218
219
220
  }

/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
Martin Reinecke's avatar
Martin Reinecke committed
221
222
unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (size_t nrings, size_t ppring, double phi0,
  ptrdiff_t stride_lon, ptrdiff_t stride_lat)
Martin Reinecke's avatar
Martin Reinecke committed
223
224
225
  {
  const double pi=3.141592653589793238462643383279502884197;

226
227
  vector<double> theta(nrings), weight(nrings), phi0_(nrings, phi0);
  vector<size_t> nph(nrings, ppring);
228
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
229
230

  weight[0]=2.;
Martin Reinecke's avatar
Martin Reinecke committed
231
  for (size_t k=1; k<=(nrings-1)/2; ++k)
Martin Reinecke's avatar
Martin Reinecke committed
232
233
234
235
236
    {
    weight[2*k-1]=2./(1.-4.*k*k)*cos((k*pi)/nrings);
    weight[2*k  ]=2./(1.-4.*k*k)*sin((k*pi)/nrings);
    }
  if ((nrings&1)==0) weight[nrings-1]=0.;
Martin Reinecke's avatar
various    
Martin Reinecke committed
237
238
  auto tmp = fmav(weight.data(),{nrings});
  mr::r2r_fftpack(tmp, tmp, {0}, false, false, 1.);
Martin Reinecke's avatar
Martin Reinecke committed
239

Martin Reinecke's avatar
Martin Reinecke committed
240
  for (size_t m=0; m<(nrings+1)/2; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
241
242
243
    {
    theta[m]=pi*(m+0.5)/nrings;
    theta[nrings-1-m]=pi-theta[m];
Martin Reinecke's avatar
Martin Reinecke committed
244
245
    ofs[m]=ptrdiff_t(m)*stride_lat;
    ofs[nrings-1-m]=ptrdiff_t(nrings-1-m)*stride_lat;
Martin Reinecke's avatar
Martin Reinecke committed
246
247
248
    weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(nrings*nph[m]);
    }

Martin Reinecke's avatar
Martin Reinecke committed
249
  return make_unique<sharp_standard_geom_info>(nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), weight.data());
Martin Reinecke's avatar
Martin Reinecke committed
250
251
252
  }

/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
Martin Reinecke's avatar
Martin Reinecke committed
253
254
unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (size_t nrings, size_t ppring, double phi0,
  ptrdiff_t stride_lon, ptrdiff_t stride_lat)
Martin Reinecke's avatar
Martin Reinecke committed
255
256
257
  {
  const double pi=3.141592653589793238462643383279502884197;

258
259
  vector<double> theta(nrings), weight(nrings,0.), phi0_(nrings, phi0);
  vector<size_t> nph(nrings, ppring);
260
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
261

Martin Reinecke's avatar
Martin Reinecke committed
262
  size_t n=nrings-1;
Martin Reinecke's avatar
Martin Reinecke committed
263
264
  double dw=-1./(n*n-1.+(n&1));
  weight[0]=2.+dw;
Martin Reinecke's avatar
Martin Reinecke committed
265
  for (size_t k=1; k<=(n/2-1); ++k)
Martin Reinecke's avatar
Martin Reinecke committed
266
267
    weight[2*k-1]=2./(1.-4.*k*k) + dw;
  weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1. -dw*((2-(n&1))*n-1);
Martin Reinecke's avatar
various    
Martin Reinecke committed
268
269
  auto tmp = fmav(weight.data(),{n});
  mr::r2r_fftpack(tmp, tmp, {0}, false, false, 1.);
Martin Reinecke's avatar
Martin Reinecke committed
270
271
  weight[n]=weight[0];

Martin Reinecke's avatar
Martin Reinecke committed
272
  for (size_t m=0; m<(nrings+1)/2; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
273
274
275
276
    {
    theta[m]=pi*m/(nrings-1.);
    if (theta[m]<1e-15) theta[m]=1e-15;
    theta[nrings-1-m]=pi-theta[m];
Martin Reinecke's avatar
Martin Reinecke committed
277
278
    ofs[m]=ptrdiff_t(m)*stride_lat;
    ofs[nrings-1-m]=ptrdiff_t(nrings-1-m)*stride_lat;
Martin Reinecke's avatar
Martin Reinecke committed
279
280
281
    weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
    }

Martin Reinecke's avatar
Martin Reinecke committed
282
  return make_unique<sharp_standard_geom_info>(nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), weight.data());
Martin Reinecke's avatar
Martin Reinecke committed
283
284
  }

285
286
287
288
289
290
291
292
293
294
295
296
297
static vector<double> get_dh_weights(size_t nrings)
  {
  vector<double> weight(nrings);

  weight[0]=2.;
  for (size_t k=1; k<=(nrings/2-1); ++k)
    weight[2*k-1]=2./(1.-4.*k*k);
  weight[2*(nrings/2)-1]=(nrings-3.)/(2*(nrings/2)-1) -1.;
  auto tmp = fmav(weight.data(),{nrings});
  mr::r2r_fftpack(tmp, tmp, {0}, false, false, 1.);
  return weight;
  }

Martin Reinecke's avatar
Martin Reinecke committed
298
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
Martin Reinecke's avatar
Martin Reinecke committed
299
300
unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (size_t nrings, size_t ppring, double phi0,
  ptrdiff_t stride_lon, ptrdiff_t stride_lat)
Martin Reinecke's avatar
Martin Reinecke committed
301
302
303
  {
  const double pi=3.141592653589793238462643383279502884197;

304
305
  vector<double> theta(nrings), weight(get_dh_weights(nrings+1)), phi0_(nrings, phi0);
  vector<size_t> nph(nrings, ppring);
306
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
307

Martin Reinecke's avatar
Martin Reinecke committed
308
  for (size_t m=0; m<nrings; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
309
310
    weight[m]=weight[m+1];

Martin Reinecke's avatar
Martin Reinecke committed
311
  for (size_t m=0; m<(nrings+1)/2; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
312
313
314
    {
    theta[m]=pi*(m+1)/(nrings+1.);
    theta[nrings-1-m]=pi-theta[m];
Martin Reinecke's avatar
Martin Reinecke committed
315
316
    ofs[m]=ptrdiff_t(m)*stride_lat;
    ofs[nrings-1-m]=ptrdiff_t(nrings-1-m)*stride_lat;
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
    weight[m]=weight[nrings-1-m]=weight[m]*2*pi/((nrings+1)*nph[m]);
    }

  return make_unique<sharp_standard_geom_info>(nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), weight.data());
  }

unique_ptr<sharp_geom_info> sharp_make_dh_geom_info (size_t nrings, size_t ppring, double phi0,
  ptrdiff_t stride_lon, ptrdiff_t stride_lat)
  {
  const double pi=3.141592653589793238462643383279502884197;

  vector<double> theta(nrings), weight(get_dh_weights(nrings)), phi0_(nrings, phi0);
  vector<size_t> nph(nrings, ppring);
  vector<ptrdiff_t> ofs(nrings);

  for (size_t m=0; m<nrings; ++m)
    {
    theta[m] = m*pi/nrings;
    ofs[m] = ptrdiff_t(m)*stride_lat;
    weight[m]*=2*pi/(nrings*ppring);
Martin Reinecke's avatar
Martin Reinecke committed
337
338
    }

Martin Reinecke's avatar
Martin Reinecke committed
339
  return make_unique<sharp_standard_geom_info>(nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), weight.data());
Martin Reinecke's avatar
Martin Reinecke committed
340
341
  }

Martin Reinecke's avatar
Martin Reinecke committed
342
343
unique_ptr<sharp_geom_info> sharp_make_mw_geom_info (size_t nrings, size_t ppring, double phi0,
  ptrdiff_t stride_lon, ptrdiff_t stride_lat)
Martin Reinecke's avatar
Martin Reinecke committed
344
345
346
  {
  const double pi=3.141592653589793238462643383279502884197;

347
348
  vector<double> theta(nrings), phi0_(nrings, phi0);
  vector<size_t> nph(nrings, ppring);
349
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
350

Martin Reinecke's avatar
Martin Reinecke committed
351
  for (size_t m=0; m<nrings; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
352
353
354
    {
    theta[m]=pi*(2.*m+1.)/(2.*nrings-1.);
    if (theta[m]>pi-1e-15) theta[m]=pi-1e-15;
Martin Reinecke's avatar
Martin Reinecke committed
355
    ofs[m]=ptrdiff_t(m)*stride_lat;
Martin Reinecke's avatar
Martin Reinecke committed
356
357
    }

Martin Reinecke's avatar
Martin Reinecke committed
358
  return make_unique<sharp_standard_geom_info>(nrings, nph.data(), ofs.data(), stride_lon, phi0_.data(), theta.data(), nullptr);
Martin Reinecke's avatar
Martin Reinecke committed
359
  }