not sure the ctl file got through (i couldn't access it through
gname), so here it is:
; compute R and T for square lattice of spheres for infinite plane
wave(s).
;size of computation x = [ -sx/2, sx/2]
; y = [ -sy/2, sy/2]
; z = [ -sz/2-pmlt, sz/2+pmlt] - direction of
propagation
(print "starting control file\n")
(define a 1) ; lattice constant
(define sx a)
(define sy a)
; size of simulation
(define-param sz 15) ; guess z direction length
(define-param auto_z? true) ; scale z according to longest
wavelength
(define-param tol_e 0.01) ; run until energy ratio decays to
this value
(define-param max_steps 500000) ; or until we do this many timesteps
(define-param rad (/ a 2)) ; radius of spheres
(define-param no-slab? false) ;is there a slab of spheres?
(define-param pc_metal? false) ; is photonic crystal metal spheres?
(define-param pc_eps (* 1.45 1.45)) ; default epsilon for spheres
; source
(define-param f_min 0.3) ;minimum frequency
(define-param f_max 1.2) ;maximum frequency
(define fcen (* 0.5 (+ f_min f_max)))
(define df (* 2 (abs (- f_max f_min))))
(define-param nfreq 601) ;number of frequency sample points
; misc
(define-param pmlt 1) ;thickness of PML
(define-param res 64) ;resolution in unit length, a
(if pc_metal?
(define pc_mat (make dielectric (epsilon 1) ; Drude
approximation for lossy metal
(polarizations (make polarizability (omega 1e-20) (gamma
0.0313167) (delta-epsilon 9.704152e40)))))
(define pc_mat (make dielectric (epsilon pc_eps))) ; otherwise
simple dielectric
)
(if auto_z? ; set z direction so about 3 of the longest wavelengths
are present before and after slab?
(if (< f_min 0.1)
(set-param! sz (/ 10 fcen))
(set-param! sz (/ 6 f_min))
)
)
(define z_tot (+ sz a (* 2 pmlt))) ; size of cell in propagation
direction = 2 * PML + sz + sphere 'layer'
; record what we are doing
(print "res=" res ", pmlt=" pmlt ", z_tot=" z_tot "\n")
(print "sx=" sx ", sy=" sy ", sz=" sz ", dg=" (/ a res) ", radius="
rad "\n")
(print "no-slab?=" no-slab? "pc_metal?=" pc_metal? ", pc_eps=" pc_eps
"\n")
(set! geometry-lattice ;make computational cell
(make lattice (size sx sy z_tot )))
(set-param! k-point (vector3 0 0 0) ) ; perioidic
(if (not no-slab?) ;make the slab or not?
(set! geometry
(list (make sphere (center 0 0 0) (radius rad) (material pc_mat)))
)
)
(print "setting PML\n")
(set! pml-layers ;PML at negative and positive Z bounds
(list
(make pml (thickness pmlt) (direction Z) (side Low))
(make pml (thickness pmlt) (direction Z) (side High))
)
)
(print "setting resolution of " res "\n")
(set-param! resolution res) ;resolution per unit length
(print "setting source\n")
(set! sources (list ; gaussian at Low Z, plane wave towards High Z
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ex)
(size sx sx 0)
(center 0 0 (- (+ 0.5 pmlt) (/ z_tot 2)))
))
)
(print "defining flux planes\n")
; define a bunch of flux planes
(define xa (/ (- (/ sz 2) pmlt 2.5) 3)) ; distances between planes
(define xb (+ pmlt 1.0)) ; starting location of planes
;transmission
(define trans1 (add-flux fcen (- f_max f_min) nfreq (make flux-region
(center 0 0 (- (/ z_tot 2) xb (* xa 3.0))) (size sx sy 0) ))) ; near
to PC
(define trans2 (add-flux fcen (- f_max f_min) nfreq (make flux-region
(center 0 0 (- (/ z_tot 2) xb (* xa 2.0))) (size sx sy 0) )))
(define trans3 (add-flux fcen (- f_max f_min) nfreq (make flux-region
(center 0 0 (- (/ z_tot 2) xb (* xa 1.0))) (size sx sy 0) )))
(define trans4 (add-flux fcen (- f_max f_min) nfreq (make flux-region
(center 0 0 (- (/ z_tot 2) xb (* xa 0.0))) (size sx sy 0) ))) ; far
from PC
(print "t1=" (- (/ z_tot 2) xb (* xa 3.0)) ", t2=" (- (/ z_tot 2) xb
(* xa 2.0)) ", t3=" (- (/ z_tot 2) xb (* xa 1.0)) ", t4=" (- (/ z_tot
2) xb (* xa 0.0)) "\n")
;reflection
(define refl1 (add-flux fcen (- f_max f_min) nfreq (make flux-region
(center 0 0 (- (+ xb (* xa 3.0)) (/ z_tot 2))) (size sx sy 0) ))) ;
near to PC
(define refl2 (add-flux fcen (- f_max f_min) nfreq (make flux-region
(center 0 0 (- (+ xb (* xa 2.0)) (/ z_tot 2))) (size sx sy 0) )))
(define refl3 (add-flux fcen (- f_max f_min) nfreq (make flux-region
(center 0 0 (- (+ xb (* xa 1.0)) (/ z_tot 2))) (size sx sy 0) )))
(define refl4 (add-flux fcen (- f_max f_min) nfreq (make flux-region
(center 0 0 (- (+ xb (* xa 0.0)) (/ z_tot 2))) (size sx sy 0) ))) ;
far from PC
(print "r1=" (- (+ xb (* xa 3.0)) (/ z_tot 2)) ", r2=" (- (+ xb (* xa
2.0)) (/ z_tot 2)) ", r3=" (- (+ xb (* xa 1.0)) (/ z_tot 2)) ",
r4=" (- (+ xb (* xa 0.0)) (/ z_tot 2)) "\n")
(print "commencing simulation...\n")
; if we are doing the slab, read in the results from not doing the slab
(if (not no-slab?)
(begin
(print "loading old flux...") (newline)
(load-minus-flux "refl1-flux" refl1) ;near to PC
(load-minus-flux "refl2-flux" refl2)
(load-minus-flux "refl3-flux" refl3)
(load-minus-flux "refl4-flux" refl4) ;far from PC
(print "done!\n") (newline)
)
)
; do the simulation
(define tim (- (/ z_tot 3) pmlt)) ;time to travel third of the z
direction
(print "time to sample energy: " tim "\n")
(print "run until energy ratio is " tol_e " or # of time steps > "
max_steps "\n")
(define tot_energy 0)
(define max_energy -1)
(define ratio 1)
(if (not no-slab?) (output-epsilon))
(run-sources+ (lambda () (or (< ratio tol_e) (> (meep-fields-t-get
fields) max_steps)))
(after-sources (lambda ()
(if (< max_energy 0)
(begin
(set! max_energy (meep-fields-field-energy-in-box fields
(volume (center 0 0 0) (size sx sy z_tot))))
(print "Energy after sources off = " max_energy "\n")
)
))
)
(at-every tim (lambda ()
(set! tot_energy (meep-fields-field-energy-in-box fields
(volume (center 0 0 0) (size sx sy z_tot)) ))
(set! max_energy (max max_energy tot_energy))
(set! ratio (/ tot_energy max_energy))
(print "energy decay (t=" (meep-time) "s, " (meep-fields-t-get
fields) " timestep) = " tot_energy " / " max_energy " = " ratio "\n"))
)
)
;if we are not doing a slab, then save the results so we can find R &
T later...
(if no-slab?
(begin
(save-flux "refl1-flux" refl1)
(save-flux "refl2-flux" refl2)
(save-flux "refl3-flux" refl3)
(save-flux "refl4-flux" refl4)
)
)
; display the fluxes- start close to PC and move out for both
transmission and reflectance
(display-fluxes trans1 trans2 trans3 trans4 refl1 refl2 refl3 refl4)
(print "done with control file\n")
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
