mbl.lisp (3623B)
1 ;; random lattice potential 2 3 ;; parameters 4 5 (defparameter *width* 4) 6 (defparameter *height* 4) 7 (defparameter *diameter* 2) 8 (defparameter *epp* 1.4d-1) 9 (defparameter *fwhm* 1d0) 10 11 ;; fixed seed 12 (defparameter *seed* 0) 13 14 ;; random seed, uncomment to get a real random potential 15 ;; (defparameter *seed* (random #x100000000 (make-random-state t))) 16 17 ;; libraries 18 19 (load (merge-pathnames "gsl.lisp" *load-truename*)) 20 (load (merge-pathnames "conv.lisp" *load-truename*)) 21 (use-package '(:gsl :conv)) 22 23 ;; main 24 25 (defun harmonic-potential (width height) 26 (loop 27 with a = (make-array (list width height)) 28 for i below width 29 do (loop 30 for j below height 31 do (setf (aref a i j) 32 (+ (expt (- i (/ (1- width) 2)) 2) 33 (expt (- j (/ (1- height) 2)) 2)))) 34 finally (return a))) 35 36 (defun half-circle-mean-energy (potential diameter) 37 (loop 38 with m = (ceiling (- (array-dimension potential 0) diameter) 2) 39 with n = (floor (- (array-dimension potential 1) diameter) 2) 40 with p = 0 41 with e = 0 42 for i below (/ diameter 2) 43 do (loop 44 for j below diameter 45 when (<= (sqrt (+ (expt (- i (/ (1- diameter) 2)) 2) 46 (expt (- j (/ (1- diameter) 2)) 2))) 47 (/ diameter 2)) 48 do 49 (incf p) 50 (incf e (aref potential (+ i m) (+ j n)))) 51 finally (return (/ e p)))) 52 53 (defun scaled-harmonic-potential (width height diameter epp) 54 (loop 55 with a = (harmonic-potential width height) 56 with c = (/ (* 24.4d0 epp) (half-circle-mean-energy a diameter)) 57 for i below (array-total-size a) 58 do (setf (row-major-aref a i) (* c (row-major-aref a i))) 59 finally (return a))) 60 61 (defparameter *gsl-rng* (gsl-rng-alloc +gsl-rng-ranlxd2+)) 62 63 (gsl-rng-set *gsl-rng* *seed*) 64 65 (defun binomial-variate (p n) 66 (gsl-ran-binomial *gsl-rng* p n)) 67 68 (defun squared-binomial-mean (p n) 69 (/ (+ p (* (expt p 2) (1- n))) n)) 70 71 (defconstant +squared-binomial-fwhm+ 1.06d-1) ;; measured value 72 73 (defun squared-binomial-potential (width height fwhm) 74 (loop 75 with c = (/ fwhm +squared-binomial-fwhm+) 76 with m = (squared-binomial-mean 5d-1 49) 77 with a = (make-array (list width height)) 78 for i below (array-total-size a) 79 do (setf (row-major-aref a i) 80 (* c (- (expt (/ (binomial-variate 5d-1 49) 49) 2) m))) 81 finally (return a))) 82 83 (defun convolved-potential (width height fwhm) 84 (loop 85 with p = (squared-binomial-potential (+ width 8) (+ height 8) fwhm) 86 with g = (normalized-gaussian 9 9 5d-1) 87 with a = (convolve p g) 88 with b = (make-array (list width height)) 89 for i below width 90 do (loop 91 for j below height 92 do (setf (aref b i j) 93 (aref a (+ i 4) (+ j 4)))) 94 finally (return b))) 95 96 (defun full-potential (width height diameter epp fwhm) 97 (loop 98 with a = (scaled-harmonic-potential width height diameter epp) 99 with b = (convolved-potential width height fwhm) 100 for i below (array-total-size a) 101 do (incf (row-major-aref a i) (row-major-aref b i)) 102 finally (return a))) 103 104 (defun save-potential (potential pathname) 105 (when (or (not (probe-file pathname)) (y-or-n-p "Overwrite potential ~S?" pathname)) 106 (with-standard-io-syntax 107 (let ((*print-pretty* t)) 108 (with-open-file (s pathname :direction :output :if-exists :supersede) 109 (format s ";; width = ~A, height = ~A, diameter = ~A, epp = ~A, fwhm = ~A, seed = ~A~%" 110 *width* *height* *diameter* *epp* *fwhm* *seed*) 111 (prin1 potential s) 112 (terpri s)))))) 113 114 (save-potential 115 (full-potential *width* *height* *diameter* *epp* *fwhm*) 116 (pathname 117 (format nil "inputs/mbl-w-~A-h-~A.pot" *width* *height*))) 118 119 (gsl-rng-free *gsl-rng*)