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