meep.scm.in

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

; Copyright (C) 2005 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

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

; 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 (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-class material-type no-parent)

(define-class polarizability no-parent
  (define-property omega no-default 'number)
  (define-property gamma no-default 'number)
  (define-property delta-epsilon no-default 'number)
  (define-property sigma 1.0 'number)
  (define-property energy-saturation 0.0 'number))

(define-class dielectric material-type
  (define-property epsilon no-default 'number)
  (define-property polarizations '() (make-list-type 'polarizability))
  (define-property chi2 0 'number)
  (define-property chi3 0 'number)
)

(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")

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

(define CYLINDRICAL -2) ; special value of dimensions for cylindrical coords

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

(define vacuum (make dielectric (epsilon 1.0)))
(define air vacuum)
(define metal (make perfect-metal))

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

(set! default-material vacuum)

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

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

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

; 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

(define structure '())

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

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

(define-param k-point false)

(define fields '())

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

; ****************************************************************
; 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))

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

(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) (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-geometric-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 default-material
		     pml-layers symmetries num-chunks Courant))))

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

(define-class src-time no-parent)
(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) 
      (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)))))
(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) 
      (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)))))))
(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) 
      (new-meep-custom-src-time 
       (object-property-value o 'src-func)
       (object-property-value o 'start-time)
       (object-property-value o 'end-time)))))

(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))

(define (add-source s f) ; add source s to fields f
  (let ((A (object-property-value s 'amp-func))
	(cen (object-property-value s 'center))
	(sz (object-property-value s 'size)))
    (if (null? A)
	(meep-fields-add-volume-source
	 fields
	 (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
	 fields
	 (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))
  (if (= dimensions CYLINDRICAL)
      (set! fields (new-meep-fields structure m))
      (set! fields (new-meep-fields structure)))
  (if verbose? (meep-fields-verbose fields))
  (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) (vector3= k-point (vector3 0))))))
      (meep-fields-use-real-fields fields))
  (if k-point (meep-fields-use-bloch fields 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 (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 '())