libmeep-mpi-default8 1.3-3build5 / usr / share / meep-mpi-default / casimir.scm

This file is indexed.

/usr/share/meep-mpi-default/casimir.scm is in libmeep-mpi-default8 1.3-3build5.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below. 

  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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
;given m1 m2, make a composit index n; the inverse is casimir-source-info below
; m1 = r-c, m2 = c  => r = m1+m2;  n - s = m2 => n = m2 + s; s = (sum_(k=1)^r k)
; => n = m2 + (sum_(k=1)^(m1+m2) k) = m2 + 1/2 (m1+m2) * (m1+m2+1)
(define (make-casimir-src-index m1 m2) (+ m2 (* (/ 2) (+ m1 m2) (+ m1 m2 1))))

; return a list (source-vol mx my mz)
; given the volume integration-vol and n, pick out the appropriate side and mx my mz to use
; sides are ordered by decreasing weight
; weights are w(side, m) = area(side)/total area * 1/(m+1)^4, a rough estimate of the 
; contribution to the stress tensor from that side and that m 
(define (casimir-source-info integration-vol n)
  (define (modround x n) (modulo (inexact->exact (round x)) n))
  (define (get-src-index n) ;given n, extract out the two values of m for 3-d 
    (let* ((s 0) ;sum of diagonals
	 (r 0) ;row intersection
	 (c 0));column intersection
      (while (< (+ s r) n) 
	     (set! r (+ r 1))
	     (set! s (+ s r)))
      (set! c (- n s))
      (list (- r c) c)))

  (let* ((min-corner (meep-volume-get-min-corner integration-vol))
	 (max-corner (meep-volume-get-max-corner integration-vol))
	 (size-vec (vector3- max-corner min-corner))
	 (center-vec (vector3+ (vector3-scale 0.5 size-vec) min-corner))
	 (sx (vector3-x size-vec))
	 (sy (vector3-y size-vec))
	 (sz (vector3-z size-vec))
	 (xshift (vector3 (/ sx 2)    0       0))
	 (yshift (vector3    0    (/ sy 2)    0))
	 (zshift (vector3    0        0   (/ sz 2))))
    (if (and (> sy 1e-15) (> sz 1e-15))  ;3d cartesian: n = 6*f(m1,m2) + s
	(let* ((s (modround n 6))
	       (nr (/ (- n s) 6))
	       (ms (get-src-index nr)) ;get (m1 m2)
	       (m1 (first ms))
	       (m2 (second ms))
	       (x-const (vector3 0 sy sz))
	       (y-const (vector3 sx 0 sz))
	       (z-const (vector3 sx sy 0))
	       (center-list
		(list (vector3- center-vec xshift) (vector3+ center-vec xshift)
		      (vector3- center-vec yshift) (vector3+ center-vec yshift)
		      (vector3- center-vec zshift) (vector3+ center-vec zshift)))
	       (m-list 
		(list (vector3 0 m1 m2) (vector3 0 m1 m2) (vector3 m1 0 m2) 
		      (vector3 m1 0 m2) (vector3 m1 m2 0) (vector3 m1 m2 0)))
	       (orientation-list (list -1 1 -1 1 -1 1))
	       (size-list
		(list x-const x-const y-const y-const z-const z-const))
	       (surface-vol
		(volume (center (list-ref center-list s)) (size (list-ref size-list s))))
	       (surface-m (list-ref m-list s)))
	  (print "Computing in 3d\n")
	  (list surface-vol
		(vector3-x surface-m) (vector3-y surface-m) (vector3-z surface-m)
		(list-ref orientation-list s) 1))
	(if (= dimensions -2) ;cylindricals - must make sure that the volume has only r >= 0
	    (let* ((3-sides? (if (<= (vector3-x min-corner) 0) true false)) ;volume passes through the origin
		   (s (if 3-sides? (modround n 3) (modround n 4)))
		   (nr (if 3-sides? (/ (- n s) 3) (/ (- n s) 4))) ;reduced index
		   (ms (get-src-index nr)) ;extract out both m-phi and m-dct
		   (m-phi (first ms))
		   (m-dct (second ms))
		   (DR (if 3-sides? (/ 1 resolution) 0)) ;cannot include r = 0!!
		   (sr (- (vector3-x max-corner) (+ (max 0 (vector3-x min-corner)) DR)))
		   (r-cen (+ (* 0.5 sr) (max 0 (vector3-x min-corner)) DR))
		   (new-center-vec (vector3 r-cen 0 (vector3-z center-vec)))
		   (r-shift (vector3 (/ sr 2) 0 0))
		   (z-shift (vector3 0 0 (/ sz 2)))
		   (r-const-size (vector3 0 0 sz))
		   (z-const-size (vector3 sr 0 0))
		   (center-list ;if 3-sides? = true, only first 3 list elements are used
		    (list (vector3- new-center-vec z-shift) (vector3+ new-center-vec z-shift)
			  (vector3+ new-center-vec r-shift) (vector3- new-center-vec r-shift)))
		   (m-list
		    (list (vector3 m-dct m-phi 0) (vector3 m-dct m-phi 0) 
			  (vector3 0 m-phi m-dct) (vector3 0 m-phi m-dct)))
		   (orientation-list (list -1 1 1 -1))
		   (size-list
		    (list z-const-size z-const-size r-const-size r-const-size))
		   (surface-vol
		    (volume (center (list-ref center-list s)) (size (list-ref size-list s))))
		   (surface-m (list-ref m-list s)))
	      (print "Computing in Cylindrical coordinates: m-phi = "m-phi", m-dct = "m-dct", 3-sides? = "3-sides?"\n")
	      (list surface-vol
		    (vector3-x surface-m) (vector3-y surface-m) (vector3-z surface-m)
		    (list-ref orientation-list s) (if (or (= s 0) (= s 1)) 1 0)))
	    (let* ((s (modround n 4)) ;2d or quasi-3d cartesian: n = 4m + s, no ambiguity in m
		   (m (/ (- n s) 4))
		   (x-const-size (vector3 0 sy )) ;sz may be non-zero for quasi-3d systems
		   (y-const-size (vector3 sx 0 ))
		   (center-list 
		    (list (vector3- center-vec xshift) (vector3+ center-vec xshift)
			  (vector3- center-vec yshift) (vector3+ center-vec yshift)))
		   (m-list
		    (list (vector3 0 m 0) (vector3 0 m 0) (vector3 m 0 0) (vector3 m 0 0)))
		   (orientation-list (list -1 1 -1 1))
		   (size-list
		    (list x-const-size x-const-size y-const-size y-const-size))
		   (surface-vol
		    (volume (center (list-ref center-list s)) (size (list-ref size-list s))))
		   (surface-m (list-ref m-list s)))
	      (print "Casimir.scm: working in 2 dimensions\n")
	      (print "  Surface center: "(list-ref center-list s)"\n")
	      (print "  Surface size: "(list-ref size-list s)"\n")
	      (list surface-vol 
		    (vector3-x surface-m) (vector3-y surface-m) (vector3-z surface-m)
		    (list-ref orientation-list s) 1))))))

