Hi everybody,
I'm fairly new to meep and have a problem with the simulation of an opal
structure.
By searching in google I found a discussion from 2006 about a code for
multilayer opals in which a code was posted that worked for a simple
opal. Still I get transmissions larger that one. Anyone here with an
idea why that happens?
Here is the code:
--------------------------------------------------------------------
; transmittance/reflectance of a fcc structure a=1.
(set! eps-averaging? false); No epsilon averaging
;-------------------EXTERNAL PARAMETERS-----------------------------------
(define-param r_sph 0.36); Oxide spheres radius
(define-param r_max (/ 1 (sqrt 6))); Thickness of the Si layer grown
inside the opal
(define-param nl 3) ; number of periods = (3 layers each)
(define-param pmlt 2) ; thickness of the pml
(define-param pad 3); distance between source and sample
(define-param no-phc 1) ; if true (1), have air, false (0), PhC
;--------------CELL DIMENSIONS-----------------------------------------
(define unitcell (sqrt 3)) ; size of the unitcell in Z
(define distz (/ 1 (sqrt 3))); distance between single layers
(define sx (/ 1 (sqrt 2))) ; size of cell in X direction
(define sy (sqrt 1.5)) ; size of cell in Y direction
(define sz (+ (* 2 pad) (* 2 pmlt) (round (+ 0.5 (* nl unitcell))))) ;
size of cell in Z direction
(set! geometry-lattice (make lattice (size sx sy sz)))
;------------------USEFUL VALUES----------------------------------------
(define sx_half (/ sx 2)); (1/sqrt(8))
(define sy_half (/ sy 2)); (sqrt(3)/sqrt(8))
(define sy_sixth (/ sy 6)); (1/sqrt(3x8))
(define nsx_half (/ sx -2)); -(1/sqrt(8))
(define nsy_half (/ sy -2)); -(sqrt(3)/sqrt(8))
(define nsy_sixth (/ sy -6)); -(1/sqrt(3x8))
(define cell_center (* unitcell (+ 0.5 (* -0.5 nl)))); z-coordinate of
the first unit cell
(define bottom_sample (+ (* unitcell (- nl 1)) cell_center distz (/ 1
(sqrt 8)))) ; Z coordinate bottom of spheres
(define substrate_centerZ (+ (/ sz 4) (/ bottom_sample 2)))
(define substrate_sizeZ (- (/ sz 2) bottom_sample))
;----------------MATERIALS-----------------------------------------------
(define silica (make dielectric (epsilon (* 1.45 1.45))))
(define silicon (make dielectric (epsilon (* 3.5 3.5))))
;--------------GEOMETRY-------------------------------------------------
(set! geometry
(if (= no-phc 1)
(list
(make sphere (center 0 0 0) (radius r_sph) (material air))
(make block (center 0 0 substrate_centerZ) (material silicon)
(size sx sy substrate_sizeZ)) ; substrate
)
(append
;OPAL OF Si SPHERES
(geometric-objects-duplicates (vector3 0 0 unitcell) 0 (- nl 1)
(list
;Layer A - thin Si
(make sphere (center nsx_half nsy_sixth (- cell_center distz))
(radius r_max)
(material silicon))
(make sphere (center sx_half nsy_sixth (- cell_center distz))
(radius r_max)
(material silicon))
(make sphere (center 0 (+ nsy_sixth nsy_half) (- cell_center
distz)) (radius r_max)
(material silicon))
(make sphere (center 0 (* 2 sy_sixth) (- cell_center distz))
(radius r_max)
(material silicon))
;Layer B - thin Si
(make sphere (center nsx_half sy_half cell_center)
(radius r_max)
(material silicon))
(make sphere (center nsx_half nsy_half cell_center)
(radius r_max)
(material silicon))
(make sphere (center sx_half sy_half cell_center) (radius r_max)
(material silica))
(make sphere (center sx_half nsy_half cell_center)
(radius r_max)
(material silicon))
(make sphere (center 0 0 cell_center) (radius r_max)
(material silicon))
;Layer C - thin Si
(make sphere (center nsx_half sy_sixth (+ cell_center distz))
(radius r_max)
(material silicon))
(make sphere (center sx_half sy_sixth (+ cell_center distz))
(radius r_max)
(material silicon))
(make sphere (center 0 (+ sy_sixth sy_half) (+ cell_center
distz)) (radius r_max)
(material silicon))
(make sphere (center 0 (* 2 nsy_sixth) (+ cell_center distz))
(radius r_max)
(material silicon))
))
;OPAL OF OXIDE SPHERES
(geometric-objects-duplicates (vector3 0 0 unitcell) 0 (- nl 1)
(list
;Layer A - oxide
(make sphere (center nsx_half nsy_sixth (- cell_center distz))
(radius r_sph)
(material silica))
(make sphere (center sx_half nsy_sixth (- cell_center distz))
(radius r_sph)
(material silica))
(make sphere (center 0 (+ nsy_sixth nsy_half) (- cell_center
distz)) (radius r_sph)
(material silica))
(make sphere (center 0 (* 2 sy_sixth) (- cell_center distz))
(radius r_sph)
(material silica))
;Layer B - oxide
(make sphere (center nsx_half sy_half cell_center)
(radius r_sph)
(material silica))
(make sphere (center nsx_half nsy_half cell_center)
(radius r_sph)(material silica))
(make sphere (center sx_half sy_half cell_center) (radius r_sph)
(material silica))
(make sphere (center sx_half nsy_half cell_center)
(radius r_sph)(material silica))
(make sphere (center 0 0 cell_center) (radius r_sph)
(material silica))
;Layer C - oxide
(make sphere (center nsx_half sy_sixth (+ cell_center distz))
(radius r_sph)
(material silica))
(make sphere (center sx_half sy_sixth (+ cell_center distz))
(radius r_sph)
(material silica))
(make sphere (center 0 (+ sy_sixth sy_half) (+ cell_center
distz)) (radius r_sph)
(material silica))
(make sphere (center 0 (* 2 nsy_sixth) (+ cell_center distz))
(radius r_sph)
(material silica))
))
;SUBSTRATE
(list (make block (center 0 0 substrate_centerZ) (material silicon)
(size sx sy substrate_sizeZ))) ; substrate
)))
;--------------SOURCES-------------------------------------------------
(define-param fcen 0.4) ; pulse center frequency
(define-param df 0.4) ; pulse width (in frequency)
(set! sources (list
(make source
(src (make gaussian-src (frequency (+ fcen 0.1)) (fwidth df)))
(component Ex)
(center 0 0 (+ pmlt (* -0.5 sz)))
(size sx sy 0))
(make source
(src (make gaussian-src (frequency (- fcen 0.1)) (fwidth df)))
(component Ex)
(center 0 0 (+ pmlt (* -0.5 sz)))
(size sx sy 0))
))
;-----------RESOLUTION, PERIODICAL BOUNDARIES, PML-----------------
(set-param! k-point (vector3 0))
(set! pml-layers (list (make pml (direction Z) (thickness pmlt))))
(set-param! resolution 20)
;--------------FLUX----------------------------------------------------
(define-param nfreq 600) ; number of frequencies at which to compute flux
(define trans ; transmitted flux
(add-flux fcen df nfreq
(make flux-region
(center 0 0 (- (/ sz 2) 0.5 pmlt)) (size sx sy 0))))
(define refl ; reflected flux
(add-flux fcen df nfreq
(make flux-region
(center 0 0 (+ 0.5 pmlt (* -0.5 sz))) (size sx sy 0))))
;-------------RUN-----------------------------------------------
; for normal run, load negated fields to subtract incident from refl. fields
(if (= no-phc 0) (load-minus-flux "refl-flux" refl))
(run-sources+
(stop-when-fields-decayed 50 Ex
(vector3 0 0 (- (/ sz 2) 0.5 pmlt))
1e-3)
(at-beginning output-epsilon))
; for normalization run, save flux fields for refl. plane
(if (= no-phc 1) (save-flux "refl-flux" refl))
(display-fluxes trans refl)
----------------------------------------------------------------------------------------------------------
the only thing I changed was to delete the silicon core leaving the
simulation with just silica particles that - according to the author -
should work.
Hoping to get some hints to where the problem is.
Regards
Marc Zöller
_______________________________________________
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
