u

src/yapsut/__init__.py

@@ -17,6 +17,9 @@ from .periodic_stats import periodic_stats, periodic_centroid
 from .stats import CumulativeOfData
 from .colored_noise import gaussian_colored_noise
 
+# planck legacy is under editing
+#from .planck_legacy import cosmological_dipole
+
 #from import .grep
 #import .intervalls
 

src/yapsut/planck_legacy/cgs_blackbody.py

+__DESCRIPTION__=""" 
+CGS       = an object returning the CGS units and constants
+BlackBody = an object returning black body fluxes and conversion facilities
+WMAP_CMB  = a specialized blackbody for CMB using WMAP parameters
+CMB       = a specialized blackbody for CMB using LFI parameters
+
+See: 
+   import blackbody 
+   help(blackbody.CGS)
+   help(blackbody.BlackBody)
+   help(blackbody.CMB)
+
+The package manages numpy.array().
+
+Version 1.1 - 2011 Set 1 - 2012 Apr 5 - 2017 Jun 29 - M.Maris
+"""
+
+class Singleton(object):
+  """
+  Class implementing the Singleton design pattern. There are several ways to implement
+  the singleton pattern in python. The most elegant way is to use the decorators but such
+  construct has been introduced only starting from python 2.6
+  """
+  def __new__(cls, *args, **kwds):
+    if not '_the_instance' in cls.__dict__:
+      cls._the_instance =  object.__new__(cls)
+      cls._the_instance.init(*args, **kwds)
+    return cls._the_instance
+  def init(self, *args, **kwds):
+    pass
+
+class __physical_parameters_cgs(Singleton) :
+   def init(self) :
+      self.K = 1.380658e-16    # erg/K
+      self.c = 2.99792458e10   # cm/sec
+      self.h = 6.6260755e-27   # erg / sec i.e. erg Hz
+      self.flux_Jansky = 1e23  # (cm2 sec cm sterad)/erg * Jansky
+      self.ScaleToMJ = 1.e17   # 1.e23 Jy/erg cm2 sec * 1e-6 MJy/Jy
+   def keys(self) :
+      return self.__dict__.keys()
+   def __call__(self,name) :
+      units = {}
+      units['K'] = 'erg/K'
+      units['c'] = 'cm/sec'
+      units['h'] = 'erg/sec'
+      units['flux_Jansky'] = '(cm2 sec cm sterad)/erg * Jansky'
+      units['ScaleToMJ'] = '1.e23 Jy/erg cm2 sec * 1e-6 MJy/Jy'
+      units['TCMB'] = 'K'
+      return "%s = %e : %s "%(name,self.__dict__[name],units[name])
+CGS=__physical_parameters_cgs()
+
+#class __physical_parameters_mks(Singleton) :
+   #def init(self) :
+      #from release_setup import RELEASE
+      #self.K = RELEASE['K_BOLTZMAN']  # J/K
+      #self.c = RELEASE['SPEED_OF_LIGHT']   # m/sec
+      #self.h = RELEASE['H_PLANCK']   # j / sec i.e. erg Hz
+      #self.flux_Jansky = 1e26  # (cm2 sec cm sterad)/erg * Jansky
+      ##self.ScaleToMJ = 1.e17   # 1.e23 Jy/erg m2 sec * 1e-6 MJy/Jy
+   #def keys(self) :
+      #return self.__dict__.keys()
+   #def __call__(self,name) :
+      #units = {}
+      #units['K'] = 'J/K'
+      #units['c'] = 'm/sec'
+      #units['h'] = 'J/sec'
+      #units['flux_Jansky'] = '(m2 sec m sterad)/J * Jansky'
+      #units['ScaleToMJ'] = '1.e26 Jy/J m2 sec * 1e-6 MJy/Jy'
+      #units['TCMB'] = 'K'
+      #return "%s = %e : %s "%(name,self.__dict__[name],units[name])
+#MKS=__physical_parameters_mks()
+
+class __black_body(Singleton) :
+   def init(self) : 
+      self.ScaleToMJ = CGS.ScaleToMJ
+      return
+   def __call__(self,FreqGHz,T) :
+      import numpy as np
+      return self.bbn_cgs(FreqGHz*1e9,T)*1e23*self.valid(FreqGHz,T)
+   def valid(self,FreqGHz,T) :
+      return (FreqGHz >= 0.) * (T >= 0.)
+   def bbl_cgs(self,Lambda,T) :
+      """
+!
+! BB(lambda,T) thermal radiance function cgs
+!
+! lambda in cm
+!
+! bbl_cgs = erg/(cm2 sec cm sterad)
+!
+"""
+      import numpy as np
+      FT = Lambda*Lambda*Lambda*Lambda*Lambda
+      FT = 2.*CGS.h*CGS.c*CGS.c/FT;
+      ET = np.exp( CGS.h*CGS.c/(Lambda*CGS.K*T) );
+      ET = ET - 1.
+      return FT / ET
+   def bbn_cgs(self,nu,T) :
+      """
+!
+! BB(nu,T) thermal radiance function cgs
+!
+! nu in Hz
+!
+! bbn_cgs = erg/(cm2 sec cm sterad)
+!
+      """
+      import numpy as np
+      FT = nu*nu*nu
+      ET = CGS.c*CGS.c
+      FT = 2.*CGS.h*FT/ET
+      ET = CGS.h*nu/(CGS.K*T)
+      ET = np.exp(ET)
+      ET = ET - 1
+      return FT / ET
+   def bbl_rj_cgs(self,Lambda,T) :
+      """
+!
+! BB thermal radiance function cgs in Rayleight-Jeans approx
+!
+! lambda in cm
+! T in K
+! rj_cgs = erg/(cm2 sec cm sterad)
+!
+      """
+      FT = Lambda*Lambda*Lambda*Lambda
+      return 2.*CGS.K*T*CGS.c/FT 
+   def bbn_rj_cgs(self,nu,T) :
+      """
+!
+! BB thermal radiance function cgs in Rayleight-Jeans approx
+!
+! nu in Hz
+! T in K
+!
+! rj_cgs = erg/(cm2 sec Hz sterad)
+!
+      """
+      FT = (nu/CGS.c)
+      FT = FT*FT
+      return 2.* CGS.K*T*FT
+   def Krj2MJysr(self,nu_ghz,Hz=False) :
+      """generates the conversion factor from K_rj to MJy/sr for a given frequency in GHz (Hz=True to pass it in Hz)"""
+      if Hz :
+         return self.bbn_rj_cgs(nu_ghz,1.)*CGS.ScaleToMJ
+      return self.bbn_rj_cgs(nu_ghz*1e9,1.)*CGS.ScaleToMJ
+   def Tb(self,nu,Bn,FreqGHz=False,MJySr=False) :
+      """returns for a given CGS brightness the correspondiing temperature"""
+      # B=2*h*n^3/c^2 /(exp(h*nu/K/T)-1)
+      # 1/(log( 1/(B/(2*h*nu^3/c^2)+1 )/(h*nu)*K)
+      import numpy as np
+      if FreqGHz : 
+         f=nu*1e9
+      else :
+         f=nu
+      if MJySr :
+         B=Bn/CGS.ScaleToMJ
+      else :
+         B=Bn
+      a=CGS.h*f/(CGS.K*np.log((2*CGS.h*f**3/CGS.c**2)/B+1))
+      return a
+   def bbn_diff(self,nu_ghz,T,Hz=False,MJySr=True) :
+      """
+%
+% [bbn_diff,bbn_diff_ratio]=bbn_diff_mks(nu_hz,T)
+%
+% bbn_diff =  derivative of the BB thermal radiance in mks
+%             for a temperature change
+%
+% bbn_diff_ratio = bbn_diff/bbn
+%
+% nu_hz in Hz
+% T in K
+%
+% bbn_diff_mks = erg/(cm2 sec Hz sterad)/K
+% bbn_diff_ratio = 1/K
+%
+% 1 Jy/sterad = 1e-23 erg/(cm2 sec Hz sterad)
+%
+% the derivative is defines as:
+%
+%     d(Bnu(T))/dT = Bnu(T) X exp(X)/(exp(X)-1) / T
+%
+% the variation is defined as
+%
+%     DeltaBnu(T) = Bnu(T) X exp(X)/(exp(X)-1) DeltaT/ T
+%
+      """
+      import numpy as np
+      if not Hz :
+         _nu=nu_ghz*1e9
+      else :
+         _nu=nu_ghz
+      if MJySr :
+         cf = CGS.ScaleToMJ
+      else :
+         cf = 1.
+      FT = 2.*CGS.h*(_nu**3)/CGS.c**2;
+      X =  CGS.h*_nu/(CGS.K*T);
+      ET = np.exp(X);
+      Lbbn = FT / (ET - 1);
+      FACT = X*ET/(ET-1)/T;
+      return {'bbn_diff':cf*Lbbn*FACT,'bbn_diff_ratio':cf*FACT,'bbn':cf*Lbbn}
+BlackBody=__black_body()
+
+class __cmb_wmap(Singleton) :
+   """CMB parameters for WMAP"""
+   def init(self) :
+      import numpy as np
+      self.Tcmb=2.72548 #Fixsen, D. J. 1999, Volume 707, Issue 2, pp. 916-920 (2009) #2.725 # K 
+      self.DeltaTDipole = 3.3e-3; #K
+      self.sqC2 = np.sqrt(211)*1e-6/np.sqrt(2*(2+1)/(2*np.pi)); #K sqrt of Cl derived from quadrupole moment = l*(l  +1)/(2pi)*Cl
+      self.sqC3 = np.sqrt(1041)*1e-6/np.sqrt(3*(3+1)/(2*np.pi)); #K
+      self.sqC100 = np.sqrt(3032.9299)*1e-6/np.sqrt(100*(100+1)/(2*np.pi)); #K
+      self.source = "Source: WMAP"
+      self.ScaleToMJ=CGS.ScaleToMJ
+   def __call__(self,FreqGHz,MJySr=True) :
+      if MJySr :
+         return BlackBody.bbn_cgs(FreqGHz*1e9,self.Tcmb)*1e23/1e6
+      return BlackBody.bbn_cgs(FreqGHz*1e9,self.Tcmb)
+      #return "Source: WMAP"
+   def fluctuations(self,FreqGHz) :
+      """Return a dictionary with the list of most important cmb components"""
+      D=BlackBody.bbn_diff(FreqGHz*1e9,self.Tcmb,Hz=True);
+      cmbf={}
+      cmbf['FreqGHz'] = FreqGHz;
+      cmbf['Inu'] = D['bbn']/1e-23/1e6;
+      cmbf['dInu_dT'] = D['bbn_diff']/1e-23/1e6;
+      cmbf['dlogInu_dT'] = D['bbn_diff_ratio'];
+      cmbf['Dipole'] = cmbf['dInu_dT']*self.DeltaTDipole;
+      cmbf['sqC2'] = cmbf['dInu_dT']*self.sqC2;
+      cmbf['sqC3'] = cmbf['dInu_dT']*self.sqC3;
+      cmbf['sqC100'] = cmbf['dInu_dT']*self.sqC100;
+      return cmbf
+   def Kcmb2MJysr(self,FreqGHz,TKCMB) :
+      """Converts Kcmb to MJy/sr"""
+      c=BlackBody.bbn_diff(FreqGHz,self.Tcmb)['bbn_diff']
+      return c*TKCMB
+   def MJysr2Kcmb(self,FreqGHz,IMJY_Sr) :
+      """Converts MJy/sr to Kcmb"""
+      c=BlackBody.bbn_diff(FreqGHz,self.Tcmb)['bbn_diff']
+      return IMJY_Sr/c
+   def etaDeltaT(self,FreqGHz) :
+      """Computes EtaDeltaT(nu) (see Eq.(34) of Zacchei et al. (2011), A&A, 536, A5) 
+      defined as (\partial B(\nu,T) / \partial T)_{T=Tcmb} / (2 K \nu^2/c^2)
+      """
+      import numpy as np
+      X=CGS.h*FreqGHz*1e9/(CGS.K*self.Tcmb)
+      eX=np.exp(X)
+      return eX*(X/(eX-1.))**2
+WMAP_CMB=__cmb_wmap()
+
+class __cmb(Singleton) :
+   """CMB parameters for LFI for the current release"""
+   def init(self) :
+      import numpy as np
+      from release_setup import RELEASE
+      self.Tcmb,self.DeltaTDipole =RELEASE.LFI_CMB_PARAMS()
+      self.sqC2 = np.sqrt(211)*1e-6/np.sqrt(2*(2+1)/(2*np.pi)); #K sqrt of Cl derived from quadrupole moment = l*(l  +1)/(2pi)*Cl
+      self.sqC3 = np.sqrt(1041)*1e-6/np.sqrt(3*(3+1)/(2*np.pi)); #K
+      self.sqC100 = np.sqrt(3032.9299)*1e-6/np.sqrt(100*(100+1)/(2*np.pi)); #K
+      self.source = "Source: LFI %s and WMAP for multipoles"%RELEASE.release
+      self.ScaleToMJ=CGS.ScaleToMJ
+   def __call__(self,FreqGHz,MJySr=True) :
+      if MJySr :
+         cf = CGS.ScaleToMJ
+      else :
+         cf = 1.
+      #if MJySr :
+         #return BlackBody.bbn_cgs(FreqGHz*1e9,self.Tcmb)*1e23/1e6
+      return BlackBody.bbn_cgs(FreqGHz*1e9,self.Tcmb)*cf
+   def fluctuations(self,FreqGHz) :
+      """Return a dictionary with the list of most important cmb components"""
+      D=BlackBody.bbn_diff(FreqGHz*1e9,self.Tcmb,Hz=True);
+      cmbf={}
+      cmbf['FreqGHz'] = FreqGHz;
+      cmbf['Inu'] = D['bbn']/1e-23/1e6;
+      cmbf['dInu_dT'] = D['bbn_diff']/1e-23/1e6;
+      cmbf['dlogInu_dT'] = D['bbn_diff_ratio'];
+      cmbf['Dipole'] = cmbf['dInu_dT']*self.DeltaTDipole;
+      cmbf['sqC2'] = cmbf['dInu_dT']*self.sqC2;
+      cmbf['sqC3'] = cmbf['dInu_dT']*self.sqC3;
+      cmbf['sqC100'] = cmbf['dInu_dT']*self.sqC100;
+      return cmbf
+   def Kcmb2MJysr(self,FreqGHz,TKCMB) :
+      """Converts Kcmb to MJy/sr"""
+      c=BlackBody.bbn_diff(FreqGHz,self.Tcmb)['bbn_diff']
+      return c*TKCMB
+   def MJysr2Kcmb(self,FreqGHz,IMJY_Sr) :
+      """Converts MJy/sr to Kcmb"""
+      c=BlackBody.bbn_diff(FreqGHz,self.Tcmb)['bbn_diff']
+      return IMJY_Sr/c
+   def etaDeltaT(self,FreqGHz) :
+      """Computes EtaDeltaT(nu) (see Eq.(34) of Zacchei et al. (2011), A&A, 536, A5) 
+      defined as (\partial B(\nu,T) / \partial T)_{T=Tcmb} / (2 K \nu^2/c^2)
+      """
+      import numpy as np
+      X=CGS.h*FreqGHz*1e9/(CGS.K*self.Tcmb)
+      eX=np.exp(X)
+      return eX*(X/(eX-1.))**2
+CMB=__cmb()

