sharp_almhelpers.cc 4.84 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_almhelpers.cc
Martin Reinecke's avatar
Martin Reinecke committed
22 23
 *  Spherical transform library
 *
24
 *  Copyright (C) 2008-2020 Max-Planck-Society
Martin Reinecke's avatar
Martin Reinecke committed
25 26 27
 *  \author Martin Reinecke
 */

28 29
#include "mr_util/infra/error_handling.h"
#include "mr_util/sharp/sharp_almhelpers.h"
Martin Reinecke's avatar
Martin Reinecke committed
30

31 32 33 34
namespace mr {

namespace detail_sharp {

35 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
using dcmplx = complex<double>;
using fcmplx = complex<float>;

sharp_standard_alm_info::sharp_standard_alm_info (size_t lmax__, size_t nm_, ptrdiff_t stride_,
  const size_t *mval__, const ptrdiff_t *mstart)
  : lmax_(lmax__), mval_(nm_), mvstart(nm_), stride(stride_)
  {
  for (size_t mi=0; mi<nm_; ++mi)
    {
    mval_[mi] = mval__[mi];
    mvstart[mi] = mstart[mi];
    }
  }

sharp_standard_alm_info::sharp_standard_alm_info (size_t lmax__, size_t mmax_, ptrdiff_t stride_,
  const ptrdiff_t *mstart)
  : lmax_(lmax__), mval_(mmax_+1), mvstart(mmax_+1), stride(stride_)
  {
  for (size_t i=0; i<=mmax_; ++i)
    {
    mval_[i]=i;
    mvstart[i] = mstart[i];
    }
  }

Martin Reinecke's avatar
Martin Reinecke committed
62
template<typename T> void sharp_standard_alm_info::tclear (T *alm) const
63 64 65
  {
  for (size_t mi=0;mi<mval_.size();++mi)
    for (size_t l=mval_[mi];l<=lmax_;++l)
66
      reinterpret_cast<T *>(alm)[mvstart[mi]+ptrdiff_t(l)*stride]=0.;
67
  }
Martin Reinecke's avatar
Martin Reinecke committed
68
void sharp_standard_alm_info::clear_alm(const any &alm) const
69
  {
Martin Reinecke's avatar
Martin Reinecke committed
70 71 72
  if (alm.type()==typeid(dcmplx *)) tclear(any_cast<dcmplx *>(alm));
  else if (alm.type()==typeid(fcmplx *)) tclear(any_cast<fcmplx *>(alm));
  else MR_fail("bad a_lm data type");
73
  }
Martin Reinecke's avatar
Martin Reinecke committed
74
template<typename T> void sharp_standard_alm_info::tget(size_t mi, const T *alm, dcmplx *almtmp, size_t nalm) const
75 76
  {
  for (auto l=mval_[mi]; l<=lmax_; ++l)
77
    almtmp[nalm*l] = alm[mvstart[mi]+ptrdiff_t(l)*stride];
78
  }
Martin Reinecke's avatar
Martin Reinecke committed
79
void sharp_standard_alm_info::get_alm(size_t mi, const any &alm, dcmplx *almtmp, size_t nalm) const
80
  {
Martin Reinecke's avatar
Martin Reinecke committed
81 82 83 84 85
  if (alm.type()==typeid(dcmplx *)) tget(mi, any_cast<dcmplx *>(alm), almtmp, nalm);
  else if (alm.type()==typeid(const dcmplx *)) tget(mi, any_cast<const dcmplx *>(alm), almtmp, nalm);
  else if (alm.type()==typeid(fcmplx *)) tget(mi, any_cast<fcmplx *>(alm), almtmp, nalm);
  else if (alm.type()==typeid(const fcmplx *)) tget(mi, any_cast<const fcmplx *>(alm), almtmp, nalm);
  else MR_fail("bad a_lm data type");
86
  }
Martin Reinecke's avatar
Martin Reinecke committed
87
template<typename T> void sharp_standard_alm_info::tadd(size_t mi, const dcmplx *almtmp, T *alm, size_t nalm) const
88 89
  {
  for (auto l=mval_[mi]; l<=lmax_; ++l)
90
    alm[mvstart[mi]+ptrdiff_t(l)*stride] += T(almtmp[nalm*l]);
91
  }
Martin Reinecke's avatar
Martin Reinecke committed
92
void sharp_standard_alm_info::add_alm(size_t mi, const dcmplx *almtmp, const any &alm, size_t nalm) const
93
  {
Martin Reinecke's avatar
Martin Reinecke committed
94 95 96
  if (alm.type()==typeid(dcmplx *)) tadd(mi, almtmp, any_cast<dcmplx *>(alm), nalm);
  else if (alm.type()==typeid(fcmplx *)) tadd(mi, almtmp, any_cast<fcmplx *>(alm), nalm);
  else MR_fail("bad a_lm data type");
97
  }
Martin Reinecke's avatar
Martin Reinecke committed
98

99
ptrdiff_t sharp_standard_alm_info::index (size_t l, size_t mi)
Martin Reinecke's avatar
Martin Reinecke committed
100
  { return mvstart[mi]+stride*ptrdiff_t(l); }
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
/* This currently requires all m values from 0 to nm-1 to be present.
   It might be worthwhile to relax this criterion such that holes in the m
   distribution are permissible. */
size_t sharp_standard_alm_info::mmax() const
  {
  //FIXME: if gaps are allowed, we have to search the maximum m in the array
  auto nm_=mval_.size();
  vector<bool> mcheck(nm_,false);
  for (auto m_cur : mval_)
    {
    MR_assert(m_cur<nm_, "not all m values are present");
    MR_assert(mcheck[m_cur]==false, "duplicate m value");
    mcheck[m_cur]=true;
    }
  return nm_-1;
  }
Martin Reinecke's avatar
Martin Reinecke committed
117

118
unique_ptr<sharp_standard_alm_info> sharp_make_triangular_alm_info (size_t lmax, size_t mmax, ptrdiff_t stride)
Martin Reinecke's avatar
Martin Reinecke committed
119
  {
Martin Reinecke's avatar
Martin Reinecke committed
120
  vector<ptrdiff_t> mvstart(mmax+1);
121 122 123
  size_t tval = 2*lmax+1;
  for (size_t m=0; m<=mmax; ++m)
    mvstart[m] = stride*ptrdiff_t((m*(tval-m))>>1);
124
  return make_unique<sharp_standard_alm_info>(lmax, mmax, stride, mvstart.data());
Martin Reinecke's avatar
Martin Reinecke committed
125 126
  }

127
unique_ptr<sharp_standard_alm_info> sharp_make_rectangular_alm_info (size_t lmax, size_t mmax, ptrdiff_t stride)
Martin Reinecke's avatar
Martin Reinecke committed
128
  {
Martin Reinecke's avatar
Martin Reinecke committed
129
  vector<ptrdiff_t> mvstart(mmax+1);
130 131
  for (size_t m=0; m<=mmax; ++m)
    mvstart[m] = stride*ptrdiff_t(m*(lmax+1));
132
  return make_unique<sharp_standard_alm_info>(lmax, mmax, stride, mvstart.data());
Martin Reinecke's avatar
Martin Reinecke committed
133
  }
134 135

}}