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