;;;;;;;;;;;;;;;;;;;;;;;;
;; defect_slab_bic.ctl
;; Alex Cerjan
;; 7.13.18
;;
;; Can you have a BIC persist when the phc slab is cladded
;; in a photonic crystal?
;;;;;;;;;;;;;;;;;;;;;;;;

;;For organizational purposes:
(define-param outputEpsBool true)
(if (not outputEpsBool)
    (set! output-epsilon (lambda () (print "skipping output-epsilon\n"))))

(define-param outputFilenamePrefix "test")
(set-param! filename-prefix (string-append outputFilenamePrefix "-"))

(define-param outDirInput "/home/alex/Data/bics_in_phcs/two_sided_system/test/")
(use-output-directory outDirInput)

;; Parameters for defect slab:
(define-param epsDefect 7) ; dielectric constant of waveguide
(define-param hDefect 0.5) ; width of waveguide
(define-param rDefect 0.22) ; radius of holes

;; Parameters for PhC:
(define-param epsPhC 12)
(define-param wUpPhC 0.18)
(define-param wDnPhC 0.18)

;; Frequency parameters:
(define-param fcen 0.359) ; pulse center frequency                            
(define-param df .05) ; pulse freq. width: large df = short impulse

;; The cell dimensions:
(define-param phcZ 11) ; # of phc unit cells in z direction (perpendicular to slab.)
(define-param dpml 6) ; PML thickness (z direction only!)

;; Simulation resolution + tolerances:
(set-param! resolution 16)
(define-param totalTime 500)

(define-param cw-solver-tol 1e-6)         ; CW solver tolerance
(define-param cw-solver-max-iters 2000)   ; CW solver max iterations
(define-param cw-solver-L 10)             ; CW solver L

;; Wavevector parameters:
(define-param single-kx 0.35)
(define-param single-ky 0.35)

;; Shift vector parameters: (offset the center of one environment from the other)
(define-param uppX 0)
(define-param lowX 0)
(define-param uppY 0)
(define-param lowY 0)

;; use symmetry?
(define-param sym-setting 0) ; -1 = TM, +1 = TE, 0 = don't use symmetry.

;;;;;;;;;;;;;;;;;;;;;;
;; Run it:
;;;;;;;;;;;;;;;;;;;;;;

;; Set system:
(define geo-latt (make lattice (size 1 1 (+ phcZ phcZ hDefect))))

;; Add sources:
(define s1 (vector3 0.2436 0.0894 0.2))
(define srcComp Ex)
(define gauss (make gaussian-src (frequency fcen) (fwidth df)))

;; Add PML:
(define pml-lay (list (make pml (direction Z) (thickness dpml))))

;; Corners of the irreducible Brillouin zone
(define Gamma (vector3 0 0 0))
(define X (vector3 0.5 0 0))
(define Y (vector3 0 0.5 0))
(define M (vector3 0.5 0.5 0))

(define kSingleChoice (vector3 single-kx single-ky 0))

;; Set k-point:
(define kList (list kSingleChoice))

;; Define sources:
(define srcList (list (make source (src gauss) (component srcComp) (center s1))))
		   
