+myio + astro

python/spherelib/__init__.py

@@ -31,3 +31,4 @@ from alm import *
 from utils import *
 from cov import *
 import myIO
+import astro

python/spherelib/astro.py

+""" astronomical units conversion tools. 
+"""
+
+## Copyright (C) 2008 APC CNRS Universite Paris Diderot <betoule@apc.univ-paris7.fr>  
+## 
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, see http://www.gnu.org/licenses/gpl.html
+
+
+from pylab import *
+import numpy
+import re
+import pyfits
+
+Jansky2SI=1.0e-26
+SI2Jansky=1.0e+26
+speedOfLight=2.99792458e8
+kBoltzmann=1.380658e-23
+sigmaStefanBoltzmann=5.67051e-8
+hPlanck=6.6260755e-34
+astronomicalUnit=1.49597893e11
+solarConstant=1368.0
+tropicalYear=3.15569259747e7
+tcmb = 2.726
+prefixes = {'n':1e-9,'u':1e-6,'m':1e-3,'k':1e3,'M':1e6,'G':1e9}
+
+def Jy(freq):
+    """ Return conversion factor from Jansky to SI at a given frequency.
+
+    Parameters
+    ----------
+    freq : scalar. Frequency value. 
+
+    Returns
+    -------
+    scalar, conversion factor value. 
+    """
+    return Jansky2SI
+
+
+def K_CMB(freq):
+    """ Return conversion factor from Kelvin CMB to SI at a given frequency.
+    
+    Parameters
+    ----------
+    freq : scalar. Frequency value. 
+
+    Returns
+    -------
+    scalar, conversion factor value. 
+    """
+    x = (hPlanck *freq)/ (kBoltzmann * tcmb)
+    ex = exp(x)
+    den = ex-1
+    den *= den
+    den = 1/den
+    fc = freq /speedOfLight 
+    return 2*kBoltzmann * fc *fc * x * x * ex * den
+
+def K_RJ(freq):
+    """ Return conversion factor from Kelvin Raleigh Jeans to SI at a given frequency.
+
+    Parameters
+    ----------
+    freq : scalar. Frequency value. 
+
+    Returns
+    -------
+    scalar, conversion factor value. 
+    """
+    return 2*kBoltzmann*freq*freq/(speedOfLight*speedOfLight)
+
+def K_KCMB(freq):
+    """ Return conversion factor from K/ KCMB to SI at a given frequency.
+
+    Parameters
+    ----------
+    freq : scalar. Frequency value. 
+
+    Returns
+    -------
+    scalar, conversion factor value. 
+    """
+    
+    return tcmb*K_CMB(freq)
+
+def y_sz(freq):
+    """ Return conversion factor from y SZ to SI at a given frequency.
+
+    Parameters
+    ----------
+    freq : scalar. Frequency value. 
+
+    Returns
+    -------
+    scalar, conversion factor value. 
+    """
+    x = (hPlanck *freq)/ (kBoltzmann * tcmb)
+    ex = exp(x)
+    den = ex-1
+    fc = freq /speedOfLight
+    den = 1/den
+    return x*ex*den*(x*(ex+1)*den - 4)*2*hPlanck*freq*den*fc*fc
+
+
+def si(freq):
+    """ Return conversion factor to SI at a given frequency.
+
+    Parameters
+    ----------
+    freq : scalar. Frequency value. 
+
+    Returns
+    -------
+    scalar, conversion factor value. 
+    """
+    return 1.0
+
+def taubeta2(freq):
+    """ Return conversion factor to taubeta2 at a given frequency.
+
+    Parameters
+    ----------
+    freq : scalar. Frequency value. 
+
+    Returns
+    -------
+    scalar, conversion factor value. 
+    """
+
+    return 1.0
+
+
+
+def convfact(freq=1.0e10, fr=r'mK_CMB',to=r'mK_CMB'):
+    """ Compute the conversion factor between two units at a given frequency in GHz.
+
+    Parameters
+    ----------
+    freq : scalar. Frequency value. 
+    fr : string as prefixunit. unit can be either 'Jy/sr', 'K_CMB', 'K_RJ', 'K/KCMB', 'y_sz', 'si', 'taubeta'
+    with optionnal prefix 'n', 'u', 'm', 'k', 'M', 'G' 
+    to : string as prefixunit. 
+
+    Returns
+    -------
+    scalar, conversion factor value. 
+    """
+
+    frpre, fru =  re.match(r'(n|u|m|k|M|G)?(Jy/sr|K_CMB|K_RJ|K/KCMB|y_sz|si|taubeta2)', fr).groups()
+    topre, tou = re.match(r'(n|u|m|k|M|G)?(Jy/sr|K_CMB|K_RJ|K/KCMB|y_sz|si|taubeta2)', to).groups()
+
+    if fru == 'Jy/sr':
+        fru = 'Jy'
+    if tou == 'Jy/sr':
+        tou = 'Jy'
+    if fru == tou:
+        fac = 1.0
+    else:
+        fac = eval(fru+'(freq)/'+tou+'(freq)')
+    
+    if not frpre is None:
+        fac *= prefixes[frpre]
+    if not topre is None:
+        fac /= prefixes[topre]
+
+    return fac
+

