📣 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 10e20127 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 24/11/15 clean up test3D

parent 42092f70
......@@ -8,7 +8,6 @@
#include "LagSHT/walltimer.h"
#include "LagSHT/lagSphericTransform.h"
#include "LagSHT/laguerre2bessel.h"
using namespace LagSHT;
......@@ -64,7 +63,6 @@ struct PARAM {
int Lmax;
int N;
r_8 R;
int Pmax;
int alpha;
string geometry;
int ntheta;
......@@ -173,7 +171,6 @@ void test3() {
//initial values
int Lmax = param.Lmax;
int Nmax = param.N;
int Pmax = param.Pmax; //JEC new
r_8 Rmax = param.R;
string geom = param.geometry;
int ntheta = param.ntheta;
......@@ -254,85 +251,7 @@ void test3() {
sphlagtrans.Analysis(fijk,flmn, almk); // <---------- the job
tstack_pop("Lag Spherical Transf.");
//JEC 9/11/15 START
tstack_push("Laguerre 2 Bessel");
BaseGeometry* sphere = sphlagtrans.GetSphericalGeometry();
LaguerreTransform* lagTrans = sphlagtrans.GetLagTransform();
Laguerre2Bessel lag2bess(sphere,lagTrans,Nmax,Pmax,Rmax);
vector< complex<r_8> > FBlmp(Nalm*Pmax);
lag2bess.Lag2BessCoeff(flmn,FBlmp);
vector< complex<r_8> > FBalmk;
vector<r_8> fFBijk;
lag2bess.Synthesis(FBlmp,FBalmk,fFBijk);
tstack_pop("Laguerre 2 Bessel");
{//Check 1
cout << "Dump FL or FB reconstructed Alm(r_k)" <<endl;
for(int n=0;n<Nmax;n++){
r_8 err_abs(0.);
r_8 err_rel(0.);
for(int l=0;l<Lmax;l++){
for (int m=0;m<=l;m++) {
int id= n*Nalm + l+m*Lmax-m*(m+1)/2; // LagSHT numbering
cout << "n,lm: " << n << " ("<<l<<","<<m<<") : " << almk[id] << ", "<<FBalmk[id] << endl;
complex<r_8> cdiff = (almk[id] - FBalmk[id])*conj(almk[id] - FBalmk[id]);
r_8 diff = sqrt(cdiff.real());
if(diff>err_abs){
err_abs = diff;
}
complex<r_8> cref= almk[id]*conj(almk[id]);
r_8 ref = sqrt(cref.real());
r_8 relatdiff = diff/ref;
if(relatdiff>err_rel) err_rel = relatdiff;
}
}
cout << " >>>>>>>>>>>>>>>>>>>>> Shell["<<n<<"] <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
cout << "Almk: Err. Max. " << err_abs << ", Err. Rel. " << err_rel << endl;
cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
}
}//check 1
{//check 2
cout << "Dump FL or FB reconstructed fijk on each shell k" <<endl;
for (int ish=0;ish<Nshell; ish++){
r_8 err_abs(0.);
r_8 err_rel(0.);
for (int ip=0; ip<Npix; ip++) {
int id = ish*Npix+ip;
r_8 diff = fabs(fijk[id]-fFBijk[id]);
if(diff>err_abs){
err_abs = diff;
}
r_8 relatdiff = diff/ fijk[id];
if(relatdiff>err_rel) err_rel = relatdiff;
}//loop on px
cout << " >>>>>>>>>>>>>>>>>>>>> Shell["<<ish<<"] <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
cout << "Fijk: Err. Max. " << err_abs << ", Err. Rel. " << err_rel << endl;
cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
}//loop on shell
}//check 2
//JEC 9/11/15 END
tstack_push("Clk and Cln OUT");
cout << "Compute back the Cl(k) on each shell to compare to input, and also the analog Cl(n) from the flmn" << endl;
......@@ -609,7 +528,6 @@ int main(int narg, char *arg[]) {
int N = 128;
r_8 R = 5.;
int Lmax = 128;
int Pmax = 128;
int test=3;
......@@ -658,16 +576,11 @@ int main(int narg, char *arg[]) {
R = atof(arg[ka+1]);
ka+=2;
}
else if (strcmp(arg[ka],"-p")==0) {
Pmax = atof(arg[ka+1]);
ka+=2;
}
else ka++;
}//eo while
param.Lmax = Lmax;
param.N = N;
param.Pmax = Pmax;
param.R = R;
param.alpha = 2;
param.geometry = geometry;
......@@ -679,7 +592,6 @@ int main(int narg, char *arg[]) {
<< " Lmax, Nmax, Pmax, Rmax , alpha: "
<< param.Lmax << ", "
<< param.N << ", "
<< param.Pmax << ", "
<< param.R << ", "
<< param.alpha << "\n"
<< "Geometry: " << param.geometry << " Ntheta, Nphi: "
......
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