u

src/yapsut/blackbody_mks.py

+__DESCRIPTION__=""" 
+MKS       = an object returning the CGS units and constants
+
+See: 
+   import blackbody_mks 
+   help(blackbody.MKS)
+
+The package manages numpy.array().
+
+Version 1.0 - 2011 Set 1 - 2012 Apr 5 - M.Maris
+"""
+
+import numpy as np
+
+from yapsut import MKS
+
+class _CnPlanck :
+   def __call__(self,n,x) :
+      return self.bx(n,x)
+   def Norm(self,n) :
+      """normalization constant of x^n/(exp(x)-1) integrated between 0 and 1"""
+      from scipy.special import gamma, zeta 
+      if n==0 :
+         return 1
+      elif n==1 : 
+         return np.pi**2/12
+      elif n==2 :
+         # gamma(3)*zeta(3)
+         return 2*zeta(3)
+      elif n==3 :
+         # gamma(4)*zeta(4)
+         return np.pi**4/15
+      else :
+         # gamma(n+1)*zeta(n+1)
+         return gamma(n+1)*zeta(n+1)
+   def bx(self,n,x) :
+      return x**n/(np.exp(x)-1)
+   def Tcumulant(self,n,xmax,nsteps) :
+      x=np.linspace(0.,xmax,nsteps)
+      dx=xmax/float(nsteps)
+      bx=np.zeros(nsteps)
+      bx[1:]=self.bx(n,x[1:])
+      if n == 0 :
+         bx[0]=+np.infty
+      elif n==1 :
+         bx[0]=1.
+      #else :
+      #   bx[0]=0
+      out=np.zeros(nsteps)
+      out[1:]=(bx[1:]+bx[:-1]).cumsum()
+      out*=dx/2/self.Norm(n)
+      return x,out
+   def Scumulant(self,n,xmax=708,nsteps=300000) :
+      """ returns a spline interpolating the cumulant of the planck function """
+      class scumspline :
+         """ spline interpolation for cumulants of x^n/(exp(2)-1) 
+         
+             with 300000 points typical accuracy some 10^-6
+             with 3000000 points typical accuracy some 10^-7
+             
+             produced through Scumulant of AdimPlanckFunction
+         """
+         @property
+         def n(self) :
+            return self._n
+         @property
+         def norm(self) :
+            return self._norm
+         @property
+         def tck(self) :
+            return self._tck
+         @property
+         def xmax(self) :
+            return self._xmax
+         @property
+         def max(self) :
+            return self._max
+         @property
+         def cumulant_accuracy(self) :
+            return 1-self._max
+         def __init__(self,x,y,norm,n) :
+            from scipy.interpolate import splrep
+            self._n=n
+            self._nsamp=len(x)
+            self._tck,self._fp,self._ierr,self._msg=splrep(x,y,full_output=True)
+            self._xmax=x.max()
+            self._norm=norm
+            self._max=self(self.xmax)
+         def __call__(self,x) :
+            """ is x > xmax or x < 0 a nan is produced 
+            x : either a number or a numpy array
+            """
+            from scipy.interpolate import splev
+            _x=np.array(x)
+            out=splev(_x,self._tck)
+            out[(_x>self.xmax)+(_x<0)]=np.nan
+            return out
+         def __len__(self) :
+            return self._nsamp
+         def copy(self) :
+            import copy
+            return copy.deepcopy(self)
+         def pdf(self,x) :
+            """ computes the pdf(n) = x^n/(np.exp(x)-1)/norm """
+            return x**self._n/(np.exp(x)-1)/self._norm
+         def __repr__(self) :
+            return "<"+"n = "+str(self.n)+"; xmax = "+str(self.xmax)+"; nsmp = "+str(len(self))+"; ierr = "+str(self._ierr)+"; accuracy = "+str(self.cumulant_accuracy)+">"
+      #
+      x,y=self.Tcumulant(n,xmax,nsteps)
+      return scumspline(x,y,self.Norm(n),n) 
+   def cumulant(self,n,xmax,nsteps=1000) :
+      x=np.linspace(0,xmax,nsteps)
+      dx=x[1]-x[0]
+      bx=np.zeros(nsteps)
+      bx[1:]=self.bx(n,x[1:])
+      if n == 0 :
+         b[0]=+np.infty
+      elif n==1 :
+         b[0]=1.
+      return dx*bx.sum()/self.Norm(n)
+   @property 
+   def xwien(self) :
+      return 4.965114231744276303
+   @property 
+   def nu_wien(self) :
+      return MKS.kB/MKS.h*self.xwien
+   @property 
+   def lambda_wien(self) :
+      return MKS.kB/(MKS.c*MKS.h)*self.xwien
+
+Wien_x      = 4.965114231744276303
+Wien_nu     = MKS.kB/MKS.h*Wien_x
+Wien_lambda = MKS.kB/(MKS.c*MKS.h)*Wien_x
+
+def MeanPhotonEnergy(T) :
+   """ returns the mean photon energy in J for a bb at temperature T """
+   #return 3.7294e-23*T
+   return MKS.meanBBPhotEn*T
+
+# Wien
+#Wien_x=4.965114231744276303
+#Wien_nu=MKS.kB/MKS.h*Wien_x
+#Wien_lambda=(MKS.c*MKS.h)/(MKS.kB*Wien_x)
+
+PlanckCn=_CnPlanck()
+#PlanckC2=PlanckCn.Scumulant(2)
+#PlanckC3=PlanckCn.Scumulant(3)
+#PlanckC4=PlanckCn.Scumulant(4)
+#PlanckC5=PlanckCn.Scumulant(5)
+      
+class _black_body_mks() :
+   def __init__(self,Type) :
+      if Type == 'bbn' :
+         self.A=2*MKS.h/MKS.c**2
+         self.NA=2/MKS.c**2/MKS.Navo
+         self.B=MKS.h/MKS.kB
+      elif Type == 'sedn' :
+         self.A=2*MKS.h/MKS.c**2*MKS.rad2sed
+         self.NA=2/MKS.c**2/MKS.Navo*MKS.rad2sed
+         self.B=MKS.h/MKS.kB
+      elif Type == 'bbl' :
+         self.A=2*MKS.h*MKS.c**2
+         self.NA=2*MKS.c/MKS.Navo
+         self.B=MKS.h*MKS.c/MKS.kB
+      elif Type == 'sedl' :
+         self.A=2*MKS.h*MKS.c**2*MKS.rad2sed
+         self.NA=2*MKS.c/MKS.Navo*MKS.rad2sed
+         self.B=MKS.h*MKS.c/MKS.kB
+      else :
+         raise Exception("Black Body Type %s is not valid"%Type,"")
+   def valid(self,FreqTHz,T) :
+      return (FreqTHz >= 0.) * (T >= 0.)
+
+class _bbl_mks(_black_body_mks) :
+   """
+!
+! BB(lambda,T) thermal radiance function mks
+!
+! lambda in m
+!
+! bbl_mks = J/(m2 s m sterad)
+!
+   """
+   def __init__(self):
+      _black_body_mks.__init__(self,'bbl')
+   def __call__(self,Lambda,T) :
+      FT = self.A/Lambda**5
+      ET = np.exp( self.B/(Lambda*T) )-1;
+      return FT / ET
+bbl_mks=_bbl_mks()      
+
+class _pi_bbl_mks(_black_body_mks) :
+   """
+!
+! pi*BB(lambda,T) thermal radiance function mks
+!
+! lambda in m
+!
+! bbl_mks = J/(m2 s m sterad)
+!
+   """
+   def __init__(self):
+      _black_body_mks.__init__(self,'bbl')
+   def __call__(self,Lambda,T) :
+      FT = np.pi*self.A/Lambda**5
+      ET = np.exp( self.B/(Lambda*T) )-1;
+      return FT / ET
+pbbl_mks=_pi_bbl_mks()      
+
+class _bbl_sed_mks(_black_body_mks) :
+   """
+!
+! BB(lambda,T) thermal sed function mks
+!
+! lambda in m
+!
+! bbl_mks = J/(m2 s m sterad)
+!
+   """
+   def __init__(self):
+      _black_body_mks.__init__(self,'sedl')
+   def __call__(self,Lambda,T) :
+      FT = self.A/Lambda**5
+      ET = np.exp( self.B/(Lambda*T) )-1;
+      return FT / ET
+bbl_sed_mks=_bbl_sed_mks()
+
+class _Npbbl_mks(_black_body_mks) :
+   """
+!
+! NpBB(lambda,T) thermal photons function mks
+!
+! lambda in m
+!
+! Nbbl_mks = photons/(m2 s m)
+!
+   """
+   def __init__(self):
+      _black_body_mks.__init__(self,'bbl')
+   def __call__(self,Lambda,T) :
+      FT = np.pi*self.NA/Lambda**4
+      ET = np.exp( self.B/(Lambda*T) )-1;
+      return FT / ET
+Npbbl_mks=_Npbbl_mks()      
+
+class _Nbbl_sed_mks(_black_body_mks) :
+   """
+!
+! NBB(lambda,T) thermal radiance function mks
+!
+! lambda in m
+!
+! Nbbl_mks = photons/(m2 s m sterad)
+!
+   """
+   def __init__(self):
+      _black_body_mks.__init__(self,'sedl')
+   def __call__(self,Lambda,T) :
+      FT = self.NA/Lambda**4
+      ET = np.exp( self.B/(Lambda*T) )-1;
+      return FT / ET
+Nbbl_sed_mks=_Nbbl_sed_mks()      
+
+class _bbn_mks(_black_body_mks) :
+   """
+!
+! BB(nu,T) thermal radiance function mks
+!
+! nu in Hz
+!
+! bbn_mks = J/(m2 s Hz sterad)
+!
+   """
+   def __init__(self):
+      _black_body_mks.__init__(self,'bbn')
+   def __call__(self,nu,T) :
+      FT = self.A*nu**3
+      ET = np.exp(self.B*nu/T)-1
+      return FT / ET
+bbn_mks=_bbn_mks()      
+
+class _bbn_sed_mks(_black_body_mks) :
+   """
+!
+! BB(nu,T) thermal sed function mks
+!
+! nu in Hz
+!
+! bbn_mks = J/(m2 s Hz sterad)
+!
+   """
+   def __init__(self):
+      _black_body_mks.__init__(self,'sedn')
+   def __call__(self,nu,T) :
+      FT = self.A*nu**3
+      ET = np.exp(self.B*nu/T)-1
+      return FT / ET
+bbn_sed_mks=_bbn_sed_mks()      
+
+class _Nbbn_mks(_black_body_mks) :
+   """
+!
+! NBB(nu,T) photons radiance mks = bbn_mks / h nu
+!
+! nu in Hz
+!
+! Nbbn_mks = photons/(m2 s Hz sterad) 
+!
+   """
+   def __init__(self):
+      _black_body_mks.__init__(self,'bbn')
+   def __call__(self,nu,T) :
+      FT = self.NA*nu**2
+      ET = np.exp(self.B*nu/T)-1
+      return FT / ET
+Nbbn_mks=_bbn_mks()      
+
+class _Nbbn_sed_mks(_black_body_mks) :
+   """
+!
+! NBB(nu,T) photons sed mks = bbn_mks / h nu
+!
+! nu in Hz
+!
+! Nbbn_mks = photons/(m2 s Hz sterad) 
+!
+   """
+   def __init__(self):
+      _black_body_mks.__init__(self,'sedn')
+   def __call__(self,nu,T) :
+      FT = self.NA*nu**2
+      ET = np.exp(self.B*nu/T)-1
+      return FT / ET
+Nbbn_sed_mks=_Nbbn_sed_mks()      
+