diff --git a/src/yapsut/__init__.py b/src/yapsut/__init__.py
index b8965bab1d40f772a123b013696b959f7ce734e3..708e8b6a827b56f94a494c8d5a877c2bcbc82ac2 100644
--- a/src/yapsut/__init__.py
+++ b/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
diff --git a/src/yapsut/planck_legacy/cgs_blackbody.py b/src/yapsut/planck_legacy/cgs_blackbody.py
new file mode 100644
index 0000000000000000000000000000000000000000..9ed7b3402d8b937c026648bf5c96f68d04abb468
--- /dev/null
+++ b/src/yapsut/planck_legacy/cgs_blackbody.py
@@ -0,0 +1,301 @@
+__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()
diff --git a/src/yapsut/planck_legacy/cosmological_dipole.py b/src/yapsut/planck_legacy/cosmological_dipole.py
new file mode 100644
index 0000000000000000000000000000000000000000..6258a9783cafa963b5de979fdbf16505e0ccd258
--- /dev/null
+++ b/src/yapsut/planck_legacy/cosmological_dipole.py
@@ -0,0 +1,370 @@
+__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()
+
+
+
diff --git a/src/yapsut/planck_legacy/time_support.py b/src/yapsut/planck_legacy/time_support.py
new file mode 100644
index 0000000000000000000000000000000000000000..88b62af6f20257ed72ef0b4c9108d47018342454
--- /dev/null
+++ b/src/yapsut/planck_legacy/time_support.py
@@ -0,0 +1,918 @@
+#!/usr/bin//python
+import time
+import datetime
+import struct
+import array
+import os.path
+import sys
+
+import numpy as np
+
+from juldate import juldate as _juldate
+from juldate import caldate as _caldate
+
+"""
+ Module to handle different time scales used by Planck, in particular:
+
+ UTC in ISO and ZULU format by using classes ISO, ZULU
+
+ UTC in seconds by using class UTCsec
+
+ TAI, PsudoOBT, SCET by using the classes TAI, OBT, SCET
+
+ it uses Python modules: time and datetime
+
+ Python ASSUMES UTCsec starts from 1970 Jan 01, 01:00:00.
+
+ Conversions are described by the following conversion ring:
+
+ ZULU --- ISO
+ \ /
+ JulianDay ------- PseudoOD
+ / \
+ TAI PseudoOBT
+ \ /
+ SCET - UTCsec
+
+ where /,\,- denotes interconversion possibility.
+
+ PseudoOBT is OBT calculate ASSUMING
+ -- onboard clock properly set at zero time
+ -- onboard clock advancing at a constant rate of 1sec/sec
+ -- onboard clodk without drifts.
+
+ Differences between PseudoOBT and true OBT are expected to be of the order of some second.
+
+ Note that since Python Datetime class assumes use of UTCsec, the only method
+ returnin a Datetime object is UTCsec.datetime() use to convert an UTCsec into
+ a Datetime.
+
+ PseudoOD is calculate ASSUMING to count the number of days after the launch epoch.
+ So PseudoOD is interconnected just to JulianDay, since PseudoOD is nothing else than
+ Julian Day with a different starting time.
+
+ WARNING OD DEPRECATED:
+ To avoid confusion with the true OD (Level1 OD) which can not be computed but which can be
+ tabulated from AHF the class OD should be avoided. The class OD and JuliandDay-OD conversion
+ methods are DEPRECATED.
+
+ USE PseudoOD instead.
+
+ Author:
+ Michele Maris, Issue 1.0, 30 Mar 2010
+ Issue 1.1, 14 May 2010 - added management of ODs
+ Issue 1.2, 25 Jun 2010 - corrected OD zero point
+ OD is now deprecated USE PseudoOD instead
+ added the possibility to use it as a command lines tool
+"""
+
+class _TIME_SCALE_BASE() :
+ """ base class to handle times: private class
+ for each base: TAI, SCET, PseudoOBT, UTC defines:
+ zero point (example ZERO_TAI) and if needed TICKS number of tickes per second
+
+ It provides methods to recover such quantities,
+ a __str__ method and a set method.
+
+ Basic representation of times is Julian Dates:
+ """
+ ZERO_MJD = _juldate((1858,11,17,12,0,0))
+ ZERO_TAI = _juldate((1958,1,1,0,0,0,0))
+ ZERO_SCET = _juldate((1958,1,1,0,0,0,0))
+ ZERO_UTC = _juldate((1970,1,1,0,0,0,0))
+ ZERO_PseudoOBT = _juldate((1958,1,1,0,0,0,0))
+
+ #ZERO_OD = _juldate((2009,5,14,10,13,00))
+
+ ZERO_OD = _juldate((2009,5,13,13,12,00))
+
+ TICKS_SCET_SEC = 1e6
+ TICKS_PseudoOBT_SEC = float(2**16)
+
+ TICKS_PseudoOBT_DAY = float(86400)*TICKS_PseudoOBT_SEC
+ TICKS_SCET_DAY = float(86400)*TICKS_SCET_SEC
+ TICKS_UTC_DAY = float(86400)
+ TICKS_TAI_DAY = float(86400)
+
+ def __init__(self,v=0.) :
+ """ creator accepts scalars, strings, and lists
+
+ scalar (float, int, long): saves the value as double
+
+ list : assumes a tuple of type: (YYYY, MM, DD, HH, MM, SS_DEC) and converts in Julian Day
+
+ string : assumes a string of type: "YYYY-MM-DDTHH:MM:SS.DEC" and converts in Julian Day
+ if empty is equal to put v = 0
+ use method set(v) to change value to a scalar, time() to return value
+ """
+ if self._isNumericScalar(v) :
+ self.v = np.float64(v)
+ return
+
+ if type(v) == type(()) or type(v) == type([]) :
+ if len(v) == 3 :
+ self.v = _juldate((v[0],v[1],v[2],0.,0.,0.))
+ return
+ if len(v) == 5 :
+ self.v = _juldate((v[0],v[1],v[2],v[3],v[4],0.))
+ return
+ else :
+ self.v = _juldate(v)
+ return
+
+ if type(v) == type('') :
+ if v == '' :
+ self.v = 0
+ return
+ else :
+ self.v = _juldate((float(v[0:4]), float(v[5:7]), float(v[8:10]), float(v[11:13]), float(v[14:16]), float(v[17:len(v)])))
+ return
+
+ print "unrecognized input type"
+ return
+
+ def __str__(self) :
+ return str(self.v)
+
+ def _isNumericScalar(self,v) :
+ tv = type(v)
+ return tv == type(0.) or tv == type(0) or tv == type(long(0)) or tv == type(np.float64(0))
+
+ def _tuple2isostring(self,tp) :
+ """ converts a tuple to an ISO time string """
+ if len(tp) == 6 :
+ return "%04d-%02d-%02dT%02d:%02d:%010.7f" % tp
+ return ""
+
+ def set(self,v) :
+ self.v = float(v)
+
+ def time(self) :
+ return self.v
+
+ def kind(self,test='') :
+ return False
+
+ def zero_tai(self) :
+ return self.__class__.ZERO_TAI
+
+ def zero_utc(self) :
+ return self.__class__.ZERO_UTC
+
+ def zero_pseudoObt(self) :
+ return self.__class__.ZERO_PseudoOBT
+
+ def zero_od(self) :
+ return self.__class__.ZERO_OD
+
+ def step_scet_sec(self) :
+ return 1./self.__class__.TICKS_SCET_SEC
+
+ def step_pseudoObt_sec(self) :
+ return 1./self.__class__.TICKS_PseudoOBT_SEC
+
+class _TIME_SCALE_PART(_TIME_SCALE_BASE) :
+ """ base class to handle times divided into two parts:
+ - a coarse part (long)
+ - a fine part (double)
+ """
+ def __init__(self,coarse=0,fine=-1.) :
+ if fine < 0 :
+ self.coarse = floor(coarse)
+ self.fine = coarse-self.coarse
+ else :
+ self.coarse = coarse
+ self.fine = fine
+
+ def tuple(self) :
+ """ returns a tuple """
+ return (self.coarse, self.fine)
+
+ def __str__(self,sept=':') :
+ return "%s%s%s" % (self.coarse,sept,self.fine)
+
+ def getTime(self,fine_ticks_sec=1.) :
+ """ getTime([ticks_fine=1.]) combines coarse and fine part into a float
+ getTime = float(coarse)+self_fine/float(fine_ticks_sec)
+ """
+ return self.coarse+self.fine/float(fine_ticks_sec)
+
+class _TIME_SCALE_STRING(_TIME_SCALE_BASE) :
+ """ base class to handle times in form of strings """
+ def __init__(self,v='') :
+ self.v=str(v)
+
+ def __str__(self) :
+ return str(self.v)
+
+ def set(self,v) :
+ self.v=str(v)
+
+
+class CalendarDate(_TIME_SCALE_BASE) :
+ def __init__(self,v) :
+ if type(v) == type(()) or type(v) == type([]) :
+ if len(v) == 3 :
+ self.v = (long(v[0]),long(v[1]),long(v[2]),long(0),long(0),float(0.))
+ return
+ if len(v) == 5 :
+ self.v = (long(v[0]),long(v[1]),long(v[2]),long(v[3]),long(v[4]),float(0.))
+ return
+ else :
+ self.v = v
+ return
+ print "unrecongnized ntuple format ",v
+ return
+
+ if type(v) == type('') :
+ if v == '' :
+ self.v = (long(0), long(0), long(0), long(0), long(0), float(0.))
+ return
+ else :
+ self.v = ((long(v[0:4]), long(v[5:7]), long(v[8:10]), long(v[11:13]), long(v[14:16]), float(v[17:len(v)])))
+ return
+ print "unrecognized string format ",v
+ return
+
+ if type(v) == type(JulianDay(0.)) :
+ self.v=_caldate(v.v)
+
+ print "unrecognized format ",v
+ return
+
+ def kind(self,test='') :
+ """ returns kind of object """
+ if test == '' :
+ return 'CalendarDate'
+ return test == 'CalendarDate'
+
+ def julianday(self) :
+ """ conversion to julian day """
+ return JulianDay(self.v)
+
+ def iso(self) :
+ """ Conversion in ISO format """
+ return ISO(self._tuple2isostring(self.v))
+
+ def zulu(self) :
+ """ Conversion in ISO format """
+ return ISO(self._tuple2isostring(self.v)).zulu()
+
+ def tuple(self) :
+ """ returns the tuple """
+ return self.v
+
+class ISO(_TIME_SCALE_STRING) :
+ """ class to handle times in form of ISO strings YYYY-MM-DDTHH:MM:SS """
+ def kind(self,test='') :
+ """ returns kind of object """
+ if test == '' :
+ return 'ISO'
+ return test == 'ISO'
+
+ def julianday(self) :
+ """ returns Julian Day """
+ return JulianDay(self.v)
+
+ def calendardate(self) :
+ """ returns the calendar date """
+ return CalendarDate(self.v)
+
+ def tuple(self) :
+ """ returns the tuple for the date"""
+ return (long(self.v[0:4]), long(self.v[5:7]), long(self.v[8:10]), long(self.v[11:13]), long(self.v[14:16]), float(self.v[17:len(self.v)]))
+
+ def zulu(self) :
+ """ returns Zulu """
+ return ZULU(self.v+'Z')
+
+class ZULU(_TIME_SCALE_STRING) :
+ """ class to handle times in form of ISO ZULU strings YYYY-MM-DDTHH:MM:SSZ """
+ def __init__(self,v='Z') :
+ _TIME_SCALE_STRING.__init__(self,v)
+ if not (self.v[len(self.v)-1] == 'Z') :
+ self.v=self.v+'Z'
+
+ def kind(self,test='') :
+ """ returns kind of object """
+ if test == '' :
+ return 'ZULU'
+ return test == 'ZULU'
+
+ def julianday(self) :
+ """ returns Julian Day """
+ return JulianDay(self.v[0:(len(self.v)-1)])
+
+ def calendardate(self) :
+ """ returns the calendar date """
+ return self.iso().calendardate()
+
+ def tuple(self) :
+ """ returns the tuple for the date"""
+ return (long(self.v[0:4]), long(self.v[5:7]), long(self.v[8:10]), long(self.v[11:13]), long(self.v[14:16]), long(self.v[17:(len(self.v)-1)]))
+
+ def iso(self) :
+ """ returns Zulu """
+ return ISO(self.v[0:(len(self.v)-1)])
+
+class JulianDay(_TIME_SCALE_BASE) :
+ """ class to handle JulianDay """
+ def kind(self,test='') :
+ if test == '' :
+ return 'JulianDay'
+ return test == 'JulianDay'
+
+ def pseudoObt(self) :
+ return PseudoOBT((self.v-self.__class__.ZERO_PseudoOBT)*self.__class__.TICKS_PseudoOBT_DAY)
+
+ def utcsec(self) :
+ return UTCsec((self.v-self.__class__.ZERO_UTC)* self.__class__.TICKS_UTC_DAY)
+
+ def scet(self) :
+ return SCET((self.v-self.__class__.ZERO_SCET)* self.__class__.TICKS_SCET_DAY)
+
+ def tai(self) :
+ return TAI((self.v-self.__class__.ZERO_TAI)* self.__class__.TICKS_TAI_DAY)
+
+ def calendardate(self) :
+ return CalendarDate(_caldate(self.v))
+
+ def tuple(self) :
+ return CalendarDate(_caldate(self.v)).tuple()
+
+ def iso(self) :
+ return ISO("%04d-%02d-%02dT%02d:%02d:%010.7f" % _caldate(self.v))
+
+ def zulu(self) :
+ return ZULU("%04d-%02d-%02dT%02d:%02d:%010.7f" % _caldate(self.v))
+
+ def od(self) :
+ return OD(self.v - self.__class__.ZERO_OD)
+
+ def pseudoOd(self) :
+ return PseudoOD(self.v - self.__class__.ZERO_OD)
+
+class OD(_TIME_SCALE_BASE) :
+ """ class to handle OD
+ (a different way to handle a Julian Day)
+ Note, the only operations allowed are: JulianDay <-> OD conversions
+
+ DEPRECATED
+ This class is deprecated, use PseudoOD instead
+ """
+ def __init__(self,v=0.) :
+ _TIME_SCALE_BASE.__init__(self,v)
+ if self.v < 0. :
+ print "Waning: negative od"
+
+ def kind(self,test='') :
+ if test == '' :
+ return 'OD'
+ return test == 'OD'
+
+ def julianday(self) :
+ return JulianDay( self.v+self.__class__.ZERO_OD)
+
+ def pseudoObt(self) :
+ return
+
+ def utcsec(self) :
+ return
+
+ def scet(self) :
+ return
+
+ def tai(self) :
+ return
+
+ def calendardate(self) :
+ return
+
+ def tuple(self) :
+ return
+
+ def iso(self) :
+ return
+
+ def zulu(self) :
+ return
+
+class PseudoOD(_TIME_SCALE_BASE) :
+ """ class to handle PseudoOD
+ (a different way to handle a Julian Day)
+ Note, the only operations allowed are: JulianDay <-> PseudoOD conversions
+ """
+ def __init__(self,v=0.) :
+ _TIME_SCALE_BASE.__init__(self,v)
+ if self.v < 0. :
+ print "Waning: negative od"
+
+ def kind(self,test='') :
+ if test == '' :
+ return 'PseudoOD'
+ return test == 'PseudoOD'
+
+ def julianday(self) :
+ return JulianDay( self.v+self.__class__.ZERO_OD)
+
+ def pseudoObt(self) :
+ return
+
+ def utcsec(self) :
+ return
+
+ def scet(self) :
+ return
+
+ def tai(self) :
+ return
+
+ def calendardate(self) :
+ return
+
+ def tuple(self) :
+ return
+
+ def iso(self) :
+ return
+
+ def zulu(self) :
+ return
+
+class PseudoOBT(_TIME_SCALE_BASE) :
+ """ class to handle PseudoOBT : in 2^-16 seconds since 1958 Jan 01, 00:00:00
+
+ PseudoOBT is OBT calculate ASSUMING
+ -- onboard clock properly set at zero time
+ -- onboard clock advancing at a constant rate of 1sec/sec
+ -- onboard clodk without drifts.
+
+ Differences between PseudoOBT and true OBT are expected to be of the order of some second.
+ """
+ def __init__(self,v=0.) :
+ _TIME_SCALE_BASE.__init__(self,v)
+ # if v is not a numeric scalar then v has been coded as a JD and it has to be converted
+ if not self._isNumericScalar(v) :
+ self.v = (self.v - self.ZERO_PseudoOBT)*self.TICKS_PseudoOBT_DAY
+
+ def kind(self,test='') :
+ if test == '' :
+ return 'PseudoOBT'
+ return test == 'PseudoOBT'
+
+ def utcsec(self) :
+ return UTCsec((self.v/self.__class__.TICKS_PseudoOBT_DAY+self.__class__.ZERO_PseudoOBT-self.__class__.ZERO_UTC)* self.__class__.TICKS_UTC_DAY)
+
+ def scet(self) :
+ return SCET((self.v/self.__class__.TICKS_PseudoOBT_DAY+self.__class__.ZERO_PseudoOBT-self.__class__.ZERO_SCET)*self.__class__.TICKS_SCET_DAY)
+
+ def tai(self) :
+ return TAI( (self.v/self.__class__.TICKS_PseudoOBT_DAY+self.__class__.ZERO_PseudoOBT-self.__class__.ZERO_TAI)*self.__class__.TICKS_TAI_DAY)
+
+ def julianday(self) :
+ return JulianDay( (self.v/self.__class__.TICKS_PseudoOBT_DAY+self.__class__.ZERO_PseudoOBT))
+
+ def iso(self) :
+ return self.julianday().iso()
+
+ def zulu(self) :
+ return self.julianday().zulu()
+
+
+class UTCsec(_TIME_SCALE_BASE) :
+ """ class to handle UTC times in seconds since 1970 Jan 01, 00:00:00 """
+ def __init__(self,v=0.) :
+ _TIME_SCALE_BASE.__init__(self,v)
+ # if v is not a numeric scalar then v has been coded as a JD and it has to be converted
+ if not self._isNumericScalar(v) :
+ self.v = (self.v - self.ZERO_UTC)*86400
+
+ def kind(self,test='') :
+ if test == '' :
+ return 'UTCsec'
+ return test == 'UTCsec'
+
+ def datetime(self) :
+ """ Returns a Datetime object from UTCsec """
+ return datetime.datetime.fromtimestamp(self.v,None)
+
+ def scet(self) :
+ return SCET( (self.v/self.__class__.TICKS_UTC_DAY+self.__class__.ZERO_UTC-self.__class__.ZERO_SCET) * self.__class__.TICKS_SCET_DAY )
+
+ def tai(self) :
+ return TAI( (self.v/self.__class__.TICKS_UTC_DAY+self.__class__.ZERO_UTC-self.__class__.ZERO_TAI)*self.__class__.TICKS_TAI_DAY )
+
+ def julianday(self) :
+ return JulianDay( (self.v/self.__class__.TICKS_UTC_DAY+self.__class__.ZERO_UTC))
+
+ def pseudoObt(self) :
+ return PseudoOBT((self.v/self.__class__.TICKS_UTC_DAY+self.__class__.ZERO_UTC-self.__class__.ZERO_PseudoOBT)*self.__class__.TICKS_PseudoOBT_DAY)
+
+ def iso(self) :
+ return self.julianday().iso()
+
+ def zulu(self) :
+ return self.julianday().zulu()
+
+
+class SCET(_TIME_SCALE_BASE) :
+ """ class to handle SCET : in microseconds since 1958 Jan 01, 00:00:00 """
+ def __init__(self,v=0.) :
+ _TIME_SCALE_BASE.__init__(self,v)
+ # if v is not a numeric scalar then v has been coded as a JD and it has to be converted
+ if not self._isNumericScalar(v) :
+ self.v = (self.v - self.__class__.ZERO_SCET)*self.__class__.TICKS_SCET_DAY
+
+ def kind(self,test='') :
+ if test == '' :
+ return 'SCET'
+ return test == 'SCET'
+
+ def tai(self) :
+ return TAI( (self.v/self.__class__.TICKS_SCET_DAY+self.__class__.ZERO_SCET-self.__class__.ZERO_TAI) * self.__class__.TICKS_TAI_DAY )
+
+ def julianday(self) :
+ return JulianDay( (self.v/self.__class__.TICKS_SCET_DAY+self.__class__.ZERO_SCET))
+
+ def pseudoObt(self) :
+ return PseudoOBT((self.v/self.__class__.TICKS_SCET_DAY+self.__class__.ZERO_SCET-self.__class__.ZERO_PseudoOBT)*self.__class__.TICKS_PseudoOBT_DAY)
+
+ def utcsec(self) :
+ return UTCsec( (self.v/self.__class__.TICKS_SCET_DAY+self.__class__.ZERO_SCET-self.__class__.ZERO_UTC)*self.__class__.TICKS_UTC_DAY )
+
+ def iso(self) :
+ return self.julianday().iso()
+
+ def zulu(self) :
+ return self.julianday().zulu()
+
+
+
+class TAI(_TIME_SCALE_BASE) :
+ """ class to handle TAI : same than UTC but from 1958 Jan 01, 00:00:00 """
+ def __init__(self,v=0.) :
+ _TIME_SCALE_BASE.__init__(self,v)
+ # if v is not a numeric scalar then v has been coded as a JD and it has to be converted
+ if not self._isNumericScalar(v) :
+ self.v = (self.v - self.__class__.ZERO_TAI)*self.__class__.TICKS_TAI_DAY
+
+ def kind(self,test='') :
+ if test == '' :
+ return 'TAI'
+ return test == 'TAI'
+
+ def julianday(self) :
+ return JulianDay( (self.v/self.__class__.TICKS_TAI_DAY+self.__class__.ZERO_TAI))
+
+ def pseudoObt(self) :
+ return PseudoOBT((self.v/self.__class__.TICKS_TAI_DAY+self.__class__.ZERO_TAI-self.__class__.ZERO_PseudoOBT)*self.__class__.TICKS_PseudoOBT_DAY)
+
+ def utcsec(self) :
+ return UTCsec( (self.v/self.__class__.TICKS_TAI_DAY+self.__class__.ZERO_TAI-self.__class__.ZERO_UTC)*self.__class__.TICKS_UTC_DAY )
+
+ def scet(self) :
+ return SCET( (self.v/self.__class__.TICKS_TAI_DAY+self.__class__.ZERO_TAI-self.__class__.ZERO_SCET) * self.__class__.TICKS_SCET_DAY )
+
+ def iso(self) :
+ return self.julianday().iso()
+
+ def zulu(self) :
+ return self.julianday().zulu()
+
+
+
+def _TEST_time_support() :
+ """ test and examples """
+
+ x = JulianDay(2400000.5)
+ print
+ print x.kind(),x, "input ad number "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.calendardate().kind()),x.calendardate()
+ print " %10s -> %10s : " % (x.kind(),type(x.tuple())),x.tuple()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+
+ x = JulianDay((1858,11,17,0,0,0))
+ print
+ print x.kind(),x,"input as tuple "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.calendardate().kind()),x.calendardate()
+ print " %10s -> %10s : " % (x.kind(),type(x.tuple())),x.tuple()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+
+ x = JulianDay("1858-11-17T00:00:00")
+ print
+ print x.kind(),x,"input as string "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.calendardate().kind()),x.calendardate()
+ print " %10s -> %10s : " % (x.kind(),type(x.tuple())),x.tuple()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+
+ x = PseudoOBT(0)
+ print
+ print x.kind(),x,"input as number "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+
+
+ x = PseudoOBT((1958,1,1,0,0,0))
+ print
+ print x.kind(),x,"input as tuple "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+
+
+ x = PseudoOBT("1958-01-01T00:00:00")
+ print
+ print x.kind(),x,"input as string "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+
+
+ x = UTCsec(0)
+ print
+ print x.kind(),x,"input as number "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+
+
+ x = UTCsec((1970,1,1,0,0,0))
+ print
+ print x.kind(),x,"input as tuple "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+
+
+ x = UTCsec("1970-01-01T00:00:00")
+ print
+ print x.kind(),x,"input as string "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+
+
+ x = SCET(0)
+ print
+ print x.kind(),x,"input as number "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+
+
+ x = SCET((1970,1,1,0,0,0))
+ print
+ print x.kind(),x,"input as tuple "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+
+
+ x = SCET("1970-01-01T00:00:00")
+ print
+ print x.kind(),x,"input as string "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.tai().kind()),x.tai()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+
+
+ x = TAI(0)
+ print
+ print x.kind(),x,"input as number "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+
+
+ x = TAI((1970,1,1,0,0,0))
+ print
+ print x.kind(),x,"input as tuple "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+
+
+ x = TAI("1970-01-01T00:00:00")
+ print
+ print x.kind(),x,"input as string "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.pseudoObt().kind()),x.pseudoObt()
+ print " %10s -> %10s : " % (x.kind(),x.utcsec().kind()),x.utcsec()
+ print " %10s -> %10s : " % (x.kind(),x.scet().kind()),x.scet()
+
+
+ x = ISO("1970-01-01T00:00:00")
+ print
+ print x.kind(),x,"input as string "
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.calendardate().kind()),x.calendardate()
+ print " %10s -> %10s : " % (x.kind(),type(x.tuple())),x.tuple()
+
+ x = ZULU("1970-01-01T00:00:00")
+ print
+ print x.kind(),x,"input as string "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.calendardate().kind()),x.calendardate()
+ print " %10s -> %10s : " % (x.kind(),type(x.tuple())),x.tuple()
+
+ x = CalendarDate((1970,1,1,0,0,0.))
+ print
+ print x.kind(),x,"input as tuple "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),type(x.tuple())),x.tuple()
+
+ x = CalendarDate((2009,5,14,10,13,0.)).julianday()
+ print
+ print x.kind(),x,"input as scalar "
+ print " %10s -> %10s : " % (x.kind(),x.iso().kind()),x.iso()
+ print " %10s -> %10s : " % (x.kind(),x.zulu().kind()),x.zulu()
+ # print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),type(x.tuple())),x.tuple()
+ print " %10s -> %10s : " % (x.kind(),x.od().kind()),x.od()
+
+ x = CalendarDate((2009,5,14,10,13,0.)).julianday().od()
+ print
+ print x.kind(),x,"input as scalar "
+ print " %10s -> %10s : " % (x.kind(),x.julianday().kind()),x.julianday()
+ print " %10s -> %10s : " % (x.kind(),x.julianday().calendardate().kind()),x.julianday().calendardate()
+
+# end
+#_TEST_time_support()
+
+if __name__ == '__main__' :
+ # creates the obj by od database
+ import os
+ import sys
+ import shutil
+ from optparse import OptionParser
+
+ class _CmdLine() :
+ def __init__(self) :
+ self.ArgList = ['in','out']
+
+ cmdln = OptionParser(usage="usage: %prog [options] "+
+ " ".join(self.ArgList)
+ +"[options]",version='1.0 - 2010 Jun 25 - M.Maris')
+ cmdln.add_option("-d","--dry-run",dest="DryRun", action = 'store_true',default=False,help='Performs a dry run')
+ cmdln.add_option("-T","--test",dest="ACTION_TEST", action = 'store_true',default=False,help='run a test')
+
+ (self.Options,self.Args) = cmdln.parse_args()
+
+ if len(self.Args) != len(self.ArgList) and not self.Options.ACTION_TEST :
+ print
+ print "TIME_SUPPORT "
+ cmdln.print_version()
+ cmdln.print_help()
+ print
+ print "Allowed codes for times (not case sensitive): "
+ print " ISO, ZULU, JD, SCET, TAI, PseudoOBT, PseudoOD"
+ print
+ print "Example:"
+ print " > time_support.py PseudoOD=0 ISO "
+ print " 2009-05-14T13:11:59.9999839"
+ print " > time_support.py PseudoOD=0 Z "
+ print " 2009-05-14T13:11:59.9999839Z"
+ print
+ sys.exit(1)
+ else :
+ if len(self.Args) > 0 :
+ self.IN = self.Args[0]
+ self.OUT = self.Args[1]
+ else :
+ self.IN = ""
+ self.OUT = ""
+
+ if len(self.IN) == 0 :
+ self.IN_TYPE = ''
+ self.IN_VALUE = ''
+ else :
+ ll = self.IN.split('=')
+ if len(ll) != 2 :
+ sys.exit(1)
+ else :
+ self.IN_TYPE = ll[0]
+ self.IN_VALUE = ll[1]
+ self.DryRun()
+
+ def DryRun(self) :
+ if self.Options.DryRun :
+ print "Dry Run\nExit\n"
+ sys.exit(1)
+
+
+ # execution of commands
+ cmdline = _CmdLine()
+ if cmdline.Options.ACTION_TEST :
+ _TEST_time_support()
+ sys.exit(1)
+
+ JD = None
+ if cmdline.IN_TYPE.lower() == 'ISO'.lower() :
+ JD = JulianDay(cmdline.IN_VALUE)
+
+ if cmdline.IN_TYPE.lower() == 'ZULU'.lower() :
+ JD = ZULU(cmdline.IN_VALUE).julianday()
+
+ if cmdline.IN_TYPE.lower() == 'JD'.lower() :
+ JD = JulianDay(float(cmdline.IN_VALUE))
+
+ if cmdline.IN_TYPE.lower() == 'PseudoOBT'.lower() :
+ JD = PseudoOBT(float(cmdline.IN_VALUE)).julianday()
+
+ if cmdline.IN_TYPE.lower() == 'TAI'.lower() :
+ JD = TAI(float(cmdline.IN_VALUE)).julianday()
+
+ if cmdline.IN_TYPE.lower() == 'SCET'.lower() :
+ JD = SCET(float(cmdline.IN_VALUE)).julianday()
+
+ if cmdline.IN_TYPE.lower() == 'PseudoOD'.lower() :
+ JD = OD(float(cmdline.IN_VALUE)).julianday()
+
+ if JD == None :
+ print "Unknown input type"
+
+ if cmdline.OUT.lower() == 'ISO'.lower() :
+ print JD.iso()
+ sys.exit(1)
+
+ if cmdline.OUT.lower() == 'PseudoOBT'.lower() :
+ print JD.pseudoObt()
+ sys.exit(1)
+
+ if cmdline.OUT.lower() == 'TAI'.lower() :
+ print JD.tai()
+ sys.exit(1)
+
+ if cmdline.OUT.lower() == 'SCET'.lower() :
+ print JD.scet()
+ sys.exit(1)
+
+ if cmdline.OUT.lower() == 'ZULU'.lower() :
+ print JD.zulu()
+ sys.exit(1)
+
+ if cmdline.OUT.lower() == 'PseudoOD'.lower() :
+ print JD.od()
+ sys.exit(1)
+
+ if cmdline.OUT.lower() == 'JD'.lower() :
+ print JD.v
+ sys.exit(1)
+
+ print "Unknown output type"
+
diff --git a/src/yapsut/planck_legacy/vec3dict.py b/src/yapsut/planck_legacy/vec3dict.py
new file mode 100644
index 0000000000000000000000000000000000000000..95596066511a7b145deae4198c19015bf896a514
--- /dev/null
+++ b/src/yapsut/planck_legacy/vec3dict.py
@@ -0,0 +1,320 @@
+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
+