; 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.

(print
 "**************************************************************************\n"
 " Test case: Square lattice of dielectric rods in air.\n"
 "**************************************************************************\n"
)

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

(run-te)
(check-freqs '((0.0 0.567309249824085 0.785322464218187 0.78533907114069 0.919065314955226 1.01180231506563 1.01180439396285 1.0986987562052) (0.089747533106419 0.565478962844045 0.770909925689672 0.786976045327448 0.911716296164207 1.00922747558406 1.01332658097629 1.12204775958112) (0.178766635242442 0.558183927765058 0.733727796380775 0.791352310664549 0.894869659413299 1.01797556826237 1.01831492275664 1.12431478868499) (0.266049753471709 0.538687892094626 0.690403576061145 0.796919728740922 0.879307985310368 1.02568655836391 1.04031515320176 1.11602169025839) (0.34960857098187 0.496962222790986 0.66067962204941 0.801560565910096 0.869495965837476 1.0350024763644 1.06884798601087 1.10467620761345) (0.413082417117357 0.446395530337967 0.65145059507318 0.803368099819892 0.866189516688412 1.04010576541796 1.09753776038 1.09827456525779) (0.424034318475705 0.450850117081375 0.647112289739006 0.807637589363962 0.86115857996382 0.99417414554491 1.0595794172637 1.11992767773851) (0.455082616824108 0.463338791124478 0.635212403292713 0.819735856354933 0.838315627232708 0.939938902646522 1.01486587270942 1.12838695905218) (0.481198597820564 0.501609229763058 0.618522175753187 0.786155073927434 0.838317450345758 0.913939091232541 0.970884228080751 1.13184861079901) (0.499280260937176 0.556434275007696 0.601659180916135 0.721570044167946 0.862218793047699 0.905816154022519 0.931039089966407 1.13347236398634) (0.507900932587 0.593526292833574 0.593596593124933 0.680566062367046 0.88387598815364 0.903945296603264 0.903949249935734 1.13397288782519) (0.4766282365696 0.551926957646813 0.607311768805225 0.746127961399139 0.852479455285177 0.903785999923616 0.949896076238274 1.13279836961847) (0.373412079642824 0.547353877862403 0.645745076852373 0.820546048524794 0.833821226258539 0.90364772169965 1.02088954489018 1.12643735967833) (0.252290144958668 0.555861520941481 0.701023986684804 0.802063555906004 0.904841846692795 0.906101988454332 1.09483334397003 1.10211421754981) (0.126868830236271 0.564051618328441 0.757943082512391 0.789566003897487 0.91058219271459 0.971875854812309 1.04852328393723 1.12559272973884) (0.0 0.567309249697206 0.785322464388339 0.785339071278791 0.919065315098791 1.01180231580376 1.0118043941002 1.09869875597474)))

(run-tm)
(check-freqs '((0.0 0.550430733979162 0.566241371578254 0.566241371590288 0.835198392541763 0.872522248550314 0.972318183103308 1.08920254179769) (0.0654861240038261 0.527080481428998 0.566799276739625 0.589165207611584 0.835576055520486 0.870980062620355 0.960846516368135 1.05902201031889) (0.128358448231762 0.495999172056421 0.568271183796953 0.619721421825485 0.833562960542279 0.867048180733246 0.930077475889248 1.04306779424378) (0.185088034665959 0.463940092597231 0.570113900868756 0.653567078126228 0.817881354346977 0.862372411985142 0.90187550142523 1.04123744478054) (0.229126816125475 0.435593887869437 0.571624367586527 0.691156829464893 0.782224998842116 0.858716310201942 0.891825123400888 1.04461063360963) (0.247300696006748 0.422810121349707 0.572206120768814 0.722519491352313 0.749366763369664 0.857346110926564 0.890391539673791 1.0468678713734) (0.250809748878834 0.429554276590146 0.565096567434357 0.720355410169051 0.758471422549885 0.858178250955784 0.890488360212771 1.03491655457074) (0.260255079137602 0.448078995516004 0.547933657276557 0.713707158688594 0.782090347346113 0.86101638088251 0.890741336275652 1.00608368845575) (0.272587744232561 0.473364092545489 0.528315799682105 0.702811847529422 0.814130266536548 0.8667463670677 0.891050068065077 0.96973788503882) (0.28329017996901 0.496884116939824 0.512907153651487 0.690607961839939 0.850483079390603 0.876011183352539 0.89126643613903 0.93077987299414) (0.287601798196358 0.507006838456059 0.507006838475394 0.684706840348946 0.883530037853841 0.883530037866573 0.889824048148961 0.898389144733936) (0.277734985850212 0.494789217195798 0.512794766069262 0.693585317484579 0.841827176583588 0.86555977664715 0.91182965298438 0.91364862207229) (0.241654191448818 0.481926694377318 0.527837662343768 0.687466863320862 0.831639003520604 0.850551381755549 0.915313305877234 0.946702579779498) (0.176222602981777 0.490886438976766 0.546175771973787 0.650073604194652 0.841125231854928 0.853522369326247 0.92744923560386 0.987456182887339) (0.0920144002040036 0.518511998133382 0.560750416146476 0.60406552818857 0.836485468004983 0.868856191329739 0.954250999892416 1.03716877067079) (0.0 0.550430733878627 0.566241371572124 0.566241371574826 0.835198392501417 0.872522248512609 0.972318183150721 1.0892025391774)))

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

; Using the targeted solver to find a defect state in a 5x5 triangular
; lattice of rods.

(print
 "**************************************************************************\n"
 " Test case: 3x3 triangular lattice of rods in air, dipole defect states.\n"
 "**************************************************************************\n"
)

(set! geometry-lattice (make lattice (size 3 3 1)
                         (basis1 (/ (sqrt 3) 2) 0.5)
                         (basis2 (/ (sqrt 3) 2) -0.5)))
(set! k-points (list (vector3 0 0.5 0))) ; K
(set! geometry (list
		(make cylinder (material (make dielectric (epsilon 12))) 
		      (center 0 0) (radius 0.2) (height infinity))))
(set! geometry (geometric-objects-lattice-duplicates geometry))
(set! geometry (append geometry 
                       (list (make cylinder (center 0 0 0) 
                                   (radius 0.33) (height infinity)
                                  (material (make dielectric (epsilon 12)))))))
(set! grid-size (vector3 (* 16 5) (* 16 5) 1))
(set! num-bands 2)
(set! target-freq 0.35)
(run-tm)

(define ct-save check-tolerance)
(set! check-tolerance (* ct-save 10))
(check-freqs '((0.335269274821527 0.337652210814253)))
(set! check-tolerance ct-save)

(set! target-freq 0)

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

(print
 "**************************************************************************\n"
 " Test case: fcc lattice of air spheres in dielectric.\n"
 "**************************************************************************\n"
)

(set! geometry-lattice (make lattice
			 (basis1 0 1 1)
			 (basis2 1 0 1)
			 (basis3 1 1 0)))
(set! k-points (interpolate 1 (list
			       (vector3 0 0.5 0.5)            ; X
			       (vector3 0 0.625 0.375)        ; U
			       (vector3 0 0.5 0)              ; L
			       (vector3 0 0 0)                ; Gamma
			       (vector3 0 0.5 0.5)            ; X
			       (vector3 0.25 0.75 0.5)        ; W
			       (vector3 0.375 0.75 0.375))))  ; K
(set! geometry (list (make sphere (center 0) (radius 0.5) (material air))))
(set! default-material (make dielectric (epsilon 11.56)))
(set! grid-size (vector3 16 16 16))
(set! mesh-size 5)
(set! num-bands 10)
(run)
(check-freqs '((0.368116356209037 0.369980241542025 0.379067426858622 0.380355825481755 0.49108300146986 0.509321674949346 0.519652077187981 0.521986848659634 0.589861118480048 0.657547105547574) (0.366168952442295 0.375372406889818 0.381873532967479 0.384010857377027 0.468647962388753 0.503195261165308 0.5211472033283 0.529177355618733 0.605687520175585 0.641010575966524) (0.355246298748627 0.378149323797633 0.390663162132607 0.398125290587968 0.435685869130532 0.491330183134628 0.525117050044954 0.539115916262814 0.632010879838287 0.636550876453458) (0.321237491869962 0.32910168786935 0.394743530898447 0.397841804199817 0.460891072117548 0.511490690292887 0.531312372506223 0.544975783308705 0.626070719045844 0.638849253793513) (0.304996789723474 0.306490030344765 0.38399757686853 0.386731926388109 0.49072085711492 0.534558361399195 0.535010673589626 0.538001547079021 0.618866648935078 0.622227387005411) (0.177728088406801 0.178577250117289 0.470894345930347 0.474296006694487 0.503036395546177 0.534075106832263 0.536910305396605 0.53914018449749 0.61872969865244 0.62126604417647) (0.0 0.0 0.517651812996759 0.521927400026579 0.521927400063738 0.542668842573521 0.54266884259157 0.545547143439792 0.606298961911071 0.608727263482667) (0.205019177848376 0.205741665417928 0.47151402965987 0.473988482367252 0.505357910726724 0.524526583531831 0.52820880373352 0.530578318673899 0.599219767232565 0.650319562639214) (0.368116356275522 0.369980241620613 0.379067426974632 0.380355825610059 0.49108300157512 0.509321674971744 0.519652077334165 0.521986848798089 0.589861118539645 0.657547147736461) (0.370158114992911 0.374979196275145 0.382883733677173 0.389256458552047 0.46024288611451 0.498652679435423 0.50235588768182 0.54800304095417 0.61691332252954 0.630865854030288) (0.371583334129105 0.38388611636557 0.384082592553739 0.407084651822841 0.431532334012398 0.487181534553358 0.488970298812729 0.56643459591284 0.619483995411123 0.64959455705408) (0.36210028608351 0.37995099453707 0.387759257999464 0.402824645197872 0.434538805815819 0.490156436339851 0.502759548908983 0.558462259582065 0.626513356662261 0.64047989670839) (0.356962911196037 0.377342903138122 0.389053322326131 0.399531555896904 0.434217203098185 0.491795282316179 0.527518754499718 0.53933002172486 0.633265155959861 0.637548652538048)))

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

(print
 "**************************************************************************\n"
 " Test case: simple cubic lattice with anisotropic dielectric.\n"
 "**************************************************************************\n"
)

(set! geometry-lattice (make lattice))
(set! default-material air)
(set! k-points (list (vector3 0) (vector3 0.5)
		     (vector3 0.5 0.5) (vector3 0.5 0.5 0.5)))
(set! grid-size (vector3 16 16 16))
(set! mesh-size 5)
(define hi-all (make dielectric (epsilon 12)))
(define hi-x (make dielectric-anisotropic (epsilon-diag 12 1 1)))
(define hi-y (make dielectric-anisotropic (epsilon-diag 1 12 1)))
(define hi-z (make dielectric-anisotropic (epsilon-diag 1 1 12)))
(set! geometry
	(list (make block (center 0) (size 0.313 0.313 1) (material hi-z))
	      (make block (center 0) (size 0.313 1 0.313) (material hi-y))
	      (make block (center 0) (size 1 0.313 0.313) (material hi-x))
	      (make block (center 0) (size 0.313 0.313 0.313) 
		    (material hi-all))))
(set! num-bands 3)
(run)
(check-freqs '((0.0 0.0 0.54604008820461) (0.250972002681914 0.250975086788192 0.428868017659483) (0.291417685034807 0.33426762303395 0.48002200778248) (0.351623530055896 0.3516235300955 0.483724678054855)))

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

(display-eigensolver-stats)
(print "Relative error ranged from " min-err " to " max-err
	      ", with a mean of " (/ sum-err num-err) "\n")
(print "PASSED all tests.\n")