diff --git a/src/yapsut/blackbody_mks.py b/src/yapsut/blackbody_mks.py new file mode 100644 index 0000000000000000000000000000000000000000..baf4f9bb46eeb3b60441d31c3b7f55fcc9deeae3 --- /dev/null +++ b/src/yapsut/blackbody_mks.py @@ -0,0 +1,335 @@ +__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() +