- sympy: I used this package to build the final complicated expression that gives the heat transfer. Unfortuately, it's immature software. For example it couldn't solve: sympy.solve((x**2)**0.5-9). Also the lambda features were confusing since there were two lambda features in the package.
- numexpr: 1.4 I used this blazingly fast expression evaluation package to evaluate the built expression.
11/8: It's published on nanohub.org.
Prototype GUI built using hubzero's Rappture GUI builder |
#Author: Majid S. al-Dosari #http://www.thermalfluidscentral.org/encyclopedia/index.php #/Near-field_radiative_transfer_between_two_semi-infinite_media #Calculate Bulk 1-D Near-Field Radiation between two Bulks. #Input: #1. permittivity of both bulks #2. distance b/w them #3. Temperature of each bulk T1,T2 #Output: #Power as a function of frequency or just total power #physics notes #the integrand is split into TE and TM components #these components are further split into propagating (omega/c<1 data-blogger-escaped-and="and" data-blogger-escaped-c="c" data-blogger-escaped-evanescent="evanescent" data-blogger-escaped-omega="omega" data-blogger-escaped-waves="waves">1) #computation procedure #1. build expression for TE(S-polarization) and TM(P-polarization) waves. #"s" in the web reference. #See 'sept' and 'sest' names below as examples #2. convert sympy expression to numexpr #3. use evalq for total power, evalqw for power as a function of frequency, #evalqwsections to get power as a function of frequency when the practical #upper limit of integration in wavevector is unknown (infinity) #they share the common argument sequence #(stringformofintegrand,beta1,beta2,T1,T2,nb=1000,omega1=0,no=1000) #beta here is normalized to omega/c. omega2 is determined as that frequency at #which 99.99...?% blackbody radiation integral is covered at the larger #of T1, and T2. #nb and no are the number of points discretized on beta and omega respectively #don't include beta=1 in the range b/w beta1 and beta2. #use beta2=.9999 as an upperlimit for integrating propagating waves and #use beta1=1.0001 as a lower limit for integrating evanescent waves. #for evalq and evalqwsections, if beta2 is in the evanescent range (>1), it will #be taken as an integration 'chunck' size and will continue integration until #asymptotes. be careful, because if you know TM waves's contribution #are negligible it may continue to integrate to ridiculous wavevectors and #therefore give unrealistic results. there is an absolute limit for beta2 # around a few interatomic distances #programdata is a dictionary that holds frequencies and wavevectors that the #integrand was evaluated over. i use them for plot axes. #example: evalq(sym2num(str(sept)),1.001,50,300,299,no=3e3,nb=3e2) #use the output from evals to make a pretty picture of the evaluated func #eg use matplotlib's contourf(swapaxis(evals output)) #coding note: #a=1+2j #numexpr.evaluate('(a*1j)+0') #gives TypeError: no ordering relation is defined for complex numbers #commented out line 263 #ordered_constants.sort() in necompiler #see 'HACK' in sym2num #evaluate new option to sub output inplace of input import numpy import sympy import numexpr #might be able to get away w/ 100sx100s grid #.0001 and .9999 , 1.00001 #Material Permittivity Input #--- def invcm2rad(ic):return ic*2*numpy.pi*2.998e10 consts={'c':2.998e8,'hbar':1.054571596e-34,'kb':1.380650277e-23 ,'omega_LO':invcm2rad(969),'omega_TO':invcm2rad(793),'Gamma':invcm2rad(4.76) ,'eps_inf':6.7} eps_inf,omega_LO ,omega_TO ,Gamma,omega=\ sympy.symbols('eps_inf omega_LO omega_TO Gamma omega') eps_SiC=eps_inf*(1+(omega_LO**2-omega_TO**2)\ /(omega_TO**2-1j*Gamma*omega-omega**2)) #eps_Au=1-(13.71e15**2)/(omega*(omega+2*numpy.pi*1j*4.05e13)) #indexed epsilon eps=[1,eps_SiC,eps_SiC] #eps=[1,eps1,eps2] beta,c=sympy.symbols('beta c') beta=beta*(omega/c) #in terms of omega/c def gama(j):return (eps[j]*((omega/c)**2)-beta**2)**.5 #def gama(j):return (eps[j]-(beta/(omega/c))**2 )**.5 #in terms of omega/c #0j to fix ValueError: negative number cannot be raised to a fractional power #not needed with sympy i think def r0_s(j):return (gama(0)-gama(j))*(gama(0)+gama(j)) def r0_p(j):return (eps[j]*gama(0)-gama(j))/(eps[j]*gama(0)+gama(j)) r0=dict({'s':r0_s,'p':r0_p}) d=sympy.Symbol('d') def sevan(sorp): return sympy.im(r0[sorp](1))*sympy.im(r0[sorp](2)) \ *beta*sympy.exp(-2*sympy.im(gama(0))*d)/ \ sympy.abs(1-r0[sorp](1)*r0[sorp](2)*sympy.exp(-2*sympy.im(gama(0))*d))**2 #the sympy on nanohub doesn't have sympy.Abs rather it has sympy.abs def sprop(sorp): return (1.0/4)*beta*(1-sympy.abs(r0[sorp](1))**2)*(1-sympy.abs(r0[sorp](2))**2)\ /sympy.abs(1-r0[sorp](1)*r0[sorp](2)*sympy.exp((2j)*gama(0)*d))**2 #the s funcs /give/ the expr. they aren't funcs themselves hbar,T,kb=sympy.symbols('hbar T kb') Thetaexpr=hbar*omega/(sympy.exp(hbar*omega/(kb*T))-1) Thetaexpr=Thetaexpr.subs(consts) def Theta(omega,T): if T==0: return 0 return Thetaexpr.subs({'T':T,'omega':omega}) def q_omega(omega,T1,T2,intofs): return ((1./sympy.pi)**2)*(Theta(omega,T1)-Theta(omega,T2))*intofs #import scipy #from scipy import integrate #useless def intofs(sfunc,beta1,beta2): #sfunc1=lambda b: sfunc(b+0j)#needed for numexpr #return scipy.integrate.quad(sfunc1,beta1,beta2)[0] #slowwwww #i made this func bc i could use scipy quad or something else return sympy.mpmath.quad(sfunc,[beta1,beta2]) #compress const array #evaluated numerically from the start in a simple way def nintofs(numstr,beta1,beta2,omega,num=100):#sfunc w/o omega :( #when i subs in an omega i get weird results [beta,omega],db,do=geninputarray(beta1,beta2,num,omega,1234,1) return db*numpy.sum(numexpr.evaluate(numstr)) def evals_(numstr,*args,**kwargs):#sfunc w/o omega :( #when i subs in an omega i get weird results [beta,omega],db,do=geninputarray(*args,**kwargs) return (numexpr.evaluate(numstr)),db,do #with zeros you get NaNs #swap(reshape()) def evals(sstr,*args,**kwargs): so,db,do=evals_(sstr,*args,**kwargs) so=numpy.reshape(so,(args[2],args[5])) #swap args for reshape other way for contourf, then swapaxis return so,db,do def evalDTs(DTstr,omegas):#sfunc w/o omega :( #when i subs in an omega i get weird results omega=omegas if omega[0]<.1:omega[0]=.1 #get's rid of NaNs return numexpr.evaluate(DTstr) #the limits of integration planckexpr=(hbar*omega**3/(4*sympy.pi**3*c**2))\ /(sympy.exp(hbar*omega/(kb*T))-1) planckexpr=planckexpr.subs(consts) def planckd(o,T):return planckexpr.subs({'omega':o,'T':T}) #dplanckexpr=planckexpr.diff(omega) #from scipy import optimize def omegalim(T): return consts['c']*2*numpy.pi/(1000e-6/T) #-.000321+1 of BB #% this might have an effect def solvebeta(omega): return omega/consts['c'] #return sympy.solve(gama(0).subs(consts).subs({'omega':omega}))[0] #it can't solve this! def betalim(T): return omegalim(T)/consts['c'] def evalqwint(sfunc,beta1,beta2,T1,T2,omega1=1e-20,no=500): os=numpy.linspace(omega1,omegalim(max(T1,T2)),no) sfb=lambda o:sfunc.subs({'omega':o})#gives func of sympy spec=[intofs( lambda b:sfb(o).evalf(subs={'beta':b}) ,beta1,beta2) for o in os] #spec=[nintofs( sym2num(str(sfunc)) ,beta1,beta2,o,num=1e3) for o in os] #np.random.seed(123) #DTs=np.random.rand(len(os)) #this is introducing a good amoutn of error in SiC for d=100e-9 DTs=(evalDTs(sym2num(str((Theta(omega,T1)-Theta(omega,T2)))) ,os)) qw= spec*DTs*(os/consts['c'])/numpy.pi**2# no 2pi ?? return qw,os[1]-os[0] #dw def evalqw(sstr,beta1,beta2,T1,T2,nb=1000,omega1=0,no=1000): #evals(sym2num(str(sppt+spst)),0,.999999,1000,0,omegalim(300),1000) so=evals(sstr,beta1,beta2,nb,omega1,omegalim(max(T1,T2)),no) #so=evals(sstr,beta1,beta2,nb,1.7e14,1.85e14,no) os=numpy.real(programdata['os']) dBetas=so[1]*(os/consts['c']) #np.random.seed(123) #DTs=np.random.rand(len(os)) #this is introducing a good amoutn of error in SiC for d=100e-9 DTs=(evalDTs(sym2num(str((Theta(omega,T1)-Theta(omega,T2)))) ,os)) qw= numpy.sum(numpy.real(so[0]),axis=0)*DTs*dBetas/numpy.pi**2# no /2pi ?? return qw,so[2] #dw import scipy from scipy import integrate def evalq(*args,**kwargs): if args[2]<=1: o=evalqw(*args,**kwargs) #for propagating waves else: o=evalqwsections(*args,**kwargs) #return sum(o[0])*o[1] return scipy.integrate.trapz(o[0],dx=o[1]) def evalqwsections(*args,**kwargs): #take beta2 as sectioning #if beta2>1: args=list(args) interval=args[2]-args[1] ac=asymptote() passn=1 xlim=(numpy.pi/.5e-9)/(omegalim(max(args[3],args[4]))/consts['c'])#~=7e3 #(2*3.14/10e-10)/ #per zhang review of NF thermal rad #say 10angstrom is local limit while 1: if args[1]>xlim:break o=evalqw(*args,**kwargs) if passn==1:cur=numpy.zeros(len(programdata['os'])) cur+=o[0] sm=numpy.sum(o[0])*o[1] if ac.didit(sm)==True:break passn+=1 #shift for next interval args[1]+=interval;args[2]+=interval; return cur, o[1] class asymptote(): #expectation of the trend is that it goes up then goes down in the end #it could vary up and down in b/w def __init__(self): self.nums=[] self.curmax=float("-inf") def didit(self,num): tol=.01 if len(self.nums)==0:self.curmax=num self.nums.append(num) if self.up()==True: self.curmax=num return False if self.up()==None: return None if self.up()==False and self.curmax*tol>num:return True else: return False def up(self): if len(self.nums)<2:return data-blogger-escaped-if="if" data-blogger-escaped-none="none" data-blogger-escaped-self.nums="self.nums">0:return True else: return False onejay=1j def sym2num(astr): """ to go from a sympy expr to a numexpr expr """ astr=astr.replace('I','1j') astr=astr.replace('1j','onejay') #A HACK!!! astr=astr.replace('im','imag') astr=astr.replace('re','real') astr=astr.replace('Abs','abs') return astr #a hold for intermediate data programdata={} def geninputarray(betamin,betamax,bnum,omegamin,omegamax,onum):#,num=None): #onum=100,bnum=100,num=None # ,betamin=1e-20,betamax=2e6,omegamin=1e-20,omegamax=2e14): """inputs in numexpr extreme:say you want just one line instead of 2D: specify a close bmin and bmax with bnum=1 betamin will always be evaluated. to evaluate at a certain beta put it in betamin """ # if num!=None:#if num specifided # bnum=onum=num #bs=numpy.add(numpy.linspace(betamin,betamax,num=num),0j) bs=numpy.zeros(bnum,dtype='c8') #could use bigger float but i dont care bs+=numpy.linspace((betamin),(betamax),num=bnum,endpoint=False) #add 0j needed to work in expr os=numpy.zeros(onum,dtype='c8') os+=numpy.linspace((omegamin),(omegamax),num=onum,endpoint=False) if os[0]==0:os[0]=1e-20 if bs[0]==0:bs[0]=1e-20 programdata.update({'os':os,'bs':bs}) return (cartesian([bs,os]).transpose() ,(float(betamax)-betamin)/len(bs),(float(omegamax)-omegamin)/len(os)) #fxy eval, dx*dy #evanescent waves test expressions sep=sevan('p').subs(consts) sept=sep.subs({'d':100e-9}) ses= sevan('s').subs(consts) sest=ses.subs({'d':100e-9}) #propagatin waves test expressions spp=sprop('p').subs(consts) sppt= spp.subs({'d':1}) sps=sprop('s').subs(consts) spst=sps.subs({'d':1}) import numpy as np #WOW! def cartesian(arrays, out=None): #http://stackoverflow.com/questions/1208118/using-numpy-to-build-an-array-of-all-combinations-of-two-arrays """ Generate a cartesian product of input arrays. Parameters ---------- arrays : list of array-like 1-D arrays to form the cartesian product of. out : ndarray Array to place the cartesian product in. Returns ------- out : ndarray 2-D array of shape (M, len(arrays)) containing cartesian products formed of input arrays. Examples -------- >>> cartesian(([1, 2, 3], [4, 5], [6, 7])) array([[1, 4, 6], [1, 4, 7], [1, 5, 6], [1, 5, 7], [2, 4, 6], [2, 4, 7], [2, 5, 6], [2, 5, 7], [3, 4, 6], [3, 4, 7], [3, 5, 6], [3, 5, 7]]) """ arrays = [np.asarray(x) for x in arrays] dtype = arrays[0].dtype n = np.prod([x.size for x in arrays]) if out is None: out = np.zeros([n, len(arrays)], dtype=dtype) m = n / arrays[0].size out[:,0] = np.repeat(arrays[0], m) if arrays[1:]: cartesian(arrays[1:], out=out[0:m,1:]) for j in xrange(1, arrays[0].size): out[j*m:(j+1)*m,1:] = out[0:m,1:] return out