;; Add geometry:
(define wDiff (* 2 (- wDnPhC wUpPhC)))
(define geo-list
  (append
   ;; defect layer:
   (list (make block (center 0 0 0) (size 1 1 hDefect)
	       (material (make dielectric (epsilon epsDefect))))
	 (make cylinder (center 0 0 0) (radius rDefect) (height hDefect) (axis (vector3 0 0 1)) (material air)))
   ;; vertical pillars above slab:
   (list (shift-geometric-object
	  (make block (center 0.5 0.5 (+ (/ hDefect 2) (/ phcZ 2)))
		(size (vector3 (* 2 wUpPhC) (* 2 wUpPhC) phcZ))
		(e1 (vector3 1 0 0)) (e2 (vector3 0 1 0)) (e3 (vector3 0 0 1))
		(material (make dielectric (epsilon epsPhC))))
	  (vector3 uppX uppY 0)))
   (list (shift-geometric-object
	  (make block (center 0.5 -0.5 (+ (/ hDefect 2) (/ phcZ 2)))
		(size (vector3 (* 2 wUpPhC) (* 2 wUpPhC) phcZ))
		(e1 (vector3 1 0 0)) (e2 (vector3 0 1 0)) (e3 (vector3 0 0 1))
		(material (make dielectric (epsilon epsPhC))))
	  (vector3 uppX uppY 0)))
   (list (shift-geometric-object
	  (make block (center -0.5 0.5 (+ (/ hDefect 2) (/ phcZ 2)))
		(size (vector3 (* 2 wUpPhC) (* 2 wUpPhC) phcZ))
		(e1 (vector3 1 0 0)) (e2 (vector3 0 1 0)) (e3 (vector3 0 0 1))
		(material (make dielectric (epsilon epsPhC))))
	  (vector3 uppX uppY 0)))
   (list (shift-geometric-object
	  (make block (center -0.5 -0.5 (+ (/ hDefect 2) (/ phcZ 2)))
		(size (vector3 (* 2 wUpPhC) (* 2 wUpPhC) phcZ))
		(e1 (vector3 1 0 0)) (e2 (vector3 0 1 0)) (e3 (vector3 0 0 1))
		(material (make dielectric (epsilon epsPhC))))
	  (vector3 uppX uppY 0)))
   ;; vertical pillars below slab:
   (list (shift-geometric-object
	  (make block (center 0.5 0.5 (+ (/ hDefect -2) (/ phcZ -2)))
		(size (vector3 (+ (* 2 wUpPhC) wDiff) (- (* 2 wUpPhC) wDiff) phcZ))
		(e1 (vector3 1 0 0)) (e2 (vector3 0 1 0)) (e3 (vector3 0 0 1))
		(material (make dielectric (epsilon epsPhC))))
	  (vector3 lowX lowY 0)))
   (list (shift-geometric-object
	  (make block (center 0.5 -0.5 (+ (/ hDefect -2) (/ phcZ -2)))
		(size (vector3 (+ (* 2 wUpPhC) wDiff) (- (* 2 wUpPhC) wDiff) phcZ))
		(e1 (vector3 1 0 0)) (e2 (vector3 0 1 0)) (e3 (vector3 0 0 1))
		(material (make dielectric (epsilon epsPhC))))
	  (vector3 lowX lowY 0)))
   (list (shift-geometric-object
	  (make block (center -0.5 0.5 (+ (/ hDefect -2) (/ phcZ -2)))
		(size (vector3 (+ (* 2 wUpPhC) wDiff) (- (* 2 wUpPhC) wDiff) phcZ))
		(e1 (vector3 1 0 0)) (e2 (vector3 0 1 0)) (e3 (vector3 0 0 1))
		(material (make dielectric (epsilon epsPhC))))
	  (vector3 lowX lowY 0)))
   (list (shift-geometric-object
	  (make block (center -0.5 -0.5 (+ (/ hDefect -2) (/ phcZ -2)))
		(size (vector3 (+ (* 2 wUpPhC) wDiff) (- (* 2 wUpPhC) wDiff) phcZ))
		(e1 (vector3 1 0 0)) (e2 (vector3 0 1 0)) (e3 (vector3 0 0 1))
		(material (make dielectric (epsilon epsPhC))))
	  (vector3 lowX lowY 0)))
   ))

;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Really run it:
;;;;;;;;;;;;;;;;;;;;;;;;;;

;; setup:
(define kIn (car kList))
(set! geometry-lattice geo-latt)
(set! pml-layers pml-lay)
(set! geometry geo-list)
(set! sources srcList)

(define (run-k-points-test T k-points)
  (define k-index 0)
  (define all-freqs '())
  (map (lambda (k)
	 (set! k-index (+ k-index 1))
	 (if (= k-index 1) (begin (init-fields) (output-epsilon)))
	 (let ((freqs (run-k-point T k)))
	   ;; (print "freqs:, " k-index
	   ;; 	  ", " (vector3-x k) ", " (vector3-y k) ", " (vector3-z k))
	   ;; (map (lambda (x) (print ", " x)) (map real-part freqs))
	   ;; (print "\n")
	   ;; (print "freqs-im:, " k-index
	   ;; 	  ", " (vector3-x k) ", " (vector3-y k) ", " (vector3-z k))
	   ;; (map (lambda (x) (print ", " x)) (map imag-part freqs))
	   ;; (print "\n")
	   (set! all-freqs (cons freqs all-freqs))
	   ))
       k-points)
  ;; (reverse all-freqs)
  )


(run-k-points-test totalTime kList)


;; (run-k-points totalTime kList)
;; (run-k-point totalTime (car kList))