maths.hpp

Go to the documentation of this file.

// Triumvirate: Three-Point Clustering Measurements in LSS
//
// Copyright (C) 2023 Mike S Wang & Naonori S Sugiyama [GPL-3.0-or-later]
//
// This file is part of the Triumvirate program. See the COPYRIGHT
// and LICENCE files at the top-level directory of this distribution
// for details of copyright and licensing.
//
// This program 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 3 of the License, or
// (at your option) any later version.
//
// This program 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 this program. If not, see <https://www.gnu.org/licenses/>.

 
#ifndef TRIUMVIRATE_INCLUDE_MATHS_HPP_INCLUDED_
#define TRIUMVIRATE_INCLUDE_MATHS_HPP_INCLUDED_
 
#include <gsl/gsl_interp.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_coupling.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_spline.h>
 
#include <cmath>
#include <complex>
#include <vector>
 
#include "monitor.hpp"

#ifdef __GNUC__
#define PURE __attribute__((pure))
#else
#define PURE
#endif
 
namespace trv {
namespace maths {
 
// ***********************************************************************
// Complex numbers
// ***********************************************************************
 
extern const std::complex<double> M_I;  

std::complex<double> eval_complex_in_polar(double r, double theta);
 
 
// ***********************************************************************
// Vectors
// ***********************************************************************

double get_vec3d_magnitude(std::vector<double> vec);

PURE double get_vec3d_magnitude(double* vec);
 
 
// ***********************************************************************
// Gamma function
// ***********************************************************************
 
// -----------------------------------------------------------------------
// Lanzcos approximation
// -----------------------------------------------------------------------

std::complex<double> eval_lanczos_approx_series(std::complex<double> z);

std::complex<double> eval_gamma_lanczos(std::complex<double> z);
 
 
// -----------------------------------------------------------------------
// Component evaluation
// -----------------------------------------------------------------------

void get_lngamma_parts(double x, double y, double& lnr, double& theta);

std::complex<double> eval_gamma_lnratio(double mu, std::complex<double> nu);
 
 
// ***********************************************************************
// Spherical harmonics
// ***********************************************************************

extern const double eps_coupling;

double wigner_3j(int j1, int j2, int j3, int m1, int m2, int m3);

class SphericalHarmonicCalculator {
 public:
  static std::complex<double> calc_reduced_spherical_harmonic(
    const int ell, const int m, double pos[3]
  );

  static void store_reduced_spherical_harmonic_in_fourier_space(
    const int ell, const int m,
    const double boxsize[3], const int ngrid[3],
    std::vector< std::complex<double> >& ylm_out
  );

  static void store_reduced_spherical_harmonic_in_config_space(
    const int ell, const int m,
    const double boxsize[3], const int ngrid[3],
    std::vector< std::complex<double> >& ylm_out
  );
};
 
 
// ***********************************************************************
// Spherical Bessel function
// ***********************************************************************

class SphericalBesselCalculator {
 public:
  int order;  

  explicit SphericalBesselCalculator(const int ell);

  SphericalBesselCalculator(const SphericalBesselCalculator& other);

  ~SphericalBesselCalculator();

  double eval(double x);
 
 private:
  // CAVEAT: This calculator is designed for the range of @f$ x = kr @f$
  // such that max(kr) > 2048π.  For @f$ x <= \max\{1000, \ell^2\} @f$,
  // cubic spline interpolation is used with Δ(kr) = 0.05;
  // for @f$ x > \max\{1000, \ell^2\} @f$, direct evaluation via
  // asymptotic expansion is used.
  double split = 1000.;    
  double step = 0.05;      
 
  gsl_interp_accel* accel;  
  gsl_spline* spline;       
};
 
}  // namespace trv::maths
}  // namespace trv
 
#endif  // !TRIUMVIRATE_INCLUDE_MATHS_HPP_INCLUDED_