-
Notifications
You must be signed in to change notification settings - Fork 592
/
perturbation-theory.ctl
104 lines (71 loc) · 3.77 KB
/
perturbation-theory.ctl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
(set-param! resolution 100) ; pixels/um
(define-param perpendicular? true)
(define src-cmpt (if perpendicular? Hz Ez))
(define fcen (if perpendicular? 0.21 0.17)) ; pulse center frequency
(define-param n 3.4) ; index of waveguide
(define-param w 1) ; width of waveguide
(define-param r 1) ; inner radius of ring
(define-param pad 4) ; padding between waveguide and edge of PML
(define-param dpml 2) ; thickness of PML
(set-param! m 5) ; angular dependence
(set! pml-layers (list (make pml (thickness dpml))))
(define sr (+ r w pad dpml)) ; radial size (cell is from 0 to sr)
(set! dimensions CYLINDRICAL)
(set! geometry-lattice (make lattice (size sr no-size no-size)))
(set! geometry (list (make block
(center (+ r (/ w 2)))
(size w infinity infinity)
(material (make dielectric (index n))))))
;; find resonant frequency of unperturbed geometry using broadband source
(define df (* 0.2 fcen)) ; pulse width (in frequency)
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component src-cmpt)
(center (+ r 0.1)))))
(run-sources+ 100 (after-sources (harminv src-cmpt (vector3 (+ r 0.1)) fcen df)))
(define frq-unperturbed (harminv-freq-re (car harminv-results)))
(reset-meep)
;; unperturbed geometry with narrowband source centered at resonant frequency
(set! pml-layers (list (make pml (thickness dpml))))
(set! geometry-lattice (make lattice (size sr no-size no-size)))
(set! geometry (list (make block
(center (+ r (/ w 2)))
(size w infinity infinity)
(material (make dielectric (index n))))))
(set! fcen frq-unperturbed)
(set! df (* 0.05 fcen))
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component src-cmpt)
(center (+ r 0.1)))))
(run-sources+ 100)
(define deps (- 1 (* n n)))
(define deps-inv (- 1 (/ (* n n))))
(define numerator-integral 0)
(if perpendicular?
(let ((para-integral (* deps 2 pi (- (* r (sqr (magnitude (get-field-point Ep (vector3 r))))) (* (+ r w) (sqr (magnitude (get-field-point Ep (vector3 (+ r w)))))))))
(perp-integral (* deps-inv 2 pi (- (* (+ r w) (sqr (magnitude (get-field-point Dr (vector3 (+ r w)))))) (* r (sqr (magnitude (get-field-point Dr (vector3 r)))))))))
(set! numerator-integral (+ para-integral perp-integral)))
(set! numerator-integral (* deps 2 pi (- (* r (sqr (magnitude (get-field-point Ez (vector3 r))))) (* (+ r w) (sqr (magnitude (get-field-point Ez (vector3 (+ r w))))))))))
(define denominator-integral (electric-energy-in-box (volume (center (* 0.5 sr)) (size sr))))
(define perturb-theory-dw-dR (/ (* -1 frq-unperturbed numerator-integral) (* 4 denominator-integral)))
(reset-meep)
;; perturbed geometry with narrowband source
(define-param dr 0.01)
(set! pml-layers (list (make pml (thickness dpml))))
(set! geometry-lattice (make lattice (size sr no-size no-size)))
(set! geometry (list (make block
(center (+ r dr (/ w 2)))
(size w infinity infinity)
(material (make dielectric (index n))))))
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component src-cmpt)
(center (+ r dr 0.1)))))
(run-sources+ 100 (after-sources (harminv src-cmpt (vector3 (+ r 0.1)) fcen df)))
(define frq-perturbed (harminv-freq-re (car harminv-results)))
(define finite-diff-dw-dR (/ (- frq-perturbed frq-unperturbed) dr))
(print "dwdR:, " perturb-theory-dw-dR " (pert. theory), " finite-diff-dw-dR " (finite diff.)\n")