File: wigner.h

package info (click to toggle)
healpix-cxx 3.50.0-3
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 3,008 kB
  • sloc: cpp: 13,882; ansic: 5,930; sh: 4,451; makefile: 216
file content (196 lines) | stat: -rw-r--r-- 5,982 bytes parent folder | download | duplicates (6)
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
31
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
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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
/*
 *  This file is part of libcxxsupport.
 *
 *  libcxxsupport 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.
 *
 *  libcxxsupport 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 libcxxsupport; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/*
 *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
 *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
 *  (DLR).
 */

/*! \file wigner.h
 *  Several C++ classes for calculating Wigner matrices
 *
 *  Copyright (C) 2009-2011 Max-Planck-Society
 *  \author Martin Reinecke and others (see individual classes)
 */

#ifndef PLANCK_WIGNER_H
#define PLANCK_WIGNER_H

#include <cmath>
#include "arr.h"

#include "sse_utils_cxx.h"

/*! Class for calculation of the Wigner matrix at pi/2, using Risbo recursion
    in a way that cannot easily be parallelised, but is fairly efficient on
    scalar machines. */
class wigner_d_halfpi_risbo_scalar
  {
  private:
    double pq;
    arr<double> sqt;
    arr2<double> d;
    int n;

    void do_line0 (double *l1, int j);
    void do_line (const double *l1, double *l2, int j, int k);

  public:
    wigner_d_halfpi_risbo_scalar(int lmax);

    const arr2<double> &recurse ();
  };

/*! Class for calculation of the Wigner matrix at arbitrary angles, using Risbo
    recursion in a way that cannot easily be parallelised, but is fairly
    efficient on scalar machines. */
class wigner_d_risbo_scalar
  {
  private:
    double p,q;
    arr<double> sqt;
    arr2<double> d;
    int n;

    void do_line0 (double *l1, int j);
    void do_line (const double *l1, double *l2, int j, int k);

  public:
    wigner_d_risbo_scalar(int lmax, double ang);

    const arr2<double> &recurse ();
  };

/*! Class for calculation of the Wigner matrix at pi/2, using Risbo recursion
    in a way that can be OpenMP-parallelised. This approach uses more memory
    and is slightly slower than wigner_d_halfpi_risbo_scalar. */
class wigner_d_halfpi_risbo_openmp
  {
  private:
    double pq;
    arr<double> sqt;
    arr2<double> d,dd;
    int n;

  public:
    wigner_d_halfpi_risbo_openmp(int lmax);

    const arr2<double> &recurse ();
  };

/*! Class for calculation of the Wigner matrix at arbitrary angles, using Risbo
    recursion in a way that can be OpenMP-parallelised. This approach uses more
    memory and is slightly slower than wigner_d_risbo_scalar. */
class wigner_d_risbo_openmp
  {
  private:
    double p,q;
    arr<double> sqt;
    arr2<double> d, dd;
    int n;

  public:
    wigner_d_risbo_openmp(int lmax, double ang);

    const arr2<double> &recurse ();
  };

/*! Class for calculation of the Wigner matrix elements by l-recursion.
    For details, see Prezeau & Reinecke 2010, http://arxiv.org/pdf/1002.1050 */
class wignergen_scalar
  {
  protected:
    typedef double dbl3[3];

    // members set in the constructor and kept fixed afterwards
    double fsmall, fbig, eps;
    int lmax;
    arr<long double> logsum, lc05, ls05;
    arr<double> flm1, flm2, cf, costh, xl;
    arr<bool> thetaflip;

    // members depending on m and m'
    int m1, m2, am1, am2, mlo, mhi, cosPow, sinPow;
    long double prefactor;
    arr<dbl3> fx;
    bool preMinus;

    // members depending on theta
    arr<double> result;

    enum { large_exponent2=90, minscale=-4, maxscale=14 };

  public:
    /*! Constructs an object that can compute Wigner matrix elements up
        to a maximum \a l value of \a lmax_, at the colatitudes provided
        in \a thetas. The generator will be allowed to regard values with
        absolute magnitudes smaller than \a epsilon as zero; a typical value
        is 1e-30. */
    wignergen_scalar (int lmax_, const arr<double> &thetas, double epsilon);

    /*! Prepares the object to produce Wigner matrix elements with \a m=m1_
        and \a m'=m2_ in subsequent calls to calc(). This operation is not cheap
        so it is recommended to use calc() for many different colatitudes after
        every call to prepare(), if possible. */
    void prepare (int m1_, int m2_);

    /*! Computes the Wigner matrix elements for the values of \a m and \a m'
        set by the preceding call to prepare(), for all \a l up to \a lmax
        (set in the constructor), and for the \a nth colatitude passed to the
        constructor. On return, \a firstl contains the index of the first
        matrix element larger than \a epsilon; all values with smaller indices
        in the result array are undefined. */
    const arr<double> &calc (int nth, int &firstl);
    void calc (int nth, int &firstl, arr<double> &resx) const;
  };

class wignergen: public wignergen_scalar
  {
#ifdef __SSE2__
  private:
    arr_align<V2df,16> result2;

  public:
    wignergen (int lmax_, const arr<double> &thetas, double epsilon)
      : wignergen_scalar (lmax_,thetas,epsilon), result2(lmax_+1) {}

    using wignergen_scalar::calc;
    const arr_align<V2df,16> &calc (int nth1, int nth2, int &firstl);
    void calc (int nth1, int nth2, int &firstl, arr_align<V2df,16> &resx) const;
#else
  public:
    wignergen (int lmax_, const arr<double> &thetas, double epsilon)
      : wignergen_scalar (lmax_,thetas,epsilon) {}
#endif
  };

class wigner_estimator
  {
  private:
    int lmax, m1, m2, mbig;
    double xlmax, epsPow, cosm1m2;

  public:
    wigner_estimator (int lmax_, double epsPow_);

    void prepare_m (int m1_, int m2_);
    bool canSkip (double theta) const;
  };

#endif