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 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

Prototype GUI built using hubzero's Rappture GUI builder
The code for the core of the program is below.

#Author: Majid S. al-Dosari

#Calculate Bulk 1-D Near-Field Radiation between two Bulks.
#1. permittivity of both bulks
#2. distance b/w them
#3. Temperature of each bulk T1,T2
#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
#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:
#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
eps_inf,omega_LO ,omega_TO ,Gamma,omega=\
sympy.symbols('eps_inf omega_LO omega_TO Gamma omega')


#indexed epsilon

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))

def sevan(sorp):
    return[sorp](1))*[sorp](2)) \
        *beta*sympy.exp(-2**d)/ \
#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)\
#the s funcs /give/ the expr. they aren't funcs themselves

hbar,T,kb=sympy.symbols('hbar T kb')
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
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
    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
    return (numexpr.evaluate(numstr)),db,do #with zeros you get NaNs
def evals(sstr,*args,**kwargs):
    #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
    if omega[0]<.1:omega[0]=.1 #get's rid of NaNs
    return numexpr.evaluate(DTstr)

#the limits of integration
def planckd(o,T):return planckexpr.subs({'omega':o,'T':T})

#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):
    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]
    #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):
    #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:
    #(2*3.14/10e-10)/ #per zhang review of NF thermal rad    
    #say 10angstrom is local limit
    while 1:
        if args[1]>xlim:break
        if passn==1:cur=numpy.zeros(len(programdata['os']))
        if ac.didit(sm)==True:break
        #shift for next 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):
    def didit(self,num):
        if len(self.nums)==0:self.curmax=num
        if self.up()==True:
            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

def sym2num(astr):
     to go from a sympy expr to a numexpr expr
    astr=astr.replace('1j','onejay') #A HACK!!!
    return astr

#a hold for intermediate data

def geninputarray(betamin,betamax,bnum,omegamin,omegamax,onum):#,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.zeros(bnum,dtype='c8') #could use bigger float but i dont care
    #add 0j needed to work in expr
    if os[0]==0:os[0]=1e-20
    if bs[0]==0:bs[0]=1e-20
    return (cartesian([bs,os]).transpose()
    #fxy eval, dx*dy

#evanescent waves test expressions
ses=  sevan('s').subs(consts)
#propagatin waves test expressions
sppt= spp.subs({'d':1})

import numpy as np
def cartesian(arrays, out=None):
    Generate a cartesian product of input arrays.

    arrays : list of array-like
        1-D arrays to form the cartesian product of.
    out : ndarray
        Array to place the cartesian product in.

    out : ndarray
        2-D array of shape (M, len(arrays)) containing cartesian products
        formed of input arrays.

    >>> 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 =[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