(define-param eps1 12)
(define-param eps2 1)
(define-param deps2 0) ; the change in one eps2 layer (default is 0)
; use quarter-wave thickness by default
(define-param d1 (/ (sqrt eps2) (+ (sqrt eps1) (sqrt eps2))))
(define d2 (- 1 d1))
(set! default-material (make dielectric (epsilon eps1)))
; this is just one block of eps2 at the center... we duplicate it below
(set! geometry (list (make block
(center 0)
(size d2 infinity infinity)
(material (make dielectric (epsilon eps2))))))
; Before we compute the defect modes, let's compute the gap by using
; a size-1 unit cell.
(set! geometry-lattice (make lattice (size 1 no-size no-size)))
(set! k-points (list (vector3 0.5))) ; edge of B.Z. is where gap edges are
(set! num-bands 2) ; only need 2 bands
(run-tm)
(define gap-min (car freqs)) ; the first band's frequency
(define gap-max (cadr freqs)) ; the second band's frequency
(print "gap is from " gap-min " to " gap-max "\n")
(define-param N 15) ; the size of the supercell
(set! geometry-lattice (make lattice (size N no-size no-size)))
; to create the geometry, we have to repeat the epsilon 1 block N times,
(set! geometry (geometric-objects-lattice-duplicates geometry))
; Finally, we have to create the defect.
(set! geometry (append
geometry
(list (make block
(center 0)
(size d2 infinity infinity)
(material (make dielectric
(epsilon (+ eps2 deps2))))))))
(set-param! resolution 32) ; 32 pixels/a
; for the defect mode, k does matter (if the supercell is big enough),
; so we just set it to zero. If k-interp is not zero, however, we will
; indeed compute the (folded) band structure as in problem 4(a).
(define-param k-interp 0)
(if (zero? k-interp)
(set! k-points (list (vector3 0)))
(set! k-points (interpolate k-interp (list (vector3 0) (vector3 0.5)))))
; because of the folding, the first band (before the gap) will be folded
; N times. So, we need to compute N bands plus some extra bands to
; get whatever defect states lie in the gap. The number of extra bands
; is set by the parameter extra-bands.
(define-param extra-bands 1)
(set! num-bands (+ N extra-bands))
; Now, instead of outputting all of the modes, we'll define a little band
; function to just output the modes in the gap. A band function is
; a function of band index, b; when passed to run, this function is
; called for every band computed. We'll also print the corresponding
; frequency so that we can easily grep for it. (See also the MPB manual).
(define (output-in-gap b)
(let ((f (list-ref freqs (- b 1)))) ; the frequency of band b
(if (and (> f gap-min) (< f gap-max))
(begin
(print "gapfreq" b ":, " deps2 ", " f "\n")
(output-efield-z b)))))
(run-tm output-in-gap) ; run, outputting Ez
; Now that we're done, let's compute the fraction of the field energy
; in the defect layer for the last band, for use in getting the slope
; domega/d(deps2) for problem 4(d).
(get-dfield num-bands) ; get the last band's D
(compute-field-energy) ; compute its epsilon |E|^2 energy density
; compute the fraction of the energy in the defect object (center block).
(define energy-in-defect (compute-energy-in-objects
(make block
(center 0) (size d2 infinity infinity)
(material air))))
(print "predicted slope:, " deps2 ", "
(* -0.5
(list-ref freqs (- num-bands 1)) ; omega of last band
(/ (+ eps2 deps2))
energy-in-defect) "\n")
; Finally, for problem 4(f), let's compute the decay rate by taking the
; relative amplitude of the field at a random point and the random point
; + a. (Note that the phase should be -1, or a phase angle of 180 degrees).
(get-efield num-bands)
(define-param fx (+ (/ N 4) 0.182341)) ; a "random" pt. far from the defect
(define field-ratio (/ (vector3-z (get-field-point (vector3 fx)))
(vector3-z (get-field-point (vector3 (- fx 1))))))
(print "decay:, "
deps2 ", "
(list-ref freqs (- num-bands 1)) ", " ; omega of last band
(- (log (magnitude field-ratio))) ", " ; decay coefficient kappa
(magnitude field-ratio) ", "
(rad->deg (angle field-ratio))
"\n")