sharp_geomhelpers.h 9.99 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.h
 *  SHARP helper function for the creation of grid geometries
 *
24
 *  Copyright (C) 2006-2020 Max-Planck-Society
Martin Reinecke's avatar
Martin Reinecke committed
25
26
27
28
29
30
 *  \author Martin Reinecke
 */

#ifndef SHARP2_GEOMHELPERS_H
#define SHARP2_GEOMHELPERS_H

Martin Reinecke's avatar
Martin Reinecke committed
31
#include <memory>
Martin Reinecke's avatar
Martin Reinecke committed
32
33
#include "libsharp2/sharp.h"

34
35
36
37
38
39
40
41
42
43
44
45
46
class sharp_standard_geom_info: public sharp_geom_info
  {
  private:
    struct Tring
      {
      double theta, phi0, weight, cth, sth;
      ptrdiff_t ofs;
      size_t nph;
      };
    std::vector<Tring> ring;
    std::vector<Tpair> pair_;
    ptrdiff_t stride;
    size_t nphmax_;
Martin Reinecke's avatar
Martin Reinecke committed
47
48
49
    template<typename T> void tclear (T *map) const;
    template<typename T> void tget (bool weighted, size_t iring, const T *map, double *ringtmp) const;
    template<typename T> void tadd (bool weighted, size_t iring, const double *ringtmp, T *map) const;
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

  public:
/*! Creates a geometry information from a set of ring descriptions.
    All arrays passed to this function must have \a nrings elements.
    \param nrings the number of rings in the map
    \param nph the number of pixels in each ring
    \param ofs the index of the first pixel in each ring in the map array
    \param stride the stride between consecutive pixels
    \param phi0 the azimuth (in radians) of the first pixel in each ring
    \param theta the colatitude (in radians) of each ring
    \param wgt the pixel weight to be used for the ring in map2alm
      and adjoint map2alm transforms.
      Pass nullptr to use 1.0 as weight for all rings.
 */
    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);
    virtual size_t nrings() const { return ring.size(); }
    virtual size_t npairs() const { return pair_.size(); }
    virtual size_t nph(size_t iring) const { return ring[iring].nph; }
    virtual size_t nphmax() const { return nphmax_; }
    virtual double theta(size_t iring) const { return ring[iring].theta; }
    virtual double cth(size_t iring) const { return ring[iring].cth; }
    virtual double sth(size_t iring) const { return ring[iring].sth; }
    virtual double phi0(size_t iring) const { return ring[iring].phi0; }
    virtual Tpair pair(size_t ipair) const { return pair_[ipair]; }
Martin Reinecke's avatar
Martin Reinecke committed
75
76
77
    virtual void clear_map(const std::any &map) const;
    virtual void get_ring(bool weighted, size_t iring, const std::any &map, double *ringtmp) const;
    virtual void add_ring(bool weighted, size_t iring, const double *ringtmp, const std::any &map) const;
78
79
  };

Martin Reinecke's avatar
Martin Reinecke committed
80
81
82
83
/*! Creates a geometry information describing a HEALPix map with an
    Nside parameter \a nside. \a weight contains the relative ring
    weights and must have \a 2*nside entries. The rings array contains
    the indices of the rings, with 1 being the first ring at the north
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
84
85
    pole; if nullptr then we take them to be sequential. Pass 4 * nside - 1
    as nrings and nullptr to rings to get the full HEALPix grid.
Martin Reinecke's avatar
Martin Reinecke committed
86
87
88
    \note if \a weight is a null pointer, all weights are assumed to be 1.
    \note if \a rings is a null pointer, take all rings
    \ingroup geominfogroup */
Martin Reinecke's avatar
Martin Reinecke committed
89
90
std::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
91
92
93
94
95
96

/*! Creates a geometry information describing a HEALPix map with an
    Nside parameter \a nside. \a weight contains the relative ring
    weights and must have \a 2*nside entries.
    \note if \a weight is a null pointer, all weights are assumed to be 1.
    \ingroup geominfogroup */
Martin Reinecke's avatar
Martin Reinecke committed
97
std::unique_ptr<sharp_geom_info> sharp_make_weighted_healpix_geom_info (size_t nside, ptrdiff_t stride,
Martin Reinecke's avatar
Martin Reinecke committed
98
  const double *weight);
Martin Reinecke's avatar
Martin Reinecke committed
99
100
101
102

/*! Creates a geometry information describing a HEALPix map with an
    Nside parameter \a nside.
    \ingroup geominfogroup */
Martin Reinecke's avatar
Martin Reinecke committed
103
static inline std::unique_ptr<sharp_geom_info> sharp_make_healpix_geom_info (size_t nside, ptrdiff_t stride)
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
104
  { return sharp_make_weighted_healpix_geom_info (nside, stride, nullptr); }
Martin Reinecke's avatar
Martin Reinecke committed
105
106
107
108
109
110
111
112

/*! Creates a geometry information describing a Gaussian map with \a nrings
    iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
    pixel in each ring is \a phi0 (in radians). The index difference between
    two adjacent pixels in an iso-latitude ring is \a stride_lon, the index
    difference between the two start pixels in consecutive iso-latitude rings
    is \a stride_lat.
    \ingroup geominfogroup */
