check.ctl

麻省理工的计算光子晶体的程序

; Test suite for MPB.  This file runs MPB for a variety of cases,
; and compares it against known results from previous versions.  If the
; answers aren't sufficiently close, it exits with an error.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Some general setup and utility routines first:

(set! tolerance 1e-9) ; use a low tolerance to get consistent results

; keep track of some error statistics:
(define min-err infinity)
(define max-err 0)
(define sum-err 0)
(define num-err 0)

; function to check if two results are sufficently close:
(define-param check-tolerance 1e-4)
(define (almost-equal? x y)
  (if (> (abs x) 1e-3)
      (let ((err (/ (abs (- x y)) (* 0.5 (+ (abs x) (abs y))))))
	(set! min-err (min min-err err))
	(set! max-err (max max-err err))
	(set! num-err (+ num-err 1))
	(set! sum-err (+ sum-err err))))
  (or 
   (< (abs (- x y)) (* 0.5 check-tolerance (+ (abs x) (abs y))))
   (and (< (abs x) 1e-3) (< (abs (- x y)) 1e-3))))

; Convert a list l into a list of indices '(1 2 ...) of the same length.
(define (indices l)
  (if (null? l)
      '()
      (cons 1 (map (lambda (x) (+ x 1)) (indices (cdr l))))))

; Check whether the freqs returned by a run (all-freqs) match correct-freqs.
(define (check-freqs correct-freqs)
  (define (check-freqs-aux fc-list f-list ik)
    (define (check-freqs-aux2 fc f ib)
      (if (not (almost-equal? fc f))
	  (error "check-freqs: k-point " ik " band " ib " is "
		 f " instead of " fc)))
    (if (= (length fc-list) (length f-list))
	(map check-freqs-aux2 fc-list f-list (indices f-list))
	(error "check-freqs: wrong number of bands at k-point " ik)))
  (if (= (length correct-freqs) (length all-freqs))
      (begin
	(map check-freqs-aux correct-freqs all-freqs (indices all-freqs))
	(print "check-freqs: PASSED\n"))
      (error "check-freqs: wrong number of k-points")))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(if (not (using-mpi?)) ; MPI code currently doesn't support 1d systems
(begin

; Use a lower tolerance for the 1d cases, since it is cheap; otherwise,
; the Bragg-sine case perennially causes problems.
(set! tolerance (/ tolerance 10000))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; First test: a simple 1d Bragg mirror:

(print
 "**************************************************************************\n"
 " Test case: 1d quarter-wave stack.\n"
 "**************************************************************************\n"
)

(set! geometry (list (make cylinder (material (make dielectric (epsilon 9.0)))
			   (center 0) (axis 1)
			   (radius infinity) (height 0.25))))
(set! k-points (interpolate 4 (list (vector3 0 0 0) (vector3 0.5 0 0))))
(set! grid-size (vector3 32 1 1))
(set! num-bands 8)

(define correct-freqs '((0.0 0.648351064758882 0.666667517098436 1.29488075758612 1.33336075391282 1.93757672906421 2.00024045502069 2.57413377730507) (0.0567106459132849 0.599851835766128 0.715264618493464 1.2533516328214 1.37508036430362 1.9023030310515 2.03577843681173 2.5447607302785) (0.111808338549745 0.544964035370552 0.770470132482679 1.19886555431273 1.43019048406127 1.84869942534301 2.09026194963917 2.49305185532527) (0.162554443030826 0.494234387214952 0.821807214979002 1.14762423868064 1.48265526473298 1.79656499475259 2.14420182927269 2.44000721223413) (0.202728586444506 0.454051807431849 0.862903053647553 1.1065252897017 1.52568848270994 1.75360874628753 2.19029794218108 2.3942414201896) (0.219409188989471 0.437366603189744 0.880190598314617 1.08923081762878 1.5443398403343 1.7349704450792 2.21143711778794 2.3731975719957)))

(run-tm)
(check-freqs correct-freqs)

(run-te)
(check-freqs correct-freqs)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Second test: a less-simple 1d Bragg mirror, consisting of a sinusoidally
; varying dielectric index (see also bragg-sine.ctl):

(print
 "**************************************************************************\n"
 " Test case: 1d sinusoidal Bragg mirrors.\n"
 "**************************************************************************\n"
)

(define pi (* 4 (atan 1))) ; 3.14159...
(define (eps-func p)
  (make dielectric (index (+ 2 (cos (* 2 pi (vector3-x p)))))))
(set! default-material (make material-function (material-func eps-func)))

(set! k-points (interpolate 9 (list (vector3 0 0 0) (vector3 0.5 0 0))))
(set! grid-size (vector3 32 1 1))
(set! num-bands 8)

(run-tm)
(check-freqs '((0.0 0.460752283717811 0.542059418927331 0.968412011064827 1.01592861336746 1.48283705488259 1.48374633575784 1.9677172977623) (0.0231436084314742 0.454366673639236 0.548556852963559 0.958177155747471 1.02618660946863 1.45880923902184 1.50778401398672 1.94923038830116) (0.046211542432863 0.43910864609908 0.56416208641797 0.938077927197344 1.0463565943812 1.43434696639083 1.53227646120221 1.92515716944815) (0.069114503567596 0.420004628065799 0.583886947263537 0.915558839464672 1.06899842076893 1.40990248972574 1.55677337160932 1.90076649748739) (0.0917307818062977 0.399446033671752 0.605415353767922 0.892293385266758 1.09244722981119 1.38548487907905 1.58126958017841 1.87630217942291) (0.113873900258811 0.378515461211418 0.627800248852586 0.868769802653355 1.11623067344131 1.36110772786128 1.60575883658692 1.85181932375803) (0.135229085992452 0.357912554145201 0.650585137964263 0.845213453792503 1.14015387025798 1.33679639228579 1.63023029383128 1.82734480769861) (0.155218781840958 0.338398194422463 0.673453381248859 0.821821692886828 1.16408502920442 1.31260712790103 1.65466024698586 1.8029089496478) (0.17271618668617 0.321193301395356 0.695971967442052 0.798962935723639 1.18782670146768 1.28869970774005 1.67898024504888 1.77858329794573) (0.185552735718465 0.308516113919293 0.716800284215745 0.777942357002287 1.21065518331675 1.26577816162819 1.70286746606745 1.75469220803656) (0.190472580448255 0.303646927430686 0.728334358202394 0.76634602375042 1.22537034565911 1.25103296039095 1.72050247023166 1.73705584292802)))

(set! default-material air) ; don't screw up later tests

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(set! tolerance (* tolerance 10000))

)) ; if (not (using-mpi?))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Square lattice of dielectric rods in air.