python/spherelib/myIO.py

+"""
+ I/O wrapper which handles both FITS and DMC format
+"""
+
+# Copyright (C) 2009 APC CNRS Universite Paris Diderot
+# <lejeune@apc.univ-paris7.fr>
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, see http://www.gnu.org/licenses/gpl.html        
+
+try:
+    import piolib as pio
+except ImportError:
+    print "Can't import piolib"
+import healpy as hp
+import fitsfunc2 as sp
+import os
+import numpy as np
+
+def __piolib(file):
+    return not(os.path.realpath(file)[-1-3:]=='fits')
+
+def save_pie_param (file, dic):
+    """ Save a python dictionnary into a pie parameter file
+
+    Parameters
+    ----------
+    file: string, pie file name
+    dic: dictionnary
+    """
+    f = open(file, "w")
+    for name, value in dic.items():
+        if isinstance(value, list):
+            f.write("number_of_"+name+" ="+str(len(value))+"\n")
+            for el in range(len(value)):
+                f.write(name+str(el+1)+' = '+str(value[el])+'\n')
+        else:
+            f.write(name+' = '+str(value)+'\n')
+    f.close()   
+    
+            
+        
+def checkgroup (groupname, PIOTYPE, *args):
+    """ Check that groupname of type PIOTYPE exists, if not create it with arguments args.
+    """
+    try:
+        exec("test_grp = pio.Open"+PIOTYPE+"Grp(groupname,\"r\")")
+        exec("pio.Close"+PIOTYPE+"Grp(test_grp)")
+    except Exception, e:
+        print "Grp " + groupname + " not found, creating it..."
+        if len(args)>0:
+            arglist = ","+",".join(args)
+        else:
+            arglist = ""
+        try:
+            exec("pio.Create"+PIOTYPE+"Grp(groupname"+arglist+")")
+        except Exception, e2:
+            print "cannot create group "+groupname+" !!!"
+                
+def ch_data_fn (grpname, filename):
+    """ Change pipelet filename into DMC object name.
+
+    Also create an empty file with DMC name into the pipelet directory
+
+    Parameters
+    ----------
+    grpname: string, DMC group name
+    filename: string, pipelet file name
+
+    Returns
+    -------
+    string, DMC file name
+    """
+    ifilename = filename
+    pipepath = os.path.dirname(filename)
+    pipepath2 = pipepath.split("/data")[-2]
+    key      = pipepath2[-1-6:]
+    filename = os.path.basename(filename)
+    filename = filename.replace(".fits", "_"+key)
+    filename = grpname+"/"+filename
+    pipename = filename.replace("/", ":")
+    ##os.system("ln -s %s %s"%(filename, ifilename))
+    os.system("touch "+pipepath+"/"+pipename)
+    filename_source = filename
+    filename_dest   = ifilename
+    return filename_source
+
+def glob_seg(glob_default, x, y):
+    """ Return the list of filename matching y in the working
+    directory of segment x.
+    
+    Parameters
+    ----------
+    x: string, segment name
+    y: string, regexp of file to glob.
+    
+    Returns
+    -------
+    list of filenames.
+    """                   
+    f = glob_default(x, y)
+    if not f:
+        z = y.replace(".fits", "")
+        z = "*:"+z
+        f = glob_default(x, z)
+        if f:
+            lst = []
+            for fi in f:
+                if len(fi.split("."))==1:
+                    lst.append(os.path.basename(fi).replace(":", "/"))
+            f = lst
+    return f
+
+def read_map(file):
+    """ Read Healpix map from file
+
+    Parameters
+    ----------
+    file: string, FITS file name or DMC object name
+
+    Returns
+    -------
+    m, array-like.
+    """
+    if __piolib(file):
+        m = pio.ReadMAPObject(file, "PIODOUBLE", "")
+    else:
+        m = hp.read_map(file)
+    return m
+
+def write_map (file, map):
+    """ Write map to file.
+
+    Parameters
+    ----------
+    file: string, FITS file name or DMC object name
+    map: array-like, Healpix map
+    
+    """
+    if __piolib(file):
+        order  = "RING"
+        coord  = "GALACTIC"
+        nside  = hp.npix2nside(len(map))
+        try:
+            grpname = pio.GetGrpName(file)
+            checkgroup (grpname, "MAP", "\""+coord+"\"", "\""+order+"\"", str(nside)) 
+            grp     = pio.OpenMAPGrp(grpname, "w")
+            val     = pio.CreateMAPObject(file,'PIODOUBLE')
+            pio.WriteMAPObject (map, file, 'PIODOUBLE','',grp)
+            pio.CloseMAPGrp(grp)
+        except pio.pioError,val: 
+            return str(val)
+    else:
+        hp.write_map(file, map)
+
+def read_alm (file):
+    """ Read Healpix alm from file
+
+    Parameters
+    ----------
+    file: string, FITS file name or DMC object name
+    
+    Returns
+    -------
+    alm, array-like.
+    """
+    if __piolib(file):
+        vec  = pio.ReadVECTObject(file, "PIODOUBLE", "")
+        lmax = np.round((np.sqrt(1+8*len(vec)/2)-3)/2);
+        alm  = np.zeros(( (lmax+1)*(lmax+2)/2, 1), dtype="complex")
+        j = 0
+        for i in range(0,len(vec),2):
+            alm[j] = complex(vec[i], vec[i+1])
+            j = j+1
+    else:
+        alm = sp.read_alm(file)
+    return np.squeeze(alm)
+    
+
+def write_alm(file, alm):
+    """ Write alm to file.
+    
+    Parameters
+    ----------
+    file: string, FITS file name or DMC object name
+    alm: array-like, Healpix alm
+    
+    """
+
+    if __piolib(file):
+        salm = len(alm)  
+        tabalm  = np.zeros(( salm*2, 1))
+        j = 0
+        for i in range(0,salm*2,2):
+            tabalm[i]   = alm[j].real
+            tabalm[i+1] = alm[j].imag
+            j = j+1   
+        try:
+            grpname = pio.GetGrpName(file)
+            checkgroup(grpname, "VECT")
+            grp     = pio.OpenVECTGrp(grpname, "w")
+            val     = pio.CreateVECTObject(file,'PIODOUBLE')
+            pio.WriteVECTObject (tabalm, file, 'PIODOUBLE','begin=0 &\n end='+str(salm*2-1),grp)
+            pio.CloseVECTGrp(grp)
+        except pio.pioError,val:
+            return str(val)
+    else:
+        sp.write_alm(file, map)     
+
+def read_cl (file):
+    """ Read cl from file
+    
+    Parameters
+    ----------
+    file: string, FITS file name or DMC object name
+    
+    Returns
+    -------
+    cl, array-like.
+    """
+    if __piolib(file):
+        cl  = pio.ReadVECTObject(file, "PIODOUBLE", "")
+    else:
+        cl = sp.read_cl(file)
+    return np.squeeze(cl)
+
+def write_cl(file, cl):
+    """ Write cl to file.
+    
+    Parameters
+    ----------
+    file: string, FITS file name or DMC object name
+    cl: array-like
+    
+    """
+    if __piolib(file):
+        try:
+            grpname = pio.GetGrpName(file)
+            checkgroup (grpname, "VECT")
+            grp     = pio.OpenVECTGrp(grpname, "w")
+            val     = pio.CreateVECTObject(file,'PIODOUBLE')
+            lmax    = len(cl)-1
+            pio.WriteVECTObject (cl, file, 'PIODOUBLE','begin=0 &\n end='+str(lmax),grp)
+            pio.CloseVECTGrp(grp)
+        except pio.pioError,val:
+            return str(val)
+    else:
+        sp.write_cl(file, cl) 
+            
+def read_cat (file):
+    """ Read cat from file
+    
+    Parameters
+    ----------
+    file: string, FITS file name or DMC object name
+    
+    Returns
+    -------
+    cat, array-like.
+    """
+    if __piolib(file):
+        cat  = pio.ReadTAB2DObject(file, "PIODOUBLE", "")
+        grpname = pio.GetGrpName(file)
+        grp  = pio.OpenTAB2DGrp(grpname, "r")   
+        ncols= pio.GetTAB2DColumnGrp(grp)
+        nrows   = len(cat)/ncols
+        cat = cat.reshape(nrows, ncols)
+    else:
+        cat = np.loadtxt(file)
+    return cat   
+
+def write_covmat(file, cov, nbmodes=None):
+    """ Write cl to file.
+    
+    Parameters
+    ----------
+    file: string, FITS file name or DMC object name
+    stats: array-like
+    nbmodes: array-like
+    """
+    if __piolib(file):
+        try:
+            grpname = pio.GetGrpName(file)
+            checkgroup (grpname, "TAB2D")
+            grp     = pio.OpenTAB2DGrp(grpname, "w")
+            val     = pio.CreateTAB2DObject(file,'PIODOUBLE')
+            nr = shape(R)[0]
+            nc = shape(R)[1]+1
+            val = pio.WriteTAB2DObject(w,file, 'PIO'+precision,"tab=0:0,0:"+str(nr-1) , grp)
+            val = pio.WriteTAB2DObject(R,file, 'PIO'+precision,"tab=1:"+str(nc-1)+",0:"+str(nr-1) , grp)
+            val     = pio.CloseTAB2DGrp(grp)
+            pio.CloseTAB2DGrp(grp)
+        except pio.pioError,val:
+            return str(val)
+        else:
+            sp.write_covmat(file, cov, nbmodes=nbmodes)    
+            
+            
+            
+