apotome.lisp

   1 (in-package :run)
   2
   3
   4 ;; calculate e field
   5 (let* ((dx .1d0)
   6              (dz .1d0)
   7              (x 100)
   8              (y x)
   9              (z 100))
  10         (multiple-value-bind (ex ey ez)
  11        (psf:electric-field-psf z y x (* dz z) (* dx x) :integrand-evaluations 200)
  12           (defparameter *kex* (fftshift3 (ft3 ex)))
  13           (defparameter *key* (fftshift3 (ft3 ey)))
  14           (defparameter *kez* (fftshift3 (ft3 ez)))
  15           nil))
  16 #+Nil
  17 (progn
  18  (save-stack-ub8 "/home/martin/tmp/kex" (normalize3-cdf/ub8-realpart *kex*))
  19  (save-stack-ub8 "/home/martin/tmp/key" (normalize3-cdf/ub8-realpart *key*))
  20  (save-stack-ub8 "/home/martin/tmp/kez" (normalize3-cdf/ub8-realpart *kez*))
  21  )
  22 #+NIL
  23 (write-pgm "/home/martin/tmp/kex.pgm"
  24            (normalize2-cdf/ub8-abs (cross-section-xz *kex*)))
  25
  26 ;; calculate volume containing a grating
  27 (destructuring-bind (z y x)
  28      (array-dimensions *kex*)
  29    (let* ((grat (make-array (list z y x) :element-type '(complex double-float))))
  30      (let ((k (floor z 2))
  31            (q 30))
  32        (do-rectangle (j i q (- y q) q (- x q))
  33          (setf (aref grat k j i) (complex (* 1d0 (mod i 8))))))
  34      (defparameter *kgrat* (fftshift3 (ft3 grat)))
  35      nil))
  36 #+nil
  37 (save-stack-ub8 "/home/martin/tmp/kgrat" (normalize3-cdf/ub8-abs *kgrat*))
  38
  39 ;; calculate image of grating
  40 (destructuring-bind (z y x)
  41     (array-dimensions *kex*)
  42  (let* ((ex (fftshift3 (ift3 (fftshift3 (.* *kgrat* *kex*)))))
  43         (ey (fftshift3 (ift3 (fftshift3 (.* *kgrat* *key*)))))
  44         (ez (fftshift3 (ift3 (fftshift3 (.* *kgrat* *kez*)))))
  45         (intens (make-array (array-dimensions ex)
  46                             :element-type '(complex double-float))))
  47    (do-box (k j i 0 z 0 y 0 x)
  48      (setf (aref intens k j i)
  49            (+ (* (conjugate (aref ex k j i)) (aref ex k j i))
  50               (* (conjugate (aref ey k j i)) (aref ey k j i))
  51               (* (conjugate (aref ez k j i)) (aref ez k j i)))))
  52    (defparameter *intens* intens)
  53    (save-stack-ub8 "/home/martin/tmp/intens-grat" (normalize3-cdf/ub8-realpart intens))
  54    (write-pgm "/home/martin/tmp/intens-grat.pgm" (normalize2-cdf/ub8-realpart (cross-section-xz intens)))))
  55
  56 ;; store laser psf
  57 (destructuring-bind (z y x)
  58     (array-dimensions *kex*)
  59  (let* ((ex (ift3 (fftshift3 *kex*)))
  60         (ey (ift3 (fftshift3 *key*)))
  61         (ez (ift3 (fftshift3 *kez*)))
  62         (intens (make-array (array-dimensions ex)
  63                             :element-type '(complex double-float))))
  64    (do-box (k j i 0 z 0 y 0 x)
  65      (setf (aref intens k j i)
  66            (+ (* (conjugate (aref ex k j i)) (aref ex k j i))
  67               (* (conjugate (aref ey k j i)) (aref ey k j i))
  68               (* (conjugate (aref ez k j i)) (aref ez k j i)))))
  69    (defparameter *intens-psf-laser* intens)
  70    (save-stack-ub8 "/home/martin/tmp/intens-psf-laser" (normalize3-cdf/ub8-realpart intens))))
  71 #+nil
  72 (write-pgm
  73  "/home/martin/tmp/intens-psf-laser.pgm"
  74  (normalize2-cdf/ub8-realpart (cross-section-xz *intens-psf-laser*)))
  75
  76 #+nil
  77 (time
  78  (save-stack-ub8 "/home/martin/tmp/kintens-grat"
  79                  (normalize3-cdf/ub8-abs (fftshift3 (ft3 *intens*)))))
  80
  81
  82 (time (destructuring-bind (z y x)
  83      (array-dimensions *intens*)
  84    (let ((n 1.515d0)
  85          (f (lens:focal-length-from-magnification 63d0))
  86          (Rap 4d0)
  87          (dz .21d0)
  88          (result (make-array (array-dimensions *intens*)
  89                              :element-type '(complex double-float))))
  90      ;; center plane doesn't have to be convolved
  91      (let ((k (floor z 2)))
  92        (do-rectangle (j i 0 y 0 x)
  93          (setf (aref result k j i) (aref *intens* k j i))))
  94
  95      ;; convolve all other planes with disks of different sizes
  96      (dotimes (k z)
  97        (unless (eq k (floor z 2))
  98          (let* ((Rd-pixels (abs (* (- k (floor z 2)) dz (/ Rap (* n f)))))
  99                 (disk (draw-unit-energy-disk-precise Rd-pixels y x))
 100                 (img (make-array (list y x) :element-type '(complex double-float))))
 101            (write-pgm (format nil "/home/martin/tmp/disk/disk~3,'0d.pgm" k)
 102                       (normalize2-cdf/ub8-realpart disk))
 103            (do-rectangle (j i 0 y 0 x)
 104              (setf (aref img j i) (aref *intens* k j i)))
 105            (let ((conv (fftshift2 (convolve2-circ img disk))))
 106              (do-rectangle (j i 0 y 0 x)
 107                (setf (aref result k j i) (aref conv j i)))))))
 108      (defparameter *intens-disk* result)
 109      (save-stack-ub8 "/home/martin/tmp/intens-grat-conv" (normalize3-cdf/ub8-abs result))
 110      (save-stack-ub8 "/home/martin/tmp/kintens-grat-conv"
 111                      (normalize3-cdf/ub8-realpart (ft3 result)))
 112      (write-pgm "/home/martin/tmp/intens-grat-conv.pgm"
 113                 (normalize2-cdf/ub8-realpart (cross-section-xz *intens-disk*))))))