📣 An issue occured with the embedded container registry on October 25 2021, between 10:30 and 12:10 (UTC+2). Any persisting issues should be reported to CC-IN2P3 Support. 🐛

Commit 321da181 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 4/2/16 change lowerBound for Jlnp computation

parent 4fef23e3
......@@ -17,7 +17,7 @@ namespace LagSHT {
/* Version 0: test de calcul des J_{lnp} via une quadrature de Laguerre
mais ca n'a pas l'air de converger (28 Sept 15)
Version 1: on fait une integration a la Clenshaw-Curtis avec startegie Localeo ou Globale
Version 1: on fait une integration a la Clenshaw-Curtis avec startegie Locale ou Globale
recursive
Version 2: Jlnp calcule via une Discrete Fourier-Bessel Transform and a Clenshaw-Curtis quadrature and global integration strategy.
Version 3: calcul via une quadrature de Clenshaw-Curtis quadrature and global integration en controlant les bornes
......@@ -86,158 +86,6 @@ Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
}//Ctor
/*!
Computation of the Jlnp = Jln(k_{lp}) = Sqrt(2/pi) tau^3/2 Int_0^infty jl(k_{lp} tau x) K_n(x) x^2 dx
with
K_n(x) = Exp[-x/2] Ln(x) Sqrt(n!/(n+alpha)!) ->cf. LaguerreFuncQuad
k_{lp}*tau = q_{lp}/r_{Nmax-1}
Apply a Discrete Spherical transform Lanusse et al A&A 540, A92 (2012) Eq.40 with l'=l
Jlnp ~ Kapa_l^(-3) tau^(3/2) Sqrt(2 Pi) Sum_{p'=0}^{PPmax-1} Kn(r_{lp'}) jl(q_{lp'}q_{lp}/q_{l,PPmax-1})/[jl+1(q_{lp})]^2
with r_{lp'} = q_{lp'}/q_{l,PPmax-1}r_{Nmax-1}
notice PPmax as nothing to do in principle with Pmax the Fourier-Bessel max value of the 3rd index
Kapa_l = q_{l,PPmax-1}/r_{Nmax-1}
JEC: 12/11/15 Add the contribution of Integral_{r_{Nmax-1}}^{r_{Nmax-1} + 3(r_{Nmax-1}-r_{Nmax-2})} by Clenshaw-Curtis quadrature
void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
r_8 tau = lagTrans_->GetTau();
r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
// cout << "(JEC) tau3sqrt= " << tau3sqrt << endl;
//JEC 12/11/15 Start
r_8 rNm1 = Rmax_ / tau; // R = r[N-1] voir acces direct = (lagTrans_->GetNodes()).back()
r_8 rNm2 = (lagTrans_->GetNodes())[Nmax_-2];
r_8 deltaR = (rNm1-rNm2);
r_8 integLowerBound = rNm1;
r_8 integUpperBound = rNm1+3.0*deltaR; //for the integral
r_8 tol = 1e-10;
typedef ClenshawCurtisQuadrature<double> Integrator_t;
string clenshawFile = clenshawDir+"/ClenshawCurtisRuleData-20.txt";
Integrator_t theQuad(20,clenshawFile,false); //false=do not recompute the weights and nodes of the quadrature
Quadrature<r_8,Integrator_t>::values_t integ_val;
//JEC 12/11/15 end
int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
Jlnp_.resize(Nmax_*Jstride);
//JEC 16/11/15 START: Todo all that might be in common somewhere (hint: laguerreTransform.h)
int alpha = lagTrans_->GetAlpha();
r_8 alphaFact = 1; // alpha!
for(int i=1;i<=alpha; i++) alphaFact *= i;
alphaFact = sqrt((r_8)alphaFact); //sqrt(alpha!)
r_8 invalphaFact = 1.0/alphaFact; // 1/sqrt(alpha)!
vector<r_8> facts(Nmax_); //sqrt(n!/(n+alpha)!) en iteratif
facts[0] = invalphaFact;
for(int n=1;n<Nmax_;n++) {
facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha)) );
}
//END
int PPmaxm1 = PPmax_ -1;
tstack_push("J Integ compute ");
stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<"-PP"<<PPmax_<<".txt";
std::ofstream ofs;
string fname(jlnpDir+"/"+ss.str());
ofs.open (fname.c_str(), std::ofstream::out);
// cout << "(JEC) Jlpn recompute..." << endl;
for(int p=0; p<Pmax_; p++){
for(int l=0; l<Lmax_; l++){
int ploff7 = p + l*Pmax_;
r_8 qlp = (*jnu_)(l,p); //ql,p
r_8 klp = qlp/rNm1; //JEC 25/1/16
r_8 qlpx = (*jnu_)(l,PPmaxm1); //ql,PPmax-1
r_8 Kapa_l = qlpx/rNm1; // Kapa_l = (ql,PPmax-1)/R
r_8 Kapa_l3 = Kapa_l * Kapa_l * Kapa_l;
r_8 coeff = sqrt(2.0 *M_PI)/Kapa_l3; // Sqrt(2 Pi)/Kapa_l^3
for (int n=0; n<Nmax_; n++){
LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
r_8 sum = 0.;
//contribution of the integral: [0,rNm1] estimated
//using Discrete Fourier-Bessel Transform
for (int pp=0; pp<PPmax_; pp++){
r_8 qlpp = (*jnu_)(l,pp); //ql,pp
r_8 jlp1 = Bessel::jn(l+1,qlpp); //j(l+1,ql,pp)
r_8 rlpp = qlpp/qlpx * rNm1; // rl,pp = (ql,pp/ql,PPmax-1) * rNm1
r_8 jlqq = Bessel::jn(l,qlpp * qlp/qlpx);
sum += jlqq/(jlp1*jlp1) * Ln->Value(rlpp);
}//pp-loop
// Jlnp_[ploff7 + n*Jstride] = (coeff * sum) * (tau3sqrt*facts[n]);
//JEC 23/11/15 save the Rmax-independant part
Jlnp_[ploff7 + n*Jstride] = coeff * facts[n] * sum;
//contribution [rNm1, rNm1+deltaR] estimated with Clenshaw-Curtis quadrature (JEC 12/11/15) Start
IntFunc1D f(l,qlp/rNm1,Ln);
theQuad.SetFuncBounded(f,integLowerBound,integUpperBound);
integ_val = theQuad.GlobalAdapStrat(tol,2,1);
// Jlnp_[ploff7 + n*Jstride] += integ_val.first * (tau3sqrt*facts[n]*sqrt(2.0/M_PI));
//JEC 23/11/15 save the Rmax-independant part
Jlnp_[ploff7 + n*Jstride] += integ_val.first * (facts[n]*sqrt(2.0/M_PI));
//Limber approx JEC 25/1/16 the sqrt(2/Pi) is canceled by the definition of SphericalBesselJ function used in this approx
r_8 limber = pow(klp,(r_8)(-3.0)) * pow((r_8)(l+0.5),(r_8)1.5) * Ln->Value(((r_8)(l)+0.5)/klp)*facts[n];
//New Approx JEC 27/1/16
theQuad.SetFuncBounded(f,0.8*r_8(l),rNm1+4.0*deltaR);
integ_val = theQuad.GlobalAdapStrat(tol,2,1);
cout <<l<<" "
<<n<<" "
<<p<< " "
<< setprecision(30) << Jlnp_[ploff7 + n*Jstride]
<< " " << limber
<< " " << integ_val.first * (facts[n]*sqrt(2.0/M_PI))
<< endl;
//JEC 23/11/15 save the Rmax-independant part
ofs <<l<<" "
<<n<<" "
<<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
Jlnp_[ploff7 + n*Jstride] *= tau3sqrt;
delete Ln;
}//n
}//l
}//p
ofs.close();
tstack_pop("J Integ compute ");
}//ComputeJlpnInteg
*/
/*!
Computation of the Jlnp = Jln(k_{lp}) = Sqrt(2/pi) tau^3/2 Int_0^infty jl(k_{lp} tau x) K_n(x) x^2 dx
......@@ -290,22 +138,45 @@ void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
}
stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
stringstream ss; ss << "NEW-jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
std::ofstream ofs;
string fname(jlnpDir+"/"+ss.str());
ofs.open (fname.c_str(), std::ofstream::out);
// //JEC TMP
// stringstream ssTMP; ssTMP << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
// string fnameTMP(jlnpDir+"/"+ss.str());
// std::ifstream ifs(fnameTMP.c_str(), std::ifstream::in);
// if(!ifs.is_open())
// throw LagSHTError("ERROR: jlnp file open failed.:"+ssTMP.str());
// //JEC TMP
// int pcur, lcur, ncur; //TMP
// r_8 val; //TMP
for(int p=0; p<Pmax_; p++){
for(int l=0; l<Lmax_; l++){
int ploff7 = p + l*Pmax_;
r_8 qlp = (*jnu_)(l,p); //ql,p
r_8 klptau = qlp/rNm1;
r_8 intLowBound = 0.8*l;
r_8 intLowBound = (l<=10)? 0.0: 0.5*l/klptau; //JEC 4/2/16 for low l just start at 0
for (int n=0; n<Nmax_; n++){
// ifs >> lcur >> ncur >> pcur >> val; //TMP read
// if( (l == lcur) || (n == ncur) || (p == pcur)){ //TMP check
// ofs <<l<<" "
// <<n<<" "
// <<p<< " " << setprecision(30) << val << endl;
// Jlnp_[ploff7 + n*Jstride] = tau3sqrt * val;
// continue;
// }
LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
IntFunc1D f(l,klptau,Ln);
......@@ -334,6 +205,7 @@ void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
}//p
// ifs.close(); //TMP
ofs.close();
tstack_pop("J Integ compute ");
......@@ -497,3 +369,157 @@ void Laguerre2Bessel::Print(ostream& os, int what){
}//Print
}//namespace
/*!
Computation of the Jlnp = Jln(k_{lp}) = Sqrt(2/pi) tau^3/2 Int_0^infty jl(k_{lp} tau x) K_n(x) x^2 dx
with
K_n(x) = Exp[-x/2] Ln(x) Sqrt(n!/(n+alpha)!) ->cf. LaguerreFuncQuad
k_{lp}*tau = q_{lp}/r_{Nmax-1}
Apply a Discrete Spherical transform Lanusse et al A&A 540, A92 (2012) Eq.40 with l'=l
Jlnp ~ Kapa_l^(-3) tau^(3/2) Sqrt(2 Pi) Sum_{p'=0}^{PPmax-1} Kn(r_{lp'}) jl(q_{lp'}q_{lp}/q_{l,PPmax-1})/[jl+1(q_{lp})]^2
with r_{lp'} = q_{lp'}/q_{l,PPmax-1}r_{Nmax-1}
notice PPmax as nothing to do in principle with Pmax the Fourier-Bessel max value of the 3rd index
Kapa_l = q_{l,PPmax-1}/r_{Nmax-1}
JEC: 12/11/15 Add the contribution of Integral_{r_{Nmax-1}}^{r_{Nmax-1} + 3(r_{Nmax-1}-r_{Nmax-2})} by Clenshaw-Curtis quadrature
void Laguerre2Bessel::ComputeJInteg(string clenshawDir, string jlnpDir){
r_8 tau = lagTrans_->GetTau();
r_8 tau3sqrt = pow(tau,(r_8)(3./2.));
// cout << "(JEC) tau3sqrt= " << tau3sqrt << endl;
//JEC 12/11/15 Start
r_8 rNm1 = Rmax_ / tau; // R = r[N-1] voir acces direct = (lagTrans_->GetNodes()).back()
r_8 rNm2 = (lagTrans_->GetNodes())[Nmax_-2];
r_8 deltaR = (rNm1-rNm2);
r_8 integLowerBound = rNm1;
r_8 integUpperBound = rNm1+3.0*deltaR; //for the integral
r_8 tol = 1e-10;
typedef ClenshawCurtisQuadrature<double> Integrator_t;
string clenshawFile = clenshawDir+"/ClenshawCurtisRuleData-20.txt";
Integrator_t theQuad(20,clenshawFile,false); //false=do not recompute the weights and nodes of the quadrature
Quadrature<r_8,Integrator_t>::values_t integ_val;
//JEC 12/11/15 end
int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
Jlnp_.resize(Nmax_*Jstride);
//JEC 16/11/15 START: Todo all that might be in common somewhere (hint: laguerreTransform.h)
int alpha = lagTrans_->GetAlpha();
r_8 alphaFact = 1; // alpha!
for(int i=1;i<=alpha; i++) alphaFact *= i;
alphaFact = sqrt((r_8)alphaFact); //sqrt(alpha!)
r_8 invalphaFact = 1.0/alphaFact; // 1/sqrt(alpha)!
vector<r_8> facts(Nmax_); //sqrt(n!/(n+alpha)!) en iteratif
facts[0] = invalphaFact;
for(int n=1;n<Nmax_;n++) {
facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha)) );
}
//END
int PPmaxm1 = PPmax_ -1;
tstack_push("J Integ compute ");
stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<"-PP"<<PPmax_<<".txt";
std::ofstream ofs;
string fname(jlnpDir+"/"+ss.str());
ofs.open (fname.c_str(), std::ofstream::out);
// cout << "(JEC) Jlpn recompute..." << endl;
for(int p=0; p<Pmax_; p++){
for(int l=0; l<Lmax_; l++){
int ploff7 = p + l*Pmax_;
r_8 qlp = (*jnu_)(l,p); //ql,p
r_8 klp = qlp/rNm1; //JEC 25/1/16
r_8 qlpx = (*jnu_)(l,PPmaxm1); //ql,PPmax-1
r_8 Kapa_l = qlpx/rNm1; // Kapa_l = (ql,PPmax-1)/R
r_8 Kapa_l3 = Kapa_l * Kapa_l * Kapa_l;
r_8 coeff = sqrt(2.0 *M_PI)/Kapa_l3; // Sqrt(2 Pi)/Kapa_l^3
for (int n=0; n<Nmax_; n++){
LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
r_8 sum = 0.;
//contribution of the integral: [0,rNm1] estimated
//using Discrete Fourier-Bessel Transform
for (int pp=0; pp<PPmax_; pp++){
r_8 qlpp = (*jnu_)(l,pp); //ql,pp
r_8 jlp1 = Bessel::jn(l+1,qlpp); //j(l+1,ql,pp)
r_8 rlpp = qlpp/qlpx * rNm1; // rl,pp = (ql,pp/ql,PPmax-1) * rNm1
r_8 jlqq = Bessel::jn(l,qlpp * qlp/qlpx);
sum += jlqq/(jlp1*jlp1) * Ln->Value(rlpp);
}//pp-loop
// Jlnp_[ploff7 + n*Jstride] = (coeff * sum) * (tau3sqrt*facts[n]);
//JEC 23/11/15 save the Rmax-independant part
Jlnp_[ploff7 + n*Jstride] = coeff * facts[n] * sum;
//contribution [rNm1, rNm1+deltaR] estimated with Clenshaw-Curtis quadrature (JEC 12/11/15) Start
IntFunc1D f(l,qlp/rNm1,Ln);
theQuad.SetFuncBounded(f,integLowerBound,integUpperBound);
integ_val = theQuad.GlobalAdapStrat(tol,2,1);
// Jlnp_[ploff7 + n*Jstride] += integ_val.first * (tau3sqrt*facts[n]*sqrt(2.0/M_PI));
//JEC 23/11/15 save the Rmax-independant part
Jlnp_[ploff7 + n*Jstride] += integ_val.first * (facts[n]*sqrt(2.0/M_PI));
//Limber approx JEC 25/1/16 the sqrt(2/Pi) is canceled by the definition of SphericalBesselJ function used in this approx
r_8 limber = pow(klp,(r_8)(-3.0)) * pow((r_8)(l+0.5),(r_8)1.5) * Ln->Value(((r_8)(l)+0.5)/klp)*facts[n];
//New Approx JEC 27/1/16
theQuad.SetFuncBounded(f,0.8*r_8(l),rNm1+4.0*deltaR);
integ_val = theQuad.GlobalAdapStrat(tol,2,1);
cout <<l<<" "
<<n<<" "
<<p<< " "
<< setprecision(30) << Jlnp_[ploff7 + n*Jstride]
<< " " << limber
<< " " << integ_val.first * (facts[n]*sqrt(2.0/M_PI))
<< endl;
//JEC 23/11/15 save the Rmax-independant part
ofs <<l<<" "
<<n<<" "
<<p<< " " << setprecision(30) << Jlnp_[ploff7 + n*Jstride] << endl;
Jlnp_[ploff7 + n*Jstride] *= tau3sqrt;
delete Ln;
}//n
}//l
}//p
ofs.close();
tstack_pop("J Integ compute ");
}//ComputeJlpnInteg
*/
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