I'm going to be doing FDTD stuff but need a baseline transport calculation using a closed-form
expression. Maple, Maxima, Sympy couldn't do the integration. So I had to do it numerically in a more manual way. So I used my trusty Python environment with the following modules that did most of the work:
-
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.
(I had to comment line 263 in necompiler.py because I was working with imaginary numbers. It's an optimization step that needs to order numbers and there is no order to imaginary numbers) See better workaround.
11/8: It's
published on
nanohub.org.
 |
| Prototype GUI built using hubzero's Rappture GUI builder |
The code for the core of the program is below.
#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
