Commit 0cc000d4 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 13/1/16 update spherical part and doc for alpha =/= 2

parent 51122e08
......@@ -29,20 +29,26 @@ using the spherical harmonic functions $Y_{l,m}$ and the genralized Laguerre fun
$$\begin{eqnarray}
Y_{l,m}(\Omega) &=& \lambda_{lm}(\theta) e^{i m\phi} = \sqrt{\frac{2 l +1}{4 \pi}\frac{(l-m)!}{(l+m)!}}P_l^m({\cos \theta}) e^{i m\phi} \\
{\cal K}_n(r,\tau) &=& \sqrt{\frac{n!}{(n+\alpha)!}} \frac{e^{-r/2\tau} L_n^{(\alpha)}(r/\tau)}{\tau^{3/2}}
{\cal K}_n(r,\tau) &=& \tau^{-3/2} \sqrt{\frac{n!}{(n+\alpha)!}} e^{-r/2\tau} \left(\frac{r}{\tau}\right)^{\frac{\alpha}{2}-1}
L_n^{(\alpha)}(r/\tau) = \tau^{-3/2} {\cal K}_n(r/\tau,1)
\end{eqnarray}$$
with $ L_n^{(\alpha)}(r)$ the generalized Laguerre polynomials (notice that by default $ \alpha = 2$). The $\tau$-scaling parameter is a free parameter but can be nicely set to $\tau = R_{max}/r_{N-1}$ with $R_{max}$ the limit in radial direction of the function $f(r,\Omega)$, and $r_{N-1}$ the largest node of the Gauss-Laguerre quadrature (see below).
with $ L_n^{(\alpha)}(r)$ the generalized Laguerre polynomials (notice that by default $ \alpha = 2$ but one is free to use real numbers). The $\tau$-scaling parameter is a free parameter but can be nicely set to $\tau = R_{max}/r_{N-1}$ with $R_{max}$ the limit in radial direction of the function $f(r,\Omega)$, and $r_{N-1}$ the largest node of the Gauss-Laguerre quadrature (see below).
For band-limited functions the discret sums are truncated with $ l <L$ and $ n < N$ (notice the L band-limit definition is not the same used for instance in `libsharp`). Focusing on the radial part only, one gets
$$\begin{eqnarray}
f(R_i) &=& f(\tau r_i) = \sum_{n=0}^{N-1} f_n\ {\cal K}_n(r_i,\tau) \\
&=& \tau^{-3/2} \sum_{n=0}^{N-1} f_n \sqrt{\frac{n!}{(n+\alpha)!}} e^{-r_i/2} L^{(\alpha)}_n(r_i) r_i^{\frac{\alpha}{2}-1}
\end{eqnarray}$$
The function is evaluated at $R_i$ values which are scaled values of the roots of the generalized Laguerre polynomial of order $N$ noted $r_i$ with $i$ in $[0,N-1]$. To get the $f_n$ coefficients one uses Gauss-Laguerre quadrature to compute
$$
f(r) = \sum_{n=0}^{N-1} f_n\ {\cal K}_n(r)
f_n = \int_{R^+} dR\ R^2\ f(R)\ {\cal K}_n(R,\tau)
$$
and the r-values are the roots of the generalized Laguerre polynomial of order $N$ noted $r_i$ with $i$ in $[0,N-1]$. From the orthogonality of the Laguerre functions, it yields the Gauss-Laguerre quadrature
It yields
$$
f_n = \sum_{i=0}^{N-1} w_i\ f(r_i)\ {\cal L}^{(\alpha)}_n (r_i) \quad\quad n\in \{ 0,\dots, N-1\}
f_n = \tau^{3/2} \sqrt{\frac{n!}{(n+\alpha)!}} \sum_{i=0}^{N-1} w_i\ f(\tau r_i)\ e^{-r_i/2} L^{(\alpha)}_n(r_i) r_i^{1-\frac{\alpha}{2}}\quad\quad n\in \{ 0,\dots, N-1\}
$$
with the weight
with the weight *(scaled with $e^{r_i}$ for numerical purpose as it is included in the computation by recurence of Laguerre polynomials)*
$$
w_i = \frac{(N + \alpha)!}{N!(N+1)^2}\frac{r_i e^{r_i}}{\left[L^{(\alpha)}_{N+1}(r_i)\right]^2}
$$
......@@ -58,9 +64,9 @@ The analysis process starts by gathering the 3D pixels real values $f_{ijk} \equ
+ for the $k$th-radial shell, one performs a Spherical Harmonic decomposition using for instance the `map2alm` routine of `libsharp` to get the set of complex coefficients $a_{lmk}=a_{lm}(r_k)$ with $l\in\{0,\dots,L-1\}$ and $m \in\{0,\dots,l\}$ (the notation remind the underlying 2D Spherical Harmonic transform and the usual way to note the resulting coefficients $a_{lm}$);
+ for each $(l,m)$, one uses the set of values $a_{lmk}$ with $k\in\{0,\dots,N-1\}$ to determine the $f_{lmn}$ complex coefficients thanks to the following simplified relation
$$
f_{lmn} = \sum_{k=0}^{N-1}w_k\ a_{lmk}\ {\cal K}_n(r_k) \quad\quad n\in\{0,\dots,N-1\}
f_{lmn} = \sum_{k=0}^{N-1}w_k\ a_{lmk}\ {\cal L}_n(r_k) \quad\quad n\in\{0,\dots,N-1\}
$$
*(here to simplify the notation I have introduced the ${\cal L}$ function but it is derived from the Gauss-Quadrature expression given above)*
+ **Synthesis**:
The synthesis also proceeds in two steps as the above detailed analysis. The algorithm receives as input the $f_{lmn}$ complex coefficients
+ to first compute for each $(l,m)$ couple of indices the intermediate $a_{lmk}$ complex values obtained from the Inverse Laguerre Transform:
......
......@@ -15,8 +15,6 @@
namespace LagSHT {
// void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
// vector<r_8>& fijk) {
void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& flmn,
vector<r_8>& fijk,
vector< complex<r_8> >& flmk) {
......@@ -83,8 +81,6 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
#endif
//Perform the Spherical Hamonic analysis first
// vector< complex<r_8> > flmk(NtotCoeff);
int off7n, off7k;
tstack_push("SHT Analysis");
......
......@@ -50,9 +50,12 @@ namespace LagSHT {
class LaguerreSphericalTransform {
public:
//JEC 13/1/16 alpha is double
// LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax,
// r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2, bool R2tau=true):
LaguerreSphericalTransform(string spheregeom, int Lmax, int Mmax, int Nmax,
r_8 R, int Nrings = -1, int Nphi = -1, int alpha=2, bool R2tau=true):
r_8 R, int Nrings = -1, int Nphi = -1, r_8 alpha=2, bool R2tau=true):
Lmax_(Lmax), Mmax_(Mmax), R_(R) {
//Factory
......@@ -108,7 +111,8 @@ class LaguerreSphericalTransform {
//! Getters
int GetOrder() const { return N_; }
int GetAlpha() const { return alpha_; }
// int GetAlpha() const { return alpha_; }
r_8 GetAlpha() const { return alpha_; } //JEC 13/1/16 r_8
int GetNPixels() const { return Npix_; }
int GetNAlm() const { return Nalm_; }
......@@ -145,9 +149,6 @@ class LaguerreSphericalTransform {
}
//! Analysis
/*! \brief Pixels -> Coeffs with Input/Output using T floating representation
\input fijk: 3D real function values
......@@ -159,8 +160,6 @@ class LaguerreSphericalTransform {
Analysis(fijk,flmn,almk);
}
//! Analysis
/*! \brief Pixels -> Coeffs with Input/Output using T floating representation
\input fijk: 3D real function values
......@@ -186,7 +185,8 @@ protected:
int N_; //!< Order of the Generalized Lagurerre polynomial
int alpha_; //!< second parameter of the generalized Laguerre polynomial
// int alpha_; //!< second parameter of the generalized Laguerre polynomial
r_8 alpha_; //!< second parameter of the generalized Laguerre polynomial
r_8 R_; //!< Radial maximal value (not yet used)
......
......@@ -176,6 +176,7 @@ void MultiLaguerreTransformSdtAlone() {
//---------------------------------------
void MultiSphericalLaguerreTransform() {
r_8 alpha = param.alpha; //JEC 13/1/16
int Nmax = param.N;
int Lmax = param.Lmax;
r_8 Rmax = param.R;
......
......@@ -31,29 +31,58 @@ LaguerreTransform::LaguerreTransform(int N, r_8 R, r_8 alpha, bool R2tau) :
if(alpha<0) throw LagSHTError("LagTransform call with alpha<0");
//reminder: LaguerreFuncQuad is called with N to compute nodes & weioghts order N
LaguerreFuncQuad lag(N_, alpha_);
lag.QuadWeightNodes(nodes_,weights_);
//JEC 12/1/16 alpha is a r_8
// int tmp = 1;
// for(int i=1;i<=alpha_; i++) tmp *= i;
// alphaFact_ = sqrt((r_8)tmp); //sqrt(alpha!)
alphaFact_ = sqrt(boost::math::tgamma(alpha_+1.0)); //<---------- usage of alpha for normalieation
{
//reminder: LaguerreFuncQuad is called with N to compute nodes & weioghts order N
LaguerreFuncQuad lag(N_, alpha_);
lag.QuadWeightNodes(nodes_,weights_);
}
r_8 tau3sqrt; // tau^(3/2)
r_8 tau3sqrtinv; //tau^(-3/2)
if(R2tau) {
tau_ = R_/nodes_[N_-1];
r_8 pow3inv2 = 3./2.;
tau3sqrt_ = pow(tau_,pow3inv2);
tau3sqrtinv_ = 1./tau3sqrt_;
tau3sqrt = pow(tau_,pow3inv2);
tau3sqrtinv = 1./tau3sqrt;
} else { //this case may be used when Rmax is not bounded
tau_ = 1.;
tau3sqrt_ = 1.;
tau3sqrtinv_ = 1.;
tau3sqrt = 1.;
tau3sqrtinv = 1.;
R_ = nodes_[N_-1];
}
//JEC 12/1/16 alpha is a r_8
// int tmp = 1;
// for(int i=1;i<=alpha_; i++) tmp *= i;
// alphaFact_ = sqrt((r_8)tmp); //sqrt(alpha!)
//JEC 12/1/16 still valid for alpha r_8
vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) the normalization
facts[0] = 1.0/sqrt(boost::math::tgamma(alpha_+1.0)); //1/sqrt(alpha)!
for(int n=1;n<N_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
{// compute LnkMtx matrices for Analysis and Synthesis
LaguerreFuncQuad lag(N_-1, alpha_); //<----- make explicit alpha here JEC 12/1/16; N-1 => Values() computes with n in [0,N-1]
vector<r_8> LnAllANA(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnAllSYN(N_); //all the values Lag_n(r) n:0...N-1
LnkMtxANA_.resize(N_*N_);
LnkMtxSYN_.resize(N_*N_);
for (int k=0; k<N_; k++){
r_8 rk = nodes_[k];
r_8 prk = pow(rk, alphap_); //< r_k^(alpha/2-1)
lag.Values(rk,LnAllANA,1./prk); //<----- use scaling to introduce r_k^(1-alpha/2) in the computation (Analysis)
lag.Values(rk,LnAllSYN,prk); //<----- use scaling to introduce r_k^(alpha/2-1) in the computation (Synthesis)
for (int n = 0; n<N_; n++ ){
LnkMtxANA_[n+N_*k] = LnAllANA[n]*facts[n]*weights_[k]*tau3sqrt;
LnkMtxSYN_[n*N_+k] = LnAllSYN[n]*facts[n]*tau3sqrtinv;
}
}
}//end scope lag
}//Ctor
int LaguerreTransform::R2Index(r_8 r) const {
......@@ -78,28 +107,27 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
std::copy(fi.begin(), fi.end(), std::ostream_iterator<char>(std::cout, " "));
#endif
//JEC 12/1/16 still valid for alpha r_8 (Todo: put facts[n] in Ctor)
r_8 invalphaFact = 1.0/alphaFact_; //1/sqrt(alpha)!
vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) the normalization
facts[0] = invalphaFact;
for(int n=1;n<N_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
LaguerreFuncQuad lag(N_-1, alpha_); //<----- make explicit alpha here JEC 12/1/16; N-1 => Values() computes with n in [0,N-1]
vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r_k) k,n:0...N-1
for (int k=0; k<N_; k++){
r_8 rk = nodes_[k];
r_8 prk = pow(rk, -alphap_); //JEC 12/1/16 mind the sign
lag.Values(rk,LnAll,prk);
for (int n = 0; n<N_; n++ ){
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k]*tau3sqrt_;
}
}
// //JEC 13/1/16 put in Ctor
// r_8 invalphaFact = 1.0/alphaFact_; //1/sqrt(alpha)!
// vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) the normalization
// facts[0] = invalphaFact;
// for(int n=1;n<N_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
// LaguerreFuncQuad lag(N_-1, alpha_); //<----- make explicit alpha here JEC 12/1/16; N-1 => Values() computes with n in [0,N-1]
// vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
// vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r_k) k,n:0...N-1
// for (int k=0; k<N_; k++){
// r_8 rk = nodes_[k];
// r_8 prk = pow(rk, -alphap_); //JEC 12/1/16 mind the sign
// lag.Values(rk,LnAll,prk);
// for (int n = 0; n<N_; n++ ){
// LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k]*tau3sqrt;
// }
// }
#ifdef CBLAS
cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (r_8*)(&fi[0]), 2*stride, &LnkMtx[0], N_, 0, (r_8 *)(&fn[0]), 2*stride);
//JEC 13/1/16 LnkMtx -> LnkMtx_
cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (r_8*)(&fi[0]), 2*stride, &LnkMtxANA_[0], N_, 0, (r_8 *)(&fn[0]), 2*stride);
// cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (r_8*)(&fi[0]), 2*stride, &LnkMtx[0], N_, 0, (r_8 *)(&fn[0]), 2*stride);
#else
vector<complex<r_8> > vtmp(N_);
for (int l=0; l<stride; l++) {
......@@ -108,8 +136,8 @@ void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi,
complex<r_8> fli = fi[l+ i*stride];
//////////JEC TO BE CHECK IF UNUSED r_8 wi = weights_[i];
for (int n=0; n<N_; n++){
// vtmp[n] += fli * wi * LnkMtx[n*N_+i];
vtmp[n] += fli * LnkMtx[n+N_*i];
// vtmp[n] += fli * LnkMtx[n+N_*i];
vtmp[n] += fli * LnkMtxANA_[n+N_*i]; //JEC 13/1/16 use LnkMtx_
}//loop on k
}//loop on n
for (int n=0; n<N_; n++) {
......@@ -133,27 +161,22 @@ void LaguerreTransform::MultiAnalysisR(const vector<r_8>& fi,
#if DEBUG >= 2
std::copy(fi.begin(), fi.end(), std::ostream_iterator<char>(std::cout, " "));
#endif
//JEC 12/1/16 still valid for alpha r_8
r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!
vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) for normalization
facts[0] = invalphaFact;
for(int n=1;n<N_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
LaguerreFuncQuad lag(N_-1, alpha_); //<----- make explicit alpha here JEC 12/1/16
vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r) n:0...N-1 TRY r_16->r_8
for (int k=0; k<N_; k++){
r_8 rk = nodes_[k];
r_8 prk = pow(rk, -alphap_); //JEC 12/1/16 mind the sign
lag.Values(rk,LnAll,prk);
for (int n = 0; n<N_; n++ ){
LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k]*tau3sqrt_;
}
}
// //JEC 13/1/16 put in Ctor
// r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!
// vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) for normalization
// facts[0] = invalphaFact;
// for(int n=1;n<N_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
// LaguerreFuncQuad lag(N_-1, alpha_); //<----- make explicit alpha here JEC 12/1/16
// vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
// vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r) n:0...N-1 TRY r_16->r_8
// for (int k=0; k<N_; k++){
// r_8 rk = nodes_[k];
// r_8 prk = pow(rk, -alphap_); //JEC 12/1/16 mind the sign
// lag.Values(rk,LnAll,prk);
// for (int n = 0; n<N_; n++ ){
// LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k]*tau3sqrt;
// }
// }
//#ifdef CBLAS
......@@ -167,8 +190,8 @@ void LaguerreTransform::MultiAnalysisR(const vector<r_8>& fi,
r_8 fli = fi[l+ i*stride];
//////////JEC TO BE CHECK IF UNUSED r_8 wi = weights_[i];
for (int n=0; n<N_; n++){
// vtmp[n] += fli * wi * LnkMtx[n*N_+i];
vtmp[n] += fli * LnkMtx[n+N_*i];
// vtmp[n] += fli * LnkMtx[n+N_*i];
vtmp[n] += fli * LnkMtxANA_[n+N_*i]; //JEC 13/1/16 use LnkMtx_
}//loop on k
}//loop on n
for (int n=0; n<N_; n++) {
......@@ -196,27 +219,30 @@ void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn,
#endif
//JEC 12/1/16 still valid for alpha r_8 (Todo: put facts[n] in Ctor)
r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!
vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
facts[0] = invalphaFact;
for(int n=1; n<N_; n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
// //JEC 13/1/16 put in Ctor
// r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!
// vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
// facts[0] = invalphaFact;
// for(int n=1; n<N_; n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
// LaguerreFuncQuad lag(N_-1, alpha_); //<----- make explicit alpha here JEC 12/1/16
// vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
// vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r_k) n,:0...N-1
// for (int k=0; k<N_; k++){
// r_8 rk = nodes_[k];
// r_8 prk = pow(rk, alphap_); //JEC 12/1/16
// lag.Values(rk,LnAll,prk); //<---- use the scaling
// for (int n = 0; n<N_; n++ ){
// LnkMtx[n*N_+k] = LnAll[n]*facts[n]*tau3sqrtinv; //JEC 12/1/16 add prk
// }
// }
LaguerreFuncQuad lag(N_-1, alpha_); //<----- make explicit alpha here JEC 12/1/16
vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r_k) n,:0...N-1
for (int k=0; k<N_; k++){
r_8 rk = nodes_[k];
r_8 prk = pow(rk, alphap_); //JEC 12/1/16
lag.Values(rk,LnAll,prk); //<---- use the scaling
for (int n = 0; n<N_; n++ ){
LnkMtx[n*N_+k] = LnAll[n]*facts[n]*tau3sqrtinv_; //JEC 12/1/16 add prk
}
}
#if CBLAS
cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (r_8*)(&fn[0]), 2*stride, &LnkMtx[0], N_, 0, (r_8 *)(&fi[0]), 2*stride);
//JEC 13/1/16 use LnkMtx_
cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (r_8*)(&fn[0]), 2*stride, &LnkMtxSYN_[0], N_, 0, (r_8 *)(&fi[0]), 2*stride);
// cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (r_8*)(&fn[0]), 2*stride, &LnkMtx[0], N_, 0, (r_8 *)(&fi[0]), 2*stride);
#else
cout << "CBLAS desactivated................" << endl;
vector<complex<r_8> > vtmp(N_);
......@@ -225,7 +251,8 @@ void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn,
for (int n = 0; n<N_; n++ ){
complex<r_8> fln = fn[l+ n*stride];
for (int k=0; k<N_; k++){
vtmp[k] += fln * LnkMtx[n*N_+k];
// vtmp[k] += fln * LnkMtx[n*N_+k];
vtmp[k] += fln * LnkMtxSYN_[n*N_+k]; //JEC 13/1/16
}//loop on k
}//loop on n
for (int k=0; k<N_; k++) fi[l+k*stride] += vtmp[k];
......
......@@ -132,15 +132,19 @@ class LaguerreTransform : public BaseLaguerreTransform {
protected:
r_8 alphaFact_; //!< sqrt(alpha!)
// r_8 alphaFact_; //!< sqrt(alpha!) NOT USED except in Ctor 13/1/16
r_8 R_; //!< maximal radius
vector<r_8> nodes_; //!< nodes of the Gauss-Laguerre quadrature used
vector<r_8> weights_; //!< weights of the Gauss-Laguerre quadrature used
r_8 tau_; //!< a scaling factor = R/nodes[N-1] JEC 22/9/15
r_8 tau3sqrt_; // tau^(3/2)
r_8 tau3sqrtinv_; // tau^(-3/2)
vector<r_8> LnkMtxANA_; //all the values Lag_n(r_k) k,n:0...N-1 for Analysis (JEC: 13/1/16)
vector<r_8> LnkMtxSYN_; //all the values Lag_n(r_k) k,n:0...N-1 for Synthesis (JEC: 13/1/16)
r_8 tau_; //!< a scaling factor = R/nodes[N-1]
// r_8 tau3sqrt_; // tau^(3/2) NOT USED except in Ctor 13/1/16
//r_8 tau3sqrtinv_; // tau^(-3/2) NOT USED except in Ctor 13/1/16
public:
class distance {
public:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment