Commit e3f2e9b3 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

simplify

parent a2311484
Pipeline #70093 passed with stages
in 8 minutes and 54 seconds
......@@ -10,6 +10,14 @@ def nalm(lmax, mmax):
return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
def random_alm(lmax, mmax):
res = np.random.uniform(-1., 1., nalm(lmax, mmax)) \
+ 1j*np.random.uniform(-1., 1., nalm(lmax, mmax))
# make a_lm with m==0 real-valued
res[0:lmax+1].imag = 0.
return res
def deltabeam(lmax,kmax):
beam=np.zeros(nalm(lmax, kmax))+0j
for l in range(lmax+1):
......@@ -24,10 +32,7 @@ kmax=2 # doesn't make any sense for the beam we are using, but just for demonst
# get random sky a_lm
# the a_lm arrays follow the same conventions as those in healpy
slmT = np.random.uniform(-1., 1., nalm(lmax, mmax)) \
+ 1j*np.random.uniform(-1., 1., nalm(lmax, mmax))
# make a_lm with m==0 real-valued
slmT[0:lmax+1].imag = 0.
slmT = random_alm(lmax, mmax)
# build beam a_lm (pencil beam for now)
blmT = deltabeam(lmax,kmax)
......
......@@ -99,25 +99,20 @@ template<typename T> class Interpolator
{
double spinsign = (k==0) ? 1. : -1.;
for (size_t m=0; m<=lmax; ++m)
{
T mfac=T((m&1) ? -1.:1.);
for (size_t l=m; l<=lmax; ++l)
{
if (l<k)
a1(l,m)=a2(l,m)=0.;
else
{
complex<T> v1=slmT(l,m)*blmT(l,k),
v2=conj(slmT(l,m))*blmT(l,k)*mfac;
a1(l,m) = (v1+conj(v2)*mfac)*T(0.5*spinsign*lnorm[l]);
a1(l,m) = slmT(l,m)*blmT(l,k).real()*T(spinsign*lnorm[l]);
if (k>0)
{
complex<T> tmp = (v1-conj(v2)*mfac)*T(-spinsign*0.5*lnorm[l]);
complex<T> tmp = slmT(l,m)*blmT(l,k).imag()*T(-spinsign*lnorm[l]);
a2(l,m) = complex<T>(-tmp.imag(), tmp.real());
}
}
}
}
size_t kidx1 = (k==0) ? 0 : 2*k-1,
kidx2 = (k==0) ? 0 : 2*k;
auto quadrant=k%4;
......
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