;compute the casimir force for a single n and single polarization
;n contains both the side number and the harmonic expansion index
(define (casimir-force-contrib force-direction integration-vol n Sigma T source-component gt . step-funcs)
  (define (cos-func X mx my mz source-vol)
    (let* 
	((min-corner (meep-volume-get-min-corner source-vol))
	 (max-corner (meep-volume-get-max-corner source-vol))
	 (size-vec (vector3- max-corner min-corner))
	 (X-start (vector3+ X (vector3-scale 0.5 size-vec)))
	 (sx (vector3-x size-vec))
	 (sy (vector3-y size-vec))
	 (sz (vector3-z size-vec))
	 (x (vector3-x X-start))
	 (y (vector3-y X-start))
	 (z (vector3-z X-start))
	 (kx (if (> sx 1e-15) (/ (* mx pi) sx) 0))
	 (ky (if (> sy 1e-15) (/ (* my pi) sy) 0))
	 (kz (if (> sz 1e-15) (/ (* mz pi) sz) 0))
	 (Nx (if (> sx 1e-15) (/ (if (= mx 0) 1 2) sx) 1))
	 (Ny (if (> sy 1e-15) (/ (if (= my 0) 1 2) sy) 1))
	 (Nz (if (> sz 1e-15) (/ (if (= mz 0) 1 2) sz) 1)))
      (* (sqrt (* Nx Ny Nz)) (cos (* kx x)) (cos (* ky y)) (cos (* kz z)))))

  (let* ((ft (meep-type source-component))
	 (source-info (casimir-source-info integration-vol n))
	 (source-vol (first source-info))
	 (mx (second source-info))
	 (my (third source-info)) ;m-phi in cylindrical coordinates
	 (mz (fourth source-info))
	 (source-orientation (fifth source-info))
	 (dt (/ Courant resolution)))
    (if (= ft E-stuff)
	(begin
	  (set! global-D-conductivity Sigma)
	  (set! global-B-conductivity 0))
	(begin
	  (set! global-B-conductivity Sigma)
	  (set! global-D-conductivity 0)))
    (if (eq? dimensions -2) 
	(begin (print "Cylindricals: m = "my" and (nr nz) = ("mx", "mz")\n")
	       (print "  surface center = "(meep-volume-center source-vol)"\n")
	       (print "  source size = "(vector3- 
					  (meep-volume-get-max-corner source-vol) 
					  (meep-volume-get-min-corner source-vol)))
	       (set! m my))) ;set exp(i m phi) field dependence
    (set! sources
	  (list (make source
		  (src (make custom-src ; delta function pulse
			 (src-func (lambda (t) (/ 1 dt)))
			 (start-time (* -0.25 dt))
			 (end-time (* 0.75 dt))
			 (is-integrated? false)))
		  (center (meep-volume-center source-vol))
		  (size (vector3- (meep-volume-get-max-corner source-vol)
				  (meep-volume-get-min-corner source-vol)))
		  (component source-component)
		  (amp-func (lambda (p) (cos-func p mx my mz source-vol))))))
    (reset-meep)
    (init-fields)
    (let* ((counter 0)
	   (force-integral 0))
      (define (integrate-function)
	(let* ((f-temp (meep-fields-casimir-stress-dct-integral
			fields
			force-direction 
			(meep-component-direction source-component)
			mx (if (eq? dimensions -2) 0 my) mz
			ft source-vol)))
	  (set! force-integral
		(+ force-integral
		   (imag-part (* (list-ref gt counter) dt 
				 source-orientation
				 (if (eq? dimensions -2)
				     (* (if (eq? my 0) 1 2)
					(real-part f-temp))
				     f-temp)))))
	  (set! counter (+ counter 1))))
      (apply run-until (cons (- T 1) (cons integrate-function step-funcs)))
      force-integral)))

