Update for ver. 1.2 of MEEP. Apparently fixed.
1 Meep 1.2 (20 July 2012) 2 3 * Fixed to work with Guile version 2.x (older versions still work); 4 requires libctl 3.2 or later.
1 Meep 1.2 (20 July 2012) 2 3 * Fixed to work with Guile version 2.x (older versions still work); 4 requires libctl 3.2 or later.
(set! geometry-lattice (make lattice (size 20 no-size no-size)))
(set! pml-layers (list (make pml (thickness 1))))
(set-param! resolution 30)
(define myAg (make dielectric (epsilon 1)
(polarizations
(make polarizability
(omega 1e-20) (gamma 0.0038715) (sigma 4.4625e+39))
; (make polarizability
; (omega 0.065815) (gamma 0.31343) (sigma 7.9247))
; (make polarizability
; (omega 0.36142) (gamma 0.036456) (sigma 0.50133))
; (make polarizability
; (omega 0.66017) (gamma 0.0052426) (sigma 0.013329))
; (make polarizability
; (omega 0.73259) (gamma 0.07388) (sigma 0.82655))
; (make polarizability
; (omega 1.6365) (gamma 0.19511) (sigma 1.1133))
)))
(set! geometry (list
(make block (center -3 0 0 0.0) (size 5 infinity infinity)
(material myAg ))
(make block (center 3 0 0 0.0) (size 5 infinity infinity)
(material myAg ))
))
(set! sources (list
(make source
(src (make gaussian-src (frequency 2.5) (fwidth 5)))
(component Ex)
(center 0 0 0))
(make source
(src (make gaussian-src (frequency 2.5) (fwidth 5)))
(component Ex)
(center -.5 0 0))
(make source
(src (make gaussian-src (frequency 2.5) (fwidth 5)))
(component Ex)
(center .5 0 0))
))
(define-param kmin 0.0)
(define-param kmax 4.0)
(define-param k-interp 50)
(define kpts (interpolate k-interp (list (vector3 0 kmin 0) (vector3 0 kmax 0))))
(define all-freqs (run-k-points 200 kpts)) ; a list of lists of frequencies
And here's a python snippet that gives a list of k-vectors their frequencies by giving it the MEEP output as a file object
def getfreqs(meepoutfileobj):
fd=[]
for aline in meepoutfileobj:
if 'freqs:' in aline:
ld=aline.split(',')
kpt=ld[2:5];freqs=ld[5:]
pts=[( [float(an) for an in kpt] ,float(afreq)) for afreq in freqs]
fd.extend(pts)
if fd==[]:return None
return fd
#Author: Majid al-Dosari
#code that takes in lorentz-drude model parameters and converts them to code
#usable by meep. just use eVLD2meepcode
#example: eVLD2meepcode('mysilver',100e-9,9,[[.845,.048,0],[.065,3.886,.816]])
#output:
#(define mysilver (make dielectric (epsilon 1)
#(polarizations
#(make polarizability
#( omega 8.06554724676e-22 ) ( gamma 0.00387146267844 ) ( sigma 6.8445e+41 ))
#( omega 0.0658148655336 ) ( gamma 0.313427166009 ) ( sigma 7.9071150519 ))
#)))
#it's an implementation of
#http://juluribk.com/2011/04/27/plasmonic-materials-in-meep/ by Bala Krishna Juluri
#normalization length a
#w omega
#p plasma
#G Gamma
def sip(stuff):#stuff in parens
return '( ' +str(stuff)+ ' )'
def wordandvalue(word,value):
return sip(word+' '+str(value))
def meepfactor(a):
return 2*3.14159265359*299792458/a #2pi*c/a
def eV2w(eV):return eV*(2*3.14159265359)/4.135666e-15
def eV2meep(a,eV):return eV2w(eV)/meepfactor(a)
def sigmaa(f,wp,w):
"""f is oscillator strength, wp is plamsa freq"""
return (f*wp**2)/w**2
def eVLD2meepL(a,wp,fGw0_list):
oGs_list=[]
smallno=1e-20;
for f,G,w0 in fGw0_list:
if w0==0:w0=smallno #hack
oGs_list.append([eV2meep(a,w0)
,eV2meep(a,G)
,sigmaa(f,eV2meep(a,wp),eV2meep(a,w0))])
return oGs_list
def eVLD2meepcode(name,a,wp,fGw0_list):
"""
input:
name: for code
a: meep normalization const
wp: plasma freq in eV
fGw0_list: nested list of
f oscillator strengths
Gamma relaxation in eV
omega0 in eV,w=0 is handled by converting it to a very small no.
"""
return writelorentz(name,oGs_list=eVLD2meepL(a,wp,fGw0_list))
def writelorentz(name,oGs_list=None,epsilon=1):#meep is in lorentz
txt=''
txt+='(define '+str(name)+' (make dielectric (epsilon '+str(epsilon)+') \n'
if oGs_list==None or len(oGs_list)==0:
txt+='))'
return txt
txt+='(polarizations \n'
for omega,Gamma,sigma in oGs_list:
txt+='(make polarizability \n'
txt+=wordandvalue('omega',omega)+' '
txt+=wordandvalue('gamma',Gamma)+' '
txt+=wordandvalue('sigma',sigma)+')\n'
txt+=')))'
return txt
#so 1. normalize then 2. cast into meep