Source code for rayoptics.raytr.raytrace

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright © 2018 Michael J. Hayford
""" Functions to support ray tracing a sequential optical model

.. Created on Thu Jan 25 11:01:04 2018

.. codeauthor: Michael J. Hayford
"""

import numpy as np
from numpy.linalg import norm
from math import sqrt, copysign

import rayoptics.optical.model_constants as mc
from .traceerror import (TraceMissedSurfaceError, TraceTIRError,
                         TraceRayBlockedError, TraceEvanescentRayError)

[docs] def bend(d_in, normal, n_in, n_out): """ refract incoming direction, d_in, about normal """ try: normal_len = norm(normal) cosI = np.dot(d_in, normal)/normal_len sinI_sqr = 1.0 - cosI*cosI n_cosIp = copysign(sqrt(n_out*n_out - n_in*n_in*sinI_sqr), cosI) alpha = n_cosIp - n_in*cosI d_out = (n_in*d_in + alpha*normal)/n_out return d_out except ValueError: raise TraceTIRError(d_in, normal, n_in, n_out)
[docs] def reflect(d_in, normal): """ reflect incoming direction, d_in, about normal """ normal_len = norm(normal) cosI = np.dot(d_in, normal)/normal_len d_out = d_in - 2.0*cosI*normal return d_out
[docs] def phase(ifc, inc_pt, d_in, normal, ifc_cntxt): """ apply phase shift to incoming direction, d_in, about normal """ try: d_out, dW = ifc.phase(inc_pt, d_in, normal, ifc_cntxt) return d_out, dW except ValueError: z_dir, wvl, n_in, n_out, interact_mode = ifc_cntxt raise TraceEvanescentRayError(ifc, inc_pt, d_in, normal, n_in, n_out)
[docs] def trace(seq_model, pt0, dir0, wvl, **kwargs): """ fundamental raytrace function Args: seq_model: the sequential model to be traced pt0: starting point in coords of first interface dir0: starting direction cosines in coords of first interface wvl: wavelength in nm eps: accuracy tolerance for surface intersection calculation Returns: (**ray**, **op_delta**, **wvl**) - **ray** is a list for each interface in **path_pkg** of these elements: [pt, after_dir, after_dst, normal] - pt: the intersection point of the ray - after_dir: the ray direction cosine following the interface - after_dst: after_dst: the geometric distance to the next interface - normal: the surface normal at the intersection point - **op_delta** - optical path wrt equally inclined chords to the optical axis - **wvl** - wavelength (in nm) that the ray was traced in """ path = seq_model.path(wvl) kwargs['first_surf'] = kwargs.get('first_surf', 1) kwargs['last_surf'] = kwargs.get('last_surf', seq_model.get_num_surfaces()-2) return trace_raw(path, pt0, dir0, wvl, **kwargs)
[docs] def trace_raw(path, pt0, dir0, wvl, eps=1.0e-12, check_apertures=False, **kwargs): """ fundamental raytrace function Args: path: an iterator containing interfaces and gaps to be traced. for each iteration, the sequence or generator should return a list containing: **Intfc, Gap, Trfm, Index, Z_Dir** pt0: starting point in coords of first interface dir0: starting direction cosines in coords of first interface wvl: wavelength in nm eps: accuracy tolerance for surface intersection calculation check_apertures: if True, do point_inside() test on inc_pt Returns: (**ray**, **op_delta**, **wvl**) - **ray** is a list for each interface in **path** of these elements: [pt, after_dir, after_dst, normal] - pt: the intersection point of the ray - after_dir: the ray direction cosine following the interface - after_dst: the geometric distance to the next interface - normal: the surface normal at the intersection point - **op_delta** - optical path wrt equally inclined chords to the optical axis - **wvl** - wavelength (in nm) that the ray was traced in """ ray = [] first_surf = kwargs.get('first_surf', 0) last_surf = kwargs.get('last_surf', None) def in_gap_range(gap_indx, include_last_surf=False): if first_surf == last_surf: return False if gap_indx < first_surf: return False if last_surf is None: return True else: return (gap_indx <= last_surf if include_last_surf else gap_indx < last_surf) def in_surface_range(s): if s < first_surf: return False if last_surf is None: return True elif s > last_surf: return False else: return True # trace object surface obj = next(path) srf_obj = obj[mc.Intfc] dst_b4, pt_obj = srf_obj.intersect(pt0, dir0, z_dir=obj[mc.Zdir]) before = obj before_pt = pt_obj before_dir = dir0 before_normal = srf_obj.normal(before_pt) tfrm_from_before = before[mc.Tfrm] z_dir_before = before[mc.Zdir] op_delta = 0.0 opl = 0.0 surf = 0 # loop of remaining surfaces in path while True: try: after = next(path) surf += 1 # transform ray data from previous ifc coords to current one rt, t = tfrm_from_before b4_pt, b4_dir = rt.dot(before_pt - t), rt.dot(before_dir) pp_dst = -b4_pt.dot(b4_dir) pp_pt_before = b4_pt + pp_dst*b4_dir ifc = after[mc.Intfc] z_dir_after = after[mc.Zdir] # intersect ray with profile pp_dst_intrsct, inc_pt = ifc.intersect(pp_pt_before, b4_dir, eps=eps, z_dir=z_dir_before) dst_b4 = pp_dst + pp_dst_intrsct # add *previous* intersection point, direction, etc., to ray ray.append([before_pt, before_dir, dst_b4, before_normal]) if in_gap_range(surf-1): opl += before[mc.Indx] * dst_b4 normal = ifc.normal(inc_pt) if check_apertures and in_surface_range(surf): if not ifc.point_inside(inc_pt[0], inc_pt[1]): raise TraceRayBlockedError(ifc, inc_pt) # if present, use the phase element to calculate after_dir if hasattr(ifc, 'phase_element'): ifc_cntxt = (z_dir_before, wvl, before[mc.Indx], after[mc.Indx], ifc.interact_mode) after_dir, phs = phase(ifc, inc_pt, b4_dir, normal, ifc_cntxt) op_delta += phs else: # refract or reflect ray at interface if ifc.interact_mode == 'reflect': after_dir = reflect(b4_dir, normal) elif ifc.interact_mode == 'transmit': after_dir = bend(b4_dir, normal, before[mc.Indx], after[mc.Indx]) elif ifc.interact_mode == 'dummy': after_dir = b4_dir else: # no action, input becomes output after_dir = b4_dir before_pt = inc_pt before_normal = normal before_dir = after_dir z_dir_before = z_dir_after before = after tfrm_from_before = before[mc.Tfrm] except TraceMissedSurfaceError as ray_miss: ray.append([before_pt, before_dir, pp_dst, before_normal]) ray_miss.surf = surf ray_miss.ifc = ifc ray_miss.prev_tfrm = before[mc.Tfrm] ray_miss.ray_pkg = ray, opl, wvl raise ray_miss except TraceTIRError as ray_tir: ray.append([inc_pt, before_dir, 0.0, normal]) ray_tir.surf = surf ray_tir.ray_pkg = ray, opl, wvl raise ray_tir except TraceRayBlockedError as ray_blocked: ray.append([inc_pt, before_dir, 0.0, normal]) ray_blocked.surf = surf ray_blocked.ray_pkg = ray, opl, wvl raise ray_blocked except TraceEvanescentRayError as ray_evn: ray.append([inc_pt, before_dir, 0.0, normal]) ray_evn.surf = surf ray_evn.ray_pkg = ray, opl, wvl raise ray_evn except StopIteration: ray.append([inc_pt, after_dir, 0.0, normal]) op_delta += opl break return ray, op_delta, wvl
[docs] def calc_optical_path(ray, path): """ computes equally inclined chords and path info for ray Args: ray: ray data for traced ray path: an iterator containing interfaces and gaps to be traced. for each iteration, the sequence or generator should return a list containing: **Intfc, Gap, Trfm, Index, Z_Dir** Returns: **ray_op** - **ray_op** - the optical path between the first and last optical surfaces """ num_items = len(ray) ray_seq_iter = zip(ray, path) next(ray_seq_iter) ray_op = 0 for i in range(1, num_items-2): after_ray_seg, surf = next(ray_seq_iter) n_after = surf[mc.Indx] dist = after_ray_seg[mc.dst] ray_op += n_after*dist return ray_op