File: perturbation-theory-2d.ctl

package info (click to toggle)
meep-openmpi 1.25.0-2
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 64,556 kB
  • sloc: cpp: 32,214; python: 27,958; lisp: 1,225; makefile: 505; sh: 249; ansic: 131; javascript: 5
file content (142 lines) | stat: -rw-r--r-- 5,080 bytes parent folder | download | duplicates (8)
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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
(set-param! resolution 30)  ; 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 sxy (* 2 (+ r w pad dpml)))
(set! geometry-lattice (make lattice (size sxy sxy no-size)))

(set! symmetries (list (make mirror-sym (direction X) (phase (if perpendicular? +1 -1)))
                       (make mirror-sym (direction Y) (phase (if perpendicular? -1 +1)))))

(set! geometry (list
		(make cylinder
                  (material (make dielectric (index n)))
                  (radius (+ r w))
                  (height infinity)
                  (center 0 0))
		(make cylinder
                  (material vacuum)
                  (radius r)
                  (height infinity)
                  (center 0 0))))

;; 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)))
               (make source
                 (src (make gaussian-src (frequency fcen) (fwidth df)))
                 (component src-cmpt)
                 (center (- (+ r 0.1)))
                 (amplitude -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

(set! pml-layers (list (make pml (thickness dpml))))

(set! geometry-lattice (make lattice (size sxy sxy no-size)))

(set! geometry (list
		(make cylinder
                  (material (make dielectric (index n)))
                  (radius (+ r w))
                  (height infinity)
                  (center 0 0))
		(make cylinder
                  (material vacuum)
                  (radius r)
                  (height infinity)
                  (center 0 0))))

(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)))
               (make source
                 (src (make gaussian-src (frequency fcen) (fwidth df)))
                 (component src-cmpt)
                 (center (- (+ r 0.1)))
                 (amplitude -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 Ey (vector3 r))))) (* (+ r w) (sqr (magnitude (get-field-point Ey (vector3 (+ r w)))))))))
          (perp-integral (* deps-inv 2 pi (- (* (+ r w) (sqr (magnitude (get-field-point Dy (vector3 0 (+ r w)))))) (* r (sqr (magnitude (get-field-point Dy (vector3 0 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 0) (size (- sxy (* 2 dpml)) (- sxy (* 2 dpml))))))
(define perturb-theory-dw-dR (/ (* -1 frq-unperturbed numerator-integral) (* 8 denominator-integral)))

(reset-meep)

;; perturbed geometry with narrowband source

(define-param dr 0.04)

(set! pml-layers (list (make pml (thickness dpml))))

(set! geometry-lattice (make lattice (size sxy sxy no-size)))

(set! geometry (list
		(make cylinder
                  (material (make dielectric (index n)))
                  (radius (+ r dr w))
                  (height infinity)
                  (center 0 0))
		(make cylinder
                  (material vacuum)
                  (radius (+ r dr))
                  (height infinity)
                  (center 0 0))))

(set! sources (list
               (make source
                 (src (make gaussian-src (frequency fcen) (fwidth df)))
                 (component src-cmpt)
                 (center (+ r dr 0.1)))
               (make source
                 (src (make gaussian-src (frequency fcen) (fwidth df)))
                 (component src-cmpt)
                 (center (- (+ r dr 0.1)))
                 (amplitude -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")