meep.scm.in

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

	(set! cur-eps meep-last-eps-filename)
	(set! output-volume gv-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 (and rm? (zero? (system (string-append convert-cmd " \"" fname "\""))))
	(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 (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)))
    (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))

; ****************************************************************
; 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-time))) ; current Meep time
	(apply run-until (cons (lambda () (>= (meep-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-time)) cond?)
		   (lambda () (and (cond?) (>= (meep-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-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-time) (+ T0 dT))
	    false ; don't stop yet
	    (let ((old-cur cur-max))
	      (set! cur-max 0)
	      (set! T0 (meep-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 (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 (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 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

(set! scm-repl-prompt "meep> ")

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