7 #:back-focal-plane-radius
8 #:focal-length-from-magnification
10 #:oil-objective-etendue
11 #:magnification-from-angles
15 #:find-inverse-ray-angle
))
17 ;; for i in `cat lens.lisp|grep "^(defun"|cut -d " " -f2`;do echo \#:$i ;done
20 (declaim (optimize (speed 2) (safety 3) (debug 3)))
22 (declaim (ftype (function (disk vec vec
)
23 (values vec
(or null vec
) &optional
))
25 (defun plane-ray (disk ray-start ray-direction
)
26 "Find the point where a ray intersects a plane. It returns two
27 vectors first the unchanged direction of the beam and then the
28 intersection point. The ray is defined by one position RAY-START and a
29 direction RAY-DIRECTION. The plane is defined by a point PLANE-CENTER
30 and its normal PLANE-NORMAL."
31 (with-slots (center normal
)
33 (let* ((hess-dist (v. center normal
)) ;; distance of plane to origin
34 (div (v. normal ray-direction
)))
35 (if (< (abs div
) 1d-13
)
36 ;; plane and ray are parallel and don't intersect
37 (values ray-direction nil
)
38 (let ((eta (/ (- hess-dist
(v. normal ray-start
))
39 ;; scaling along ray to hit exactly on plane
42 (v+ ray-start
(v* ray-direction eta
))))))))
44 (plane-ray (make-disk) (v 0d0
.1d0 -
1d0
) (v 0d0
0d0
1d0
))
48 ;; | -----\ intersection
52 ;; | ----+/ r' | ----\
53 ;; f |--/ c | ----\ dir
59 ;; -------+----------<-------+-------------------------------------------
68 (declaim (ftype (function (lens vec vec
69 &key
(:normalize boolean
))
70 (values (or null vec
) vec
&optional
))
73 ray-start ray-direction
75 "Return new ray-direction and intersection after thin lens."
76 (with-slots (center normal focal-length radius
)
78 (unless (< (abs (- 1 (norm ray-direction
))) 1d-12
)
79 (error "ray-direction doesn't have unit length"))
80 (unless (< (abs (- 1 (norm normal
))) 1e-12)
81 (error "lens-normal doesn't have unit length"))
82 (multiple-value-bind (dir intersection
)
84 ray-start ray-direction
)
85 (declare (ignore dir
))
86 (let* ((rho (v- intersection center
)))
87 #+nil
(format t
"~a~%" (list 'intersection intersection
'center center
88 'rho rho
'norm
(norm rho
)))
89 (unless (< (norm rho
) radius
)
90 (return-from lens-ray
(values nil intersection
)))
91 (let* ((cosphi (v. normal ray-direction
))
92 (r (v- (v* ray-direction
(/ focal-length cosphi
))
105 ;; --------------- s |
109 ;; | ----/ | \-------
110 ;; / -----/ \ | \------
111 ;; / ----/ ru \ | \-------
112 ;; /-/ \|rho \-------
115 ;; ------------------------------+-----------------------------------------
118 ;; +--------------------
121 (declaim (ftype (function (thin-objective
123 (values (or null vec
) vec
&optional
))
125 ;; 2008 Hwang Simulation of an oil immersion objective lens...
126 (defun thin-objective-ray (objective
127 ray-start ray-direction
)
128 "Returns direction of outgoing ray and intersection with principal
129 sphere. If this sphere isn't hit inside the maximum angle (given by
130 NA), then nil as a direction and the intersection with the
131 lens-surface is returned instead."
132 (with-slots (center normal focal-length
134 immersion-index numerical-aperture
)
136 (unless (< (abs (- 1 (norm normal
))) 1d-12
)
137 (error "lens normal doesn't have unit length"))
138 (unless (< (abs (- 1 (norm ray-direction
))) 1d-12
)
139 (error "ray-direction doesn't have unit length"))
141 (multiple-value-bind (r intersection
)
143 ray-start ray-direction
:normalize nil
)
145 (return-from thin-objective-ray
(values nil intersection
)))
146 (let* ((a (v* normal
(* focal-length
(- immersion-index
1))))
148 (rho (v- intersection center
))
150 (nf (* immersion-index focal-length
))
154 (return-from thin-objective-ray
(values nil intersection
))
155 #+nil
(error "ray can't pass through objective"))
156 (let* ((s (v* ray-direction
(- nf
(sqrt rat
))))
158 (cosu (v. ro normal
))
159 (sinu2 (- 1 (* cosu cosu
)))
160 (sinu-max (/ numerical-aperture immersion-index
)))
161 (unless (<= sinu2
(* sinu-max sinu-max
))
162 (return-from thin-objective-ray
(values nil intersection
)))
163 (values ro
(v+ s intersection
)))))))
166 (thin-objective-ray (v 0d0
0d0
0d0
)
170 (normalize (v 0d0
0d0 -
1d0
)))
177 ;; X-----+-------------------------------
184 ;; ---------------------+----------------------------------------
190 (declaim (ftype (function (double-float double-float
)
191 (values double-float
&optional
))
192 back-focal-plane-radius
))
193 (defun back-focal-plane-radius (focal-length numerical-aperture
)
194 (* focal-length numerical-aperture
))
196 (back-focal-plane-radius 2.61d0
1.4d0
)
198 (declaim (ftype (function (double-float)
199 (values double-float
&optional
))
200 focal-length-from-magnification
))
201 (defun focal-length-from-magnification (mag)
204 (focal-length-from-magnification 63d0
)
206 (declaim (ftype (function (double-float double-float
207 &optional double-float
208 double-float double-float
)
209 (values double-float
&optional
))
211 (defun etendue (chief-height marginal-angle
213 (refractive-index 1d0
)
214 (marginal-height 0d0
) (chief-angle 0d0
))
215 (let ((q (- (* chief-height refractive-index marginal-angle
)
216 (* marginal-height refractive-index chief-angle
))))
219 (etendue .07d0
(* 67d0
(/ pi
180)) 1.515d0
)
222 (declaim (ftype (function (double-float &optional double-float
224 (values double-float
&optional
))
225 oil-objective-etendue
))
226 (defun oil-objective-etendue (field-radius &optional
(numerical-aperture 1.4d0
)
227 (refractive-index 1.515d0
))
228 (let ((rat (/ numerical-aperture refractive-index
)))
229 (unless (<= (abs rat
) 1)
230 (error "impossible angle, check numerical aperture and refractive index."))
231 (let* ((marginal-angle (asin rat
)))
232 (etendue field-radius marginal-angle refractive-index
))))
234 (oil-objective-etendue .07d0
)
236 (declaim (ftype (function (double-float double-float
237 &optional double-float double-float
)
238 (values double-float
&optional
))
239 magnification-from-angles
))
240 (defun magnification-from-angles (u uu
&optional
(n 1d0
) (nn 1d0
))
241 "u is the objects marginal ray angle and uu for the image."
245 (declaim (ftype (function (mirror vec vec
)
246 (values (or null vec
) vec
&optional
))
248 ;; sketch of the mirror for incoming parallel light
249 ;; --------------------+-----------------------
252 ;; / n| \ N=n*(- (p.n))
256 ;; / | \ q=p-2(p.n)*n
262 ;; p .. ray-direction
263 ;; N .. mirror-normal
265 (defun mirror-ray (mirror ray-start ray-direction
)
266 "Mirror center C, mirror normal N.
267 Return reflected direction and intersection. If the ray isn't inside
268 of the radius return nil and the intersection."
269 (with-slots (center normal radius
)
271 (unless (< (abs (- 1 (norm normal
))) 1d-12
)
272 (error "mirror normal doesn't have unit length"))
273 (unless (< (abs (- 1 (norm ray-direction
))) 1d-12
)
274 (error "ray-direction doesn't have unit length"))
275 (multiple-value-bind (dir intersection
)
276 (plane-ray mirror ray-start ray-direction
)
277 (declare (ignore dir
))
278 (unless (< (norm (v- intersection center
)) radius
)
279 (return-from mirror-ray
(values nil intersection
)))
280 (let ((dir (v+ ray-direction
(v* normal
281 (* -
2d0
(v. ray-direction normal
))))))
282 (values dir intersection
)))))
284 (deftype optical-component
()
285 `(or disk mirror lens thin-objective
))
289 (start (v) :type vec
)
290 (direction (v) :type vec
))
293 (normal (v 0d0
0d0
1d0
) :type vec
)
294 (center (v) :type vec
)
295 (radius 1d0
:type double-float
))
297 (defstruct (mirror (:include disk
)))
299 (defstruct (lens (:include disk
))
300 (focal-length 1d0
:type double-float
))
302 (defstruct (thin-objective (:include lens
))
303 (immersion-index 1.515d0
:type double-float
)
304 (numerical-aperture 1.38d0
:type double-float
))
306 (declaim (ftype (function (double-float vec vec
)
307 (values vec
&optional
))
309 (defun rotate-vector (angle axis vector
)
310 (let ((m (rotation-matrix angle axis
)))
318 ;; -----+-------+------------o------------\------------------------
320 (declaim (ftype (function (double-float double-float thin-objective
)
321 (values double-float
&optional
))
323 ;; 2005wolf p. 180 formula (8) sine condition for one image in infinity
324 (defun find-inverse-ray-angle (point-x point-y objective
)
325 "Find the angle in the back focal plane for a ray that originates
326 from a POINT in the object. The coordinates of the point are given
327 relative to the center of the front focal plane."
328 (with-slots (focal-length)
330 (let* ((h (sqrt (+ (* point-x point-x
) (* point-y point-y
))))
331 (phi (asin (/ h focal-length
))))
335 (let* ((f (focal-length-from-magnification 63d0
))
338 (bfp-radius (back-focal-plane-radius f na
))
339 (obj (make-thin-objective :normal
(v 0d0
0d0 -
1d0
)
340 :center
(v 0d0
0d0
(- f
))
342 :numerical-aperture na
343 :immersion-index ri
))
344 (field (v .0d0
0d0
(- (+ f
(* ri f
)))))
345 (pupil (v bfp-radius
0d0
0d0
)))
346 (list bfp-radius f
(find-inverse-ray-angle 0 .07 obj
)))