meep.scm.in

; Copyright (C) 2005-2022 Massachusetts Institute of Technology
;
; This program is free software; you can redistribute it and/or modify
; it under the terms of the GNU General Public License as published by
; the Free Software Foundation; either version 2 of the License, or
; (at your option) any later version.
;
; This program is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
; GNU General Public License for more details.
;
; You should have received a copy of the GNU General Public License
; along with this program; if not, write to the Free Software
; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

; ****************************************************************

; for define* (keyword args) in guile 1.8
(if (< (string->number (major-version)) 2)
    (use-modules (ice-9 optargs)))

; Get the number of arguments to a function p.  However, some
; older versions of Guile (e.g. 1.2) do not support the 'arity
; property, and procedure-property just returns false.  In
; this case, we assume that the procedure returns 1 argument,
; as this is the most useful default for our purposes.  Sigh.

(define (procedure-num-args p)
  (let ((arity (if (defined? 'procedure-minimum-arity)
                   (procedure-minimum-arity p)
                   (procedure-property p 'arity))))
    (if arity (car arity) 1)))

; ****************************************************************

; Set print-ok? to whether or not we are the MPI master process.
; Also, MPI doesn't support interactive mode.  However, don't try this
; if we are running within gen-ctl-io, as it won't work.
(if (not (defined? 'output-source)) ; (a function defined by gen-ctl-io)
    (begin
      (set-param! print-ok? (zero? (meep-my-rank)))
      (set-param! interactive? (= 1 (meep-count-processors)))))

; ****************************************************************

(define-input-var epsilon-input-file "" 'string)

(define-class material-type no-parent)

(define-class susceptibility no-parent
  (define-property sigma-offdiag (vector3 0 0 0) 'vector3)
  (define-property sigma-diag no-default 'vector3))
(define-class lorentzian-susceptibility susceptibility
  (define-property frequency no-default 'number)
  (define-property gamma no-default 'number))
(define-class drude-susceptibility susceptibility
  (define-property frequency no-default 'number)
  (define-property gamma no-default 'number))
(define-class noisy-lorentzian-susceptibility lorentzian-susceptibility
  (define-property noise-amp no-default 'number))
(define-class noisy-drude-susceptibility drude-susceptibility
  (define-property noise-amp no-default 'number))
(define-class gyrotropic-lorentzian-susceptibility lorentzian-susceptibility
  (define-property bias no-default 'vector3))
(define-class gyrotropic-drude-susceptibility drude-susceptibility
  (define-property bias no-default 'vector3))

(define-class gyrotropic-saturated-susceptibility susceptibility
  (define-property frequency no-default 'number)
  (define-property gamma no-default 'number)
  (define-property bias no-default 'vector3)
  (define-property alpha no-default 'number))

(define polarizability lorentzian-susceptibility) ; backwards compat
(define omega frequency) ; backwards compat
(define (sigma x) (sigma-diag x x x))

(define-class medium material-type
  (define-property epsilon-diag (vector3 1 1 1) 'vector3)
  (define-property epsilon-offdiag (vector3 0 0 0) 'vector3)
  (define-property mu-diag (vector3 1 1 1) 'vector3)
  (define-property mu-offdiag (vector3 0 0 0) 'vector3)
  (define-property E-susceptibilities '() (make-list-type 'susceptibility))
  (define-property H-susceptibilities '() (make-list-type 'susceptibility))
  (define-property E-chi2-diag (vector3 0 0 0) 'vector3)
  (define-property E-chi3-diag (vector3 0 0 0) 'vector3)
  (define-property H-chi2-diag (vector3 0 0 0) 'vector3)
  (define-property H-chi3-diag (vector3 0 0 0) 'vector3)
  (define-property D-conductivity-diag (vector3 0 0 0) 'vector3)
  (define-property B-conductivity-diag (vector3 0 0 0) 'vector3)
)

; backwards compatibility:
(define E-polarizations E-susceptibilities)
(define H-polarizations H-susceptibilities)

(define (epsilon eps) (epsilon-diag eps eps eps))
(define (mu m) (mu-diag m m m))

(define dielectric medium) ; old name for backwards compatibility
(define polarizations E-polarizations) ; backwards compatibility

; useful shortcuts for isotropic conductivity
(define (D-conductivity c) (D-conductivity-diag c c c))
(define (B-conductivity c) (B-conductivity-diag c c c))

; shortcuts for isotropic nonlinearities
(define (E-chi2 d) (E-chi2-diag d d d))
(define (E-chi3 d) (E-chi3-diag d d d))
(define (H-chi2 d) (H-chi2-diag d d d))
(define (H-chi3 d) (H-chi3-diag d d d))
(define chi2 E-chi2) (define chi3 E-chi3) ; backwards compatibility

(define-class perfect-metal material-type)

; arbitrary material(x)
(define-class material-function material-type
  (define-property material-func no-default 'function
    (lambda (p) (= 1 (procedure-num-args p)))))
(define (epsilon-func f) ; convenience wrapper
  (material-func (lambda (p) (make dielectric (epsilon (f p))))))

(define (index n) (epsilon (* n n))) ; convenient substitute for epsilon

; use the solid geometry classes, variables, etcetera in libgeom:
; (one specifications file can include another specifications file)
(include "@LIBCTL_DIR@/utils/geom.scm")

; Include materials library
(include "materials.scm")

(if (defined? 'meep-component-Ex) (include "meep-enums.scm"))

(define CYLINDRICAL -2) ; special value of dimensions for cylindrical coords
(define AUTOMATIC -1) ; special value for directions, when auto-determined

; ****************************************************************
; Multilevel-atom nonlinear susceptibilities

(define-class transition no-parent
  (define-property from-level no-default 'integer non-negative?)
  (define-property to-level no-default 'integer non-negative?)
  (define-property transition-rate 0 'number) ; nonradiative rate (0 if none)
  (define-property frequency 0 'number) ; radiative frequency (0 if none)
  (define-property sigma-diag (vector3 1 1 1) 'vector3) ; per-transition sigma
  (define-property gamma 0 'number) ; optical damping rate
  (define-property pumping-rate 0 'number)) ; pumping rate (0 if none)

(define (transition-time t) (transition-rate (/ t)))

(define-class multilevel-atom susceptibility
  (define-property initial-populations '() (make-list-type 'number))
  (define-property transitions '() (make-list-type 'transition)))

; ****************************************************************
; Add some predefined variables, for convenience:

(define vacuum (make dielectric (epsilon 1.0)))
(define air vacuum)
(define metal (make perfect-metal))
(define perfect-electric-conductor metal)
(define perfect-magnetic-conductor (make medium (mu (/ -0.0)))) ; -infinity

(define infinity 1.0e20) ; big number for infinite dimensions of objects

(set! default-material vacuum)

(define pi (* 4.0 (atan 1.0)))

; ****************************************************************
; some utilities

; round x to dig digits after the decimal place
(define (round-dig dig x) (/ (round (* x (expt 10 dig))) (expt 10 dig)))

; display a comma-delimited list of values, prefixed by the data-name
; and the run index, with data a list of values.
(define-param run-index 0)
(define (display-run-data data-name data)
  (print data-name run-index ":")
  (map (lambda (v) (print ", " v)) data)
  (print "\n"))

; display a list of data arrays (all required to be of same length)
; as comma-separated values, prefixed with name (and the run index).
(define (display-csv name . ds)
  (if (not (null? ds))
      (apply
       map
       (cons
	(lambda (. vs) (display-run-data name vs))
	ds))))

; ****************************************************************
; More input variables

; list of extra materials not explicit in the geometry, e.g. those
; used in material-functions, so that Meep doesn't miss them
(define-param extra-materials '())

(define structure '())

; list of sources added by init-fields (below)
(define-param sources '())

(define-param m 0) ; angular dependence exp(i m phi) in cylindrical

; If false (default), Meep forces certain field components for |m|>1 to
; be zero for |m| pixels from r=0.  This is something of a hack
; which seems to ensure stability for Courant = 0.5 regardless of m,
; at the expense of some accuracy near r=0.  If true, in order to
; remain stable, the Courant factor should be ~ min[0.5, 1 / (|m| + 0.5)]
;  or so assuming the computational cell includes the r=0 origin.
(define-param accurate-fields-near-cylorigin? false)

(define-param force-complex-fields? false)

(define-param k-point false)

; whether to use the "beta" trick to handle kz in 2d
(define-param kz-2d "complex")

(define fields '())

(define-param subpixel-tol 1e-4)
(define-param subpixel-maxeval 100000)

; a "global" conductivity to add to all materials, mostly
; for the convenience of Casimir calculations
(define-param global-D-conductivity 0)
(define-param global-B-conductivity 0)

; ****************************************************************
; Setting up the structure

(define-class symmetry no-parent
  (define-property direction no-default 'integer)
  (define-property phase 1.0 'cnumber))
(define-class rotate2-sym symmetry)
(define-class rotate4-sym symmetry)
(define-class mirror-sym symmetry)

(define ALL -1) ; special value for directions, when all values
(define-class pml no-parent
  (define-property thickness no-default 'number)
  (define-property direction ALL 'integer)
  (define-property side ALL 'integer)
  (define-property strength 1.0 'number) ; obsolete: R -> R^strength
  (define-property R-asymptotic 1e-15 'number positive?)
  (define-property mean-stretch 1.0 'number (lambda (x) (>= x 1.0)))
  (define-property pml-profile (lambda (u) (* u u)) 'function))

(define-class absorber pml) ; just a scalar conductivity ramp

(export-type (make-list-type 'symmetry))
(export-type (make-list-type 'pml))
(export-type (make-list-type 'material-type))

(define-param symmetries '())
(define-param pml-layers '())
(define-param num-chunks 0)
(define-param Courant 0.5)

(define (infer-dimensions k)
  (if (and (not (null? k)) (= dimensions 3)) ; infer dimensions
      (if (and (= (vector3-z (object-property-value geometry-lattice 'size))
		  no-size)
	       (or (not k) (not (equal? kz-2d "3d")) (zero? (vector3-z k))))
	  2
	  3)
      dimensions))
(define (require-dimensions!)
  (if (null? structure) (set-dimensions (infer-dimensions k-point))))

(define-class volume-class no-parent
  (define-property center no-default 'vector3)
  (define-property size (vector3 0 0 0) 'vector3))

(define (volume . args)
  (require-dimensions!)
  (let ((v (apply make (cons volume-class args))))
    (let ((cen (object-property-value v 'center))
	  (sz (object-property-value v 'size)))
      (new-meep-volume
       (vector3- cen (vector3-scale 0.5 sz))
       (vector3+ cen (vector3-scale 0.5 sz))))))

(define-param eps-averaging? true) ; 10% slower, but huge accuracy gains

(define (init-structure . k_)
  (let ((k (if (null? k_) '() (car k_)))
	(s (object-property-value geometry-lattice 'size)))
    (set! structure (make-structure
		     (infer-dimensions k)
		     s geometry-center
		     resolution
		     eps-averaging? subpixel-tol subpixel-maxeval
		     (and ensure-periodicity (not (not k)))
		     geometry extra-materials
		     default-material epsilon-input-file
		     pml-layers symmetries num-chunks Courant
		     global-D-conductivity global-B-conductivity))))

; ****************************************************************
; Adding sources

(define-class src-time no-parent
  (define-property is-integrated? false 'boolean))
(define-class continuous-src src-time
  (define-property frequency no-default 'number)
  (define-property start-time 0 'number)
  (define-property end-time infinity 'number)
  (define-property width 0 'number)
  (define-property cutoff 3.0 'number)
  (define-derived-property swigval 'SCM
    (lambda (o)
      (let ((s (new-meep-continuous-src-time
		(* 1.0 (object-property-value o 'frequency))
		(object-property-value o 'width)
		(object-property-value o 'start-time)
		(object-property-value o 'end-time)
		(object-property-value o 'cutoff))))
	(meep-src-time-is-integrated-set s (object-property-value o 'is-integrated?))
	s))))
(define-class gaussian-src src-time
  (define-property frequency no-default 'number)
  (define-property width no-default 'number)
  (define-property start-time 0 'number)
  (define-property cutoff 5.0 'number)
  (define-derived-property swigval 'SCM
    (lambda (o)
      (let ((s (new-meep-gaussian-src-time
		(* 1.0 (object-property-value o 'frequency))
		(object-property-value o 'width)
		(object-property-value o 'start-time)
		(+ (object-property-value o 'start-time)
		   (* 2 (object-property-value o 'width)
		      (object-property-value o 'cutoff))))))
	(meep-src-time-is-integrated-set s (object-property-value o 'is-integrated?))
	s))))
(define-class custom-src src-time
  (define-property src-func no-default 'function)
  (define-property start-time (- infinity) 'number)
  (define-property end-time infinity 'number)
  (define-derived-property swigval 'SCM
    (lambda (o)
      (let ((s (new-meep-custom-src-time
		(object-property-value o 'src-func)
		(object-property-value o 'start-time)
		(object-property-value o 'end-time))))
	(meep-src-time-is-integrated-set s (object-property-value o 'is-integrated?))
	s))))

(define (fwidth df) (width (/ df))) ; to specify frequency width instead
(define (wavelength lam) (frequency (/ lam)))
(define (period T) (frequency (/ T)))

(define-class source no-parent
  (define-property src no-default 'src-time)
  (define-property component no-default 'integer)
  (define-property center no-default 'vector3)
  (define-property size (vector3 0 0 0) 'vector3)
  (define-property amplitude 1.0 'cnumber)
  (define-property amp-func '() 'SCM))

; the following definitions are taken from MPB
(define NO-PARITY 0)
(define EVEN-Z 1)
(define ODD-Z 2)
(define EVEN-Y 4)
(define ODD-Y 8)
(define TE EVEN-Z)
(define TM ODD-Z)

; special component for eigenmode-source (not defined during gen-ctl-io)
(define ALL-COMPONENTS (if (defined? 'Dielectric) Dielectric 99))

(define-class eigenmode-source source
  ; the following two properties have special default values -- if we
  ; detect these values, we will replace the values by those of the
  ; size and center of the parent class (a bit hackish, but...)
  (define-property eig-lattice-size (vector3 -1 -1 -1) 'vector3)
  (define-property eig-lattice-center (vector3 infinity infinity infinity) 'vector3)

  (define-property component ALL-COMPONENTS 'integer) ; new default val
  (define-property direction AUTOMATIC 'integer)
  (define-property eig-band 1 'integer positive?)
  (define-property eig-kpoint (vector3 0) 'vector3)
  (define-property eig-match-freq? true 'boolean)
  (define-property eig-parity NO-PARITY 'integer)
  (define-property eig-resolution 0 'integer non-negative?)
  (define-property eig-tolerance 1e-12 'number positive?))

(define (add-source s f) ; add source s to fields f
  (define (dflt v v0 vd) (if (vector3= v v0) vd v))
  (let ((A (object-property-value s 'amp-func))
	(cen (object-property-value s 'center))
	(sz (object-property-value s 'size)))
    (if (object-member? 'eigenmode-source s)
	(let* ((ecen (dflt (object-property-value s 'eig-lattice-center)
			   (vector3 infinity infinity infinity) cen))
	       (esz (dflt (object-property-value s 'eig-lattice-size)
			  (vector3 -1 -1 -1) sz))
	       (v (volume (center cen) (size sz)))
	       (d0 (object-property-value s 'direction))
	       (d (if (negative? d0)
		      (meep-fields-normal-direction fields v)
		      d0)))
	  (if (null? A)
	      (meep-fields-add-eigenmode-source f
	       (object-property-value s 'component)
	       (object-property-value (object-property-value s 'src) 'swigval)
	       d v
	       (volume (center ecen) (size esz))
	       (object-property-value s 'eig-band)
	       (object-property-value s 'eig-kpoint)
	       (object-property-value s 'eig-match-freq?)
	       (object-property-value s 'eig-parity)
	       (object-property-value s 'eig-resolution)
	       (object-property-value s 'eig-tolerance)
	       (object-property-value s 'amplitude))
	      (meep-fields-add-eigenmode-source f
	       (object-property-value s 'component)
	       (object-property-value (object-property-value s 'src) 'swigval)
	       d v
	       (volume (center ecen) (size esz))
	       (object-property-value s 'eig-band)
	       (object-property-value s 'eig-kpoint)
	       (object-property-value s 'eig-match-freq?)
	       (object-property-value s 'eig-parity)
	       (object-property-value s 'eig-resolution)
	       (object-property-value s 'eig-tolerance)
	       (object-property-value s 'amplitude) A)))
	(if (null? A)
	    (meep-fields-add-volume-source f
	     (object-property-value s 'component)
	     (object-property-value (object-property-value s 'src) 'swigval)
	     (volume (center cen) (size sz))
	     (* 1.0 (object-property-value s 'amplitude)))
	    (meep-fields-add-volume-source f
	     (object-property-value s 'component)
	     (object-property-value (object-property-value s 'src) 'swigval)
	     (volume (center cen) (size sz))
	     A (* 1.0 (object-property-value s 'amplitude)))))))

; ****************************************************************
; Setting up the fields

(define init-fields-hooks '()) ; list of thunks to execute after init-fields
(define (init-fields)
  (if (null? structure) (init-structure k-point))
  (let ((kz2d (if (= (infer-dimensions k-point) 2) kz-2d "3d")))
    (set! fields (new-meep-fields structure
				  (if (= dimensions CYLINDRICAL) m 0)
				  (if (and (not (equal? kz2d "3d")) k-point)
				      (vector3-z k-point) 0.0)
				  (not accurate-fields-near-cylorigin?)))
    (if (not (or force-complex-fields?
	         (and (= dimensions CYLINDRICAL) (not (zero? m)))
	         (not (for-all? symmetries
			        (lambda (s)
				  (zero? (imag-part
					  (object-property-value s 'phase))))))
	         (not (or (not k-point)
			  (and (equal? kz2d "real/imag")
			       (= (vector3-x k-point) 0)
			       (= (vector3-y k-point) 0))
			  (vector3= k-point (vector3 0))))))
        (meep-fields-use-real-fields fields)
        (print "Meep: using complex fields.\n"))
    (if k-point (meep-fields-use-bloch fields
				       (if (not (equal? kz2d "3d"))
					   (vector3 (vector3-x k-point)
						    (vector3-y k-point))
					   k-point))))
  (map (lambda (s) (add-source s fields)) sources)
  (map (lambda (thunk) (thunk)) init-fields-hooks))

(define (meep-time)
  (if (null? fields) (init-fields))
  (meep-fields-time fields))

(define (meep-round-time)
  (if (null? fields) (init-fields))
  (meep-fields-round-time fields))

(define (get-field-point c pt)
  (meep-fields-get-field fields c pt))

(define (get-epsilon-point pt)
  (meep-fields-get-eps fields pt))

; ****************************************************************
; Various ways to restart all or part of the simulation.

(define (change-k-point! k)
  (set! k-point k)
  (if (not (null? fields))
      (if (and (not (or (not k-point) (vector3= k-point (vector3 0))))
	       (not (zero? (meep-fields-is-real-get fields))))
	  (begin (delete-meep-fields fields)
		 (set! fields '())
		 (init-fields))
	  (if k-point (meep-fields-use-bloch fields k-point)))))

(define (change-sources! new-sources)
  (set! sources new-sources)
  (if (not (null? fields))
      (begin
	(meep-fields-remove-sources fields)
	(map (lambda (s) (add-source s fields)) sources))))

(define (reset-meep)
  (delete-meep-fields fields) (set! fields '())
  (delete-meep-structure structure) (set! structure '()))

(define (restart-fields)
  (if (not (null? fields))
      (begin
	(meep-fields-t-set fields 0)
	(meep-fields-zero-fields fields))
      (init-fields)))

; ****************************************************************
; Flux spectra

(define-class flux-region no-parent
  (define-property center no-default 'vector3)
  (define-property size (vector3 0 0 0) 'vector3)
  (define-property direction AUTOMATIC 'integer)
  (define-property weight 1.0 'cnumber))

(define (fields-add-fluxish-stuff add-dft-stuff fields fcen df nfreq stufflist . nperiods)
  (define vl '()) ; volume_list of flux regions
  (map (lambda (f)
	 (let* ((v (volume (center (object-property-value f 'center))
			    (size (object-property-value f 'size))))
		(d0 (object-property-value f 'direction))
		(d (if (negative? d0)
		       (meep-fields-normal-direction fields v)
		       d0))
		(c (meep-direction-component Sx d)))
	   (set! vl
		 (make-volume-list
		  (volume (center (object-property-value f 'center))
			  (size (object-property-value f 'size)))
		  c
		  (object-property-value f 'weight)
		  vl))))
       stufflist)
  (let ((stuff
         (if (eq? add-dft-stuff meep-fields-add-dft-near2far)
             (add-dft-stuff fields vl
                            (- fcen (/ df 2)) (+ fcen (/ df 2)) nfreq
                            (if (not (null? nperiods)) (car nperiods)))
             (add-dft-stuff fields vl
                            (- fcen (/ df 2)) (+ fcen (/ df 2)) nfreq))))
    (delete-meep-volume-list vl)
    stuff))

(define (fields-add-flux fields fcen df nfreq . fluxes)
  (fields-add-fluxish-stuff meep-fields-add-dft-flux
			    fields fcen df nfreq fluxes))

(define (add-flux fcen df nfreq . fluxes)
  (if (null? fields) (init-fields))
  (apply fields-add-flux (append (list fields fcen df nfreq) fluxes)))

(define (scale-flux-fields s f)
  (meep-dft-flux-scale-dfts f s))

(define (get-flux-freqs f)
  (std-vector-double-to-scm (meep-dft-flux-freq-get f)))

(export-type (make-list-type 'number))
(define (get-fluxes f)
  (dft-flux-flux f))

(define (display-fluxes . fluxes)
  (if (not (null? fluxes))
      (apply display-csv
	     (append (list "flux"
			   (get-flux-freqs (car fluxes)))
		     (map get-fluxes fluxes)))))

(define (load-flux fname flux)
  (if (null? fields) (init-fields))
  (meep-dft-flux-load-hdf5 flux fields fname "" (get-filename-prefix)))

(define (save-flux fname flux)
  (if (null? fields) (init-fields))
  (meep-dft-flux-save-hdf5 flux fields fname "" (get-filename-prefix)))

(define (load-minus-flux fname flux)
  (load-flux fname flux)
  (meep-dft-flux-scale-dfts flux -1.0))

; ****************************************************************
; Mode monitor

(define mode-region flux-region)

(define (fields-add-mode-monitor fields fcen df nfreq . fluxes)
  (if (not (= (length fluxes) 1)) (error "add-mode-monitor expected just one mode-region."))
  (let* ((region (car fluxes))
         (v (volume (center (object-property-value region 'center))
                    (size (object-property-value region 'size))))
         (d0 (object-property-value region 'direction))
         (d (if (< d0 0)
                (meep-fields-normal-direction fields v)
                d0)))
    (meep-fields-add-mode-monitor fields d v (- fcen (/ df 2)) (+ fcen (/ df 2)) nfreq)))

(define (add-mode-monitor fcen df nfreq . fluxes)
  (if (null? fields) (init-fields))
  (apply fields-add-mode-monitor (append (list fields fcen df nfreq) fluxes)))

; ****************************************************************
; Energy spectra

(define-class energy-region no-parent
  (define-property center no-default 'vector3)
  (define-property size (vector3 0 0 0) 'vector3)
  (define-property direction 0 'integer)
  (define-property weight 1.0 'cnumber))

(define (fields-add-energy fields fcen df nfreq . energys)
  (fields-add-fluxish-stuff meep-fields-add-dft-energy
			    fields fcen df nfreq energys))

(define (add-energy fcen df nfreq . energys)
  (if (null? fields) (init-fields))
  (apply fields-add-energy (append (list fields fcen df nfreq) energys)))

(define (scale-energy-fields s f)
  (meep-dft-energy-scale-dfts f s))

(define (get-energy-freqs f)
  (std-vector-double-to-scm (meep-dft-energy-freq-get f)))

(define (get-electric-energy f)
  (dft-energy-electric f))

(define (get-magnetic-energy f)
  (dft-energy-magnetic f))

(define (get-total-energy f)
  (dft-energy-total f))

(define (display-electric-energy . energys)
  (if (not (null? energys))
      (apply display-csv
	     (append (list "electric-energy"
			   (get-energy-freqs (car energys)))
		     (map get-electric-energy energys)))))

(define (display-magnetic-energy . energys)
  (if (not (null? energys))
      (apply display-csv
	     (append (list "magnetic-energy"
			   (get-energy-freqs (car energys)))
		     (map get-magnetic-energy energys)))))

(define (display-total-energy . energys)
  (if (not (null? energys))
      (apply display-csv
	     (append (list "total-energy"
			   (get-energy-freqs (car energys)))
		     (map get-total-energy energys)))))

(define (load-energy fname energy)
  (if (null? fields) (init-fields))
  (meep-dft-energy-load-hdf5 energy fields fname "" (get-filename-prefix)))

(define (save-energy fname energy)
  (if (null? fields) (init-fields))
  (meep-dft-energy-save-hdf5 energy fields fname "" (get-filename-prefix)))

(define (load-minus-energy fname energy)
  (load-energy fname energy)
  (meep-dft-energy-scale-dfts energy -1.0))

; ****************************************************************
; Force spectra (from stress tensor) - very similar interface to flux spectra

(define-class force-region no-parent
  (define-property center no-default 'vector3)
  (define-property size (vector3 0 0 0) 'vector3)
  (define-property direction no-default 'integer)
  (define-property weight 1.0 'cnumber))

(define (fields-add-force fields fcen df nfreq . forcees)
  (fields-add-fluxish-stuff meep-fields-add-dft-force
			    fields fcen df nfreq forcees))

(define (add-force fcen df nfreq . forcees)
  (if (null? fields) (init-fields))
  (apply fields-add-force (append (list fields fcen df nfreq) forcees)))

(define (scale-force-fields s f)
  (meep-dft-force-scale-dfts f s))

(define (get-force-freqs f)
  (std-vector-double-to-scm (meep-dft-force-freq-get f)))

(define (get-forces f)
  (dft-force-force f))

(define (display-forces . forcees)
  (if (not (null? forcees))
      (apply display-csv
	     (append (list "force"
			   (get-force-freqs (car forcees)))
		     (map get-forces forcees)))))

(define (load-force fname force)
  (if (null? fields) (init-fields))
  (meep-dft-force-load-hdf5 force fields fname "" (get-filename-prefix)))

(define (save-force fname force)
  (if (null? fields) (init-fields))
  (meep-dft-force-save-hdf5 force fields fname "" (get-filename-prefix)))

(define (load-minus-force fname force)
  (load-force fname force)
  (meep-dft-force-scale-dfts force -1.0))

; ****************************************************************
; Near-to-far-field transformations (again similar to dft-foobar)

(define-class near2far-region no-parent
  (define-property center no-default 'vector3)
  (define-property size (vector3 0 0 0) 'vector3)
  (define-property direction AUTOMATIC 'integer)
  (define-property weight 1.0 'cnumber))

(define (fields-add-near2far fields fcen df nfreq nperiods . near2fars)
  (fields-add-fluxish-stuff meep-fields-add-dft-near2far
			    fields fcen df nfreq near2fars nperiods))

(define* (add-near2far fcen df nfreq #:key (nperiods 1) . near2fars)
  (if (null? fields) (init-fields))
  ;; Keyword arguments automatically get added to the "rest" list in guile,
  ;; so we have to remove the nperiods key and value from near2fars
  (apply fields-add-near2far (append (list fields fcen df nfreq nperiods)
                                     (delq #:nperiods (delq nperiods near2fars)))))

(define (scale-near2far-fields s f)
  (meep-dft-near2far-scale-dfts f s))

(define (get-near2far-freqs f)
  (std-vector-double-to-scm (meep-dft-near2far-freq-get f)))

(define (get-farfield f x)
  (dft-near2far-farfield f x))

(define (output-farfields near2far fname where resolution)
  (meep-dft-near2far-save-farfields near2far fname (get-filename-prefix)
                                    where resolution))

(define (flux near2far df where resolution)
  (dft-near2far-flux near2far df where resolution))

(define (load-near2far fname near2far)
  (if (null? fields) (init-fields))
  (meep-dft-near2far-load-hdf5 near2far fields fname "" (get-filename-prefix)))

(define (save-near2far fname near2far)
  (if (null? fields) (init-fields))
  (meep-dft-near2far-save-hdf5 near2far fields fname "" (get-filename-prefix)))

(define (load-minus-near2far fname near2far)
  (load-near2far fname near2far)
  (meep-dft-near2far-scale-dfts near2far -1.0))

; ****************************************************************
; dft fields

(define* (add-dft-fields components freq_min freq_max nfreq #:key (yee-grid false) . where)
  (if (null? fields) (init-fields))
  (let ((num-components (length components))
        (where-only (delq #:yee-grid (delq yee-grid where))))
    (let ((vol (if (null? where-only) (meep-fields-total-volume fields) (car where-only))))
      (meep-fields-add-dft-fields fields components vol freq_min freq_max nfreq (not yee-grid)))))

(define (output-dft dft_fields fname)
  (if (null? fields) (init-fields))
  (meep-fields-output-dft fields dft_fields fname))

; ****************************************************************
; Generic step functions: these are functions which are called
; (potentially) at every time step.  They can either be a thunk
; or they can take one argument, to-do.  to-do is either 'step
; or 'finish, where 'step means to output (or whatever)
; normally, and 'finish is passed once at the end of the run
; (and is used to close files, print summary output, etcetera).

; step functions can be either thunks (the common case), or
; can take a "to-do" argument that is currently either 'step
; or 'finish (so that they can clean up at the end of a run).
(define (eval-step-func func to-do)
  (if (= 0 (procedure-num-args func))
      (if (eq? to-do 'step) (func))
      (func to-do)))

; Some convenient wrappers for step functions passed to run.  e.g., these
; can be used to only output at certain times, instead of ata every time step.

(define (combine-step-funcs . step-funcs)
  (lambda (to-do)
    (map (lambda (f) (eval-step-func f to-do)) step-funcs)))

; generic wrapper
(define (when-true-funcs cond? step-funcs)
  (lambda (to-do)
    (if (or (eq? to-do 'finish) (cond?))
	(map (lambda (f) (eval-step-func f to-do)) step-funcs))))

; evaluate step-funcs whenever (cond?) is true/false.
(define (when-true cond? . step-funcs) (when-true-funcs cond? step-funcs))
(define (when-false cond? . step-funcs)
  (when-true-funcs (lambda () (not (cond?))) step-funcs))

; output at an interval of dT (in meep/simulation time).
(define (at-every dT . step-funcs)
  (if (null? fields) (init-fields))
  (let ((Tlast (meep-round-time)))
    (lambda (to-do)
      (let ((T (meep-round-time)))
	(if (or (eq? to-do 'finish)
		(>= T (+ Tlast dT (* -0.5 (meep-fields-dt-get fields)))))
	    (begin
	      (map (lambda (f) (eval-step-func f to-do)) step-funcs)
	      (set! Tlast T)))))))

(define (after-time T . step-funcs)
  (if (null? fields) (init-fields))
  (let ((T0 (meep-round-time)))
    (when-true-funcs (lambda () (>= (meep-round-time) (+ T0 T))) step-funcs)))

(define (before-time T . step-funcs)
  (if (null? fields) (init-fields))
  (let ((T0 (meep-round-time)))
    (when-true-funcs (lambda () (< (meep-round-time) (+ T0 T))) step-funcs)))

(define (at-time T . step-funcs)
  (let ((done? false))
    (after-time T (lambda (to-do)
		    (if (or (not done?) (eq? to-do 'finish))
			(map (lambda (f) (eval-step-func f to-do)) step-funcs))
		    (set! done? (or done? (eq? to-do 'step)))))))

(define (after-sources . step-funcs)
  (if (null? fields) (init-fields))
  (apply after-time
	 (cons (- (meep-fields-last-source-time fields) (meep-round-time))
	       step-funcs)))

; after sources plus a time T.
(define (after-sources+ T . step-funcs)
  (if (null? fields) (init-fields))
  (apply after-time
	 (cons (- (+ (meep-fields-last-source-time fields) T) (meep-round-time))
	       step-funcs)))

(define (during-sources . step-funcs)
  (if (null? fields) (init-fields))
  (apply before-time
	 (cons (- (meep-fields-last-source-time fields) (meep-round-time))
	       step-funcs)))

; the user could just call functions, but this functions saves the user
; from having to manually call init-fields
(define (at-beginning . step-funcs)
  (let ((done? false))
    (lambda (to-do)
      (if (not done?)
	  (begin
	    (map (lambda (f) (eval-step-func f to-do)) step-funcs)
	    (set! done? true))))))

; for completeness (although the user could just do this after running):
(define (at-end . step-funcs)
  (lambda (to-do)
    (if (eq? to-do 'finish)
	(begin
	  (map (lambda (f) (eval-step-func f 'step)) step-funcs)
	  (map (lambda (f) (eval-step-func f 'finish)) step-funcs)))))

; run the step-funcs with the magnetic fields synchronized in time
; with the electric fields
(define (synchronized-magnetic . step-funcs)
  (lambda (to-do)
    (meep-fields-synchronize-magnetic-fields fields)
    (map (lambda (f) (eval-step-func f to-do)) step-funcs)
    (meep-fields-restore-magnetic-fields fields)))

; ****************************************************************
; File output functions (can only be called after init-fields).

(define-param filename-prefix "")
(define (get-filename-prefix)
  (if (eq? filename-prefix false)
      ""
      (if (and (not (null? include-files))
	       (string-null? filename-prefix))
	  (string-append
	   (strip-suffix ".scm"
	     (strip-suffix ".ctl" (cdr (split-pathname (car include-files))))))
	  filename-prefix)))

; Use output directory instead of outputting in same directory;
; uses init-fields-hooks to handle fields not yet initted.
(define (use-output-directory . dname_)
  (let ((dname (if (null? dname_)
		   (string-append (get-filename-prefix) "-out")
		   (car dname_))))
    (let ((hook
	   (let ((trashed? false)) ; only trash output directory once per run
	     (lambda ()
	       (print "Meep: using output directory \"" dname "\"\n")
	       (meep-fields-set-output-directory fields dname)
	       (if (not trashed?) (meep-trash-output-directory dname))
	       (set! trashed? true)))))
      (set! init-fields-hooks (cons hook init-fields-hooks))
      (if (not (null? fields)) (hook))
      (set! filename-prefix false)
      dname)))

(define-param output-volume '()) ; region to output; NULL for everywhere

(define output-append-h5 '()) ; h5 file to append data to (NULL if none)

; hook function called with the filename after every HDF5 files is created;
; this can be used to convert the file into other formats, etcetera.
(define output-h5-hook (lambda (fname) false)) ; default is no-op

(define output-single-precision? false) ; output single-prec to save space

(define meep-last-eps-filename "") ; most recent epsilon file outputted

(define (output-component c . h5file)
  (if (null? fields) (error "init-fields is required before output-component"))
  (meep-fields-output-hdf5
   fields c
   (if (null? output-volume)
       (meep-fields-total-volume fields)
       output-volume)
   (if (null? h5file) output-append-h5 (car h5file))
   (and (null? h5file) (not (null? output-append-h5)))
   output-single-precision?
   (get-filename-prefix))
  (if (null? h5file)
      (let ((nm (meep-fields-h5file-name fields (meep-component-name c)
					 (get-filename-prefix) true)))
	(if (eq? c Dielectric) (set! meep-last-eps-filename nm))
	(output-h5-hook nm))))

; cs = list of components, and func is function of position & these components
(define (output-field-function-helper name cs func real-only? h5file)
  (if (null? fields) (error "init-fields is required before output-field-function"))
  (meep-fields-output-hdf5
   fields name
   (cons cs func)
   (if (null? output-volume)
       (meep-fields-total-volume fields)
       output-volume)
   (if (null? h5file) output-append-h5 (car h5file))
   (and (null? h5file) (not (null? output-append-h5)))
   output-single-precision?
   (get-filename-prefix)
   real-only?)
  (if (null? h5file)
      (output-h5-hook (meep-fields-h5file-name fields
                       name (get-filename-prefix) true))))

(define (output-field-function name cs func . h5file)
  (output-field-function-helper name cs func false h5file))

(define (output-real-field-function name cs func . h5file)
  (output-field-function-helper name cs func true h5file))

(define (output-components fname . cs)
  (if (null? fields) (error "init-fields is required before output-component"))
  (let ((f
	 (if (null? output-append-h5)
	     (list (meep-fields-open-h5file fields fname (meep-h5file-WRITE)
					    (get-filename-prefix) true))
	     '())))
    (map (lambda (c)
	   (apply output-component (cons c f))
	   (if (null? output-append-h5)
	       (meep-h5file-prevent-deadlock (car f))))
	 cs)
    (if (null? output-append-h5) (delete-meep-h5file (car f))))
  (if (null? output-append-h5)
      (output-h5-hook (meep-fields-h5file-name
		       fields fname (get-filename-prefix) true))))

; convenience functions, similar to MPB:
(define (output-epsilon) (output-component Dielectric))
(define (output-mu) (output-component Permeability))
(define (output-hpwr) (output-component H-EnergyDensity))
(define (output-dpwr) (output-component D-EnergyDensity))
(define (output-tot-pwr) (output-component EnergyDensity))
(defmacro-public define-output-field (name cp CP)
  `(begin
     (define (,(symbol-append 'output- cp))
       (output-components ,name
			  ,(symbol-append CP 'x)
			  ,(symbol-append CP 'y)
			  ,(symbol-append CP 'z)
			  ,(symbol-append CP 'r) ,(symbol-append CP 'p)))
     (define (,(symbol-append 'output- cp '-x))
       (output-component ,(symbol-append CP 'x)))
     (define (,(symbol-append 'output- cp '-y))
       (output-component ,(symbol-append CP 'y)))
     (define (,(symbol-append 'output- cp '-z))
       (output-component ,(symbol-append CP 'z)))
     (define (,(symbol-append 'output- cp '-r))
       (output-component ,(symbol-append CP 'r)))
     (define (,(symbol-append 'output- cp '-p))
       (output-component ,(symbol-append CP 'p)))))
(define-output-field "h" hfield H)
(define-output-field "b" bfield B)
(define-output-field "e" efield E)
(define-output-field "d" dfield D)
(define-output-field "s" poynting S) ; compat. with MPB.
(define-output-field "s" sfield S)

(define (with-prefix pre . step-funcs)
  (lambda (to-do)
    (let ((pre-save filename-prefix))
      (set! filename-prefix (string-append pre (get-filename-prefix)))
      (map (lambda (f) (eval-step-func f to-do)) step-funcs)
      (set! filename-prefix pre-save))))

; change output-volume for a few step-funcs to v
(define (in-volume v . step-funcs)
  (let ((cur-eps "")) ; allow per-volume eps filenames
    (lambda (to-do)
      (let ((v-save output-volume)
	    (eps-save meep-last-eps-filename))
	(set! output-volume v)
	(if (not (string-null? cur-eps))
	    (set! meep-last-eps-filename cur-eps))
	(map (lambda (f) (eval-step-func f to-do)) step-funcs)
	(set! cur-eps meep-last-eps-filename)
	(set! output-volume v-save)
	(if (not (string-null? eps-save))
	    (set! meep-last-eps-filename eps-save))))))
(define (in-point pt . step-funcs)
  (apply in-volume (cons (volume (center pt)) step-funcs)))

; Meep supports outputting d+1 dimensional HDF5 files where the last
; dimension is time.
(define (to-appended fname . step-funcs)
  (if (null? fields) (init-fields))
  (let ((h5 (meep-fields-open-h5file fields fname (meep-h5file-WRITE)
				     (get-filename-prefix))))
    (lambda (to-do)
      (let ((h5save output-append-h5))
	(set! output-append-h5 h5)
	(map (lambda (f) (eval-step-func f to-do)) step-funcs)
	(if (eq? to-do 'finish)
	    (begin
	      (delete-meep-h5file h5)
	      (output-h5-hook (meep-fields-h5file-name
			       fields fname (get-filename-prefix)))))
	(set! output-append-h5 h5save)))))

(define (convert-h5 rm? convert-cmd . step-funcs)
  (define (convert fname)
    (if (zero? (meep-my-rank))
	(if (and (zero? (system (string-append convert-cmd " \"" fname "\"")))
		 rm?)
	    (system (string-append "rm \"" fname "\"")))))
  (lambda (to-do)
    (let ((hooksave output-h5-hook))
      (set! output-h5-hook convert)
      (map (lambda (f) (eval-step-func f to-do)) step-funcs)
      (set! output-h5-hook hooksave))))

(define (h5topng rm? options . step-funcs)
  (apply convert-h5 (cons rm? (cons (string-append
				     "EPS=\"" meep-last-eps-filename "\"; "
				     "h5topng " options)
				    step-funcs))))

(define (output-png-rm? rm? c options)
  (let ((maxabs 0.0)) ; keep track of amplitude for image scaling
    (lambda (to-do)
      (if (eq? to-do 'step)
	  (begin
	    (set! maxabs
		  (max maxabs
		       (meep-fields-max-abs
			fields c
			(if (null? output-volume)
			    (meep-fields-total-volume fields)
			    output-volume))))
	    ((h5topng rm?
		      (string-append "-M " (number->string maxabs) " " options)
		      (lambda () (output-component c))) to-do))))))

(define (output-png c options) (output-png-rm? true c options))
(define (output-png+h5 c options) (output-png-rm? false c options))

; ****************************************************************
; harminv functions for extracting bands, etcetera

 ; for do-harminv
(export-type (make-list-type 'cnumber))
(export-type (make-list-type 'cvector3))

; generic data-collection function
(defmacro-public collect-harminv! (data data-dt)
  `(lambda (c pt)
     (set! ,data '())
     (let ((t0 0))
       (lambda ()
	 (set! ,data-dt (- (meep-time) t0))
	 (set! t0 (meep-time))
	 (set! ,data (cons (get-field-point c pt) ,data))))))

; do-harminv returns a (freq, amp, err) vector3; define accessor functions:
(define harminv-freq vector3-x)
(define (harminv-freq-re b) (real-part (vector3-x b)))
(define (harminv-freq-im b) (imag-part (vector3-x b)))
(define (harminv-Q b) (/ (harminv-freq-re b)
			 (* -2 (harminv-freq-im b))))
(define harminv-amp vector3-y)
(define harminv-err vector3-z)

(define-param harminv-spectral-density 1.1)
(define-param harminv-Q-thresh 50.0)
(define-param harminv-rel-err-thresh 1e20)
(define-param harminv-err-thresh 0.01)
(define-param harminv-rel-amp-thresh -1.0)
(define-param harminv-amp-thresh -1.0)

(define (analyze-harminv data fcen df maxbands . dt)
  (display-run-data
   "harminv"
   (list "frequency" "imag. freq." "Q" "|amp|" "amplitude" "error"))
  (let ((bands (do-harminv data
			   (if (null? dt)
			       (meep-fields-dt-get fields)
			       (car dt))
			   (- fcen (/ df 2)) (+ fcen (/ df 2)) maxbands
                           harminv-spectral-density harminv-Q-thresh
                           harminv-rel-err-thresh harminv-err-thresh
                           harminv-rel-amp-thresh harminv-amp-thresh)))
    (map (lambda (b) ; b = vector of (freq, amp, error)
	   (display-run-data
	    "harminv"
	    (list
	     (harminv-freq-re b)
	     (harminv-freq-im b)
	     (harminv-Q b)
	     (magnitude (harminv-amp b))
	     (harminv-amp b)
	     (harminv-err b))))
	 bands)
    bands))

(defmacro-public harminv! (data dt results c pt fcen df maxbands)
  `(let ((data' '()) (dt' 0)
	 (c' ,c) (pt' ,pt) (fcen' ,fcen) (df' ,df) (maxbands' ,maxbands))
     (combine-step-funcs
      (at-end
       (lambda ()
	 (set! ,data (reverse data')) ; put in correct order
	 (set! ,dt dt')
	 (set! ,results
	       (analyze-harminv ,data fcen' df'
				(if (list? maxbands')
				    (if (null? maxbands') 100 (car maxbands'))
				    (if (zero? maxbands') 100 maxbands'))
				dt'))))
      ((collect-harminv! data' dt') c' pt'))))

; collect in harminv-data + analyze and store in harminv-results
(define harminv-data '())
(define harminv-data-dt 0)
(define harminv-results '())
(define (harminv c pt fcen df . mxbands)
  (harminv! harminv-data harminv-data-dt harminv-results c pt fcen df mxbands))

; ****************************************************************
; verbosity

(define* (quiet #:optional (quiet? true))
  (verbosity (if quiet? 0 1)))
(if (and (defined? 'verbose) verbose?) (verbosity 2))

; ****************************************************************
; get-eigenmode-coefficients

(define (get-keyword-value args keyword default)
  (let ((kv (memq keyword args)))
    (if (and kv (>= (length kv) 2))
        (cadr kv)
        default)))

(define (get-eigenmode-coefficients flux bands . args)
  (if (null? fields)
      (error "init-fields is required before using get-eigenmode-coefficients"))
  (let ((eig-parity (get-keyword-value args #:eig-parity NO-PARITY))
        (eig-vol (get-keyword-value args #:eig-vol (meep-dft-flux-where-get flux)))
        (eig-resolution (get-keyword-value args #:eig-resolution 0))
        (eig-tolerance (get-keyword-value args #:eig-tolerance 1e-12))
        (direction (get-keyword-value args #:direction AUTOMATIC))
        (num-bands (length bands))
        (kpoint-func (get-keyword-value args #:kpoint-func '())))
    (let ((nfreq (length (get-flux-freqs flux))))
      (define coeffs (make-typed-array 'c64 '0 num-bands nfreq 2))
      (define vgrp (make-typed-array 'f64 '0 num-bands nfreq))
      (define kpoints-and-kdom (do-get-eigenmode-coefficients fields flux eig-vol bands eig-parity coeffs vgrp
                                                             eig-resolution eig-tolerance kpoint-func direction))
      (list coeffs vgrp (first kpoints-and-kdom) (second kpoints-and-kdom)))))

; ****************************************************************
; dft-ldos step function

(define dft-ldos-data '())
(define dft-ldos-Fdata '())
(define dft-ldos-Jdata '())

(define (get-ldos-freqs f)
  (std-vector-double-to-scm (meep-dft-ldos-freq-get f)))

(define (dft-ldos fcen df nfreq)
  (let ((ldos (new-meep-dft-ldos
	       (- fcen (/ df 2)) (+ fcen (/ df 2)) nfreq)))
    (lambda (to-do)
      (if (eq? to-do 'step)
	  (meep-dft-ldos-update ldos fields)
	  (begin
	    (set! dft-ldos-data (dft-ldos-ldos ldos))
	    (set! dft-ldos-Fdata (dft-ldos-F ldos))
	    (set! dft-ldos-Jdata (dft-ldos-J ldos))
	    (display-csv "ldos" (get-ldos-freqs ldos) dft-ldos-data))))))

; ****************************************************************
; run functions

; default time interval (seconds) between progress printouts
(define-param progress-interval 4)

; display progress from T0 until T, every dt seconds (wall time)
(define (display-progress T0 T dt)
  (let ((t0 (meep-wall-time)) (tlast (meep-wall-time)))
    (lambda ()
      (let ((t (meep-wall-time)))
	(if (>= (- t tlast) dt)
	    (begin
	      (print "Meep progress: "
		     (- (meep-time) T0) "/" T " = "
		     (round-dig 1 (/ (- (meep-time) T0) (* 0.01 T))) "% done"
		     " in " (round-dig 1 (- t t0)) "s, "
		     (round-dig 1 ( - (* (- t t0) (/ T (- (meep-time) T0)))
				      (- t t0))) "s to go\n")
	      (set! tlast t)))))))

; run until (cond?) is true or, if cond? is a number, until time cond?
; (in Meep units) has elapsed, calling step-funcs at every time step.
(define (run-until cond? . step-funcs)
  (set! interactive? false)
  (if (null? fields) (init-fields))
  (if (number? cond?) ; cond? is a time to run for
      (let ((T0 (meep-round-time))) ; current Meep time
	(apply run-until (cons (lambda () (>= (meep-round-time)
					      (+ T0 cond?)))
			       (cons (display-progress T0 (+ T0 cond?)
						       progress-interval)
				     step-funcs))))
      (begin ; otherwise, cond? is a boolean thunk
	(map (lambda (f) (eval-step-func f 'step)) step-funcs)
	(if (cond?)
	    (begin
	      (map (lambda (f) (eval-step-func f 'finish)) step-funcs)
	      (print "run " run-index " finished at t = " (meep-time)
		     " (" (meep-fields-t-get fields) " timesteps)\n")
	      (set! run-index (+ run-index 1)))
	    (begin
	      (meep-fields-step fields)
	      (apply run-until (cons cond? step-funcs)))))))

; run until all sources are finished and cond? is true.  If cond? is a number
; T, run until all sources are finished + a time T.
(define (run-sources+ cond? . step-funcs)
  (if (null? fields) (init-fields))
  (let ((Ts (meep-fields-last-source-time fields)))
  (apply run-until
	 (cons (if (number? cond?)
		   (+ (- Ts (meep-round-time)) cond?)
		   (lambda () (and (cond?) (>= (meep-round-time) Ts))))
	       step-funcs))))

; run until all sources are finished
(define (run-sources . step-funcs)
  (apply run-sources+ (cons 0 step-funcs)))

; condition function, designed to be used in conjunction with run-sources+,
; that returns true when |field|^2 at a given point has decayed more than
; a certain amount, always running for at least steps of dT.
(define (stop-when-fields-decayed dT c pt decay-by)
  (if (null? fields) (init-fields))
  (let ((T0 (meep-round-time))
	(max-abs (sqr (magnitude (meep-fields-get-field fields c pt))))
	(cur-max 0))
    (lambda ()
      (let ((fabs (sqr (magnitude (meep-fields-get-field fields c pt)))))
	(set! cur-max (max cur-max fabs))
	(if (<= (meep-round-time) (+ T0 dT))
	    false ; don't stop yet
	    (let ((old-cur cur-max))
	      (set! cur-max 0)
	      (set! T0 (meep-round-time))
	      (set! max-abs (max max-abs old-cur))
	      (if (not (zero? max-abs))
		  (print "field decay(t = " (meep-time)"): "
			 old-cur " / " max-abs " = " (/ old-cur max-abs) "\n"))
	      (<= old-cur (* max-abs decay-by))))))))


; ****************************************************************
; band diagrams

(define (run-k-point T k)
  (define cs (map (lambda (o) (object-property-value o 'component)) sources))
  (define pts (map (lambda (o) (object-property-value o 'center)) sources))
  (define As (map (lambda (o) (object-property-value o 'amplitude)) sources))
  (define fmin
    (max 0 (apply min (map (lambda (o)
			     (let ((t (object-property-value o 'src)))
			       (if (object-member? 'gaussian-src t)
				   (- (object-property-value t 'frequency)
				      (/ 1 (object-property-value t 'width) 2))
				   infinity)))
			   sources))))
  (define fmax
    (apply max (map (lambda (o)
		      (let ((t (object-property-value o 'src)))
			(if (object-member? 'gaussian-src t)
			    (+ (object-property-value t 'frequency)
			       (/ 1 (object-property-value t 'width) 2))
			    0)))
		    sources)))
  (if (or (null? cs) (> fmin fmax))
      (error "run-k-point requires a gaussian-src source"))
  ; TODO: apply harminv to multiple points and only accept freqs
  ; with correct relative amplitudes?
  (change-k-point! k)
  (restart-fields)
  (run-sources+ T (after-sources
		   (harminv (car cs) (car pts)
			    (* 0.5 (+ fmin fmax)) (- fmax fmin))))
  (map harminv-freq harminv-results))

(define (run-k-points T k-points)
  (define k-index 0)
  (define all-freqs '())
  (map (lambda (k)
	 (set! k-index (+ k-index 1))
	 (if (= k-index 1) (begin (change-k-point! k) (if (null? fields) (init-fields)) (output-epsilon)))
	 (let ((freqs (run-k-point T k)))
	   (print "freqs:, " k-index
		  ", " (vector3-x k) ", " (vector3-y k) ", " (vector3-z k))
	   (map (lambda (x) (print ", " x)) (map real-part freqs))
	   (print "\n")
	   (print "freqs-im:, " k-index
		  ", " (vector3-x k) ", " (vector3-y k) ", " (vector3-z k))
	   (map (lambda (x) (print ", " x)) (map imag-part freqs))
	   (print "\n")
	   (set! all-freqs (cons freqs all-freqs))))
       k-points)
  (reverse all-freqs))


; ****************************************************************
; field integration

(define (get-where-and-fields where-and-fields)
  (let ((f (if (= 2 (length where-and-fields))
	       (cadr where-and-fields)
	       fields)))
    (if (null? f)
	(error "init-fields is required before using field functions"))
    (let ((where (if (null? where-and-fields)
		     (meep-fields-total-volume f)
		     (car where-and-fields))))
      (cons where f))))

(define (integrate-field-function cs func . where-and-fields)
  (let ((waf (get-where-and-fields where-and-fields)))
    (meep-fields-integrate (cdr waf) (cons cs func) (car waf))))

(define (integrate2-field-function fields2 cs1 cs2 func . where-and-fields)
  (let ((waf (get-where-and-fields where-and-fields)))
    (meep-fields-integrate2 (cdr waf) fields2
			    (cons (cons cs1 cs2) func) (car waf))))

(define (max-abs-field-function cs func . where-and-fields)
  (let ((waf (get-where-and-fields where-and-fields)))
    (meep-fields-max-abs (cdr waf) (cons cs func) (car waf))))

(define (flux-in-box dir box)
  (if (null? fields)
      (error "init-fields is required before using flux-in-box"))
  (meep-fields-flux-in-box fields dir box))
(define (electric-energy-in-box box)
  (if (null? fields)
      (error "init-fields is required before using electric-energy-in-box"))
  (meep-fields-electric-energy-in-box fields box))
(define (magnetic-energy-in-box box)
  (if (null? fields)
      (error "init-fields is required before using magnetic-energy-in-box"))
  (meep-fields-magnetic-energy-in-box fields box))
(define (field-energy-in-box box)
  (if (null? fields)
      (error "init-fields is required before using field-energy-in-box"))
  (meep-fields-field-energy-in-box fields box))

; ****************************************************************
; Load helper functions for Casimir-force calculations
(if (defined? 'meep-component-Ex) (include "casimir.scm"))

; ****************************************************************

; Load GNU Readline support, for easier command-line editing support.
; This is not loaded in by default in Guile 1.3.2+ because readline is
; licensed under the GPL, which would have caused Guile to effectively
; be under the GPL itself.  However, since Meep is under the GPL too,
; we can load Readline by default with no problems.

@ACTIVATE_READLINE@  ; command to activate readline is determined by configure

(ctl-set-prompt! "meep> ")

; ****************************************************************