File: polarization_grating.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 (125 lines) | stat: -rw-r--r-- 5,209 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
;; polarization grating from C. Oh and M.J. Escuti, Optics Letters, Vol. 33, No. 20, pp. 2287-9, 2008
;; note: reference uses z as the propagation direction and y as the out-of-plane direction; this script uses x and z, respectively

(set-param! resolution 50)    ; pixels/μm

(define-param dpml 1.0)       ; PML thickness
(define-param dsub 1.0)       ; substrate thickness
(define-param dpad 1.0)       ; padding thickness

(set! k-point (vector3 0 0 0))

(define boundary-layers (list (make pml (thickness dpml) (direction X))))
(set! pml-layers boundary-layers)

(define n0 1.55)
(define delta-n 0.159)
(define eps-diag (matrix3x3 (vector3 (sqr n0) 0 0)
                            (vector3 0 (sqr n0) 0)
                            (vector3 0 0 (sqr (+ n0 delta-n)))))

(define-param wvl 0.54)      ; center wavelength
(define fcen (/ wvl))        ; center frequency

(define-param d 1.7)         ; chiral layer thickness
(define-param ph 70)         ; chiral layer twist angle
(define-param gp 6.5)        ; grating period
(define-param nmode 5)       ; number of mode coefficients to compute

(set! ph (deg->rad ph))

(define sx (+ dpml dsub d d dpad dpml))
(define sy gp)

(define cell (make lattice (size sx sy no-size)))
(set! geometry-lattice cell)

; linear-polarized planewave pulse source
(define src-pt (vector3 (+ (* -0.5 sx) dpml (* 0.3 dsub)) 0 0))
(define lp-src (list (make source
                       (src (make gaussian-src (frequency fcen) (fwidth (* 0.05 fcen))))
                       (component Ez)
                       (center src-pt)
                       (size 0 sy 0))
                     (make source
                       (src (make gaussian-src (frequency fcen) (fwidth (* 0.05 fcen))))
                       (component Ey)
                       (center src-pt)
                       (size 0 sy 0))))

(set! sources lp-src)

(set! default-material (make medium (index n0)))

(define tran-pt (vector3 (- (* 0.5 sx) dpml (* 0.5 dpad)) 0 0))
(define tran-flux1 (add-flux fcen 0 1 (make flux-region (center tran-pt) (size 0 sy 0))))

(run-sources+ 100)

(define input-flux (get-fluxes tran-flux1))

(reset-meep)

(set! geometry-lattice cell)

(set! pml-layers boundary-layers)

(set! sources lp-src)

(set! k-point (vector3 0 0 0))

(set! default-material air)

; twist angle of nematic director; from equation 1b
(define phi (lambda (p)
              (let ((xx (- (vector3-x p) (+ (* -0.5 sx) dpml dsub))))
                (if (and (>= xx 0) (<= xx d))
                    (+ (* pi (vector3-y p) (/ gp)) (* ph xx (/ d)))
                    (+ (* pi (vector3-y p) (/ gp)) (- (* ph xx (/ d))) (* 2 ph))))))

(define lc-epsilon-diag (vector3 0 0 0))
(define lc-epsilon-offdiag (vector3 0 0 0))
(define lc-epsilon (matrix3x3 (vector3 0 0 0) (vector3 0 0 0) (vector3 0 0 0)))

; return the anisotropic permittivity tensor for a uniaxial, twisted nematic liquid crystal
(define lc-mat (lambda (p)
                 (let
                     ; rotation matrix for rotation around x axis
                     ((Rx (matrix3x3 (vector3 1 0 0)
                                     (vector3 0 (cos (phi p)) (sin (phi p)))
                                     (vector3 0 (- (sin (phi p))) (cos (phi p))))))
                   (set! lc-epsilon (matrix3x3* Rx (matrix3x3* eps-diag (matrix3x3-transpose Rx))))
                   (set! lc-epsilon-diag (vector3 (vector3-x (vector3-x lc-epsilon))
                                                  (vector3-y (vector3-y lc-epsilon))
                                                  (vector3-z (vector3-z lc-epsilon))))
                   (set! lc-epsilon-offdiag (vector3 (vector3-x (vector3-y lc-epsilon))
                                                     (vector3-x (vector3-z lc-epsilon))
                                                     (vector3-y (vector3-z lc-epsilon))))
                   (make medium (epsilon-diag lc-epsilon-diag) (epsilon-offdiag lc-epsilon-offdiag)))))

(set! geometry (list
                (make block
                  (center (+ (* -0.5 sx) (* 0.5 (+ dpml dsub))) 0 0)
                  (size (+ dpml dsub) infinity infinity)
                  (material (make medium (index n0))))
                (make block
                  (center (+ (* -0.5 sx) dpml dsub d) 0 0)
                  (size (* 2 d) infinity infinity)
                  (material (make material-function (material-func lc-mat))))))

(define tran-flux2 (add-flux fcen 0 1 (make flux-region (center tran-pt) (size 0 sy 0))))

(run-sources+ 300)

(define res1 (get-eigenmode-coefficients tran-flux2 (arith-sequence 1 1 nmode) #:eig-parity (+ ODD-Z EVEN-Y)))
(define res2 (get-eigenmode-coefficients tran-flux2 (arith-sequence 1 1 nmode) #:eig-parity (+ EVEN-Z ODD-Y)))

(define t-coeffs1 (list-ref res1 0))
(define t-coeffs2 (list-ref res2 0))
(define kdom (list-ref res1 3))

(map (lambda (nm)
       (let ((mode-angle (acos (/ (vector3-x (list-ref kdom nm)) fcen)))
             (mode-tran (/ (+ (sqr (magnitude (array-ref t-coeffs1 nm 0 0))) (sqr (magnitude (array-ref t-coeffs2 nm 0 0)))) (list-ref input-flux 0))))
         (print "tran:, " nm ", " (rad->deg mode-angle) ", " mode-tran "\n")))
     (arith-sequence 0 1 nmode))