sharp_almhelpers.cc 4.78 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
24
25
26
27
/*
 *  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_almhelpers.c
 *  Spherical transform library
 *
 *  Copyright (C) 2008-2019 Max-Planck-Society
 *  \author Martin Reinecke
 */

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

Martin Reinecke's avatar
Martin Reinecke committed
31
using namespace std;
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
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
57
template<typename T> void sharp_standard_alm_info::tclear (T *alm) const
58
59
60
  {
  for (size_t mi=0;mi<mval_.size();++mi)
    for (size_t l=mval_[mi];l<=lmax_;++l)
61
      reinterpret_cast<T *>(alm)[mvstart[mi]+ptrdiff_t(l)*stride]=0.;
62
  }
Martin Reinecke's avatar
Martin Reinecke committed
63
void sharp_standard_alm_info::clear_alm(const any &alm) const
64
  {
Martin Reinecke's avatar
Martin Reinecke committed
65
66
67
  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");
68
  }
Martin Reinecke's avatar
Martin Reinecke committed
69
template<typename T> void sharp_standard_alm_info::tget(size_t mi, const T *alm, dcmplx *almtmp, size_t nalm) const
70
71
  {
  for (auto l=mval_[mi]; l<=lmax_; ++l)
72
    almtmp[nalm*l] = alm[mvstart[mi]+ptrdiff_t(l)*stride];
73
  }
Martin Reinecke's avatar
Martin Reinecke committed
74
void sharp_standard_alm_info::get_alm(size_t mi, const any &alm, dcmplx *almtmp, size_t nalm) const
75
  {
Martin Reinecke's avatar
Martin Reinecke committed
76
77
78
79
80
  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");
81
  }
Martin Reinecke's avatar
Martin Reinecke committed
82
template<typename T> void sharp_standard_alm_info::tadd(size_t mi, const dcmplx *almtmp, T *alm, size_t nalm) const
83
84
  {
  for (auto l=mval_[mi]; l<=lmax_; ++l)
85
    alm[mvstart[mi]+ptrdiff_t(l)*stride] += T(almtmp[nalm*l]);
86
  }
Martin Reinecke's avatar
Martin Reinecke committed
87
void sharp_standard_alm_info::add_alm(size_t mi, const dcmplx *almtmp, const any &alm, size_t nalm) const
88
  {
Martin Reinecke's avatar
Martin Reinecke committed
89
90
91
  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");
92
  }
Martin Reinecke's avatar
Martin Reinecke committed
93

94
ptrdiff_t sharp_standard_alm_info::index (size_t l, size_t mi)
Martin Reinecke's avatar
Martin Reinecke committed
95
  { return mvstart[mi]+stride*ptrdiff_t(l); }
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
/* 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
112

113
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
114
  {
Martin Reinecke's avatar
Martin Reinecke committed
115
  vector<ptrdiff_t> mvstart(mmax+1);
116
117
118
  size_t tval = 2*lmax+1;
  for (size_t m=0; m<=mmax; ++m)
    mvstart[m] = stride*ptrdiff_t((m*(tval-m))>>1);
119
  return make_unique<sharp_standard_alm_info>(lmax, mmax, stride, mvstart.data());
Martin Reinecke's avatar
Martin Reinecke committed
120
121
  }

122
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
123
  {
Martin Reinecke's avatar
Martin Reinecke committed
124
  vector<ptrdiff_t> mvstart(mmax+1);
125
126
  for (size_t m=0; m<=mmax; ++m)
    mvstart[m] = stride*ptrdiff_t(m*(lmax+1));
127
  return make_unique<sharp_standard_alm_info>(lmax, mmax, stride, mvstart.data());
Martin Reinecke's avatar
Martin Reinecke committed
128
  }