tebdol

Simulation of ultracold atoms in optical lattices
git clone https://miroslavurbanek.com/git/tebdol.git
Log | Files | Refs | README

conv.lisp (1573B)


      1 (defpackage :conv
      2   (:use :common-lisp)
      3   (:export :normalized-gaussian
      4 	   :convolve))
      5 
      6 (in-package :conv)
      7 
      8 (defun gaussian (width height sigma)
      9   (loop
     10      with a = (make-array (list width height))
     11      for i below width
     12      do (loop
     13 	   for j below height
     14 	   do (setf (aref a i j)
     15 		    (exp (/ (+ (expt (- i (/ (1- width) 2)) 2)
     16 			       (expt (- j (/ (1- height) 2)) 2))
     17 			    (* -2 (expt sigma 2))))))
     18      finally (return a)))
     19 
     20 (defun sum-array (array)
     21   (loop
     22      for i below (array-total-size array)
     23      sum (row-major-aref array i)))
     24 
     25 (defun normalize-array (array)
     26   (loop
     27      with a = (make-array (array-dimensions array))
     28      with s = (sum-array array)
     29      for i below (array-total-size a)
     30      do (setf (row-major-aref a i) (/ (row-major-aref array i) s))
     31      finally (return a)))
     32 
     33 (defun normalized-gaussian (width height sigma)
     34   (normalize-array (gaussian width height sigma)))
     35 
     36 (defun convolve (array kernel)
     37   (loop
     38      with a = (make-array (array-dimensions array))
     39      with aw = (array-dimension array 0)
     40      with ah = (array-dimension array 1)
     41      with kw = (array-dimension kernel 0)
     42      with kh = (array-dimension kernel 1)
     43      for i below aw
     44      do (loop
     45 	   for j below ah
     46 	   do (setf (aref a i j)
     47 		    (loop
     48 		       for k below kw
     49 		       for m = (+ i (- k (floor kw 2)))
     50 		       when (and (>= m 0) (< m aw))
     51 		       sum (loop
     52 			      for l below kh
     53 			      for n = (+ j (- l (floor kh 2)))
     54 			      when (and (>= n 0) (< n ah))
     55 			      sum (* (aref array m n)
     56 				     (aref kernel k l))))))
     57      finally (return a)))