src/yapsut/planck_legacy/cosmological_dipole.py

+__DESCRIPTION__="""
+Simulates a cosmological dipole
+"""
+
+from numpy import pi as PI
+from numpy import zeros, cos, sin, array, sqrt, arange
+
+from .cgs_blackbody import CMB
+from .vec3dict import vec3dict
+
+try :
+  from healpy.pixelfunc import pix2ang
+except :
+  print("Warning: no healpy support")
+
+class CosmologicalDipole(object) :
+  def __init__(self,TK,lat_deg,lon_deg,Tcmb=None) :
+    self.Tcmb=2.725 if Tcmb == None else Tcmb*1
+    try :
+       self.TK = TK[0]*1.
+       self.Err_TK = TK[1]*1.
+    except :
+       self.TK = TK*1.
+       self.Err_TK = 0.
+    try :
+       self.lat_deg = lat_deg[0]*1.
+       self.Err_lat_deg = lat_deg[1]*1.
+    except :
+       self.lat_deg = lat_deg*1.
+       self.Err_lat_deg = 0.
+    try :
+       self.long_deg = lon_deg[0]*1.
+       self.Err_long_deg = lon_deg[1]*1.
+    except :
+       self.long_deg = lon_deg*1.
+       self.Err_long_deg = 0.
+    self.long = self.long_deg/180.*PI
+    self.Err_long = self.Err_long_deg/180.*PI
+    self.lat = self.lat_deg/180.*PI
+    self.Err_lat = self.Err_lat_deg/180.*PI
+    self.vers=vec3dict()
+    self.vers['x'] = cos(self.long) * cos(self.lat)
+    self.vers['y'] = sin(self.long) * cos(self.lat)
+    self.vers['z'] = sin(self.lat)
+    vv = 0.
+    for k in ['x','y','z'] :
+      vv += self.vers[k]**2
+    vv = sqrt(vv)
+    for k in ['x','y','z'] :
+      self.vers[k] = self.vers[k]/vv
+    self.Dip = vec3dict()
+    for k in ['x','y','z'] :
+       self.Dip[k] = self.TK*self.vers[k]
+    self.Vsun = vec3dict()
+    for k in ['x','y','z'] :
+       self.Vsun[k] = self.Dip[k]/self.Tcmb
+  def __call__(self,PList) :
+    return self.toi(PList)
+  def amplitude(self,SpinAxis) :
+     import numpy as np
+     return self.TK*np.sin(np.arccos(self.vers['x']*SpinAxis['x']+self.vers['y']*SpinAxis['y']+self.vers['z']*SpinAxis['z']))
+  def TemperatureCMB_K(self) :
+     return self.Tcmb
+  def verErr(self,dx,dy,asArray=False) :
+    from vec3dict import vec3dict
+    vers=vec3dict()
+    vers['x'] = cos(self.long+self.Err_long*dx) * cos(self.lat+self.Err_lat*dy)
+    vers['y'] = sin(self.long+self.Err_long*dx) * cos(self.lat+self.Err_lat*dy)
+    vers['z'] = sin(self.lat+self.Err_lat*dy)
+    if asArray :
+       return array([vers['x'],vers['y'],vers['z']])
+    return vers
+  def cordErr(self,dx,dy,asArray=False,deg=True) :
+    if deg :
+      if asArray :
+         return array([self.long_deg+self.Err_long_deg*dx,self.lat_deg+self.Err_lat_deg*dy])
+      return {'long_deg':self.long_deg+self.Err_long_deg*dx,'lat_deg':self.lat_deg+self.Err_lat_deg*dy}
+    if asArray :
+       return array([self.long+self.Err_long*dx,self.lat+self.Err_lat*dy])
+    return {'long_rad':self.long+self.Err_long*dx,'lat_rad':self.lat+self.Err_lat*dy}
+  def toi(self,PList,Vsun=False) :
+     if Vsun :
+        return self.Vsun['x']*PList['x']+self.Vsun['y']*PList['y']+self.Vsun['z']*PList['z']
+     return self.Dip['x']*PList['x']+self.Dip['y']*PList['y']+self.Dip['z']*PList['z']
+  def quadrupole(self,PListBeamRF) :
+     """Needs pointing list in beam reference frame, with Z the beam center direction
+      This is just a bare approx, do not use for serious calculations
+     """
+     toi=self.toi(PListBeamRF,Vsun=True)
+     return self.TemperatureCMB()*0.5*toi**2
+  def convolution_effect(self,phase,dipole_spin_angle,boresight,dump,coscan_shift,croscan_shift,deg=True) :
+     """computes the effect of convolution causing a shift and a dump of the dipole """
+     import numpy as np
+     out = {'phase':phase,'boresight':boresight,'dipole_spin_angle':dipole_spin_angle,'dump':dump,'coscan_shift':coscan_shift,'croscan_shift':croscan_shift}
+     fact=np.pi/180. if deg else 1.
+     cb=np.cos(boresight*fact)
+     sb=np.sin(boresight*fact)
+     ct=np.cos(dipole_spin_angle*fact)
+     st=np.sin(dipole_spin_angle*fact)
+     cphi=np.cos(phase*fact)
+     out['Aver']=self.TK*ct*cb
+     out['Amp']=self.TK*st*sb
+     out['Var']=out['Amp']*cphi
+     out['Dipole']=out['Aver']+out['Var']
+     cbp=np.cos((boresight+croscan_shift)*fact)
+     sbp=np.sin((boresight+croscan_shift)*fact)
+     cphiP=np.cos((phase+coscan_shift)*fact)
+     out['Aver-cross-scan']=self.TK*ct*cbp
+     out['Amp-cross-scan']=self.TK*st*sbp
+     out['DeltaAver-cross-scan']=self.TK*ct*cbp-out['Aver']
+     out['DeltaAmp-cross-scan']=self.TK*st*sbp-out['Amp']
+     out['Var-coscan']=out['Amp']*cphiP
+     out['DeltaVar-coscan']=out['Var-coscan']-out['Var']
+     out['DipoleConvolved']=dump*(out['Aver-cross-scan']+out['Amp-cross-scan']*cphiP)
+     out['DeltaDipole']=out['DipoleConvolved']-out['Dipole']
+     return out
+  def Map(self,Nside) : 
+    npix = 12*Nside**2
+    ipix = arange(npix,dtype=int)
+    theta,phi = pix2ang(Nside,ipix)
+    x = cos(phi) * sin(theta)
+    y = sin(phi) * sin(theta)
+    z = cos(theta)
+    vv = sqrt(x**2+y**2+z**2)
+    x = x/vv
+    y = y/vv
+    z = z/vv
+    return self.TK*(x*self.vers['x']+y*self.vers['y']+z*self.vers['z'])
+    
+class CosmologicalDipoleGalactic(CosmologicalDipole) :
+  def __init__(self,TK = [3.355e-3,8e-6],lat_deg = [48.26,0.03],lon_deg = [263.99,0.14]) :
+    #CosmologicalDipole.__init__(self,TK = [3.355e-3,8e-6],lat_deg = [48.26,0.03],lon_deg = [263.99,0.14])
+    CosmologicalDipole.__init__(self,TK = TK,lat_deg = lat_deg,lon_deg = lon_deg)
+    self.Source =     """WMAP5 solar system dipole parameters, from: http://arxiv.org/abs/0803.0732"""
+    
+class CosmologicalDipoleEcliptic(CosmologicalDipole) :
+  def __init__(self,TK = 3.355e-3,lat_deg = -11.14128,lon_deg = 171.64904) :
+    CosmologicalDipole.__init__(self,TK = TK,lat_deg = lat_deg,lon_deg = lon_deg)
+    self.Source =     """WMAP5 solar system dipole parameters, from: http://arxiv.org/abs/0803.0732"""
+    
+def wmap5_parameters():
+    """WMAP5 solar system dipole parameters, 
+    from: http://arxiv.org/abs/0803.0732"""
+    # 369.0 +- .9 Km/s
+    SOLSYSSPEED = 369e3
+    ## direction in galactic coordinates
+    ##(d, l, b) = (3.355 +- 0.008 mK,263.99 +- 0.14,48.26deg +- 0.03)
+    SOLSYSDIR_GAL_THETA = np.deg2rad( 90 - 48.26 )
+    SOLSYSDIR_GAL_PHI = np.deg2rad( 263.99 )
+    SOLSYSSPEED_GAL_U = healpy.ang2vec(SOLSYSDIR_GAL_THETA,SOLSYSDIR_GAL_PHI)
+    SOLSYSSPEED_GAL_V = SOLSYSSPEED * SOLSYSSPEED_GAL_U
+    SOLSYSSPEED_ECL_U = gal2ecl(SOLSYSSPEED_GAL_U)
+    SOLSYSDIR_ECL_THETA, SOLSYSDIR_ECL_PHI = healpy.vec2ang(SOLSYSSPEED_ECL_U)
+    ########## /WMAP5
+    return SOLSYSSPEED, SOLSYSDIR_ECL_THETA, SOLSYSDIR_ECL_PHI
+
+
+class DynamicalDipole :
+   def __init__(self,ObsEphemeridsFile,AU_DAY=True,Tcmb=None) :
+      """ephmerids must be either in AU_DAY or KM_SEC
+         set AU_DAY = True if in AU_DAY (default) or False if in KM_SEC
+      """ 
+      from readcol import readcol 
+      import numpy as np
+      import copy
+      from scipy import interpolate
+      import pickle
+      #from matplotlib import pyplot as plt
+      self.Tcmb=2.725 if Tcmb == None else Tcmb*1
+      self.ObsEphemeridsFile=ObsEphemeridsFile
+      try :
+         f=open(ObsEphemeridsFile,'r') 
+      except :
+         print("Error: %s file not foud." % ObsEphemeridsFile)
+         return
+      try :
+         self.eph=pickle.load(f)
+         f.close()
+      except :
+         f.close()
+         eph=readcol(ObsEphemeridsFile,fsep=',',asStruct=True).__dict__
+         self.eph={}
+         for k in ['VY', 'VX', 'seconds', 'Month', 'hours', 'Y', 'Year', 'JD', 'X', 'VZ', 'Z', 'minutes', 'day'] :
+            self.eph[k]=copy.deepcopy(eph[k])
+      self.JD0=self.eph['JD'][0]
+      self.JDMax=self.eph['JD'].max()
+      self.itp={}
+      self.tscale = self.JDMax-self.JD0
+      self.AU_km=149597870.691
+      self.DAY_sec=86400.0  
+      lU=self.AU_km if AU_DAY else 1.
+      lT=self.DAY_sec if AU_DAY else 1.
+      for k in ['X','Y','Z'] :
+         self.itp[k]=interpolate.interp1d((self.eph['JD']-self.eph['JD'][0])/self.tscale,lU*self.eph[k])
+      for k in ['VX','VY','VZ'] :
+         self.itp[k]=interpolate.interp1d((self.eph['JD']-self.eph['JD'][0])/self.tscale,lU/lT*self.eph[k])
+      for k in ['X','Y','Z'] :
+         self.itp['Dip'+k]=interpolate.interp1d((self.eph['JD']-self.eph['JD'][0])/self.tscale,lU/lT*self.TemperatureCMB_K()/self.speed_of_light_kmsec()*self.eph['V'+k])
+   def test(self) :
+      "returns the value for dipole parallel to V, result must be order 30/3e5*2.75 Kcmb"
+      px=self.itp['VX'](0.)
+      py=self.itp['VY'](0.)
+      pz=self.itp['VZ'](0.)
+      v=(px**2+py**2+pz**2)**0.5
+      px*=1/v
+      py*=1/v
+      pz*=1/v
+      return self.Dipole(0.,px,py,pz)
+   def __len__(self) :
+      return len(self.eph('JD'))
+   def speed_of_light_kmsec(self) :
+      return 2.99792458e5
+   def TemperatureCMB_K(self) :
+      return self.Tcmb
+   def pickle(self,filename) :
+      import pickle
+      try :
+         f=open(filename,'w')
+      except :
+         print("Error: %s file not foud." % filename)
+         return
+      pickle.dump(self.eph,f)
+      f.close()
+   def Vel(self,t) :
+      from vec3dict import vec3dict
+      vx=self.itp['DipX'](t)
+      vy=self.itp['DipY'](t)
+      vz=self.itp['DipZ'](t)
+      return vec3dict(vx,vy,vz)
+   def DotVel(self,t,px,py,pz) :
+      vx = self.itp['DipX'](t)
+      vy=self.itp['DipX'](t)
+      vz=self.itp['DipZ'](t)
+      acc = vx*px
+      acc += vy*py
+      acc += vz*pz
+      return acc/((vx**2+vy**2+vz**2)*(px**2+py**2+pz**2))**0.5
+   def VelVersor(self,t) :
+      return self.Vel(t).versor()
+   def DotVelVersor(self,t,px,py,pz) :
+      return self.DotVel(t,px,py,pz)/self.Vel(t).norm()
+   def Dipole(self,t,px,py,pz) :
+      """Returns the DD dipole for given t and p where p = p['x'],p['y'],p['z'] at times t (t = days since JD0)"""
+      acc = self.itp['DipX'](t)*px
+      acc += self.itp['DipY'](t)*py
+      acc += self.itp['DipZ'](t)*pz
+      return acc
+   def __call__(self,JD,p,versor=False) :
+      if versor :
+         return self.Vel((JD-self.JD0)/self.tscale).versor().dot(p)
+      return self.Dipole((JD-self.JD0)/self.tscale,p['x'],p['y'],p['z'])
+
+if __name__=="__main__" :
+   from readcol import readcol
+   from sky_scan import ScanCircle
+   import numpy as np
+   import time
+   DD=DynamicalDipole('../DB/planck_2009-04-15T00:00_2012-08-04T00:00_ssb_1hour.csv')
+   
+   LFIEffectiveBoresightAngles=readcol('../DB/lfi_effective_boresight_angles.csv',asStruct=True,fsep=',').__dict__
+   PointingList=readcol('../DB/ahf_pointings.csv',asStruct=True,fsep=',').__dict__
+   
+   PLL = {}
+   idx = np.where(PointingList['pointID_PSO'] < 90000000)
+   PLL['PID'] = PointingList['pointID_unique'][idx]
+   PLL['JD'] = PointingList['midJD'][idx]
+   PLL['spin_ecl_lat_rad'] = PointingList['spin_ecl_lat'][idx]/180.*np.pi
+   PLL['spin_ecl_lon_rad'] = PointingList['spin_ecl_lon'][idx]/180.*np.pi
+   PLL['spin_x'] = np.cos(PLL['spin_ecl_lon_rad'])*np.cos(PLL['spin_ecl_lat_rad'])
+   PLL['spin_y'] = np.sin(PLL['spin_ecl_lon_rad'])*np.cos(PLL['spin_ecl_lat_rad'])
+   PLL['spin_z'] = np.sin(PLL['spin_ecl_lat_rad'])
+   
+   fh=27
+   beta_rad = LFIEffectiveBoresightAngles['beta_fh'][fh-18]*np.pi/180.
+
+
+   #for iipid in range(10) :
+      #ipid = 1000*iipid
+      
+
+
+   #stats = {}
+   #for k in ['t','d','c'] :
+      #for arg in ['min','max','amp'] :
+         
+   
+   ddmax=np.zeros(len(PLL['JD']))
+   ddmin=np.zeros(len(PLL['JD']))
+   cdmax=np.zeros(len(PLL['JD']))
+   cdmin=np.zeros(len(PLL['JD']))
+   tdmax=np.zeros(len(PLL['JD']))
+   tdmin=np.zeros(len(PLL['JD']))
+   tic = time.time()    
+   SC = ScanCircle(PLL['spin_ecl_lon_rad'][0],PLL['spin_ecl_lat_rad'][0],beta_rad,NSMP=360)
+   SC.generate_circle()
+   Cphi = np.cos(SC.phi)
+   Sphi = np.sin(SC.phi)
+   SumCC=(Cphi*Cphi).sum()
+   SumCS=(Cphi*Sphi).sum()
+   SumSS=(Sphi*Sphi).sum()
+   ArgVSpin = np.zeros(len(PLL['JD']))
+   ArgRSpin = np.zeros(len(PLL['JD']))
+   ArgVR = np.zeros(len(PLL['JD']))
+   ArgVPmax = np.zeros(len(PLL['JD']))
+   ArgVPmin = np.zeros(len(PLL['JD']))
+   for ipid in range(len(PLL['JD'])) :
+      SC = ScanCircle(PLL['spin_ecl_lon_rad'][ipid],PLL['spin_ecl_lat_rad'][ipid],beta_rad,NSMP=360)
+      SC.generate_circle()
+      dd=DD(PLL['JD'][ipid],SC.r)
+
+      t=(PLL['JD'][ipid]-DD.JD0)/DD.tscale
+      
+      #ArgVSpin[ipid]=DD.DotVel((PLL['JD'][ipid]-DD.JD0)/DD.tscale,PLL['spin_x'][ipid],PLL['spin_y'][ipid],PLL['spin_z'][ipid])
+      
+      sx=PLL['spin_x'][ipid]
+      sy=PLL['spin_y'][ipid]
+      sz=PLL['spin_z'][ipid]
+      vx=DD.itp['VX'](t)
+      vy=DD.itp['VY'](t)
+      vz=DD.itp['VZ'](t)
+      rx=DD.itp['X'](t)
+      ry=DD.itp['Y'](t)
+      rz=DD.itp['Z'](t)
+      
+      ArgVSpin[ipid] = 180./np.pi*np.arccos((sx*vx+sy*vy+sz*vz)/(((sx*sx+sy*sy+sz*sz)*(vx*vx+vy*vy+vz*vz))**0.5))
+      ArgRSpin[ipid] = 180./np.pi*np.arccos((sx*rx+sy*ry+sz*rz)/(((sx*sx+sy*sy+sz*sz)*(rx*rx+ry*ry+rz*rz))**0.5))
+      ArgVR[ipid] = 180./np.pi*np.arccos((vx*rx+vy*ry+vz*rz)/(((vx*vx+vy*vy+vz*vz)*(rx*rx+ry*ry+rz*rz))**0.5))
+      
+      dp=180./np.pi*np.arccos((SC.r['x']*vx+SC.r['y']*vy+SC.r['z']*vz)/(((SC.r['x']**2+SC.r['y']**2+SC.r['z']**2)*(vx**2+vy**2+vz**2))**0.5))
+      ArgVPmax[ipid]=dp.max()
+      ArgVPmin[ipid]=dp.min()
+
+
+      #ArgRSpin[ipid]=DD.DotVel((PLL['JD'][ipid]-DD.JD0)/DD.tscale,PLL['spin_x'][ipid],PLL['spin_y'][ipid],PLL['spin_z'][ipid])
+      cd = CosmologicalDipoleEcliptic().toi(SC.r)
+      td=cd+dd
+      ddmax[ipid]=dd.max()
+      ddmin[ipid]=dd.min()
+      cdmax[ipid]=cd.max()
+      cdmin[ipid]=cd.min()
+      tdmax[ipid]=td.max()
+      tdmin[ipid]=td.min()
+   print("Elapsed %e sec",time.time()-tic)
+   
+   import matplotlib.pyplot as plt
+   
+   plt.close('all')
+   plt.figure(1)
+   plt.plot(PLL['PID'],ArgVSpin,'.',label='Spin,Velocity')
+   plt.plot(PLL['PID'],ArgRSpin,'.',label='Spin,Radius')
+   plt.plot(PLL['PID'],ArgVR,'.',label='Radius,Velocity')
+   plt.legend(loc=5)
+   plt.title('Sol.Syst.Bar. Angles')
+   plt.xlabel('PID')
+   plt.ylabel('Angle [deg]')
+   plt.savefig('geometry_angles.png',dpi=300)
+   plt.show()
+      
+   plt.figure(2)
+   plt.plot(PLL['PID'],ArgVPmin,'.',label='min Pointing,Velocity angle')
+   plt.plot(PLL['PID'],ArgVPmax,'.',label='max Pointing,Velocity angle')
+   plt.legend(loc=5)
+   plt.title('Angles between Pointing direction and Velocity')
+   plt.xlabel('PID')
+   plt.ylabel('Angle [deg]')
+   plt.savefig('pointing_velocity_angles_minimax.png',dpi=300)
+   plt.show()
+   
+   
+   

