Friday, July 6, 2012

Task 6: Basline Near-Field Calculation

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