sharp_geomhelpers.cc 11.8 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
/*
 *  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 */

Martin Reinecke's avatar
Martin Reinecke committed
21
/*! \file sharp_geomhelpers.cc
Martin Reinecke's avatar
Martin Reinecke committed
22 23
 *  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
#include <algorithm>
Martin Reinecke's avatar
Martin Reinecke committed
29 30
#include <cmath>
#include <vector>
31 32 33 34 35 36
#include "mr_util/sharp/sharp_geomhelpers.h"
#include "mr_util/math/gl_integrator.h"
#include "mr_util/math/constants.h"
#include "mr_util/math/fft1d.h"
#include "mr_util/infra/error_handling.h"
#include "mr_util/math/math_utils.h"
Martin Reinecke's avatar
Martin Reinecke committed
37

38 39 40
namespace mr {

namespace detail_sharp {
41

42 43
using namespace std;

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 89 90 91 92 93
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
94

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

Martin Reinecke's avatar
Martin Reinecke committed
107
void sharp_standard_geom_info::clear_map (const any &map) const
108
  {
Martin Reinecke's avatar
Martin Reinecke committed
109 110
  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
111
  else MR_fail("bad map data type");
112 113
  }

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

Martin Reinecke's avatar
Martin Reinecke committed
145 146
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
147
  {
Martin Reinecke's avatar
Martin Reinecke committed
148 149
  size_t npix=nside*nside*12;
  size_t ncap=2*nside*(nside-1);
Martin Reinecke's avatar
Martin Reinecke committed
150

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

Martin Reinecke's avatar
Martin Reinecke committed
192
  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
193 194
  }

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

Martin Reinecke's avatar
Martin Reinecke committed
201 202
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
203
  {
204 205
  vector<size_t> nph(nrings, nphi);
  vector<double> phi0_(nrings, phi0);
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
206
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
207 208 209 210

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

Martin Reinecke's avatar
Martin Reinecke committed
218
  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
219 220 221
  }

/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
Martin Reinecke's avatar
Martin Reinecke committed
222 223
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
224
  {
225 226
  vector<double> theta(nrings), weight(nrings), phi0_(nrings, phi0);
  vector<size_t> nph(nrings, ppring);
227
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
228 229

  weight[0]=2.;
Martin Reinecke's avatar
Martin Reinecke committed
230
  for (size_t k=1; k<=(nrings-1)/2; ++k)
Martin Reinecke's avatar
Martin Reinecke committed
231 232 233 234 235
    {
    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
Martin Reinecke committed
236 237 238 239
  {
  pocketfft_r<double> plan(nrings);
  plan.exec(weight.data(), 1., false);
  }
Martin Reinecke's avatar
Martin Reinecke committed
240

Martin Reinecke's avatar
Martin Reinecke committed
241
  for (size_t m=0; m<(nrings+1)/2; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
242 243 244
    {
    theta[m]=pi*(m+0.5)/nrings;
    theta[nrings-1-m]=pi-theta[m];
Martin Reinecke's avatar
Martin Reinecke committed
245 246
    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
247 248 249
    weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(nrings*nph[m]);
    }

Martin Reinecke's avatar
Martin Reinecke committed
250
  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
251 252 253
  }

/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
Martin Reinecke's avatar
Martin Reinecke committed
254 255
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
256
  {
257 258
  vector<double> theta(nrings), weight(nrings,0.), phi0_(nrings, phi0);
  vector<size_t> nph(nrings, ppring);
259
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
260

Martin Reinecke's avatar
Martin Reinecke committed
261
  size_t n=nrings-1;
Martin Reinecke's avatar
Martin Reinecke committed
262 263
  double dw=-1./(n*n-1.+(n&1));
  weight[0]=2.+dw;
Martin Reinecke's avatar
Martin Reinecke committed
264
  for (size_t k=1; k<=(n/2-1); ++k)
Martin Reinecke's avatar
Martin Reinecke committed
265 266
    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
Martin Reinecke committed
267 268 269 270
  {
  pocketfft_r<double> plan(n);
  plan.exec(weight.data(), 1., false);
  }
Martin Reinecke's avatar
Martin Reinecke committed
271 272
  weight[n]=weight[0];

Martin Reinecke's avatar
Martin Reinecke committed
273
  for (size_t m=0; m<(nrings+1)/2; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
274 275 276 277
    {
    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
278 279
    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
280 281 282
    weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
    }

Martin Reinecke's avatar
Martin Reinecke committed
283
  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
284 285
  }

286 287 288 289 290 291 292 293
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.;
Martin Reinecke's avatar
Martin Reinecke committed
294 295 296 297
  {
  pocketfft_r<double> plan(nrings);
  plan.exec(weight.data(), 1., false);
  }
298 299 300
  return weight;
  }

Martin Reinecke's avatar
Martin Reinecke committed
301
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
Martin Reinecke's avatar
Martin Reinecke committed
302 303
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
304
  {
305 306
  vector<double> theta(nrings), weight(get_dh_weights(nrings+1)), phi0_(nrings, phi0);
  vector<size_t> nph(nrings, ppring);
307
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
308

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

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

Martin Reinecke's avatar
Martin Reinecke committed
338
  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
339 340
  }

Martin Reinecke's avatar
Martin Reinecke committed
341 342
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
343
  {
344 345
  vector<double> theta(nrings), phi0_(nrings, phi0);
  vector<size_t> nph(nrings, ppring);
346
  vector<ptrdiff_t> ofs(nrings);
Martin Reinecke's avatar
Martin Reinecke committed
347

Martin Reinecke's avatar
Martin Reinecke committed
348
  for (size_t m=0; m<nrings; ++m)
Martin Reinecke's avatar
Martin Reinecke committed
349 350 351
    {
    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
352
    ofs[m]=ptrdiff_t(m)*stride_lat;
Martin Reinecke's avatar
Martin Reinecke committed
353 354
    }

Martin Reinecke's avatar
Martin Reinecke committed
355
  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
356
  }
357 358

}}