Martin Reinecke's avatar
Martin Reinecke committed
113
114
std::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
115
116
117
118
119
120
121
122
123
124
125
126
127
128

/*! Creates a geometry information describing an ECP map with \a nrings
    iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
    pixel in each ring is \a phi0 (in radians). The index difference between
    two adjacent pixels in an iso-latitude ring is \a stride_lon, the index
    difference between the two start pixels in consecutive iso-latitude rings
    is \a stride_lat.
    \note The spacing of pixel centers is equidistant in colatitude and
      longitude.
    \note The sphere is pixelized in a way that the colatitude of the first ring
      is \a 0.5*(pi/nrings) and the colatitude of the last ring is
      \a pi-0.5*(pi/nrings). There are no pixel centers at the poles.
    \note This grid corresponds to Fejer's first rule.
    \ingroup geominfogroup */
Martin Reinecke's avatar
Martin Reinecke committed
129
130
std::unique_ptr<sharp_geom_info> sharp_make_fejer1_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
131
132
133

/*! Old name for sharp_make_fejer1_geom_info()
    \ingroup geominfogroup */
Martin Reinecke's avatar
Martin Reinecke committed
134
135
static inline std::unique_ptr<sharp_geom_info> sharp_make_ecp_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
136
  {
Martin Reinecke's avatar
Martin Reinecke committed
137
  return sharp_make_fejer1_geom_info (nrings, nphi, phi0, stride_lon, stride_lat);
Martin Reinecke's avatar
Martin Reinecke committed
138
139
140
141
142
143
144
145
146
147
148
149
150
151
  }

/*! Creates a geometry information describing an ECP map with \a nrings
    iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
    pixel in each ring is \a phi0 (in radians). The index difference between
    two adjacent pixels in an iso-latitude ring is \a stride_lon, the index
    difference between the two start pixels in consecutive iso-latitude rings
    is \a stride_lat.
    \note The spacing of pixel centers is equidistant in colatitude and
      longitude.
    \note The sphere is pixelized in a way that the colatitude of the first ring
      is \a 0 and that of the last ring is \a pi.
    \note This grid corresponds to Clenshaw-Curtis integration.
    \ingroup geominfogroup */
Martin Reinecke's avatar
Martin Reinecke committed
152
153
std::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
154
155
156
157
158
159
160
161
162
163
164
165
166

/*! Creates a geometry information describing an ECP map with \a nrings
    iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
    pixel in each ring is \a phi0 (in radians). The index difference between
    two adjacent pixels in an iso-latitude ring is \a stride_lon, the index
    difference between the two start pixels in consecutive iso-latitude rings
    is \a stride_lat.
    \note The spacing of pixel centers is equidistant in colatitude and
      longitude.
    \note The sphere is pixelized in a way that the colatitude of the first ring
      is \a pi/(nrings+1) and that of the last ring is \a pi-pi/(nrings+1).
    \note This grid corresponds to Fejer's second rule.
    \ingroup geominfogroup */
Martin Reinecke's avatar
Martin Reinecke committed
167
std::unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (size_t nrings, size_t ppring, double phi0,
168
169
170
171
172
173
174
175
176
177
178
179
180
181
  ptrdiff_t stride_lon, ptrdiff_t stride_lat);

/*! Creates a geometry information describing a Driscoll-Healy map with \a nrings
    iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
    pixel in each ring is \a phi0 (in radians). The index difference between
    two adjacent pixels in an iso-latitude ring is \a stride_lon, the index
    difference between the two start pixels in consecutive iso-latitude rings
    is \a stride_lat.
    \note The spacing of pixel centers is equidistant in colatitude and
      longitude.
    \note The sphere is pixelized in a way that the colatitude of the first ring
      is 0 and that of the last ring is \a pi-pi/nrings.
    \ingroup geominfogroup */
std::unique_ptr<sharp_geom_info> sharp_make_dh_geom_info (size_t nrings, size_t ppring, double phi0,
Martin Reinecke's avatar
Martin Reinecke committed
182
  ptrdiff_t stride_lon, ptrdiff_t stride_lat);
Martin Reinecke's avatar
Martin Reinecke committed
183
184
185
186
187
188
189
190
191
192
193
194
195
196

/*! Creates a geometry information describing a map with \a nrings
    iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
    pixel in each ring is \a phi0 (in radians). The index difference between
    two adjacent pixels in an iso-latitude ring is \a stride_lon, the index
    difference between the two start pixels in consecutive iso-latitude rings
    is \a stride_lat.
    \note The spacing of pixel centers is equidistant in colatitude and
      longitude.
    \note The sphere is pixelized in a way that the colatitude of the first ring
      is \a pi/(2*nrings-1) and that of the last ring is \a pi.
    \note This is the grid introduced by McEwen & Wiaux 2011.
    \note This function does \e not define any quadrature weights.
    \ingroup geominfogroup */
Martin Reinecke's avatar
Martin Reinecke committed
197
198
std::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
199
200

#endif