;%%%%%%%%%%%%%%%%%%%%% BLOCH PBCS %%%%%%%%%%%%%%%%%%%%%%
;here the source is specified in 
;the form exp(i k x), k = pi/L (m + k_red), 
;m = (mx,my,mz) are integers (reciprocal lattice vectors
;k_red = (kx,ky,kz) is in the 1st BZ, m an integer
;source-vol is assumed to occupy one entire plane intersecting
;the computational cell, so we don't need to extract out
;its information - there is only one side to it

;pass the vector (mx my mz) and (kx ky kz) ready-made, since
;this surface consists of only one face
(define (casimir-force-contrib-bloch force-direction source-vol k-vec Sigma T source-component gt . step-funcs)
  ;sources of the form exp(i g x); surface integration in
  ;casimir.cpp integrates against exp(-i g x)
  (define (casimir-bloch-func X gx gy gz source-vol)
    (let*
	((min-corner (meep-volume-get-min-corner source-vol))
	 (max-corner (meep-volume-get-max-corner source-vol))
	 (size-vec (vector3- max-corner min-corner)) ;crossection of computational cell
	 (sx (vector3-x size-vec))
	 (sy (vector3-y size-vec))
	 (sz (vector3-z size-vec))
	 (x (vector3-x X))
	 (y (vector3-y X))
	 (z (vector3-z X))
	 (Kx (if (> sx 1e-15) (/ (* gx pi) 1) 0)) ;phase winding is independent of unit cell size
	 (Ky (if (> sy 1e-15) (/ (* gy pi) 1) 0))
	 (Kz (if (> sz 1e-15) (/ (* gz pi) 1) 0))
	 (Nx (if (> sx 1e-15) (/ sx) 1))
	 (Ny (if (> sy 1e-15) (/ sy) 1))
	 (Nz (if (> sz 1e-15) (/ sz) 1)))
      (* (sqrt (* Nx Ny Nz)) (exp (* (sqrt -1) (+ (* Kx x) (* Ky y) (* Kz z)))))))

  (let* ((ft (meep-type source-component))
	 (dt (/ Courant resolution))
	 ;Bloch phases - exp( i * (2*pi*m + pi*k) x/L)
	 (gx (vector3-x k-vec))
	 (gy (vector3-y k-vec))
	 (gz (vector3-z k-vec)))
    (set! force-complex-fields? true)
    (set! k-point (vector3-scale 0.5 k-vec))
    (if (= ft E-stuff)
	(begin
	  (set! global-D-conductivity Sigma)
	  (set! global-B-conductivity 0))
	(begin
	  (set! global-D-conductivity 0)
	  (set! global-B-conductivity Sigma)))
    (set! sources
	  (list (make source
		  (src (make custom-src
			 (src-func (lambda (t) (/ 1 dt)))
			 (start-time (* -0.25 dt))
			 (end-time (* 0.75 dt))
			 (is-integrated? false)))
		  (center (meep-volume-center source-vol))
		  (size (vector3- (meep-volume-get-max-corner source-vol)
				  (meep-volume-get-min-corner source-vol)))
		  (component source-component)
		  (amp-func (lambda (p) (casimir-bloch-func p gx gy gz source-vol))))))
    (reset-meep)
    (init-fields)
    (let* ((counter 0)
	   (force-integral 0))
      (define (integrate-function)
	(set! force-integral
	      (+ force-integral
		 (imag-part (* (list-ref gt counter) dt
			       (meep-fields-casimir-stress-dct-integral
				fields
				force-direction
				(meep-component-direction source-component)
				gx gy gz
				ft source-vol true)))))
	(set! counter (+ counter 1)))
      (apply run-until (cons (- T 1) (cons integrate-function step-funcs)))
      force-integral)))

APT Browse - Built by Thomas Orozco.