; Compute the transmission from a normal-incident planewave through
; a 2d photonic crystal composed of a square lattice of dielectric
; rods in air.
;
; The cell has nx rods in the x direction, with absorbing boundary
; conditions.  In the y direction, we have periodic boundary
; conditions.  The planewave is incident in the x direction, and is
; separated from the crystal by a distance "pad".  We compute the
; transmitted power through a flux plane on the other side of the
; crystal.

; some parameters describing the structure
(define-param eps 8.9) ; dielectric constant of the rods
(define-param r 0.2) ; rod radius
(define-param nx 5) ; # rods in x direction
(define-param pad 4) ; distance between rods and source/flux plane
(define-param dpml 1) ; thickness of PML absorbing boundary layer

(define sx (+ nx (* 2 (+ pad dpml)))) ; size of cell in x direction
(set! geometry-lattice (make lattice (size sx 1 no-size)))

(set! k-point (vector3 0 0 0)) ; periodic boundaries
(set! pml-layers (list (make pml (direction X) (thickness dpml))))

(set! geometry
      (geometric-object-duplicates
       (vector3 1 0) (/ (- 1 nx) 2) (/ (- nx 1) 2) ; nx copies in x
       (make cylinder (center 0) (radius r) (height infinity)
	     (material (make dielectric (epsilon eps))))))

(set-param! resolution 20) ; # pixels/a

(define-param fcen 0.375) ; pulse center frequency
(define-param df 0.35)    ; pulse width (in frequency)
(set! sources (list
	       (make source
		 (src (make gaussian-src (frequency fcen) (fwidth df)))
		 (component Ez)
		 (center (+ dpml (* -0.5 sx)) 0)
		 (size 0 1))))

(set! symmetries (list (make mirror-sym (direction Y)))) ; exploit y=0 mirror

(define-param nfreq 500) ; number of frequencies at which to compute flux
(define trans ; transmitted flux                                          
  (add-flux fcen df nfreq
	    (make flux-region
	      (center (- (* 0.5 sx) dpml 0.5) 0) (size 0 1))))

; Run the simulation (stopping when Ez at flux plane grows small)
(run-sources+ (stop-when-fields-decayed
	       50 Ez
	       (vector3 (- (* 0.5 sx) dpml 0.5) 0)
	       1e-3)
	      (at-beginning output-epsilon))

(display-fluxes trans) ; print out the flux spectrum (power vs. frequency)