src/yapsut/planck_legacy/vec3dict.py

+import numpy as _np
+class vec3dict :
+   """ handles a 3d vector both as a dictionary and an indexed array """
+   def __init__(self,*karg) :
+      """ vec3dict() : empty vector 3d (Values are None)
+          vec3dict(0) : vector 3d with all the elements [int(0), int(0), int(0)]
+          vec3dict(5) : vector 3d with all the elements int(5)
+          vec3dict(1,3,5) : vector 3d with all elements [int(1), int(3), int(5)]
+          vec3dict('x') : vector 3d [1.,0.,0.]
+      """
+      import numpy as np
+      self.v=[None,None,None]
+      if len(karg) == 0 :
+         return 
+      if type(karg[0]) == type('') :
+         self.v=[0.,0.,0.]
+         if karg[0] == 'x' or karg[0] == 'X' :
+            self.v[0]=1.
+         elif karg[0] == 'y' or karg[0] == 'Y' :
+            self.v[1]=1.
+         elif karg[0] == 'z' or karg[0] == 'Z' :
+            self.v[2]=1.
+         else :
+            raise Exception('invalid versor name')
+         self.v=np.array(self.v)
+         return
+      if type(karg[0]) == type(vec3dict()) :
+         self.v[0] = karg[0][0]
+         self.v[1] = karg[0][1]
+         self.v[2] = karg[0][2]
+         self.v=np.array(self.v)
+         return
+      if len(karg) < 3 :
+         for k in range(3) :
+            self.v[k]=karg[0]
+         self.v=np.array(self.v)
+      else :
+         self.v[0]=karg[0]
+         self.v[1]=karg[1]
+         self.v[2]=karg[2]
+         self.v=np.array(self.v)
+   def name2index(self,key):
+      "return the index for a given component name"
+      if key=='x' :
+         return 0
+      elif key=='y' :
+         return 1
+      elif key=='z' :
+         return 2
+      else :
+         raise Exception('invalid keyword')
+   def index2name(self,idx):
+      "return the name for a given component index"
+      try :
+         return ['x','y','z'][idx]
+      except :
+         raise Exception('invalid index')
+   def ismybrother(self,that) :
+      return that.__class__.__name__ == self.__class__.__name__ 
+   def __len__(self) :
+      try :
+         return self.v.shape[1]
+      except :
+         return 1
+   def linterp(self,t,left=-_np.inf,right=_np.inf,noTest=False) :
+      import numpy as np
+      from scipy import interpolate
+      if noTest :
+         it=np.array(np.floor(t),dtype='int')
+         dt=t-it
+         x=(self['x'][it+1]-self['x'][it])*dt+self['x'][it]
+         y=(self['y'][it+1]-self['y'][it])*dt+self['y'][it]
+         z=(self['z'][it+1]-self['z'][it])*dt+self['z'][it]
+         return vec3dict(x,y,z)
+      x=np.zeros(t.shape)
+      y=np.zeros(t.shape)
+      z=np.zeros(t.shape)
+      it=np.array(np.floor(t),dtype='int')
+      idx=np.where(it < 0)[0]
+      if len(idx) > 0 :
+         x[idx]=left
+         y[idx]=left
+         z[idx]=left
+      idx=np.where(it > len(self)-1)[0]
+      if len(idx) > 0 :
+         x[idx]=right
+         y[idx]=right
+         z[idx]=right
+      idx=np.where(it == len(self)-1)[0]
+      if len(idx) > 0 :
+         x[idx]=self['x'][-1]
+         y[idx]=self['x'][-1]
+         z[idx]=self['x'][-1]
+      idx=np.where((0<=it)*(it<len(self)-1))[0]
+      if len(idx) > 0 :
+         dt=t[idx]-it[idx]
+         it=it[idx]
+         x[idx]=(self['x'][it+1]-self['x'][it])*dt+self['x'][it]
+         y[idx]=(self['y'][it+1]-self['y'][it])*dt+self['y'][it]
+         z[idx]=(self['z'][it+1]-self['z'][it])*dt+self['z'][it]
+      return vec3dict(x,y,z)
+   def argslice(self,idx) :
+      out=vec3dict()
+      out['x']=self['x'][idx]
+      out['y']=self['y'][idx]
+      out['z']=self['z'][idx]
+      return out
+   def __getitem__(self,idx) :
+      import numpy as np
+      if (idx.__class__ == np.zeros(2).__class__) :
+         return self.argslice(idx)
+      try :
+         idx = int(idx)
+         try :
+            return self.v[idx]
+         except :
+            raise Exception('out of bounds')
+      except :
+         if idx=='x' :
+            i=0
+         elif idx=='y' :
+            i=1
+         elif idx=='z' :
+            i=2
+         else :
+            raise Exception('invalid keyword')
+         return self.v[i]
+   def __setitem__(self,idx,that) :
+      try :
+         idx = int(idx)
+         try :
+            self.v[idx] = that
+         except :
+            raise Exception('out of bounds')
+      except :
+         if idx=='x' :
+            i=0
+         elif idx=='y' :
+            i=1
+         elif idx=='z' :
+            i=2
+         else :
+            raise Exception('invalid keyword')
+         self.v[i] = that
+   def __add__(self,that) :
+      new = vec3dict(self)
+      if self.ismybrother(that):
+         for k in range(3) :
+            new.v[k]+=that[k]
+      else :
+         for k in range(3) :
+            new.v[k]+=that
+      return new
+   def __sub__(self,that) :
+      new = vec3dict(self)
+      if self.ismybrother(that) :
+         for k in range(3) :
+            new.v[k]-=that[k]
+      else :
+         for k in range(3) :
+            new.v[k]-=that
+      return new
+   def __div__(self,that) :
+      new = vec3dict(self)
+      if not self.ismybrother(that) :
+         for k in range(3) :
+            new.v[k]=new.v[k]/that
+      else :
+         raise Exception('right side not a scalar')
+      return new
+   def __mul__(self,that) :
+      if not self.ismybrother(that) :
+         new = vec3dict(self)
+         for k in range(3) :
+            new.v[k]*=that
+         return new
+      else :
+         raise Exception('right side not a scalar')
+   def __rmul__(self,that) :
+      if not self.ismybrother(that) :
+         new = vec3dict(self)
+         for k in range(3) :
+            new.v[k]*=that
+         return new
+      else :
+         raise Exception('left side not a scalar')
+   def __neg__(self) :
+      return vec3dict(-self.v[0],-self.v[1],-self.v[2])
+   def __str__(self) :
+      return 'x : '+str(self.v[0])+', y : '+str(self.v[1])+', z : '+str(self.v[2])
+   def mul_by_array(self,that) :
+      """
+         multiplies by an array, component by component
+         it is used as an example to generate a multicomponent vec3dict
+         Example:
+            arange(10)*vec3dict(1,0,0) 
+         will result in a list of vec3dict objects with one component, while
+            vec3dict(1,0,0).by_array(arange(10)) 
+         will result in a single vec3dict objects with 10 elements in each column
+         BeWare: the resulting Vec3Dict is not granted to work with all the methods
+      """
+      return vec3dict(self.v[0]*that,self.v[1]*that,self.v[2]*that)
+   def add_array(self,that) :
+      """
+         add an array, component by component
+         it is used as an example to generate a multicomponent vec3dict
+         Example:
+            arange(10)+vec3dict(1,0,0) 
+         will result in an error
+            vec3dict(1,0,0).add_array(arange(10)) 
+         will result in a single vec3dict objects with 10 elements in each column
+         which are the sum of 
+            vec3dict(1,0,0)['x'], vec3dict(1,0,0)['y'], vec3dict(1,0,0)['z']
+         with the array
+      """
+      return vec3dict(self.v[0]+that,self.v[1]+that,self.v[2]+that)
+   def dot(self,that) :
+      """ dot product """
+      if not self.ismybrother(that) :
+         raise Exception('right side not a vector')
+      else :
+         return self.v[0]*that.v[0]+self.v[1]*that.v[1]+self.v[2]*that.v[2]
+   def norm(self) :
+      """ returns the norm """
+      return (self.v[0]*self.v[0]+self.v[1]*self.v[1]+self.v[2]*self.v[2])**0.5
+   def keys(self) :
+      """ return names of elements """
+      return ['x','y','z']
+   def array(self,dtype=None) :
+      """ returns an array """
+      from numpy import array
+      if dtype == None :
+         return array(self.v)
+      return array(self.v,dtype=dtype)
+   def dict(self) :
+      """ returns a dictionary """
+      return {'x':self.v[0],'y':self.v[1],'z':self.v[2]}
+   def ext(self,that) :
+      if not self.ismybrother(that) :
+         raise Exception('right side not a vector')
+      else :
+         new = vec3dict(0)
+         new[0]=self.v[1]*that.v[2]-self.v[2]*that.v[1]
+         new[1]=-self.v[0]*that.v[2]+self.v[2]*that.v[0]
+         new[2]=self.v[0]*that.v[1]-self.v[1]*that.v[0]
+         return new
+   def norm_ext(self,that) :
+      "norm of an external product"
+      if not self.ismybrother(that) :
+         raise Exception('right side not a vector')
+      else :
+         new=(self.v[1]*that.v[2]-self.v[2]*that.v[1])**2
+         new+=(self.v[0]*that.v[2]-self.v[2]*that.v[0])**2
+         new+=(self.v[1]*that.v[0]-self.v[0]*that.v[1])**2
+         return new**0.5
+   def angle(self,that) :
+      "angle between two vectors, radiants"
+      import numpy as np
+      if not self.ismybrother(that) :
+         raise Exception('right side not a vector')
+      else :
+         Dot=self.dot(that)
+         Ext=self.norm_ext(that)
+         return np.arctan2(Ext,Dot)
+   def copy(self) :
+      import copy
+      return copy.deepcopy(self)
+   def versor(self) :
+      "returns the versor"
+      new=self.copy()
+      nn=new.norm()
+      for k in new.keys() : new[k]=new[k]/nn
+      return new
+   def mean(self) :
+      "returns the averaged vector"
+      return vec3dict(self.v[0].mean(),self.v[1].mean(),self.v[2].mean())
+         
+#class matr3dict :
+   #""" handles a 3x3 matrix vector both as a dictionary and an indexed array , note rows are 3d matrices"""
+
+if __name__=='__main__' :
+   import numpy as np
+   v1=vec3dict()
+   print v1
+   
+   v2=vec3dict(1)
+   print v2
+
+   v3=vec3dict(1,2,3)
+   print v3
+   
+   v4=vec3dict(v3)
+   print v4
+   
+   print -v4
+   
+   print 2.*v4
+   print v4+3.
+   print v4.norm()
+   print v4 - v4
+   print v4.dot(v4)
+   print v4.ext(v4+vec3dict(0.,0.,1.))
+   print v4+vec3dict(0.,0.,1.)
+   
+   vx=vec3dict(1.,0.,0.)
+   vy=vec3dict(0.,1.,0.)
+   vz=vec3dict(0.,0.,1.)
+   
+   print
+   print "||VX x VY|| ",vx.norm_ext(vy)
+   print "||VX x V45degXY|| ",vx.norm_ext(vec3dict(np.cos(np.pi/4),np.cos(np.pi/4),0.))
+   print "||VX x V45degXZ|| ",vx.norm_ext(vec3dict(np.cos(np.pi/4),0.,np.cos(np.pi/4)))
+   print
+   print "angle(vx,vy) = ",vx.angle(vy)*180./np.pi
+   print "angle(vx,vz) = ",vx.angle(vz)*180./np.pi
+   print "angle(vy,vz) = ",vy.angle(vz)*180./np.pi
+   print
+   print "angle(vx,V45degXY) = ",vx.angle(vec3dict(np.cos(np.pi/4),np.cos(np.pi/4),0.))*180./np.pi
+   print "angle(vx,V45degXZ) = ",vx.angle(vec3dict(np.cos(np.pi/4),0.,np.cos(np.pi/4)))*180./np.pi
+