meepefc.txt

计算光子晶体的等频线,研究光子晶体MEEP软件使用到

;       -------------- calculation of equai-frequency concours--------------
;
; call ctl like this
;
; mpb band=3 f=0.4 k-find.ctl | grep iso | awk -F, '{print $2,$3}' >test.dat
; and itwill compute the 3rd band at frequency 0.4. 
;
; mpb band=2 f=0.4 k-find.ctl | grep iso | awk -F, '{print $2,$3}' >borders.dat
; will compute the 2nd band at frequency 0.4. , but there is no solution, so the
; borders of the brillouin-zone are plotted
;
; xmgrace test.dat borders.dat
; will plot the efc
;
 
;-------------------------------------------------
;-------------------------------------------------
; Define some parameters
;-------------------------------------------------
;-------------------------------------------------
(define-param r 0.42)
(define-param band 3)
(define-param f 0.4)
 

;-------------------------------------------------
;-------------------------------------------------
; Setup hexagonal lattice
;-------------------------------------------------
;-------------------------------------------------
(set! geometry-lattice (make lattice (size 1 1 no-size)
                         (basis1 (/ (sqrt 3) 2) 0.5)
                         (basis2 (/ (sqrt 3) 2) -0.5)))
 
;-------------------------------------------------
;-------------------------------------------------
; Define high symmetry points for hexagonal lattice
;-------------------------------------------------
;-------------------------------------------------
(define M (vector3 0.0 0.5 0));
(define K (vector3 (/ -3) (/ 3) 0));
 
;-------------------------------------------------
;-------------------------------------------------
; Set up host material
;-------------------------------------------------
;-------------------------------------------------
(set! default-material (make dielectric (epsilon 11.6)))
;-------------------------------------------------
;-------------------------------------------------
; Set up the object of our system
;-------------------------------------------------
;-------------------------------------------------
(set! geometry (list (make cylinder
                     (center 0 0 0) (radius r) (height infinity)
                     (material (make dielectric (epsilon 1))))))

;-------------------------------------------------
;-------------------------------------------------
; Choose the resolution
;-------------------------------------------------
;-------------------------------------------------
(set! resolution 32)
(set! mesh-size 15)
;-------------------------------------------------
;-------------------------------------------------
; Define a new find-k function
; returns iso-k-vector in cartesian coordinates
;-------------------------------------------------
;-------------------------------------------------
(define ((k-find pol kdir kmax b target mag) dw)
(vector3-scale (list-ref (find-k pol (+ target dw) b b kdir 1e-4 mag 0 kmax ) 0) (unit-vector3 (reciprocal->cartesian kdir)))
)
;-------------------------------------------------
;-------------------------------------------------
; Define a function for plotting iso-k-vectors
; if a stop gap occurs  (v=#f) the vector v2 should be printed
; (gives the edge of the brillouin-zone)
;-------------------------------------------------
;-------------------------------------------------
(define (print-iso v v2)
(if v
(print "iso, " (vector3-x v) ", " (vector3-y v) ", " (vector3-z v) "\n")
(print "iso, " (vector3-x v2) ", " (vector3-y v2) ", " (vector3-z v2) "\n")
)
)
;-------------------------------------------------
;-------------------------------------------------
; Define the borders of the brillouin-zone to seek for EFC
; To plot the whole EFC, one has to mirrow and rotate about n*60?°
;-------------------------------------------------
;-------------------------------------------------
(define k-points2 (interpolate 19 (list M K)))

;-------------------------------------------------
;-------------------------------------------------
; Print for every k-point in k-points2-list the result
; If an error occurs the result will be #f (e.g a stop direction has been found)
; so we give (print-iso) the borders, too
; TM= Ez-Pol
; TE= Hz-Pol
;
;-------------------------------------------------
;-------------------------------------------------
(for-each (lambda (k)
(print-iso (false-if-exception ((k-find TM k (vector3-norm  (reciprocal->cartesian k)) band f 0.3) 0)) (reciprocal->cartesian k))
) k-points2)