Skip to content

GitLab

Commit 3325d8d4 authored by Syl's avatar Syl
Browse files

Correction bug TB, EB, TE. Add auto spectrum

parent 58be8805
Hide whitespace changes
Inline Side-by-side
Showing with 187 additions and 194 deletions
TestEsti.py
View file @ 3325d8d4
......@@ -25,22 +25,22 @@ nside = 4
lmax = 2 * nside - 1
Slmax = 2 * nside - 1
deltal = 1
nsimu = 10000
nsimu = 1000
clth = np.array(hp.read_cl('planck_base_planck_2015_TTlowP.fits'))
Clthshape = zeros(((6,)+shape(clth)[1:]))
Clthshape[:4] = clth
clth = Clthshape
EB = 0.5
EB = 0.1
clth[4] = EB*sqrt(clth[1]*clth[2])
TB = 0.5
TB = 0.0
clth[5] = TB*sqrt(clth[0]*clth[2])
lth = arange(2, lmax+1)
spec = ['EB', 'TE', 'TB']
spec = None #['EB']
temp = True
polar = True
polar = False
corr = False
pixwin = False
pixwin = True
ellbins = arange(2, lmax + 2, deltal)
ellbins[-1] = lmax + 1
......@@ -51,7 +51,7 @@ nbins = len(ellbins) - 1
# Create mask
t, p = hp.pix2ang(nside, range(hp.nside2npix(nside)))
mask = np.ones(hp.nside2npix(nside), bool)
mask[abs(90 - rad2deg(t)) < 60] = False
# mask[abs(90 - rad2deg(t)) < 30] = False
npix = sum(mask)
fwhm = 0.5
......@@ -70,58 +70,67 @@ ipok = ip[mask]
Pl, S = libcov.compute_ds_dcb(ellbins, nside, ipok, bl, clth, Slmax, spec=spec,
pixwin=pixwin, timing=True, MC=False)
# ############## Compute spectra ###############
# ############## Initialise xqml class ###############
muKarcmin = 1.0
pixvar = muKarcmin2var(muKarcmin, nside)
varmap = ones((nstoke * npix)) * pixvar
NoiseVar = np.diag(varmap)
cmb = np.array(hp.synfast(clth, nside, fwhm=deg2rad(fwhm), pixwin=pixwin,
new=True, verbose=False, lmax=Slmax))
noise = (randn(len(varmap)) * varmap**0.5).reshape(nstoke, -1)
dm = cmb[istokes][:, mask] + noise
esti = xqml.xQML(mask, ellbins, clth, Pl=Pl, S=S, fwhm=fwhm,
spec=spec, temp=temp, polar=polar, corr=corr)
esti.construct_esti(NoiseVar, NoiseVar)
cl = esti.get_spectra(dm, dm)
V = esti.get_covariance()
esti = xqml.xQML(mask, ellbins, clth, NA=NoiseVar, NB=NoiseVar, Pl=Pl,
S=S, fwhm=fwhm, spec=spec, temp=temp, polar=polar, corr=corr)
# esti.construct_esti(NoiseVar, NoiseVar)
V = esti.get_covariance(cross=True)
Va = esti.get_covariance(cross=False)
# ############## Construct MC ###############
allcla = []
allcl = []
# allcmb = []
# esti = xqml.xQML(mask, ellbins, clth, Pl=Pl, fwhm=fwhm, spec=spec, temp=temp,
# polar=polar, corr=corr)
esti.construct_esti(NoiseVar, NoiseVar)
fpixw = extrapolpixwin(nside, lmax+2, pixwin=pixwin)
start = timeit.default_timer()
for n in np.arange(nsimu):
progress_bar(n, nsimu, timeit.default_timer() - start)
cmb = np.array(hp.synfast(clth[:, :len(fpixw)]*(fpixw*bl)**2, nside,
lmax=Slmax, fwhm=deg2rad(fwhm), new=True, verbose=False))
pixwin=False, lmax=Slmax, fwhm=deg2rad(fwhm), new=True,
verbose=False))
cmbm = cmb[istokes][:, mask]
# allcmb.append(cmbm)
dmA = cmbm + (randn(nstoke * npix) * sqrt(varmap)).reshape(nstoke, npix)
dmB = cmbm + (randn(nstoke * npix) * sqrt(varmap)).reshape(nstoke, npix)
# allcmb.append(cmbm)
allcl.append(esti.get_spectra(dmA, dmB))
allcla.append(esti.get_spectra(dmA))
figure(figsize=[10, 8])
clf()
subplot(3, 1, 1)
plot(lth, clth[ispecs][:, lth].T, '--k')
# myS = cov(np.array(allcmb).reshape(nsimu, -1), rowvar=False)
hcl = mean(allcl, 0)
scl = std(allcl, 0)
[plot(ellval, hcl[s], 'o', color='C%i' % s, label=r"$%s$" % spec[s])
hcla = mean(allcla, 0)
scla = std(allcla, 0)
figure(figsize=[12, 8])
clf()
subplot(3, 2, 1)
title("Cross")
plot(lth, clth[ispecs][:, lth].T, '--k')
[errorbar(ellval, hcl[s], scl[s], fmt='o', color='C%i' % s, label=r"$%s$" % spec[s])
for s in np.arange(nspec)]
[fill_between(ellval, (hcl - scl/sqrt(nsimu))[s], (hcl + scl/sqrt(nsimu))[s],
color='C%i' % s, alpha=0.2) for s in np.arange(nspec)]
semilogy()
ylabel(r"$C_\ell$")
legend(loc=4, frameon=True)
subplot(3, 2, 2)
title("Auto")
plot(lth, clth[ispecs][:, lth].T, '--k')
[errorbar(ellval, hcla[s], scla[s], fmt='o', color='C%i' % s, label=r"$%s$" %
spec[s]) for s in np.arange(nspec)]
semilogy()
legend(loc=4)
subplot(3, 1, 2)
cosmic = sqrt(2./(2 * lth + 1)) / mean(mask) * clth[ispecs][:, lth]
subplot(3, 2, 3)
# cosmic = sqrt(2./(2 * lth + 1)) / mean(mask) * clth[ispecs][:, lth]
# plot(lth, cosmic.transpose(), '--k')
[plot(ellval, scl[s], '--', color='C%i' % s, label=r"$\sigma^{%s}_{\rm MC}$" %
spec[s]) for s in np.arange(nspec)]
......@@ -129,18 +138,33 @@ cosmic = sqrt(2./(2 * lth + 1)) / mean(mask) * clth[ispecs][:, lth]
for s in np.arange(nspec)]
ylabel(r"$\sigma(C_\ell)$")
semilogy()
legend(loc=4)
# legend(loc=4, frameon=True)
subplot(3, 2, 4)
[plot(ellval, scla[s], ':', color='C%i' % s, label=r"$\sigma^{%s}_{\rm MC}$" %
spec[s]) for s in np.arange(nspec)]
[plot(ellval, sqrt(diag(Va)).reshape(nspec, -1)[s], 'o', color='C%i' % s)
for s in np.arange(nspec)]
semilogy()
subplot(3, 1, 3)
subplot(3, 2, 5)
[plot(ellval, (hcl[s]-clth[ispecs[s]][lth])/(scl[s]/sqrt(nsimu)), '--o',
color='C%i' % s) for s in np.arange(nspec)]
ylabel(r"$R[C_\ell]$")
xlabel(r"$\ell$")
ylim(-3, 3)
grid()
subplot(3, 2, 6)
[plot(ellval, (hcla[s]-clth[ispecs[s]][lth])/(scla[s]/sqrt(nsimu)), '--o',
color='C%i' % s) for s in np.arange(nspec)]
xlabel(r"$\ell$")
ylim(-3, 3)
grid()
show()
# savefig("../Plots/Git/"+"test0.png")
savefig("../Plots/Git/"+"test0.png")
if __name__ == "__main__":
"""
......
estimators.py
View file @ 3325d8d4
......@@ -336,21 +336,21 @@ def ClQuadEstimator(invW, y):
return Cl
# def biasQuadEstimator(NoiseN, El):
# """
# Compute bias term bl such that Cl = Fll^-1 . ( yl + bias)
def biasQuadEstimator(NoiseN, El):
"""
Compute bias term bl such that Cl = Fll^-1 . ( yl + bias)
# Parameters
# ----------
# NoiseN : ???
# ???
Parameters
----------
NoiseN : ???
???
# Returns
# ----------
# ???
# """
# lrange = np.arange((len(El)))
# return np.array([np.sum(NoiseN * El[l]) for l in lrange])
Returns
----------
???
"""
lrange = np.arange((len(El)))
return np.array([np.sum(NoiseN * El[l]) for l in lrange])
# def blEstimatorFlat(NoiseN, El):
......
libcov.py
View file @ 3325d8d4
......@@ -201,7 +201,6 @@ def compute_PlS(
fpixwin = extrapolpixwin(nside, Slmax+1, pixwin)
norm = (2*ll[2:]+1)/(4.*np.pi)*(fpixwin[2:]**2)*(bl[2:Slmax+1]**2)
clthn = clth[:, 2: Slmax+1]
masks = []
for i in np.arange(nbins):
masks.append((ll[2:] >= minell[i]) & (ll[2:] <= maxell[i]))
......@@ -234,38 +233,39 @@ def compute_PlS(
cos_chi = allcosang[i, j]
# Tegmark version
Q22 = norm*F1l2(cos_chi, Slmax)[2:]
Q22 = norm * F1l2(cos_chi, Slmax)
# # /!\ signe - !
R22 = -norm*F2l2(cos_chi, Slmax)[2:]
R22 = -norm * F2l2(cos_chi, Slmax)
# # Matt version
# d20 = dlss(cos_chi, 2, 0, Slmax+1)
# d2p2 = dlss(cos_chi, 2, 2, Slmax+1)
# d2m2 = dlss(cos_chi, 2, -2, Slmax+1)
# P02 = -d20[2:]
# Q22 = ( d2p2 + d2m2 )[2:]/2.
# R22 = ( d2p2 - d2m2 )[2:]/2.
if TE or TB:
P02 = -norm*F1l0(cos_chi, Slmax)[2:]
# # Matt version
# d20 = -dlss(cos_chi, 2, 0, Slmax)[2:]
# P02 = norm * d20
P02 = -norm * F1l0(cos_chi, Slmax)
elemA = 0
elemB = 0
elemC = 0
elemD = 0
if TE:
elemTE = P02*clthn[3]
elemA += np.sum(cji*elemTE)
elemB -= np.sum(sji*elemTE)
elemC += np.sum(cij*elemTE)
elemD -= np.sum(sij*elemTE)
elemTE = np.sum(P02*clthn[3])
elemA += cji*elemTE
elemB -= sji*elemTE
elemC += cij*elemTE
elemD -= sij*elemTE
if TB:
elemTB = P02*clthn[5]
elemA += np.sum(sji*elemTB)
elemB += np.sum(cji*elemTB)
elemC += np.sum(sij*elemTB)
elemD += np.sum(cij*elemTB)
elemTB = np.sum(P02*clthn[5])
elemA += sji*elemTB
elemB += cji*elemTB
elemC += sij*elemTB
elemD += cij*elemTB
S[i, jj] = elemA
S[i, jj+npi] = elemB
......@@ -464,7 +464,7 @@ def compute_S(
progress_bar(i, npi, -0.5 * (start-timeit.default_timer()))
for j in np.arange(i, npi):
if temp:
pl = norm*pl0(allcosang[i, j], Slmax)[2:]
pl = norm*pl0(allcosang[i, j], Slmax)
elem = np.sum((pl * clthn[0]))
S[i, j] = elem
S[j, i] = elem
......@@ -477,9 +477,9 @@ def compute_S(
cos_chi = allcosang[i, j]
# Tegmark version
Q22 = norm*F1l2(cos_chi, Slmax)[2:]
Q22 = norm*F1l2(cos_chi, Slmax)
# # /!\ signe - !
R22 = -norm*F2l2(cos_chi, Slmax)[2:]
R22 = -norm*F2l2(cos_chi, Slmax)
# # Matt version
# d20 = dlss(cos_chi, 2, 0, Slmax+1)
......@@ -490,7 +490,7 @@ def compute_S(
# R22 = ( d2p2 - d2m2 )[2:]/2.
if TE or TB:
P02 = -norm*F1l0(cos_chi, Slmax)[2:]
P02 = -norm*F1l0(cos_chi, Slmax)
elemA = 0
elemB = 0
elemC = 0
......@@ -1010,7 +1010,7 @@ def F1l0(z, lmax):
0.20392
"""
if abs(z) == 1.0:
return(np.zeros(lmax + 1))
return(np.zeros(lmax - 1))
else:
ell = np.arange(2, lmax + 1)
thepl = pl0(z, lmax)
......@@ -1020,7 +1020,7 @@ def F1l0(z, lmax):
a0 = 2.0 / np.sqrt((ell - 1) * ell * (ell + 1) * (ell + 2))
a1 = ell * z * theplm1 / (1 - z**2)
a2 = (ell / (1 - z**2) + ell * (ell - 1) / 2) * thepl
bla = np.append([0, 0], a0 * (a1 - a2))
bla = a0 * (a1 - a2)
return bla
......@@ -1047,10 +1047,10 @@ def F1l2(z, lmax):
0.58396
"""
if z == 1.0:
return np.append(np.zeros(2), np.ones(lmax - 1) * 0.5)
return np.ones(lmax - 1) * 0.5
elif z == -1.0:
ell = np.arange(lmax + 1)
return np.append(np.zeros(2), 0.5 * (-1)**ell[2:])
return 0.5 * (-1)**ell[2:]
else:
ell = np.arange(2, lmax + 1)
thepl2 = pl2(z, lmax)
......@@ -1060,7 +1060,7 @@ def F1l2(z, lmax):
a0 = 2.0 / ((ell - 1) * ell * (ell + 1) * (ell + 2))
a1 = (ell + 2) * z * theplm1_2 / (1 - z**2)
a2 = ((ell - 4) / (1 - z**2) + ell * (ell - 1) / 2) * thepl2
bla = np.append([0, 0], a0 * (a1 - a2))
bla = a0 * (a1 - a2)
return bla
......@@ -1086,10 +1086,10 @@ def F2l2(z, lmax):
0.34045
"""
if z == 1.0:
return np.append(np.zeros(2), -0.5 * np.ones(lmax - 1))
return -0.5 * np.ones(lmax - 1)
elif z == -1.0:
ell = np.arange(lmax + 1)
return np.append(np.zeros(2), 0.5 * (-1)**ell[2:])
return 0.5 * (-1)**ell[2:]
else:
ell = np.arange(2, lmax + 1)
thepl2 = pl2(z, lmax)
......@@ -1099,34 +1099,55 @@ def F2l2(z, lmax):
a0 = 4.0 / ((ell - 1) * ell * (ell + 1) * (ell + 2) * (1 - z**2))
a1 = (ell + 2) * theplm1_2
a2 = (ell - 1) * z * thepl2
bla = np.append([0, 0], a0 * (a1 - a2))
bla = a0 * (a1 - a2)
return bla
def dlss(z, s1, s2, lmax):
'''
Matt version
'''
d = np.zeros((lmax+1))
if s1 < abs(s2):
print("error spins, s1<|s2|")
'''
Matt version
'''
d = np.zeros((lmax+1))
if s1 < abs(s2):
print("error spins, s1<|s2|")
return
# sign = -1 if (s1 + s2) and 1 else 1
sign = (-1)**(s1-s2)
d[s1] = sign / 2.**s1 * math.sqrt(math.factorial(2.*s1)/math.factorial(
1.*s1 + s2)/math.factorial(1.*s1-s2)) * (1.+z)**(.5*(s1+s2)) * (
1.-z)**(.5 * (s1-s2))
l1 = s1+1.
rhoSSL1 = math.sqrt((l1*l1 - s1*s1) * (l1*l1 - s2*s2)) / l1
d[s1+1] = (2*s1+1.)*(z-s2/(s1+1.)) * d[s1] / rhoSSL1
for l in np.arange(s1+1, lmax, 1):
l1 = l+1.
rhoSSL = math.sqrt((l*l*1. - s1*s1) * (l*l*1. - s2*s2)) / (l*1.)
rhoSSL1 = math.sqrt((l1*l1 - s1*s1) * (l1*l1 - s2*s2)) / l1
d[l+1] = ((2.*l+1.)*(
z-s1*s2/(l*1.)/l1)*d[l] - rhoSSL*d[l-1]) / rhoSSL1
return d
# sign = -1 if (s1 + s2) and 1 else 1
sign = (-1)**(s1-s2)
d[s1] = sign / 2.**s1 * math.sqrt(math.factorial(2.*s1)/math.factorial(1.*s1+s2)/math.factorial(1.*s1-s2)) * (1.+z)**(.5*(s1+s2)) * (1.-z)**(.5*(s1-s2))
l1 = s1+1.
rhoSSL1 = math.sqrt( (l1*l1 - s1*s1) * (l1*l1 - s2*s2) ) / l1
d[s1+1] = (2*s1+1.)*(z-s2/(s1+1.)) * d[s1] / rhoSSL1
for l in np.arange(s1+1, lmax, 1) : #( l=s1+1; l<lmax; l++) {
l1 = l+1.
rhoSSL = math.sqrt( (l*l*1.- s1*s1) * (l*l*1.- s2*s2) ) / (l*1.)
rhoSSL1 = math.sqrt( (l1*l1 - s1*s1) * (l1*l1 - s2*s2) ) / l1
d[l+1] = ( (2.*l+1.)*(z-s1*s2/(l*1.)/l1)*d[l] - rhoSSL*d[l-1] ) / rhoSSL1
return d
# def dlss(z, s1, s2, lmax):
# '''
# Matt version
# '''
# d = np.zeros((lmax+1))
# if s1 < abs(s2):
# print("error spins, s1<|s2|")
# return
# # sign = -1 if (s1 + s2) and 1 else 1
# sign = (-1)**(s1-s2)
# d[s1] = sign / 2.**s1 * math.sqrt(math.factorial(2.*s1)/math.factorial(
# 1.*s1 + s2)/math.factorial(1.*s1-s2)) * (1.+z)**(.5*(s1+s2)) * (
# 1.-z)**(.5 * (s1-s2))
# l1 = s1+1.
# rhoSSL1 = math.sqrt((l1*l1 - s1*s1) * (l1*l1 - s2*s2)) / l1
# d[s1+1] = (2*s1+1.)*(z-s2/(s1+1.)) * d[s1] / rhoSSL1
# for l in np.arange(s1+1, lmax, 1):
# l1 = l+1.
# rhoSSL = math.sqrt((l*l*1. - s1*s1) * (l*l*1. - s2*s2)) / (l*1.)
# rhoSSL1 = math.sqrt((l1*l1 - s1*s1) * (l1*l1 - s2*s2)) / l1
# d[l+1] = ((2.*l+1.)*(
# z-s1*s2/(l*1.)/l1)*d[l] - rhoSSL*d[l-1]) / rhoSSL1
# return d
if __name__ == "__main__":
......
xqml.py
View file @ 3325d8d4
......@@ -17,19 +17,20 @@ import healpy as hp
import random as rd
from libcov import compute_ds_dcb
from simulation import getstokes
from estimators import El
from estimators import CovAB
from estimators import CrossGisherMatrix
from estimators import CrossWindowFunction
from estimators import yQuadEstimator, ClQuadEstimator
from estimators import biasQuadEstimator
class xQML(object):
""" Main class to handle the spectrum estimation """
def __init__(
self, mask, bins, clth, lmax=None, Pl=None, S=None,
fwhm=0., spec=None, temp=False, polar=True, corr=False,
self, mask, bins, clth, NA=None, NB=None, lmax=None, Pl=None,
S=None, fwhm=0., spec=None, temp=False, polar=True, corr=False,
pixwin=True):
"""
Parameters
......@@ -56,6 +57,10 @@ class xQML(object):
If True, applies pixel window function to spectra. Default: True
"""
self.bias = None
self.cross = NB is not None
self.NA = NA
self.NB = NB if self.cross else NA
# Number of pixels in the mask
# For example that would be good to have an assertion
# on the mask size, just to check that it corresponds to a valid nside.
......@@ -63,7 +68,6 @@ class xQML(object):
# Map resolution (healpix)
self.nside = hp.npix2nside(npixtot)
# ipok are pixel indexes outside the mask
self.ipok = np.arange(npixtot)[np.array(mask, bool)]
self.mask = mask
......@@ -73,7 +77,7 @@ class xQML(object):
self.ellbins = bins
# Maximum multipole based on nside (rule of thumb to avoid aliasing)
self.Slmax = 3 * self.nside - 1 if lmax is None else lmax
self.Slmax = 2 * self.nside - 1 if lmax is None else lmax
# Beam 2pt function (Gaussian)
self.bl = hp.gauss_beam(np.deg2rad(fwhm), lmax=self.Slmax+1)
......@@ -82,7 +86,7 @@ class xQML(object):
# For example that would be good to assert that the user
# set at least polar or temp to True.
self.stokes, self.spec, self.istokes, self.ispecs = self._getstokes(
self.stokes, self.spec, self.istokes, self.ispecs = getstokes(
spec, temp, polar, corr)
self.nstokes = len(self.stokes)
self.nspec = len(self.spec)
......@@ -100,13 +104,15 @@ class xQML(object):
else:
self.Pl = Pl
if S is None:
self.S = self._CorrelationMatrix(clth)
self.S = self._SignalCovMatrix(clth)
else:
self.S = S
if S is not None:
self.S = S
def construct_esti(self, NA, NB):
self.construct_esti(NA=NA, NB=NB)
def construct_esti(self, NA, NB=None):
"""
Compute the inverse of the datasets pixel covariance matrices C,
the quadratic matrix parameter E, and inverse of the window
......@@ -120,19 +126,24 @@ class xQML(object):
Noise covariance matrix of dataset B
"""
self.cross = NB is not None
self.NA = NA
self.NB = NB if self.cross else NA
# Invert (signalA + noise) matrix
self.invCa = linalg.inv(self.S + NA)
self.invCa = linalg.inv(self.S + self.NA)
# Invert (signalB + noise) matrix
self.invCb = linalg.inv(self.S + NB)
self.invCb = linalg.inv(self.S + self.NB)
# Compute E using Eq...
self.E = El(self.invCa, self.invCb, self.Pl)
if not self.cross:
self.bias = biasQuadEstimator(self.NA, self.E)
# Finally compute invW by inverting...
self.invW = linalg.inv(CrossWindowFunction(self.E, self.Pl))
def get_spectra(self, map1, map2):
def get_spectra(self, mapA, mapB=None):
"""
Return the unbiased spectra
......@@ -150,92 +161,25 @@ class xQML(object):
"""
# Should compute auto-spectra if map2 == None ?
# Define conditions based on the map size
cond_size1 = np.size(map1) == self.nstokes * self.npix
cond_size2 = np.size(map2) == self.nstokes * self.npix
d1 = map1 if cond_size1 else map1[self.istokes, self.mask]
d2 = map2 if cond_size2 else map2[self.istokes, self.mask]
# yl is...
yl = yQuadEstimator(d1.ravel(), d2.ravel(), self.E)
# cl is obtained using...
cl = ClQuadEstimator(self.invW, yl)
self.cross = mapB is not None
cond_sizeA = np.size(mapA) == self.nstokes * self.npix
dA = mapA if cond_sizeA else mapA[self.istokes, self.mask]
if self.cross:
cond_sizeB = np.size(mapB) == self.nstokes * self.npix
dB = mapB if cond_sizeB else mapB[self.istokes, self.mask]
yl = yQuadEstimator(dA.ravel(), dB.ravel(), self.E)
cl = ClQuadEstimator(self.invW, yl)
else:
if self.bias is None:
self.bias = biasQuadEstimator(self.NA, self.E)
yl = yQuadEstimator(dA.ravel(), dA.ravel(), self.E) - self.bias
cl = ClQuadEstimator(self.invW, yl)
# Return the reshaped set of cls
return cl.reshape(self.nspec, -1)
def _getstokes(self, spec=None, temp=False, polar=False, corr=False):
"""
Get the Stokes parameters number and name(s)
Parameters
----------
spec : bool
If True, get Stokes parameters for polar (default: True)
polar : bool
If True, get Stokes parameters for polar (default: True)
temp : bool
If True, get Stokes parameters for temperature (default: False)
corr : bool
If True, get Stokes parameters for EB and TB (default: False)
Returns
----------
stokes : list of string
Stokes variables names
spec : int
Spectra names
istokes : list
Indexes of power spectra
Example
----------
>>> getstokes(polar=True, temp=False, corr=False)
(['Q', 'U'], ['EE', 'BB'], [1, 2])
>>> getstokes(polar=True, temp=True, corr=False)
(['I', 'Q', 'U'], ['TT', 'EE', 'BB', 'TE'], [0, 1, 2, 3])
>>> getstokes(polar=True, temp=True, corr=True)
(['I', 'Q', 'U'], ['TT', 'EE', 'BB', 'TE', 'EB', 'TB'], [0, 1, 2, 3, 4,
5])
"""
if spec is not None:
_temp = "TT" in spec or "TE" in spec or "TB" in spec or temp
_polar = "EE" in spec or "BB" in spec or "TE" in spec or "TB" in \
spec or "EB" in spec or polar
_corr = "TE" in spec or "TB" in spec or "EB" in spec or corr
else:
_temp = temp
_polar = polar
_corr = corr
speclist = []
if _temp or (spec is None and corr):