+ bin spline, prolate

python/spherelib/__init__.py

@@ -25,6 +25,7 @@
 from map import *
 from spherelib_map import regress, fastfill, cat2mask, apodize_mask
 from spherelib_alm import _alm2cov
+from spherelib_bin import _make_prolate
 from fitsfunc2 import *
 from bin import *
 from alm import *

python/spherelib/alm.py

@@ -23,6 +23,54 @@ import utils
 import spherelib_alm
 import bin
 
+
+## from http://www.phys.uu.nl/~haque/computing/WPark_recipes_in_python.html
+## using numerical recipes
+def nextpow2(i): 
+    n = 2 
+    while n < i: 
+        n = n * 2 
+    return n 
+
+
+def lmax2nside(lmax, level=2):
+    """ Return a suggestion for nside given an lmax. 
+    
+    Parameters
+    ----------
+    lmax: integer; maximum order of analysis
+    level: integer, (from most "agressive" to safest) 1, 2 or 3 (default value is 2).
+
+    Returns
+    -------
+    integer, power of 2
+
+    """
+
+    ## most "agressive" (and inexact...) solution : as if the healpix was minimal (npix=dim(harmonic spaces))
+    ## nside >= (lmax+1)/sqrt(12)
+    ## npix >= (lmax+1)^2
+    if level==1:
+        nside = 2 ** nextpow2( (lmax+1)/sqrt(12.0) )
+
+    ## Healpix's advice
+    ## nside >= lmax/2 
+    ## i.e. npix >= 3 * lmax^2
+    elif level==2:
+         nside = 2 ** nextpow2( lmax/2.0 )
+
+    ## More conservative
+    ## nside >= lmax
+    elif level==3:
+        nside = 2 ** nextpow2( lmax )
+
+    else:
+        raise Exception ("level is 1, 2 or 3")
+
+    nside = max( nside, 4 ) ## because map2alm not happy with nside = 1 or 2
+
+
+
 def alm2cov (list_alm):
     """ Compute second order statistics of spherical harmonic coefficients.
     

python/spherelib/bin.py

@@ -20,6 +20,8 @@ spectral windows computations.
 
 
 from numpy import *
+import spherelib_bin
+
 
 def bin_stat (covmat, bin, nmode=None):
     """ Perform the averaging of spectral covariance matrices over flat bins. 
@@ -276,3 +278,72 @@ def unbin_cl (cq, bin):
     return squeeze(cl)
     
     
+def make_prolate (bands, rule=1, thetas=[]):
+    """ Build prolate window functions. 
+    
+    Parameters
+    ----------
+    bands: list of integers, windows limits
+    rule: integer, prolate rule between 1 (acos ( 1 - solid /(2*pi))), 2( 1 / sqrt(lmean)) ,3(1/lmean>1/(2*delta) ? 1/lmean : 1/(2*delta))
+    thetas: list of float: caps opening values (rule set to 0 in this case). 
+    
+    Returns
+    -------
+    array-like (lmax+1, nwin)
+    """
+    nwin = len(bands)-2
+    if thetas:
+        rule = 0
+        thetas = array(thetas, dtype=float32)
+    else:    
+        thetas = ones(nwin, dtype=float32)
+    bands = array(bands, dtype=int32)
+    lmax = bands[-1]
+    bins = zeros((nwin, lmax+1), dtype=float64)
+    spherelib_bin._make_prolate (bands, bins, thetas, rule)
+    return bins.transpose()
+
+
+def make_spline (bands, order=3):
+    """ Build spline window functions. 
+
+    Parameters
+    ----------
+    bands: list of integers, windows limits
+    order: integer, spline order, default is 3
+
+    Returns
+    -------
+    array-like (lmax+1, nwin)
+    """
+    order = 3
+    lmax = bands[-1]
+    nwin = len(bands)-2
+    bands = array(bands, dtype=int32)
+    bins = zeros((nwin, lmax+1), dtype=float64)
+    spherelib_bin._make_spline (bands, bins, order)
+    return bins.transpose()
+
+
+def norm_cl (cl, type="sphere"):
+    """ Return re-normalized spectra. 
+    
+    cl: array-like, cl values (lmax+1)
+    type: string. 'sphere': normalised in energy w.r.t. the nb of spherical modes
+                  'line'  : normalized in energy as an ell**2(N) function
+                  'max'   : normalized as an ell**infinity(N) function (with maximum value = 1)
+    """
+
+    lmax = len(cl)-1
+    ELL  = arange(lmax+1)
+
+    if type=='sphere':
+        norm = sum( (cl**2) * (2*ELL+1) ) / (4*pi)
+    elif type=='line':
+        norm = sum( cl**2 )
+    elif type=="max":
+        norm = max(abs(cl))**2
+    else:
+        raise Exception ("norm_cl: norm must be one of the following: \"sphere\", \"line\", \"max\" ");
+    
+    return cl*( norm**(-0.5) )

python/spherelib/myIO.py

@@ -215,7 +215,7 @@ def write_alm(file, alm):
         except pio.pioError,val:
             return str(val)
     else:
-        sp.write_alm(file, map)     
+        sp.write_alm(file, alm)     
 
 def read_cl (file):
     """ Read cl from file

python/wscript

@@ -107,10 +107,26 @@ def build(bld):
 		libpath      = [os.path.abspath('../healpix/'+target+'/lib')],
 		staticlib    = ['cxxsupport', 'cfitsio', 'fftpack', 'healpix_cxx'],
 		uselib_local = 'spherelib')
+
+	swig_spherelib_bin = bld(    
+		features = 'cxx cshlib pyext',
+		includes = ['../lib/src', '../healpix/'+target+'/include','spherelib', 
+'../include', '../lib/src'],
+		source = ['spherelib/spherelib_bin.i','spherelib/spherelib_bin.cpp'],
+		target = '_spherelib_bin',
+		defines      = ['HEALPIXDATA=\''+bld.root.abspath()+'../healpix/data/\''], 
+		swig_flags = '-c++ -python -Wall',
+		lib          = libs,
+		libpath      = [os.path.abspath('../healpix/'+target+'/lib')],
+		staticlib    = ['cxxsupport', 'cfitsio', 'fftpack', 'healpix_cxx'],
+		uselib_local = 'spherelib')
 	
 	swig_spherelib_map.install_path = '${PYTHONDIR}/spherelib'
 	swig_spherelib_alm.install_path = '${PYTHONDIR}/spherelib'
+	swig_spherelib_bin.install_path = '${PYTHONDIR}/spherelib'
 
 	obj = bld(features = 'py')
+	print out
 	obj.find_sources_in_dirs( ['./spherelib', out+'/default/spherelib'], exts=['.py'] )
+	obj.find_sources_in_dirs( [out+'/default/spherelib'], exts=['.py'] )
 	obj.install_path = '${PYTHONDIR}/spherelib'