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

(JEC) 18/1/16 spin =/= 0 test

parent eea2e6db
......@@ -95,8 +95,7 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
#endif
sphere_->map2alm(&(fijk.data()[off7k]),reinterpret_cast<r_8*>(&(almk.data()[off7n])),false); //the false is to forbidd accumulation
}//loop on shell
#if DEBUG >= 1
......@@ -129,4 +128,144 @@ void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk,
void LaguerreSphericalTransform::Synthesis(const vector< complex<r_8> >& Elmn, const vector<complex<r_8> >& Blmn,
int spin,
vector<r_8>& fijk_re, vector<r_8>& fijk_im,
vector< complex<r_8> >& Elmk, vector<complex<r_8> >& Blmk
) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
if(Elmn.size() != (size_t)NtotCoeff || Blmn.size() != (size_t)NtotCoeff) throw LagSHTError("SphLagTransform::Synthesis size error");
fijk_re.resize(Ntot3DPix);
fijk_im.resize(Ntot3DPix);
//Perform the Inverse Laguerre Transform first
Elmk.resize(NtotCoeff);
Blmk.resize(NtotCoeff);
#if DEBUG >= 1
cout << "SphLagTransform::Synthesis type <r_8> I/O .... Go..." << endl;
cout << "Start Inverse Laguerre...." << endl;
#endif
tstack_push("Laguerre MultiSynthesis");
lagTrans_->MultiSynthesis(Elmn, Elmk, Nalm_);
lagTrans_->MultiSynthesis(Blmn, Blmk, Nalm_);
tstack_pop("Laguerre MultiSynthesis");
// cout << "Elmk, Blmk" <<endl;
// for(int i=0; i<100; i++){
// cout << "Elmk["<<i<<"]: " << Elmk[i] << "; "
// << "Blmk["<<i<<"]: " << Blmk[i] << endl;
// }
#if DEBUG >= 1
cout << "End Inverse Laguerre...." << endl;
cout << "Start Inverse H.Spherical Transform " << endl;
#endif
tstack_push("SHT Synthesis");
int off7n, off7k;
//perform the synthesis per shell
for(int n =0; n<N_; n++){
off7n = n*Nalm_;
off7k = n*Npix_;
#if DEBUG >= 2
cout << "Synthesis... at shell num: " << n << endl;
#endif
sphere_->alm2map_spin(reinterpret_cast<const r_8*>(&(Elmk.data()[off7n])),
reinterpret_cast<const r_8*>(&(Blmk.data()[off7n])),
&fijk_re.data()[off7k],
&fijk_im.data()[off7k],
spin,
false); //the false is to forbidd accumulation
}//loop on shell
tstack_pop("SHT Synthesis");
#if DEBUG >= 1
cout << "SphLagTransform::Synthesis.... Done" << endl;
#endif
}//Synthesis
void LaguerreSphericalTransform::Analysis(const vector<r_8>& fijk_re, const vector<r_8>& fijk_im,
int spin,
vector< complex<r_8> >& Elmn, vector<complex<r_8> >& Blmn,
vector< complex<r_8> >& Elmk, vector<complex<r_8> >& Blmk
) {
int Ntot3DPix = Npix_*N_;
int NtotCoeff = Nalm_*N_;
if(fijk_re.size() != (size_t)Ntot3DPix || fijk_im.size() != (size_t)Ntot3DPix) throw LagSHTError("SphLagTransform::Analysis size error");
Elmn.resize(NtotCoeff);
Blmn.resize(NtotCoeff);
//Perform the Inverse Laguerre Transform first
Elmk.resize(NtotCoeff);
Blmk.resize(NtotCoeff);
#if DEBUG >= 1
cout << "SphLagTransform::Analysis.... Go..." << endl;
#endif
//Perform the Spherical Hamonic analysis first
int off7n, off7k;
tstack_push("SHT Analysis");
//Ylm analysis per shell
for(int n =0; n<N_; n++){
off7n = n*Nalm_;
off7k = n*Npix_;
#if DEBUG >= 1
cout << "Start Harmonic Spheric... at n = "<< n << endl;
#endif
sphere_->map2alm_spin(&(fijk_re.data()[off7k]),
&(fijk_im.data()[off7k]),
reinterpret_cast<r_8*>(&(Elmk.data()[off7n])),
reinterpret_cast<r_8*>(&(Blmk.data()[off7n])),
spin,
false); //the false is to forbidd accumulation
}//loop on shell
#if DEBUG >= 1
cout << "End Harmonic Spheric..." << endl;
cout << "Start Laguerre..." << endl;
#endif
tstack_pop("SHT Analysis");
tstack_push("Laguerre MultiAnalysis");
//Laguerre Transform
lagTrans_->MultiAnalysis(Elmk, Elmn, Nalm_);
lagTrans_->MultiAnalysis(Blmk, Blmn, Nalm_);
#if DEBUG >= 1
cout << "End Laguerre..." << endl;
#endif
tstack_pop("Laguerre MultiAnalysis");
#if DEBUG >= 2
for (int i=0;i<NtotCoeff; i++) {
cout << "flmn("<<i<<"): " << Elmn[i] << ", " << Blmn[i] << endl;
}
#endif
#if DEBUG >= 1
cout << "SphLagTransform::Analysis.... Done" << endl;
#endif
}//Analysis
}// Fin de namespace
......@@ -171,6 +171,26 @@ class LaguerreSphericalTransform {
vector< complex<r_8> >& almk
);
/*!
Test of spin-weighted decomposition
*/
void Synthesis(const vector< complex<r_8> >& Elmn, const vector<complex<r_8> >& Blmn,
int spin,
vector<r_8>& fijk_re, vector<r_8>& fijk_im,
vector< complex<r_8> >& Elmk, vector<complex<r_8> >& Blmk
);
void Analysis(const vector<r_8>& fijk_re, const vector<r_8>& fijk_im,
int spin,
vector< complex<r_8> >& Elmn, vector<complex<r_8> >& Blmn,
vector< complex<r_8> >& Elmk, vector<complex<r_8> >& Blmk
);
protected:
BaseGeometry* sphere_; //<! the 2D pixelization of the shells
......
......@@ -112,6 +112,18 @@ class BaseGeometry {
sharpjob_.map2alm(map,alm,add);
}
virtual void alm2map_spin(const r_8 *alm1, const r_8 *alm2, r_8 *map1, r_8 *map2,
int spin, bool add=false){
sharpjob_.alm2map_spin(alm1,alm2,map1,map2,spin,add);
}
virtual void map2alm_spin(const r_8 *map1, const r_8 *map2, r_8 *alm1, r_8 *alm2,
int spin, bool add){
sharpjob_.map2alm_spin(map1,map2,alm1,alm2,spin,add);
}
protected:
//! Set the vector of pixel positions
......
......@@ -30,6 +30,7 @@ struct PARAM {
int N;
r_8 R;
int Pmax;
int spin; //JEC 18/1/16
// int alpha;
double alpha;
string geometry;
......@@ -197,7 +198,8 @@ void MultiSphericalLaguerreTransform() {
Nmax,
Rmax,
ntheta,
nphi
nphi,
alpha
);
......@@ -289,6 +291,140 @@ void MultiSphericalLaguerreTransform() {
}//MultiSphericalLaguerreTransform
//---------------------------------------
void SpinMultiSphericalLaguerreTransform() {
r_8 alpha = param.alpha; //JEC 13/1/16
int Nmax = param.N;
int Lmax = param.Lmax;
r_8 Rmax = param.R;
string geometry = param.geometry;
int ntheta = param.ntheta;
int nphi = param.nphi;
int spin = param.spin;
cout << " ___________ MultiSphericalLaguerreTransform TEST _____________ " << endl;
tstack_push("data input");
int state = 1234567 + 8912 ; //random seed
LaguerreSphericalTransform sphlagtrans(geometry,
Lmax,
-1,
Nmax,
Rmax,
ntheta,
nphi,
alpha
);
int Nalm = sphlagtrans.GetSphericalGeometry()->Nalm();
int Nshell = sphlagtrans.GetOrder();
int Ntot = Nshell * Nalm; //total number of Coefficients of the Spherical Laguerre transform
int Npix = sphlagtrans.GetSphericalGeometry()->Npix();
int NpTot= Nshell * Npix; //totoal number of 3D-pixels
cout << "Verif: Npix, Nptot, Nalm, Nshell, Ntot: "
<< Npix << " "
<< NpTot << " "
<< Nalm << " "
<< Nshell << " "
<< Ntot << endl;
vector< complex<r_8> > ElmnOrig(Ntot);
vector< complex<r_8> > BlmnOrig(Ntot);
int id=-1;
for(int n=0;n<Nmax;n++){ //on each shell
for(int m=0;m<Lmax;m++){ //Warning the filling of alm is adapted for libsharp memory
for(int l=m;l<Lmax;l++){
id++;
r_8 rv = drand(-1,1,&state);
r_8 iv = (m==0) ? 0.0 : drand(-1,1,&state);
ElmnOrig[id] = complex<r_8>(rv,iv);
rv = drand(-1,1,&state);
iv = (m==0) ? 0.0 : drand(-1,1,&state);
BlmnOrig[id] = complex<r_8>(rv,iv);
}//end l-loop
}//end m-loop
}//end loop on shell
#if DEBUG >=2
std::copy(ElmnOrig.begin(), ElmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
std::copy(BlmnOrig.begin(), BlmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
#endif
tstack_pop("data input");
tstack_push("processing");
tstack_push("processing part Synthesis");
#if DEBUG >= 1
cout << "Main:Synthesis complex function..." << endl;
#endif
vector<r_8> fijk_re; //NpTot
vector<r_8> fijk_im;
vector< complex<r_8> > Elmk; //calcul intermediaire
vector< complex<r_8> > Blmk;
sphlagtrans.Synthesis(ElmnOrig, BlmnOrig, spin, fijk_re, fijk_im, Elmk, Blmk);
//#if DEBUG >= 2
for (int i=0; i<NpTot; i++){
cout << "fijk("<<i<<"): " << fijk_re[i] << ", " << fijk_im[i] << endl;
}
//#endif
tstack_pop("processing part Synthesis");
tstack_push("processing part Analysis");
#if DEBUG >= 1
cout << "Main:Analysis complex function..." << endl;
#endif
vector< complex<r_8> > Elmn;
vector< complex<r_8> > Blmn;
vector< complex<r_8> > Elmk_ana;
vector< complex<r_8> > Blmk_ana;
sphlagtrans.Analysis(fijk_re, fijk_im, spin, Elmn, Blmn, Elmk_ana, Blmk_ana);
tstack_pop("processing part Analysis");
tstack_pop("processing");
cout << "Error analysis..." << endl;
r_8 err_abs(0.);
r_8 err_rel(0.);
int imax = -1;
for(int i=0;i<Ntot;i++){
if(i>(Ntot-9)) {
cout << "(" << i << ") : Elnm Orig " << ElmnOrig[i] << " <-> Elmn Rec " << Elmn[i] << endl;
cout << "............ Blnm Orig " << BlmnOrig[i] << " <-> Blmn Rec " << Blmn[i] << endl;
}
complex<r_8> cdiff = ((ElmnOrig[i] - Elmn[i]) * conj(ElmnOrig[i] - Elmn[i])
+ (BlmnOrig[i] - Blmn[i]) * conj(BlmnOrig[i] - Blmn[i]))/2. ;
r_16 diff = sqrt(cdiff.real());
if(diff>err_abs){
err_abs = diff;
imax = i;
}
complex<r_8> foriCAbs = (ElmnOrig[i]*conj(ElmnOrig[i]) + BlmnOrig[i]*conj(BlmnOrig[i]))/2. ;
r_8 foriAbs = sqrt(foriCAbs.real());
r_8 relatdiff = diff/foriAbs;
if((relatdiff)>err_rel) err_rel = relatdiff;
}
cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
cout << "Err. Max. " << err_abs << " [" << imax << "], Err. Rel. " << err_rel << endl;
cout << " >>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<<" << endl;
}//SpinMultiSphericalLaguerreTransform
//-------------------------------------------------------------
void TestLaguerre2Bessel() {
......@@ -618,6 +754,7 @@ int main(int narg, char *arg[]) {
r_8 R = 1.;
int Lmax = 128;
int Pmax = 0;
int spin = 0;
int test=1;
......@@ -645,7 +782,8 @@ int main(int narg, char *arg[]) {
<< " [-ntheta <number of theta rings> [determined by the geometry]\n in case of Healpix gives the nside parameter]\n"
<< " [-nphi <number of pixel per rings> [determined by the geometry]]\n"
<< " [-p <Pmax for Fourier-Bessel> [Nmax] ]\n"
<< " [-r <Rmax> [1.]]"
<< " [-r <Rmax> [1.]]\n"
<< " [-spin <spin value> [0]"
<< endl;
return 0;
}
......@@ -653,6 +791,10 @@ int main(int narg, char *arg[]) {
alpha = atof(arg[ka+1]);
ka+=2;
}
else if (strcmp(arg[ka],"-spin")==0) {
spin = atoi(arg[ka+1]);
ka+=2;
}
else if (strcmp(arg[ka],"-n")==0) {
N=atoi(arg[ka+1]);
ka+=2;
......@@ -704,6 +846,7 @@ int main(int narg, char *arg[]) {
param.N = N;
param.Pmax = Pmax;
param.R = R;
param.spin = spin;
param.alpha = alpha;
param.geometry = geometry;
param.ntheta = ntheta;
......@@ -720,6 +863,7 @@ int main(int narg, char *arg[]) {
<< param.Lmax << ", "
<< param.N << ", "
<< param.Pmax << ", "
<< param.spin << ", "
<< param.R << ", "
<< param.alpha << "\n"
<< "Geometry: " << param.geometry << " Ntheta, Nphi: "
......@@ -778,6 +922,16 @@ int main(int narg, char *arg[]) {
}
break;
case 6:
{
tstack_push("SpinMultiSphericalLaguerreTransform");
SpinMultiSphericalLaguerreTransform();
tstack_pop("SpinMultiSphericalLaguerreTransform");
tstack_report("SpinMultiSphericalLaguerreTransform");
}
break;
default:
throw LagSHTError("EROOR Test type unkwown");
}//end of switch
......
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