Skip to content

Instantly share code, notes, and snippets.

@evanbiederstedt
Created July 16, 2015 15:24
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save evanbiederstedt/64207101f1cea0ab796a to your computer and use it in GitHub Desktop.
Save evanbiederstedt/64207101f1cea0ab796a to your computer and use it in GitHub Desktop.
Hogg test, array 3072, mean = 0, var = 1
CODE:
""""
#
# Create an array shaped (3000,), mean = 0.0, variance = 1.0, and compute a_lm values.
# use np.random.normal(mean, std, size)
# http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.normal.html
#
hoggarray = np.random.normal(0.0, 1.0, 3072) # mean = 0, std = 1 = var = 1
print hoggarray
print "array shape is"
print hoggarray.shape # array shape
print "median is"
print np.median(hoggarray) # median
print "mean is"
print np.mean(hoggarray) #mean
print "variance is"
print np.var(hoggarray) #variance
""""
OUTPUT:
[-0.55676395 0.48483198 -0.59881687 ..., 0.8116213 -0.10942984
-0.22426777]
array shape is
(3072,)
median is
0.00747173452011
mean is
-0.00217478806788
variance is
0.986227505318
CODE:
""""
hoggalmarr = hp.map2alm(hoggarray) # This is an array of a_lm values
print "The array of spherical harmonic coefficients a_lm is"
print hoggalmarr
print "The arr.shape is " + str(hoggalmarr.shape)
print "The length of a_lm array is " + str(len(almarr))
""""
OUTPUT:
The array of spherical harmonic coefficients a_lm is
[-0.00837408+0.j 0.04134131+0.j 0.06053477+0.j ...,
0.02215815-0.01595516j 0.06886747-0.02858377j -0.02954801+0.04785066j]
The arr.shape is (1176,)
The length of a_lm array is 1176
CODE:
""""
#
# Create an array with only the values a_3m, i.e. a_30, a_31, a_32, a_33
#
# First convert the array of alm coefficients into a real
#
realalm = hoggalmarr.real
#
print realalm[:36]
""""
OUTPUT:
[-0.00837408 0.04134131 0.06053477 -0.0705577 -0.05237246 -0.03644615
0.03338565 -0.04845627 -0.16817491 -0.00438878 0.08074669 0.05330128
0.07296344 0.03266376 0.02259408 0.03067941 -0.00937243 -0.01763272
-0.02879025 0.02725608 -0.04301382 -0.12658109 0.08010855 0.0023713
-0.02216241 -0.03100067 -0.02929653 0.02698051 0.05829383 0.03782057
0.11589424 -0.00777469 -0.17969654 -0.02823551 0.0165002 0.00345596]
CODE:
""""
empty_almlist = []
#
a30 = realalm[3]
a31 = realalm[35]
a32 = realalm[66]
a33 = realalm[96]
#
print "a30 is " + str(a30)
print "a31 is " + str(a31)
print "a32 is " + str(a32)
print "a33 is " + str(a33)
#
print str(pairs[3]) # Check with our output above
print str(pairs[35])
print str(pairs[66])
print str(pairs[96])
#
empty_almlist.append(a30)
empty_almlist.append(a31)
empty_almlist.append(a32)
empty_almlist.append(a33)
#
print empty_almlist
""""
OUTPUT:
a30 is -0.0705577002094
a31 is 0.00345595650839
a32 is -0.0167120913546
a33 is 0.0560785204204
(3, 0)
(3, 1)
(3, 2)
(3, 3)
[-0.070557700209365054, 0.0034559565083948527, -0.016712091354625973, 0.05607852042035022]
CODE:
""""
# create array of real-valued alm coefficients, a30 a31 a32 a33
realalm3 = np.asarray(empty_almlist) # np.asarray() converts input into an array
print realalm3
# Calculate (abs(alm))**2 i.e. |alm|^2
abs_alm3 = np.absolute(realalm3)
print abs_alm3
# Now calculate the squares element-wise, x**2
alm3_squared = abs_alm3**2
print alm3_squared
""""
OUTPUT:
[-0.0705577 0.00345596 -0.01671209 0.05607852]
[ 0.0705577 0.00345596 0.01671209 0.05607852]
[ 4.97838906e-03 1.19436354e-05 2.79293997e-04 3.14480045e-03]
@evanbiederstedt
Copy link
Author

Add evaluation for empirical C3
CODE
"""
anafastCl = hp.anafast(hoggarray, lmax=32)

len(anafastCl) = 33

remove monopole and dipole values, l=0, l=1

hatCl = anafastCl[2:] #len() = 31, type() = np.ndarray
hatC3 = hatCl[1] # index 0 = C2, 1 = C3, etc.
hatC4 = hatCl[2]

print "Empirical C_3 from Healpix anafast is"
print hatC3
print "The value 7_C3 is "
print 7_hatC3
"""
OUTPUT
Empirical C_3 from Healpix anafast is
0.00186721021751
The value 7*C3 is
0.0130704715226

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment