sharp.h 7.25 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
28
29
30
/*
 *  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.h
 *  Portable interface for the spherical transform library.
 *
 *  Copyright (C) 2012-2019 Max-Planck-Society
 *  \author Martin Reinecke \author Dag Sverre Seljebotn
 */

#ifndef SHARP_SHARP_H
#define SHARP_SHARP_H

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
31
32
#include <cstddef>
#include <vector>
Martin Reinecke's avatar
Martin Reinecke committed
33
#include <memory>
Martin Reinecke's avatar
Martin Reinecke committed
34
35
36

/*! \internal
    Type holding all required information about a map geometry. */
Martin Reinecke's avatar
Martin Reinecke committed
37
struct sharp_geom_info
Martin Reinecke's avatar
Martin Reinecke committed
38
  {
39
40
41
42
43
44
45
46
47
48
49
  struct Tring
    {
    double theta, phi0, weight, cth, sth;
    ptrdiff_t ofs;
    int nph, stride;
    };
  struct Tpair
    {
    Tring r1,r2;
    };
  std::vector<Tpair> pair;
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
50
  int nphmax;
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
/*! 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_geom_info(int nrings, const int *nph, const ptrdiff_t *ofs,
    const int *stride, const double *phi0, const double *theta,
    const double *wgt);
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
66
67
  void clear_map(double *map) const;
  void clear_map(float *map) const;
Martin Reinecke's avatar
Martin Reinecke committed
68
  };
Martin Reinecke's avatar
Martin Reinecke committed
69
70
71
72
73
74

/*! \defgroup almgroup Helpers for dealing with a_lm */
/*! \{ */

/*! \internal
    Helper type for index calculation in a_lm arrays. */
Martin Reinecke's avatar
Martin Reinecke committed
75
struct sharp_alm_info
Martin Reinecke's avatar
Martin Reinecke committed
76
77
78
79
80
81
  {
  /*! Maximum \a l index of the array */
  int lmax;
  /*! Number of different \a m values in this object */
  int nm;
  /*! Array with \a nm entries containing the individual m values */
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
82
  std::vector<int> mval;
Martin Reinecke's avatar
Martin Reinecke committed
83
84
  /*! Array with \a nm entries containing the (hypothetical) indices of
      the coefficients with quantum numbers 0,\a mval[i] */
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
85
  std::vector<ptrdiff_t> mvstart;
Martin Reinecke's avatar
Martin Reinecke committed
86
87
88
89
90
91
92
93
94
95
96
  /*! Stride between a_lm and a_(l+1),m */
  ptrdiff_t stride;

/*! Creates an a_lm data structure from the following parameters:
    \param lmax maximum \a l quantum number (>=0)
    \param mmax maximum \a m quantum number (0<= \a mmax <= \a lmax)
    \param stride the stride between entries with identical \a m, and \a l
      differing by 1.
    \param mstart the index of the (hypothetical) coefficient with the
      quantum numbers 0,\a m. Must have \a mmax+1 entries.
 */
Martin Reinecke's avatar
Martin Reinecke committed
97
98
  sharp_alm_info(int lmax, int mmax, int stride, const ptrdiff_t *mstart);

Martin Reinecke's avatar
Martin Reinecke committed
99
100
101
102
103
104
105
106
107
/*! Creates an a_lm data structure which from the following parameters:
    \param lmax maximum \a l quantum number (\a >=0)
    \param nm number of different \a m (\a 0<=nm<=lmax+1)
    \param stride the stride between entries with identical \a m, and \a l
      differing by 1.
    \param mval array with \a nm entries containing the individual m values
    \param mvstart array with \a nm entries containing the (hypothetical)
      indices of the coefficients with the quantum numbers 0,\a mval[i]
 */
Martin Reinecke's avatar
Martin Reinecke committed
108
  sharp_alm_info (int lmax, int nm, int stride, const int *mval,
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
109
    const ptrdiff_t *mvstart);
Martin Reinecke's avatar
Martin Reinecke committed
110
111
112
113
/*! Returns the index of the coefficient with quantum numbers \a l,
    \a mval[mi].
    \note for a \a sharp_alm_info generated by sharp_make_alm_info() this is
    the index for the coefficient with the quantum numbers \a l, \a mi. */
Martin Reinecke's avatar
Martin Reinecke committed
114
115
116
  ptrdiff_t index (int l, int mi);
  };

Martin Reinecke's avatar
Martin Reinecke committed
117
118
119
120
121
122
123
124
125
126
127
/*! \} */

/*! \defgroup geominfogroup Functions for dealing with geometry information */
/*! \{ */

/*! \} */

/*! \defgroup lowlevelgroup Low-level libsharp2 SHT interface */
/*! \{ */

/*! Enumeration of SHARP job types. */
Martin Reinecke's avatar
Martin Reinecke committed
128
enum sharp_jobtype { SHARP_YtW=0,               /*!< analysis */
Martin Reinecke's avatar
Martin Reinecke committed
129
130
131
132
133
134
               SHARP_MAP2ALM=SHARP_YtW,   /*!< analysis */
               SHARP_Y=1,                 /*!< synthesis */
               SHARP_ALM2MAP=SHARP_Y,     /*!< synthesis */
               SHARP_Yt=2,                /*!< adjoint synthesis */
               SHARP_WY=3,                /*!< adjoint analysis */
               SHARP_ALM2MAP_DERIV1=4     /*!< synthesis of first derivatives */
Martin Reinecke's avatar
Martin Reinecke committed
135
             };
Martin Reinecke's avatar
Martin Reinecke committed
136
137

/*! Job flags */
Martin Reinecke's avatar
Martin Reinecke committed
138
enum sharp_jobflags { SHARP_DP              = 1<<4,
Martin Reinecke's avatar
Martin Reinecke committed
139
140
141
142
143
               /*!< map and a_lm are in double precision */
               SHARP_ADD             = 1<<5,
               /*!< results are added to the output arrays, instead of
                    overwriting them */
               SHARP_USE_WEIGHTS     = 1<<20,    /* internal use only */
Martin Reinecke's avatar
Martin Reinecke committed
144
             };
Martin Reinecke's avatar
Martin Reinecke committed
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166

/*! Performs a libsharp2 SHT job. The interface deliberately does not use
  the C99 "complex" data type, in order to be callable from C89 and C++.
  \param type the type of SHT
  \param spin the spin of the quantities to be transformed
  \param alm contains pointers to the a_lm coefficients. If \a spin==0,
    alm[0] points to the a_lm of the SHT. If \a spin>0, alm[0] and alm[1]
    point to the two a_lm sets of the SHT. The exact data type of \a alm
    depends on whether the SHARP_DP flag is set.
  \param map contains pointers to the maps. If \a spin==0,
    map[0] points to the map of the SHT. If \a spin>0, or \a type is
    SHARP_ALM2MAP_DERIV1, map[0] and map[1] point to the two maps of the SHT.
    The exact data type of \a map depends on whether the SHARP_DP flag is set.
  \param geom_info A \c sharp_geom_info object compatible with the provided
    \a map arrays.
  \param alm_info A \c sharp_alm_info object compatible with the provided
    \a alm arrays. All \c m values from 0 to some \c mmax<=lmax must be present
    exactly once.
  \param flags See sharp_jobflags. In particular, if SHARP_DP is set, then
    \a alm is expected to have the type "complex double **" and \a map is
    expected to have the type "double **"; otherwise, the expected
    types are "complex float **" and "float **", respectively.
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
167
  \param time If not nullptr, the wall clock time required for this SHT
Martin Reinecke's avatar
Martin Reinecke committed
168
    (in seconds) will be written here.
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
169
  \param opcnt If not nullptr, a conservative estimate of the total floating point
Martin Reinecke's avatar
Martin Reinecke committed
170
171
    operation count for this SHT will be written here. */
void sharp_execute (sharp_jobtype type, int spin, void *alm, void *map,
Martin Reinecke's avatar
Martin Reinecke committed
172
  const sharp_geom_info &geom_info, const sharp_alm_info &alm_info,
Martin Reinecke's avatar
Martin Reinecke committed
173
174
175
176
177
178
179
180
181
182
183
184
  int flags, double *time, unsigned long long *opcnt);

void sharp_set_chunksize_min(int new_chunksize_min);
void sharp_set_nchunks_max(int new_nchunks_max);

/*! \} */

int sharp_get_mlim (int lmax, int spin, double sth, double cth);
int sharp_veclen(void);
const char *sharp_architecture(void);

#endif