diff --git a/docs/examples.rst b/docs/examples.rst index 010c828..5e5a074 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -65,7 +65,7 @@ We now turn the combination of X-ray source and first mirror around the y-axis ( ... rotmat[:3, :3] = euler.euler2mat(angle, 0, 0, 'szxy') ... light.position = np.dot(rotmat, light_pos) ... light.dir = np.dot(rotmat, light_dir) - ... m1.pos4d = np.dot(rotmat, m1pos4d) + ... m1.geometry.pos4d = np.dot(rotmat, m1pos4d) ... rays = light(100) ... pos.data = [] ... pos(rays) diff --git a/docs/pyplots/vis_pol.py b/docs/pyplots/vis_pol.py index 90f1bb8..7ae68e6 100644 --- a/docs/pyplots/vis_pol.py +++ b/docs/pyplots/vis_pol.py @@ -33,7 +33,7 @@ rotmat[:3, :3] = euler.euler2mat(angle, 0, 0, 'szxy') light.position = np.dot(rotmat, light_pos) light.dir = np.dot(rotmat, light_dir) - m1.pos4d = np.dot(rotmat, m1pos4d) + m1.geometry.pos4d = np.dot(rotmat, m1pos4d) rays = light(10000) rays = experiment(rays) n_detected[i] = rays['probability'].sum() diff --git a/marxs/analysis/analysis.py b/marxs/analysis/analysis.py index 88f2f39..31ab1c1 100644 --- a/marxs/analysis/analysis.py +++ b/marxs/analysis/analysis.py @@ -75,8 +75,8 @@ def width(x, photons): mdet = FlatDetector(position=np.dot(orientation, np.array([x, 0, 0])), orientation=orientation, zoom=1e5, pixsize=1.) - photons = mdet(photons) - return objective_func(photons, **objective_func_args) + p = mdet(photons.copy()) + return objective_func(p, **objective_func_args) return scipy.optimize.minimize_scalar(width, args=(photons,), options={'maxiter': 20, 'disp': False}, **kwargs) diff --git a/marxs/analysis/coords.py b/marxs/analysis/coords.py index bddc874..00595b9 100644 --- a/marxs/analysis/coords.py +++ b/marxs/analysis/coords.py @@ -35,8 +35,10 @@ class ProjectOntoPlane(FlatOpticalElement): loc_coos_name = ['proj_x', 'proj_y'] '''name for output columns of the projected position in plane coordinates.''' + display = {'shape': 'None'} + def __call__(self, photons): - vec_center_inter = - h2e(self.geometry('center')) + h2e(photons['pos']) - photons[self.loc_coos_name[0]] = np.dot(vec_center_inter, h2e(self.geometry('e_y'))) - photons[self.loc_coos_name[1]] = np.dot(vec_center_inter, h2e(self.geometry('e_z'))) + vec_center_inter = - h2e(self.geometry['center']) + h2e(photons['pos']) + photons[self.loc_coos_name[0]] = np.dot(vec_center_inter, h2e(self.geometry['e_y'])) + photons[self.loc_coos_name[1]] = np.dot(vec_center_inter, h2e(self.geometry['e_z'])) return photons diff --git a/marxs/analysis/gratings.py b/marxs/analysis/gratings.py index a4b31b1..324037c 100644 --- a/marxs/analysis/gratings.py +++ b/marxs/analysis/gratings.py @@ -8,6 +8,7 @@ from ..design import RowlandTorus from . import (find_best_detector_position) from .analysis import sigma_clipped_std +from ..math.geometry import Cylinder class AnalysisError(Exception): @@ -77,7 +78,9 @@ def resolvingpower_per_order(gratings, photons, orders, detector=None, col = 'det_x' # 0 at center, detpix_x is 0 in corner. zeropos = 0. elif isinstance(detector, RowlandTorus): - det = CircularDetector.from_rowland(detector, width=1e5) + det = CircularDetector(geometry=Cylinder.from_rowland(detector, width=1e5, + rotation=np.pi, + kwargs={'phi_lim':[-np.pi/2, np.pi/2]})) info['method'] = 'Circular detector on Rowland circle' col = 'detpix_x' # Find position of order 0. diff --git a/marxs/analysis/tests/test_analysis.py b/marxs/analysis/tests/test_analysis.py index c68c268..e7627d3 100644 --- a/marxs/analysis/tests/test_analysis.py +++ b/marxs/analysis/tests/test_analysis.py @@ -1,16 +1,17 @@ import numpy as np from ...utils import generate_test_photons from ..analysis import find_best_detector_position, mean_width_2d - +from ...simulator import propagate def test_best_detector_position(): photons = generate_test_photons(100) photons['dir'].data[:, 1] = np.random.rand(100) photons['pos'].data[:, 1] = 0.5 + photons = propagate(photons, -10) # det_x is along the second axis out = find_best_detector_position(photons, objective_func_args={'colname': 'det_x'}) - assert np.isclose(out.x, 1.0) + assert np.isclose(out.x, 1.) def test_best_detector_position_rotated(): @@ -19,6 +20,7 @@ def test_best_detector_position_rotated(): photons['dir'].data[:, 0] = np.random.rand(100) photons['dir'].data[:, 2] = - 1 photons['pos'].data[:, 2] = 0.5 + photons = propagate(photons, -10) # det_x is along the second axis out = find_best_detector_position(photons, objective_func_args={'colname': 'det_x'}, @@ -33,6 +35,7 @@ def test_best_detector_position_2d(): photons = generate_test_photons(100) photons['dir'].data[:, 2] = np.random.rand(100) photons['pos'].data[:, 1] = 0.5 + photons = propagate(photons, -10) # det_x is along the second axis out = find_best_detector_position(photons, mean_width_2d, {}) assert np.isclose(out.x, 1.0) diff --git a/marxs/base/base.py b/marxs/base/base.py index acc57e8..1f2b4b5 100644 --- a/marxs/base/base.py +++ b/marxs/base/base.py @@ -2,6 +2,7 @@ from collections import OrderedDict import inspect import warnings +from copy import deepcopy import numpy as np from transforms3d import affines @@ -9,6 +10,7 @@ from astropy.table import Column from astropy.extern.six import with_metaclass +from ..visualization.utils import DisplayDict class GeometryError(Exception): pass @@ -37,13 +39,17 @@ def __new__(mcs, name, bases, dict): for k in dict: # for items that have a docstring (i.e. methods) that is empty if hasattr(dict[k], '__doc__') and dict[k].__doc__ is None: - for b in mro: - # look if this b defines the method - # (it might not in case of multiple inheritance) - if hasattr(b, k) and (getattr(b, k).__doc__ is not None): - # copy docstring if super method has one - dict[k].__doc__ = getattr(b, k).__doc__ - break + # Some __doc__ strings cannot be rewitten in Python2, e.g. doc of a type + try: + for b in mro: + # look if this b defines the method + # (it might not in case of multiple inheritance) + if hasattr(b, k) and (getattr(b, k).__doc__ is not None): + # copy docstring if super method has one + dict[k].__doc__ = getattr(b, k).__doc__ + break + except AttributeError: + pass return type.__new__(mcs, name, bases, dict) @@ -56,7 +62,7 @@ class MarxsElement(with_metaclass(DocMeta, object)): representation such as a "Rowland Torus". ''' - display = {} + display = {'shape': 'None'} 'Dictionary for display specifications, e.g. color' def __init__(self, **kwargs): @@ -68,6 +74,7 @@ def __init__(self, **kwargs): if len(kwargs) > 0: raise ValueError('Initialization arguments {0} not understood'.format(', '.join(kwargs.keys()))) + self.display = DisplayDict(self, deepcopy(self.display)) def describe(self): return OrderedDict(element=self.name) diff --git a/marxs/design/rowland.py b/marxs/design/rowland.py index 16ed03b..1918c39 100644 --- a/marxs/design/rowland.py +++ b/marxs/design/rowland.py @@ -18,18 +18,18 @@ from __future__ import division -import warnings import numpy as np from scipy import optimize import transforms3d from transforms3d.utils import normalized_vector from ..optics.base import OpticalElement -from ..base import _parse_position_keywords, MarxsElement +from ..base import MarxsElement from ..optics import FlatDetector from ..math.rotations import ex2vec_fix from ..math.utils import e2h, h2e, anglediff from ..simulator import ParallelCalculated +from ..math.geometry import Geometry # Python 2 vs 3 (basestring does not exist in Python 3) try: @@ -68,7 +68,7 @@ def find_radius_of_photon_shell(photons, mirror_shell, x, percentile=[1,99]): wing of the PSF. ''' - p = photons[:] + p = photons.copy() mdet = FlatDetector(position=np.array([x, 0, 0]), zoom=1e8, pixsize=1.) p = mdet(p) ind = (p['probability'] > 0) & (p['mirror_shell'] == mirror_shell) @@ -76,7 +76,7 @@ def find_radius_of_photon_shell(photons, mirror_shell, x, percentile=[1,99]): return np.percentile(r, percentile) -class RowlandTorus(MarxsElement): +class RowlandTorus(MarxsElement, Geometry): '''Torus with y axis as symmetry axis. Note that the origin of the torus is the focal point, which is @@ -101,8 +101,9 @@ class RowlandTorus(MarxsElement): def __init__(self, R, r, **kwargs): self.R = R self.r = r - self.pos4d = _parse_position_keywords(kwargs) - super(RowlandTorus, self).__init__(**kwargs) + # super does not work, because the various classes have different signature + Geometry.__init__(self, kwargs) + MarxsElement.__init__(self, **kwargs) def quartic(self, xyz, transform=True): '''Quartic torus equation. @@ -186,7 +187,7 @@ def f(val_in): raise Exception('Intersection with torus not found.') return val_out - def parametric_surface(self, theta, phi): + def parametric_surface(self, theta, phi, display): '''Parametric representation of surface of torus. In contrast to `parametric` the input parameters here are 1-d arrays @@ -254,7 +255,7 @@ def xyzw2parametric(self, xyzw, transform=True, intersectvalid=True): intersectvalid : bool When ``r >=R`` the torus can intersect with itself. At these points, phi is not unique. If ``intersectvalid`` is true, those points will - be filled with on arbitrarily chosen valid value (0.), otherwise + be filled with an arbitrarily chosen valid value (0.), otherwise they will be nan. Returns @@ -275,9 +276,11 @@ def xyzw2parametric(self, xyzw, transform=True, intersectvalid=True): theta = np.arcsin(s * xyz[:, 1] / self.r) + (s < 0) * np.pi factor = self.R + self.r * np.cos(theta) - with warnings.catch_warnings(): - phi = np.arctan2(xyz[:, 2] / factor, xyz[:, 0] / factor) - phi[factor == 0] = 0 if intersectvalid else np.nan + phi = np.zeros_like(factor) + ind = factor != 0 + phi[ind] = np.arctan2(xyz[ind, 2] / factor[ind], xyz[ind, 0] / factor[ind]) + if not intersectvalid: + phi[~factor] = np.nan return theta, phi diff --git a/marxs/design/tests/test_design.py b/marxs/design/tests/test_design.py index 557076c..415dd33 100644 --- a/marxs/design/tests/test_design.py +++ b/marxs/design/tests/test_design.py @@ -25,9 +25,6 @@ def test_radius_of_photon_shell(): photons = mypointing.process_photons(photons) marxm = MarxMirror('./marxs/optics/hrma.par', position=np.array([0., 0, 0])) photons = marxm(photons) - r1, r2 = find_radius_of_photon_shell(photons, 0, 9e4) - assert abs(r1 - 5380.) < 10. - assert abs(r2 - 5495.) < 10. r1, r2 = find_radius_of_photon_shell(photons, 1, 9e3) assert abs(r1 - 433.) < 1. assert abs(r2 - 442.) < 1. @@ -280,7 +277,7 @@ def test_facet_rotation_via_facetargs(): mygas = GratingArrayStructure(mytorus, d_element=60., x_range=[5e3,1e4], radius=[538., 550.], elem_class=FlatGrating, elem_args={'zoom': 30, 'd':0.0002, 'order_selector': gratingeff}) blaze = transforms3d.axangles.axangle2mat(np.array([0,1,0]), np.deg2rad(15.)) mygascat = GratingArrayStructure(mytorus, d_element=60., x_range=[5e3,1e4], radius=[538., 550.], elem_class=FlatGrating, elem_args={'zoom': 30, 'orientation': blaze, 'd':0.0002, 'order_selector': gratingeff}) - assert np.allclose(np.rad2deg(np.arccos(np.dot(mygas.elements[0].geometry('e_x')[:3], mygascat.elements[0].geometry('e_x')[:3]))), 15.) + assert np.allclose(np.rad2deg(np.arccos(np.dot(mygas.elements[0].geometry['e_x'][:3], mygascat.elements[0].geometry['e_x'][:3]))), 15.) def test_persistent_facetargs(): '''Make sure that facet_args is still intact after generating facets. @@ -350,11 +347,11 @@ def test_LinearCCDArray(): assert len(ccds.elements) == 7 for e in ccds.elements: # For this test we don't care if z and e_z are parallel or antiparallel - assert np.isclose(np.abs(np.dot([0, 0, 1, 0], e.geometry('e_z'))), 1.) + assert np.isclose(np.abs(np.dot([0, 0, 1, 0], e.geometry['e_z'])), 1.) assert (e.pos4d[0, 3] >= 0) and (e.pos4d[0, 3] < 1.) # center ccd is on the optical axis - assert np.allclose(ccds.elements[3].geometry('e_y'), [0, 1, 0, 0]) + assert np.allclose(ccds.elements[3].geometry['e_y'], [0, 1, 0, 0]) def test_LinearCCDArray_rotatated(): '''Test an array with different rotations. @@ -370,7 +367,7 @@ def test_LinearCCDArray_rotatated(): radius=[-100., 100.], phi=0., elem_class=mock_facet) assert len(ccds.elements) == 7 for e in ccds.elements: - assert np.isclose(np.dot([0, -0.8660254, 0.5], e.geometry('e_z')[:3]), 0, atol=1e-4) + assert np.isclose(np.dot([0, -0.8660254, 0.5], e.geometry['e_z'][:3]), 0, atol=1e-4) def test_impossible_LinearCCDArray(): '''The rotation is chosen such that all requested detector positions are @@ -423,6 +420,6 @@ def test_nojumpsinCCDorientation(): theta=[np.pi - 0.2, np.pi - 0.1]) # If all is right, this will vary very smoothly along the circle. - xcomponent = np.array([e.geometry('e_y')[0] for e in det.elements]) + xcomponent = np.array([e.geometry['e_y'][0] for e in det.elements]) xcompdiff = np.diff(xcomponent) assert (np.max(np.abs(xcompdiff)) / np.median(np.abs(xcompdiff))) < 1.1 diff --git a/marxs/math/geometry.py b/marxs/math/geometry.py new file mode 100644 index 0000000..0a5623a --- /dev/null +++ b/marxs/math/geometry.py @@ -0,0 +1,602 @@ +# Licensed under GPL version 3 - see LICENSE.rst +from copy import copy, deepcopy + +import numpy as np +from transforms3d.affines import decompose44, compose +from transforms3d.axangles import axangle2mat + +from .utils import e2h, h2e, anglediff, angle_between, norm_vector +from .pluecker import point_dir2plane +from ..base import _parse_position_keywords +from ..visualization.utils import plane_with_hole, combine_disjoint_triangulations + + +def xyz_square(geometry, r_factor=1): + '''Generate Eukledian positions for the corners of a square. + + The square is centered on the center of the object and the edges are + given by ``v_y`` and ``v_z``. + + Parameters + ---------- + r_factor : float + Scaling factor for the square. + + Returns + ------- + box : np.array of shape (4, 3) + Eukledian coordinates of the corners of the square in 3d space. + ''' + g = geometry + box = h2e(g['center']) + r_factor * np.vstack([h2e( g['v_y']) + h2e(g['v_z']), + h2e(-g['v_y']) + h2e(g['v_z']), + h2e(-g['v_y']) - h2e(g['v_z']), + h2e( g['v_y']) - h2e(g['v_z']) + ]) + return box + +def xyz_circle(geometry, r_factor=1, philim=[0, 2 * np.pi], n_vertices=90): + '''Generate Eukledian positions along an ellipse. + + The circle is centered on the center of the object and the semi-major + and minor axes are given by ``v_y`` and ``v_z``. Note that this function + is usually used to generate circle position, although ellipses are possible, + thus the name. + + The circle (or ellipse) is approximated by a polygon with ``n_vertices`` + vertices, where the value of ``n_vertices`` is taken from the ``self.display`` + dictionary. + + Parameters + ---------- + r_factor : float + Scaling factor for the square. + phi_lim : list + Lower and upper limit for the angle phi to restrict the circle to a wedge. + + Returns + ------- + circle : np.array of shape (n, 3) + Eukledian coordinates of the corners of the square in 3d space. + ''' + n = n_vertices + phi = np.linspace(0.5 * np.pi, 2.5 * np.pi, n, endpoint=False) + v_y = r_factor * geometry['v_y'] + v_z = r_factor * geometry['v_z'] + + x = np.cos(phi) + y = np.sin(phi) + # phi could be less then full circle + # wrap phi at lower bound (which could be negative). + # For the default [0, 2 pi] this is a no-op + phi = (phi - philim[0]) % (2 * np.pi) + ind = phi < anglediff(philim) + x[~ind] = 0 + y[~ind] = 0 + + return h2e(geometry['center'] + x.reshape((-1, 1)) * v_y + y.reshape((-1, 1)) * v_z) + + +class NoGeometry(object): + _geometry = {} + shape = 'none' + + def __init__(self, kwargs={}): + self.geometry = self + + +class Geometry(NoGeometry): + + n_points = 50 + + def __init__(self, kwargs={}): + self.pos4d = _parse_position_keywords(kwargs) + #copy class attribute to instance attribute + self._geometry = copy(self._geometry) + super(Geometry, self).__init__(kwargs=kwargs) + + def __getitem__(self, key): + '''This function wraps access to the pos4d matrix. + + This is mostly a convenience method that gives access to vectors from the + ``pos4d`` matrix in familiar terms with string labels: + + - ``center``: The ``center`` is the origin of the local coordiante system + of the optical elemement. Typically, if will be the center of the + active plane, e.g. the center of the mirror surface. + - :math:`\hat v_{y,z}`: The box stretches in the y and z direction for + :math:`\pm \hat v_y` and :math:`\pm \hat v_z`. In optics, the relevant + interaction often happens on a surface, e.g. the surface of a mirror or + of a reflection grating. In the defaul configuration, this surface is in + the yz-plane. + - :math:`\hat v_x`: The thickness of the box is not as important in many + cases, + but is useful to e.g. render the geometry or to test if elements collide. + The box reaches from the yz-plane down to :math:`- \hat v_x`. Note, that + this definition means the size parallel to the y and z axis is + :math:`2 |\hat v_{y,z}|`, but only :math:`1 |\hat v_{x}|`. + - :math:`\hat e_{x,y,z}`: When an object is initialized, it automatically + adds unit vectors in the + direction of :math:`\hat v_{x,y,z}` called :math:`\hat e_{x,y,z}`. + - ``plane``: It also adds the homogeneous coordinates of the active plane + as ``plane``. + - Other labels get looked up in ``self._geometry`` and if the resulting + value is a 4-d vector, it gets transformed with ``pos4d``. + + Not all these labels make sense for every optical element (e.g. a curved + mirror does not really have a "plane"). + Access through this method is slower than direct indexing of ``self.pos4d``. + ''' + + if key == 'center': + return self.pos4d[:, 3] + elif key in ['v_x', 'e_x']: + val = self.pos4d[:, 0] + elif key in ['v_y', 'e_y']: + val = self.pos4d[:, 1] + elif key in ['v_z', 'e_z']: + val = self.pos4d[:, 2] + elif key == 'plane': + return point_dir2plane(self['center'], self['e_x']) + else: + val = self._geometry[key] + if isinstance(val, np.ndarray) and (val.shape[-1] == 4): + val = np.dot(self.pos4d, val) + if key[:2] == 'e_': + return val / np.linalg.norm(val) + else: + return val + + def intersect(self, dir, pos): + '''Calculate the intersection point between a ray and the element + + Parameters + ---------- + dir : `numpy.ndarray` of shape (N, 4) + homogeneous coordinates of the direction of the ray + pos : `numpy.ndarray` of shape (N, 4) + homogeneous coordinates of a point on the ray + + Returns + ------- + intersect : boolean array of length N + ``True`` if an intersection point is found. + interpos : `numpy.ndarray` of shape (N, 4) + homogeneous coordinates of the intersection point. Values are set + to ``np.nan`` if no intersection point is found. + interpos_local : `numpy.ndarray` of shape (N, 2) + y and z coordinates in the coordiante system of the active plane. + ''' + raise NotImplementedError + + def get_local_euklid_bases(self, interpos_local): + '''Obtain a local eukledian base at a set of positions. + + Parameters + ---------- + interpos_local : `numpy.ndarray` of shape (N, 2) + coordinates in the coordiante system of the geometry (e.g. (x, y), or + (r, phi)). + + Returns + ------- + e_1, e_2, n : `numpy.ndarray` of shape (N, 4) + Vectors pointing in direction 1, 2, and normal to the surface. + ''' + raise NotImplementedError + +class FinitePlane(Geometry): + '''Base class for geometrically flat optical elements. + + ''' + + shape = 'box' + + loc_coos_name = ['y', 'z'] + '''name for output columns that contain the interaction point in local coordinates.''' + + def intersect(self, dir, pos): + '''Calculate the intersection point between a ray and the element + + Parameters + ---------- + dir : `numpy.ndarray` of shape (N, 4) + homogeneous coordinates of the direction of the ray + pos : `numpy.ndarray` of shape (N, 4) + homogeneous coordinates of a point on the ray + + Returns + ------- + intersect : boolean array of length N + ``True`` if an intersection point is found. + interpos : `numpy.ndarray` of shape (N, 4) + homogeneous coordinates of the intersection point. Values are set + to ``np.nan`` if no intersection point is found. + interpos_local : `numpy.ndarray` of shape (N, 2) + y and z coordinates in the coordiante system of the active plane + (not normalized to the dimensions of the element in question, but + in absolute units). + ''' + k_nominator = np.dot(self['center'] - pos, self['e_x']) + k_denom = np.dot(dir, self['e_x']) + + is_parallel = k_denom == 0 + # To avoid warning for parallel rays, which is handeled expicitly below + with np.errstate(divide='ignore'): + k = k_nominator / k_denom + + forward = k >= 0 # with k < 0, rays would have to move backwards. + + if dir.ndim == 2: + interpos = pos + k[:, None] * dir # broadcasting array + else: + interpos = pos + k * dir # dir is a scalar + vec_center_inter = interpos - self['center'] + interpos_local = np.vstack([np.dot(vec_center_inter, self['e_y']), + np.dot(vec_center_inter, self['e_z'])]).T + + intersect = (~is_parallel & forward & + (np.abs(interpos_local[:, 0]) <= np.linalg.norm(self['v_y'])) & + (np.abs(interpos_local[:, 1]) <= np.linalg.norm(self['v_z']))) + for i in [interpos, interpos_local]: + if dir.ndim == 2: + i[~intersect, :3] = np.nan + # input of single photon. + elif (dir.ndim == 1) and not intersect: + i[:3] = np.nan + + return intersect, interpos, interpos_local + + def get_local_euklid_bases(self, interpos_local): + '''Obtain a local eukledian base at a set of positions. + + Parameters + ---------- + interpos_local : `numpy.ndarray` of shape (N, 2) + coordinates in the coordiante system of the geometry (e.g. (x, y), or + (r, phi)). + + Returns + ------- + e_1, e_2, n : `numpy.ndarray` of shape (N, 4) + Vectors pointing in direction 1, 2, and normal to the surface. + ''' + + n = interpos_local.shape[0] + x = np.tile(self['e_x'], (n, 1)) + y = np.tile(self['e_y'], (n, 1)) + z = np.tile(self['e_z'], (n, 1)) + return y, z, x + + +class PlaneWithHole(FinitePlane): + + shape='triangulation' + outer_factor = 3 + inner_factor = 0 + + def __init__(self, kwargs): + self._geometry['r_inner'] = kwargs.pop('r_inner', 0.) + super(PlaneWithHole, self).__init__(kwargs) + + def triangulate(self, display={}): + '''Return a triangulation of the aperture hole embedded in a square. + + The size of the outer square is determined by the ``'outer_factor'`` element + in ``self.display``. + + Returns + ------- + xyz : np.array + Numpy array of vertex positions in Eukeldian space + triangles : np.array + Array of index numbers that define triangles + ''' + outer_disp = self.outer_display(display) + outer_shape = self.outer_shape(display) + + xyz, triangles = plane_with_hole(outer_disp, outer_shape) + + if self['r_inner'] > 0.: + inner_shape = self.inner_shape(display) + # Inner edge of the display. If we have several stacked apertures, + # we don't want to fill is all up to r=0. + inner_disp = self.inner_display(display) + new_xyz, new_tri = plane_with_hole(inner_shape, inner_disp) + xyz, triangles = combine_disjoint_triangulations([xyz, new_xyz], + [triangles, new_tri]) + + return xyz, triangles + + def outer_shape(self, diplay): + raise NotImplementedError + + def outer_display(self, diplay): + raise NotImplementedError + + def inner_shape(self, diplay): + raise NotImplementedError + + def inner_display(self, diplay): + raise NotImplementedError + + +class RectangleHole(PlaneWithHole): + def outer_shape(self, display): + return xyz_square(self) + + def outer_display(self, display): + return xyz_square(self, r_factor=display['outer_factor']) + + +class CircularHole(PlaneWithHole): + + def __init__(self, kwargs): + phi = kwargs.pop('phi', [0, 2. * np.pi]) + if np.max(np.abs(phi)) > 10: + raise ValueError('Input angles >> 2 pi. Did you use degrees (radian expected)?') + if phi[0] > phi[1]: + raise ValueError('phi[1] must be greater than phi[0].') + if (phi[1] - phi[0]) > (2 * np.pi + 1e-6): + raise ValueError('phi[1] - phi[0] must be less than 2 pi.') + self.phi = phi + + super(CircularHole, self).__init__(kwargs) + + def outer_display(self, display): + '''Return values in Eukledian space.''' + return xyz_circle(self, r_factor=display['outer_factor']) + + def outer_shape(self, display): + return xyz_circle(self, r_factor=1, philim=self.phi) + + def inner_shape(self, display): + return xyz_circle(self, + r_factor=self['r_inner']/np.linalg.norm(self['v_y']), + philim=self.phi) + def inner_display(self, display): + # Inner edge of the display. If we have several stacked apertures, + # we don't want to fill is all up to r=0. + return xyz_circle(self, + r_factor=display['inner_factor'] * self['r_inner'] / np.linalg.norm(self['v_y']), + philim=self.phi) + +class Cylinder(Geometry): + '''A Geometry shaped like a ring or tube. + + This object is shaped like a tube. The form is a circle in the xy plane + and and flat along the z-direction. This can be used, for example to + simulate a setup that can follow the Rowland circle + geometry exactly which is useful, e.g. to study the resolution of a + spectrograph without worrying about the details of the detector geometry. + + Parameters + ---------- + position, orientation, zoom, pos4d : see description of `pos4d` + The radius of the tube is given by the ``zoom`` keyword, see `pos4d`. + Use ``zoom[0] == zoom[1]`` to make a circular tube. + ``zoom[0] != zoom[1]`` gives an elliptical profile. + ``zoom[2]`` sets the extension in the z direction. + phi_lim : list + If a cylinder does not cover the full circle, set ``phi_lim`` to the + limits, e.g. ``[-np.pi / 2, np.pi / 2]`` makes a "half-pipe". + ''' + loc_coos_name = ['phi', 'z'] + + shape = 'surface' + coos_limits = [np.array([-np.pi, np.pi]), np.array([-1, 1])] + + # This verification code does not fit here right now, but when I convert to traits + # later it might come in useful again and thus I keep it here as a comment. + # @property + # def phi_lim(self): + # return self._phi_lim + + # @phi_lim.setter + # def phi_lim(self, value): + # if (len(value) != 2) or (value[0] > value[1]): + # raise ValueError('phi_lim has the format [lower limit, upper limit]') + # for v in value: + # if (v < -np.pi) or (v > np.pi): + # raise ValueError('phi_lim must be in range -pi to +pi.') + # self.coos_limits[0] = value + + def __init__(self, kwargs={}): + self.coos_limits = deepcopy(self.coos_limits) + self.coos_limits[0] = np.asanyarray(kwargs.pop('phi_lim', [-np.pi, np.pi])) + super(Cylinder, self).__init__(kwargs) + + def __getitem__(self, value): + if value == 'R': + trans, rot, zoom, shear = decompose44(self.pos4d) + return zoom[0] + else: + return super(Cylinder, self).__getitem__(value) + + @classmethod + def from_rowland(cls, rowland, width, rotation=0., kwargs={}): + '''Generate a `Cylinder` from a `RowlandTorus`. + + According to the definition of the `marxs.design.rowland.RowlandTorus` + the origin phi=0 is at the "top". When this class method is used to + make a detector that catches all dispersed grating signal on the + Rowland torus, a ``rotation=np.pi`` places the center of the Cylinder + close to the center of the torus (the location of the focal point in + the standard Rowland geometry). + + Parameters + ---------- + rowland : `~marxs.design.RowlandTorus` + The circular detector is constructed to fit exactly into the + Rowland Circle defined by ``rowland``. + width : float + Half-width of the tube in the flat direction (z-axis) in mm + rotation : float + Rotation angle of the Cylinder around its z-axis compared to the + phi=0 position of the Rowland torus. + + ''' + # Step 1: Rotate around z axis + rot = axangle2mat(np.array([0, 0, 1.]), rotation) + # Step 2: Get position and size from Rowland torus + pos4d_circ = compose([rowland.R, 0, 0], rot, [rowland.r, rowland.r, width]) + # Step 3: Transform to global coordinate system + pos4d_circ = np.dot(rowland.pos4d, pos4d_circ) + # Step 4: Make detector + det_kwargs = {'pos4d': pos4d_circ} + det_kwargs.update(kwargs) + return cls(det_kwargs) + + + def intersect(self, dir, pos, transform=True): + '''Calculate the intersection point between a ray and the element + + Parameters + ---------- + dir : `numpy.ndarray` of shape (N, 4) + homogeneous coordinates of the direction of the ray + pos : `numpy.ndarray` of shape (N, 4) + homogeneous coordinates of a point on the ray + transform : bool + If ``True``, input is in global coordinates and needs to be transformed + here for the calculations; if ``False`` input is in local coordinates. + + Returns + ------- + intersect : boolean array of length N + ``True`` if an intersection point is found. + interpos : `numpy.ndarray` of shape (N, 4) + homogeneous coordinates of the intersection point. Values are set + to ``np.nan`` is no intersecton point is found. + interpos_local : `numpy.ndarray` of shape (N, 2) + phi, z coordiantes (in the local frame) for one of the intersection points. + If both intersection points are required, reset ``self.inner`` and call this + function again. + ''' + # This could be moved to a general function + if not np.all(dir[:, 3] == 0): + raise ValueError('First input must be direction vectors.') + # Could test pos, too... + if transform: + invpos4d = np.linalg.inv(self.pos4d) + dir = np.dot(invpos4d, dir.T).T + pos = np.dot(invpos4d, pos.T).T + + xyz = h2e(pos) + dir_e = h2e(dir) + + # Solve quadratic equation in steps. a12 = (-xr +- sqrt(xr - r**2(x**2 - R**2))) + xy = xyz[:, :2] + r = dir[:, :2] + c = np.sum(xy**2, axis=1) - 1. + b = 2 * np.sum(xy * r, axis=1) + a = np.sum(r**2, axis=1) + underroot = b**2 - 4 * a * c + # List of intersect in xy plane. + intersect = (underroot >= 0) + i = intersect # just a shorthand because it's used so much below + + interpos_local = np.ones((pos.shape[0], 2)) + interpos_local[:] = np.nan + interpos = np.ones_like(pos) + interpos[:] = np.nan + + if intersect.sum() > 0: + i_ind = intersect.nonzero()[0] + denom = 2 * a[i] + a1 = (- b[i] + np.sqrt(underroot[i])) / denom + a2 = (- b[i] - np.sqrt(underroot[i])) / denom + xy_1 = xy[i, :] + a1[:, np.newaxis] * r[i, :] + phi_1 = np.arctan2(xy_1[:, 1], xy_1[:, 0]) + xy_2 = xy[i, :] + a2[:, np.newaxis] * r[i, :] + phi_2 = np.arctan2(xy_2[:, 1], xy_2[:, 0]) + # 1, 2 look like hits in x,y but might still miss in z + z_1 = xyz[i, 2] + a1 * dir[i, 2] + z_2 = xyz[i, 2] + a2 * dir[i, 2] + hit_1 = ((a1 >= 0) & (np.abs(z_1) <= 1.) & + angle_between(phi_1, self.coos_limits[0][0], self.coos_limits[0][1])) + hit_2 = ((a2 >= 0) & (np.abs(z_2) <= 1.) & + angle_between(phi_2, self.coos_limits[0][0], self.coos_limits[0][1])) + # If both 1 and 2 are hits, use the closer one + hit_1[hit_2 & (a2 < a1)] = False + hit_2[hit_1 & (a2 >= a1)] = False + intersect[i_ind] = hit_1 | hit_2 + # Set values into array from either point 1 or 2 + interpos_local[i_ind[hit_1], 0] = phi_1[hit_1] + interpos_local[i_ind[hit_1], 1] = z_1[hit_1] + + interpos_local[i_ind[hit_2], 0] = phi_2[hit_2] + interpos_local[i_ind[hit_2], 1] = z_2[hit_2] + # Calculate pos for point 1 or 2 in local xyz coord system + interpos[i_ind[hit_1], :] = e2h(xyz[i_ind, :] + a1[:, None] * dir_e[i_ind, :], 1)[hit_1, :] + interpos[i_ind[hit_2], :] = e2h(xyz[i_ind, :] + a2[:, None] * dir_e[i_ind, :], 1)[hit_2, :] + + + trans, rot, zoom, shear = decompose44(self.pos4d) + # interpos_local in z direction is in local coordinates, i.e. + # the x coordiante is 0..1, but we want that in units of the + # global coordinate system. + interpos_local[:, 1] = interpos_local[:, 1] * zoom[2] + interpos = np.dot(self.pos4d, interpos.T).T + + return intersect, interpos, interpos_local + + + def parametric_surface(self, phi=None, z=None, display={}): + '''Parametric description of the tube. + + This is just another way to obtain the shape of the tube, e.g. + for visualization. + + Parameters + ---------- + phi : np.array + ``phi`` is the angle around the tube profile. Set to ``None`` to use the + extend of the element itself. + z : np.array + The coordiantes along the radius coordinate. Set to ``None`` to use the + extend of the element itself. + + Returns + ------- + xyzw : np.array + Ring coordinates in global homogeneous coordinate system. + ''' + phi = np.linspace(self.coos_limits[0][0], self.coos_limits[0][1], self.n_points) \ + if phi is None else np.asanyarray(phi) + trans, rot, zoom, shear = decompose44(self.pos4d) + z = self.coos_limits[1] if z is None else np.asanyarray(z) / zoom[2] + if (phi.ndim != 1) or (z.ndim != 1): + raise ValueError('input parameters have 1-dim shape.') + phi, z = np.meshgrid(phi, z) + x = np.cos(phi) + y = np.sin(phi) + w = np.ones_like(z) + coos = np.array([x, y, z, w]).T + return np.einsum('...ij,...j', self.pos4d, coos) + + def get_local_euklid_bases(self, interpos_local): + '''Obtain a local eukledian base at a set of positions. + + Parameters + ---------- + interpos_local : `numpy.ndarray` of shape (N, 2) + coordinates in the coordiante system of the geometry (e.g. (x, y), or + (r, phi)). + + Returns + ------- + e_1, e_2, n : `numpy.ndarray` of shape (N, 4) + Vectors pointing in direction 1, 2, and normal to the surface. + ''' + + phi = interpos_local[:, 0] + zeros = np.zeros_like(phi) + e_phi = np.array([-np.sin(phi), np.cos(phi), zeros, zeros]).T + e_z = np.tile(self['e_z'], (interpos_local.shape[0], 1)) + e_n = np.array([np.cos(phi), np.sin(phi), zeros, zeros]).T + + e_phi = np.einsum('...ij,...j', self.pos4d, e_phi) + e_n = np.einsum('...ij,...j', self.pos4d, e_n) + # e_z is already in the local coordinate system + return norm_vector(e_phi), e_z, norm_vector(e_n) diff --git a/marxs/math/pluecker.py b/marxs/math/pluecker.py index 3b9c0dd..302448f 100644 --- a/marxs/math/pluecker.py +++ b/marxs/math/pluecker.py @@ -74,7 +74,7 @@ def intersect_line_plane(p_line, h_plane): '''Intersect a line in Pluecker coordinates with a plane in homogeneous coordinates. This function returns the position of the intersection. - If the line is parallel (but not idendical) to the plane that intersection + If the line is parallel (but not identical) to the plane that intersection happens "at infinity" and thus the last component of the intersection coordinate is 0. If the line is in the plane, then all components of the return vector are 0. diff --git a/marxs/math/polarization.py b/marxs/math/polarization.py index 429dc62..e6ce571 100644 --- a/marxs/math/polarization.py +++ b/marxs/math/polarization.py @@ -1,6 +1,7 @@ # Licensed under GPL version 3 - see LICENSE.rst import numpy as np import astropy.units as u +from numpy.core.umath_tests import inner1d from .utils import norm_vector, e2h @@ -110,9 +111,10 @@ def paralleltransport_matrix(dir1, dir2, jones=np.eye(2), replace_nans=True): replace_nans : bool If ``True`` return an identity matrix for those rays with ``dir1=dir2``. In those cases, the local coordinate system is not well - defines and thus no Jones matrix can be applied. In MARXS ``dir1=dir2`` + defined and thus no Jones matrix can be applied. In MARXS ``dir1=dir2`` often happens if some photons in a list miss the optical element in - question - these photons just pass through. + question - these photons just pass through and their polarization vector + should be unchanged. Returns ------- @@ -120,21 +122,31 @@ def paralleltransport_matrix(dir1, dir2, jones=np.eye(2), replace_nans=True): ''' dir1 = norm_vector(dir1) dir2 = norm_vector(dir2) + jones_3 = np.eye(3) jones_3[:2, :2] = jones + + pmat = np.zeros((dir1.shape[0], 3, 3)) s = np.cross(dir1, dir2) - s = s / np.linalg.norm(s, axis=1)[:, None] - p_in = np.cross(dir1, s) - p_out = np.cross(dir2, s) - - Oininv = np.array([s, p_in, dir1]).swapaxes(1, 0) - Oout = np.array([s, p_out, dir2]).swapaxes(1, 2).T - temp = np.einsum('...ij,kjl->kil', jones_3, Oininv) - pmat = np.einsum('ijk,ikl->ijl', Oout, temp) - - if replace_nans: - ind = np.isnan(s[:, 0]) - pmat[ind, :, :] = np.eye(3)[None, :, :] + s_norm = np.linalg.norm(s, axis=1) + # Find dir values that remain unchanged + # For these the cross prodict will by 0 + # and a numerical error is raised in s / norm(s) + # Expected output value for these depends on "replace_nans" + ind = np.isclose(s_norm, 0) + + if (~ind).sum() > 0: + + s = s[~ind, :] / s_norm[~ind][:, None] + p_in = np.cross(dir1[~ind, :], s) + p_out = np.cross(dir2[~ind, :], s) + Oininv = np.array([s, p_in, dir1[~ind, :]]).swapaxes(1, 0) + Oout = np.array([s, p_out, dir2[~ind, :]]).swapaxes(1, 2).T + temp = np.einsum('...ij,kjl->kil', jones_3, Oininv) + pmat[~ind, :, :] = np.einsum('ijk,ikl->ijl', Oout, temp) + + factor = 1 if replace_nans else np.nan + pmat[ind, :, :] = factor * np.eye(3)[None, :, :] return pmat def parallel_transport(dir_old, dir_new, pol_old, **kwargs): diff --git a/marxs/math/tests/test_geometry.py b/marxs/math/tests/test_geometry.py new file mode 100644 index 0000000..b5544dd --- /dev/null +++ b/marxs/math/tests/test_geometry.py @@ -0,0 +1,244 @@ +# Licensed under GPL version 3 - see LICENSE.rst +import numpy as np +from numpy.core.umath_tests import inner1d +import transforms3d +import transforms3d.euler +import pytest + +from ..geometry import Cylinder, FinitePlane +from ...math.utils import h2e + + +def test_intersect_plane_forwards(): + + pos = np.array([[20., 0., 0.8, 1.], # start above + [0., .0, 0.4, 1.], # start inside + [-5., .1, -.2, 1]]) # start below + dir = np.array([[-.1, 0., 0., 0], + [-1., 0, 0.02, 0.], + [-2., -0.02, 0.023, 0.]]) + + plane = FinitePlane() + intersect, interpos, inter_local = plane.intersect(dir, pos) + assert np.all(intersect == [True, True, False]) + intersect, interpos, inter_local = plane.intersect(-dir, pos) + assert np.all(intersect == [False, True, True]) + +def test_intersect_plane_zoomed(): + pos = np.array([[1., 4., 0.8, 1.], + [1., -4., 8., 1.], + [1., -6., -8., 1]]) + dir = np.array([[-.1, 0., 0., 0], + [-1., 0, 0.02, 0.], + [-2., -0.02, 0.023, 0.]]) + + plane = FinitePlane({'zoom': [1, 5, 10]}) + intersect, interpos, inter_local = plane.intersect(dir, pos) + assert np.all(intersect == [True, True, False]) + assert np.all(interpos[0, :] == [0, 4, .8, 1]) + assert np.all(inter_local[0, ] == [4, .8]) + +def test_intersect_plane_moved(): + + pos = np.array([[4.1, 5.5, 6.8, 1.], + [5., 4.1, 4.99, 1.], + [4.1, 3., 6., 1]]) + dir = np.array([[-.1, 0., 0., 0], + [-1., 0, 0.02, 0.], + [-2., -0.02, 0.023, 0.]]) + + plane = FinitePlane({'position': [4,5,6]}) + intersect, interpos, inter_local = plane.intersect(dir, pos) + assert np.all(intersect == [True, True, False]) + assert np.all(interpos[1, :] == [4, 4.1, 5.01, 1]) + assert np.allclose(inter_local[1, :], [-0.9, -.99]) + +def test_intersect_tube_miss(): + '''Rays passing at larger radius or through the center miss the tube.''' + circ = Cylinder() + # Passing at larger radius + intersect, interpos, inter_local = circ.intersect(np.array([[1., 0., .1, 0], + [-1., -1., 0., 0.]]), + np.array([[0, -1.1, 0., 1.], + [2., 0., 0., 1.]])) + assert np.all(intersect == False) + assert np.all(np.isnan(interpos)) + assert np.all(np.isnan(inter_local)) + + # Through the center almost parallel to z-axis + intersect, interpos, inter_local = circ.intersect(np.array([[0., 0.1, 1., 0]]), + np.array([[0.5, 0.5, 0., 1.]])) + assert np.all(intersect == False) + assert np.all(np.isnan(interpos)) + assert np.all(np.isnan(inter_local)) + + # Repeat with a tube that's moved to make sure we did not mess up local and + # global coordinates in the implementation + circ = Cylinder({'position': [0.8, 0.8, 1.2]}) + intersect, interpos, inter_local = circ.intersect(np.array([[1., 0., .0, 0], + [1., 0., 0., 0.]]), + np.array([[0., 0., 0., 1.], + [0., -0.3, 1., 1.]])) + assert np.all(intersect == False) + assert np.all(np.isnan(interpos)) + assert np.all(np.isnan(inter_local)) + +def test_intersect_tube_rays_forward(): + '''Rays only intersect if propagating them forward works''' + pos = np.array([[20., 0., 0.8, 1.], # start above + [0., .0, 0.4, 1.], # start inside + [-5., .1, -.2, 1]]) # start below + dir = np.array([[-.1, 0., 0., 0], + [-1., 0, 0.02, 0.], + [-2., -0.02, 0.023, 0.]]) + + circ = Cylinder({'zoom': [1, 1, 1]}) + intersect, interpos, inter_local = circ.intersect(dir, pos) + assert np.all(intersect == [True, True, False]) + intersect, interpos, inter_local = circ.intersect(-dir, pos) + assert np.all(intersect == [False, True, True]) + + +def test_intersect_tube_hitmiss_zoomed(): + '''Rays that hit one tube, but miss a shorter one''' + dir = np.array([[-.1, 0., 0., 0], [-1., 0, 0.2, 0.]]) + pos = np.array([[20., 0., 0.8, 1.], [2., .0, 0.4, 1.]]) + + circ = Cylinder({'zoom': [1, 1, 1]}) + intersect, interpos, inter_local = circ.intersect(dir, pos) + assert np.all(intersect == True) + + circ = Cylinder({'zoom': [1, 1, .5]}) + intersect, interpos, inter_local = circ.intersect(dir, pos) + assert np.all(intersect == False) + assert np.all(np.isnan(interpos)) + assert np.all(np.isnan(inter_local)) + +def test_intersect_tube_2points(): + dir = np.array([[-.1, 0., 0., 0], [-0.5, 0, 0., 0.], [-1., 0., 0., 0.]]) + pos = np.array([[50., 0., 0., 1.], [10., .5, 0., 1.], [2, 0, .3, 1.]]) + + circ = Cylinder() + intersect, interpos, inter_local = circ.intersect(dir, pos) + assert np.all(intersect == True) + assert np.allclose(h2e(interpos), np.array([[1., 0., 0.], + [np.sqrt(0.75), 0.5, 0.], + [1, 0., .3]])) + assert np.allclose(inter_local, np.array([[0., np.arcsin(0.5), 0.], + [0., 0., 0.3]]).T) + +def test_intersect_tube_2points_transformed(): + '''inter_local should be the same as interpos, + even if we make the tube longer in z''' + pos = np.array([[50., 0., 0., 1.], [10., .5, 0., 1.], [2, 0, .3, 1.]]) + dir = np.array([[-.1, 0., 0., 0], [-0.5, 0, 0., 0.], [-1., 0., 0., 0.]]) + + exp_inter_loc = np.array([[0., np.arcsin(0.5), 0.], + [0., 0., 0.3]]).T + exp_interpos = np.array([[1., 0., 0.], + [np.sqrt(0.75), 0.5, 0.], + [1, 0., .3]]) + circ = Cylinder({'zoom': [1, 1., 5.]}) + intersect, interpos, inter_local = circ.intersect(dir, pos) + assert np.all(intersect == True) + assert np.allclose(h2e(interpos), exp_interpos) + assert np.allclose(inter_local, exp_inter_loc) + # Rotate: interpos should be the same as in the unrotated case. + # Beware of rotation direction! + orient = transforms3d.axangles.axangle2mat([0, 0, 1], -.3) + circ = Cylinder({'orientation': orient}) + intersect, interpos, inter_local = circ.intersect(dir, pos) + assert np.all(intersect == True) + assert np.allclose(h2e(interpos), exp_interpos) + assert np.allclose(inter_local, np.array([[0.3, 0.3 + np.arcsin(0.5), 0.3], + [0., 0., 0.3]]).T) + + # Shift: inter_local is the same, but interpos changes + circ = Cylinder({'position': [-.3, 0, 0]}) + intersect, interpos, inter_local = circ.intersect(dir, pos) + assert np.all(intersect == True) + assert np.allclose(h2e(interpos), np.array([[.7, 0, 0], + [np.sqrt(0.75) -.3, 0.5, 0.], + [.7, 0., .3]])) + assert np.allclose(inter_local, exp_inter_loc) + +def test_intersect_tube_2points_outside(): + '''Setting phi_lim + + These are the same numbers, but with the phi_lim set, all intersection are + at the "bottom" of the cylinder.''' + pos = np.array([[50., 0., 0., 1.], [10., .5, 0., 1.], [2, 0, .3, 1.]]) + dir = np.array([[-.1, 0., 0., 0], [-0.5, 0, 0., 0.], [-1., 0., 0., 0.]]) + + circ = Cylinder({'phi_lim': [.6, np.pi]}) + intersect, interpos, inter_local = circ.intersect(dir, pos) + assert np.all(intersect == True) + assert np.allclose(h2e(interpos), np.array([[-1., 0., 0.], + [-np.sqrt(0.75), 0.5, 0.], + [-1, 0., .3]])) + assert np.allclose(inter_local, np.array([[np.pi, np.pi - np.arcsin(0.5), np.pi], + [0., 0., 0.3]]).T) + # phi_lims that guarantee a miss + circ = Cylinder({'phi_lim': [-.4, -.3]}) + intersect, interpos, inter_local = circ.intersect(dir, pos) + assert np.all(intersect == False) + assert np.all(np.isnan(interpos)) + assert np.all(np.isnan(inter_local)) + +@pytest.mark.xfail +def test_phi_lim_verification(): + '''Check Errors for wrong phi_lim format.''' + with pytest.raises(ValueError) as e: + circ = Cylinder({'phi_lim': [-.4, -.3, .4]}) + assert '[lower limit, upper limit]' in str(e.value) + + with pytest.raises(ValueError) as e: + circ = Cylinder({'phi_lim': [-.4, -.5]}) + assert '[lower limit, upper limit]' in str(e.value) + + with pytest.raises(ValueError) as e: + circ = Cylinder({'phi_lim': [-.4, 5]}) + assert 'range -pi to +pi' in str(e.value) + + +def test_intersect_tube_2points_translation(): + '''Repeat test above with a tube that's moved and zoomed to make sure we + did not mess up local and global coordinates in the implementation + ''' + circ = Cylinder({'position': [0.8, 0.8, 1.2], 'zoom': [1, 1, 2]}) + intersect, interpos, inter_local = circ.intersect(np.array([[1., 0., .0, 0], + [1., 0., 0., 0.]]), + np.array([[0., 0., 0., 1.], + [0., 0., 1., 1.]])) + assert np.all(intersect == True) + assert np.allclose(h2e(interpos), np.array([[.2, 0., 0.], + [.2, 0., 1.]])) + +def test_tube_parametric(): + '''Generate points on surface using parametric. Then make rays for intersection + which should return those points as intersection points.''' + circ = Cylinder({'zoom': [2., 3., 4.]}) + parametric = circ.parametric_surface(phi=[-.1, 0., 1.]) + parametric = parametric.reshape((6, 4)) + # select dir so that it's in the right direction to recover these points. + dir = np.tile([1, 0, 0, 0], (6, 1)) + intersect, interpos, inter_local = circ.intersect(dir, parametric) + assert np.allclose(parametric, interpos) + +@pytest.mark.parametrize("geom", [FinitePlane, Cylinder]) +def test_local_coordsys(geom): + '''Ensure local coordiante systems are orthonormal''' + rot = transforms3d.euler.euler2mat(*(np.pi * 2 * np.random.rand(3))) + g = geom({'rotation': rot, + 'position': np.random.rand(3)}) + + x, y, z = g.get_local_euklid_bases(np.random.rand(5, 2)) + + assert np.allclose(inner1d(x, y), 0) + assert np.allclose(inner1d(x, z), 0) + + for vec in [x, y, z]: + assert np.allclose(np.linalg.norm(vec, axis=1), 1.) + + # Check it's a right-handed coordinage system + assert np.allclose(np.cross(x[:, :3], y[:, :3]), z[:, :3]) diff --git a/marxs/math/tests/test_utils.py b/marxs/math/tests/test_utils.py index 5174100..42f1ba2 100644 --- a/marxs/math/tests/test_utils.py +++ b/marxs/math/tests/test_utils.py @@ -42,3 +42,11 @@ def test_anglediff(): assert np.isclose(utils.anglediff([0, 3.]), 3.) assert np.isclose(utils.anglediff([-1, 1]), 2.) assert np.isclose(utils.anglediff([1., -1.]), 2 * np.pi - 2.) + +def test_angle_between(): + assert utils.angle_between(.1, 0, 2.) + assert np.all(utils.angle_between(np.array([-.1, 0., 6.1, .2, 6.5]), 5., 1.)) + assert utils.angle_between(2., 6., 1.) == False + assert utils.angle_between(4., 1., 2.) == False + assert utils.angle_between(-.5, .5, -3.1) == False + assert utils.angle_between(0, -np.pi, np.pi) == True diff --git a/marxs/math/utils.py b/marxs/math/utils.py index 15da884..50d7a5b 100644 --- a/marxs/math/utils.py +++ b/marxs/math/utils.py @@ -156,7 +156,6 @@ def norm_vector(vec): length2 = np.sum(vec * vec, axis=-1) return vec / np.sqrt(length2)[:, None] - def anglediff(phi): '''Angle range covered by phi, accounting for 2 pi properly @@ -170,3 +169,40 @@ def anglediff(phi): # If anglediff == 2 pi exactly, presumably the user want to cover the full circle. anglediff = anglediff % (2. * np.pi) return anglediff + +def angle_between(angle, border1, border2): + '''Test if an angle is between two others + + Since angles are cyclic, a simple numerical comparison is not + sufficient, e.g. 355 deg is in the interval [350 deg, 10 deg] even though + numerically 355 is not less then 10. + + Parameters + ---------- + angle : float or np.array + Angle array (in rad) + border1 : float + Lower border of interval (inclusive) in radian + border2 : float + Higher border of interval (inclusive) in radian + + Returns + ------- + comparison : bool of same shape as ``angle`` + Result of the comparison + + Example + ------- + + >>> from marxs.math.utils import angle_between + >>> angle_between(-0.1, -0.2, 6) + True + ''' + twopi = 2 * np.pi + b1 = (twopi + (border1 % twopi)) % twopi + b2 = (twopi + (border2 % twopi)) % twopi + ang = (twopi + (angle % twopi)) % twopi + if b1 < b2: + return (b1 <= ang) & (ang <= b2) + else: + return (b1 <= ang) | (ang <= b2) diff --git a/marxs/missions/chandra/acis.py b/marxs/missions/chandra/acis.py index b3667f0..ec49838 100644 --- a/marxs/missions/chandra/acis.py +++ b/marxs/missions/chandra/acis.py @@ -22,6 +22,7 @@ def __init__(self, **kwargs): self.TDET = TDET['ACIS'] self.ODET = ODET['ACIS'] self.pixsize_in_rad = np.deg2rad(PIXSIZE['ACIS']) + kwargs['ignore_pixel_warning'] = True super(ACISChip, self).__init__(**kwargs) @property diff --git a/marxs/missions/chandra/tests/test_hetg.py b/marxs/missions/chandra/tests/test_hetg.py index 67c6bbb..e871ac4 100644 --- a/marxs/missions/chandra/tests/test_hetg.py +++ b/marxs/missions/chandra/tests/test_hetg.py @@ -7,9 +7,9 @@ def test_orient(): '''Check orientation of gratings.''' hetg = HETG() for i in range(len(hetg.elements)): - ex = hetg.elements[i].geometry('e_x')[:3] - ey = hetg.elements[i].geometry('e_y')[:3] - ez = hetg.elements[i].geometry('e_z')[:3] + ex = hetg.elements[i].geometry['e_x'][:3] + ey = hetg.elements[i].geometry['e_y'][:3] + ez = hetg.elements[i].geometry['e_z'][:3] x = np.array([hetg.hess['xu'][i], hetg.hess['yu'][i], hetg.hess['zu'][i]]) y = np.array([hetg.hess['xuxf'][i], hetg.hess['yuxf'][i], hetg.hess['zuxf'][i]]) z = np.array([hetg.hess['xuyf'][i], hetg.hess['yuyf'][i], hetg.hess['zuyf'][i]]) @@ -26,10 +26,9 @@ def test_groove_dir(): '''Regression test: Groove direction was not set correctly from input table.''' hetg = HETG() for i in range(len(hetg.elements)): - a = hetg.elements[i].geometry('e_groove')[:3] - b = np.array([hetg.hess['xul'][i], hetg.hess['yul'][i], hetg.hess['zul'][i]]) - c = hetg.elements[i].geometry('e_perp_groove')[:3] - d = np.array([hetg.hess['xud'][i], hetg.hess['yud'][i], hetg.hess['zud'][i]]) + a, c, n = hetg.elements[i].e_groove_coos(np.zeros((1, 2))) + b = np.array([hetg.hess['xul'][i], hetg.hess['yul'][i], hetg.hess['zul'][i], 0]) + d = np.array([hetg.hess['xud'][i], hetg.hess['yud'][i], hetg.hess['zud'][i], 0]) # all a normalized for x in [a, b, c, d]: assert np.allclose(np.linalg.norm(x), 1., atol=1e-4) diff --git a/marxs/optics/aperture.py b/marxs/optics/aperture.py index 378f434..3d99125 100644 --- a/marxs/optics/aperture.py +++ b/marxs/optics/aperture.py @@ -6,8 +6,7 @@ from .base import FlatOpticalElement from ..base import GeometryError -from ..visualization.utils import plane_with_hole, combine_disjoint_triangulations -from ..math.utils import anglediff, h2e +from ..math.geometry import RectangleHole, CircularHole from ..simulator import BaseContainer from .. import utils @@ -55,123 +54,38 @@ def generate_local_xy(self, n): return NotImplementedError def __call__(self, photons): - self.add_output_cols(photons, self.loc_coos_name) + self.add_output_cols(photons, self.geometry.loc_coos_name) # Add ID number to ID col, if requested if self.id_col is not None: photons[self.id_col] = self.id_num # Set position in different coordinate systems x, y = self.generate_local_xy(len(photons)) - if self.loc_coos_name is not None: - photons[self.loc_coos_name[0]] = x - photons[self.loc_coos_name[1]] = y - photons['pos'] = self.geometry('center') + x.reshape((-1, 1)) * self.geometry('v_y') + y.reshape((-1, 1)) * self.geometry('v_z') - projected_area = np.dot(photons['dir'].data, - self.geometry('e_x')) + if self.geometry.loc_coos_name is not None: + photons[self.geometry.loc_coos_name[0]][:] = x + photons[self.geometry.loc_coos_name[1]][:] = y + photons['pos'] = self.geometry['center'] + x.reshape((-1, 1)) * self.geometry['v_y'] + y.reshape((-1, 1)) * self.geometry['v_z'] + projected_area = np.dot(photons['dir'].data, - self.geometry['e_x']) # Photons coming in "through the back" would have negative probabilities. # Unlikely to ever come up, but just in case we clip to 0. - photons['probability'] *= np.clip(projected_area, 0, 1.) + photons['probability'][:] *= np.clip(projected_area, 0, 1.) return photons def process_photons(self, photons): raise NotImplementedError('You probably want to use __call__.') - def inner_shape(self): - '''Return values in Eukledean space''' - raise NotImplementedError - - def xyz_square(self, r_factor): - '''Generate Eukledian positions for the corners of a square. - - The square is centered on the center of the object and the edges are - given by ``v_y`` and ``v_z``. - - Parameters - ---------- - r_factor : float - Scaling factor for the square. - - Returns - ------- - box : np.array of shape (4, 3) - Eukledian coordinates of the corners of the square in 3d space. - ''' - g = self.geometry - box = h2e(g('center')) + r_factor * np.vstack([h2e( g('v_y')) + h2e(g('v_z')), - h2e(-g('v_y')) + h2e(g('v_z')), - h2e(-g('v_y')) - h2e(g('v_z')), - h2e( g('v_y')) - h2e(g('v_z')) - ]) - return box - - def xyz_circle(self, r_factor, philim=[0, 2 * np.pi]): - '''Generate Eukledian positions along an ellipse. - - The circle is centered on the center of the object and the semi-major - and minor axes are given by ``v_y`` and ``v_z``. Note that this function - is usually used to generate circle position, although ellipses are possible, - thus the name. - - The circle (or ellipse) is approximated by a polygon with ``n_vertices`` - vertices, where the value of ``n_vertices`` is taken from the ``self.display`` - dictionary. - - Parameters - ---------- - r_factor : float - Scaling factor for the square. - phi_lim : list - Lower and upper limit for the angle phi to restrict the circle to a wedge. - - Returns - ------- - circle : np.array of shape (n, 3) - Eukledian coordinates of the corners of the square in 3d space. - ''' - - n = self.display.get('n_vertices', 90) - phi = np.linspace(0.5 * np.pi, 2.5 * np.pi, n, endpoint=False) - v_y = r_factor * self.geometry('v_y') - v_z = r_factor * self.geometry('v_z') - - x = np.cos(phi) - y = np.sin(phi) - # phi could be less then full circle - # wrap phi at lower bound (which could be negative). - # For the default [0, 2 pi] this is a no-op - phi = (phi - philim[0]) % (2 * np.pi) - ind = phi < anglediff(philim) - x[~ind] = 0 - y[~ind] = 0 - - return h2e(self.geometry('center') + x.reshape((-1, 1)) * v_y + y.reshape((-1, 1)) * v_z) - - def outer_shape(self): '''Return values in Eukledian space.''' return self.xyz_square(self.display.get('outer_factor', 3)) - def triangulate(self): - '''Return a triangulation of the aperture hole embedded in a square. - - The size of the outer square is determined by the ``'outer_factor'`` element - in ``self.display``. - - Returns - ------- - xyz : np.array - Numpy array of vertex positions in Eukeldian space - triangles : np.array - Array of index numbers that define triangles - ''' - outer = self.outer_shape() - inner = self.inner_shape() - return plane_with_hole(outer, inner) - class RectangleAperture(FlatAperture): '''Select the position where a parallel ray from an astrophysical source starts the simulation. ''' + + default_geometry = RectangleHole + def generate_local_xy(self, n): x = np.random.random(n) * 2. - 1. y = np.random.random(n) * 2. - 1. @@ -180,15 +94,7 @@ def generate_local_xy(self, n): @property def area(self): '''Area covered by the aperture''' - return 4 * np.linalg.norm(self.geometry('v_y')) * np.linalg.norm(self.geometry('v_z')) * u.mm**2 - - def inner_shape(self): - g = self.geometry - return h2e(g('center')) + np.vstack([h2e( g('v_y')) + h2e(g('v_z')), - h2e(-g('v_y')) + h2e(g('v_z')), - h2e(-g('v_y')) - h2e(g('v_z')), - h2e( g('v_y')) - h2e(g('v_z'))]) - + return 4 * np.linalg.norm(self.geometry['v_y']) * np.linalg.norm(self.geometry['v_z']) * u.mm**2 class CircleAperture(FlatAperture): '''Select the position where a parallel ray from an astrophysical source starts the simulation. @@ -212,28 +118,22 @@ class CircleAperture(FlatAperture): ``self.display['inner_factor']`` can be used to restrict the radius range where the inner disk is displayed in a plot. ''' + + default_geometry = CircularHole + def __init__(self, **kwargs): - phi = kwargs.pop('phi', [0, 2. * np.pi]) - if np.max(np.abs(phi)) > 10: - raise ValueError('Input angles >> 2 pi. Did you use degrees (radian expected)?') - if phi[0] > phi[1]: - raise ValueError('phi[1] must be greater than phi[0].') - if (phi[1] - phi[0]) > (2 * np.pi + 1e-6): - raise ValueError('phi[1] - phi[0] must be less than 2 pi.') - self.phi = phi - self.r_inner = kwargs.pop('r_inner', 0.) super(CircleAperture, self).__init__(**kwargs) - if self.r_inner > np.linalg.norm(self.geometry('v_y')): + if self.geometry['r_inner'] > np.linalg.norm(self.geometry['v_y']): raise ValueError('r_inner must be less than size of full aperture.') - if not np.isclose(np.linalg.norm(self.geometry('v_y')), - np.linalg.norm(self.geometry('v_z'))): + if not np.isclose(np.linalg.norm(self.geometry['v_y']), + np.linalg.norm(self.geometry['v_z'])): raise GeometryError('Aperture does not have the same size in y, z direction.') def generate_local_xy(self, n): - phi = np.random.uniform(self.phi[0], self.phi[1], n) + phi = np.random.uniform(self.geometry.phi[0], self.geometry.phi[1], n) # normalize r_inner - r_inner = self.r_inner / np.linalg.norm(self.geometry('v_y')) + r_inner = self.geometry['r_inner'] / np.linalg.norm(self.geometry['v_y']) r = np.sqrt(np.random.uniform(r_inner**2, 1., n)) x = r * np.cos(phi) @@ -243,32 +143,8 @@ def generate_local_xy(self, n): @property def area(self): '''Area covered by the aperture''' - A_circ = np.pi * (np.linalg.norm(self.geometry('v_y'))**2 - self.r_inner**2) - return (self.phi[1] - self.phi[0]) / (2 * np.pi) * A_circ * u.mm**2 - - def outer_shape(self): - '''Return values in Eukledian space.''' - return self.xyz_circle(self.display.get('outer_factor', 3)) - - def inner_shape(self): - return self.xyz_circle(1, self.phi) - - def triangulate(self): - xyz, triangles = super(CircleAperture, self).triangulate() - if (self.r_inner > 0): - inneredge = self.xyz_circle(self.r_inner / - np.linalg.norm(self.geometry('v_y')), - self.phi) - # Inner edge of the display. If we have several stacked apertures, - # we don't want to fill is all up to r=0. - innerdisplay = self.xyz_circle(self.display.get('inner_factor', 0) * - self.r_inner / - np.linalg.norm(self.geometry('v_y')), - self.phi) - new_xyz, new_tri = plane_with_hole(inneredge, innerdisplay) - xyz, triangles = combine_disjoint_triangulations([xyz, new_xyz], - [triangles, new_tri]) - return xyz, triangles + A_circ = np.pi * (np.linalg.norm(self.geometry['v_y'])**2 - self.geometry['r_inner']**2) + return (self.geometry.phi[1] - self.geometry.phi[0]) / (2 * np.pi) * A_circ * u.mm**2 class MultiAperture(BaseAperture, BaseContainer): diff --git a/marxs/optics/base.py b/marxs/optics/base.py index ee2c2e6..d83a100 100644 --- a/marxs/optics/base.py +++ b/marxs/optics/base.py @@ -1,12 +1,13 @@ # Licensed under GPL version 3 - see LICENSE.rst from functools import wraps -from copy import copy +from copy import copy, deepcopy import numpy as np from astropy.table import Table, Row from ..math.pluecker import point_dir2plane, dir_point2line, intersect_line_plane from ..math.utils import e2h, h2e +from ..math.geometry import Geometry, FinitePlane from ..base import SimulationSequenceElement, _parse_position_keywords class OpticalElement(SimulationSequenceElement): @@ -25,86 +26,31 @@ class OpticalElement(SimulationSequenceElement): `process_photon`. Marxs will call `process_photons`, which (if not overwritten) contains a simple for-loop to loop over all photons in the array and call `process_photon` on each of them. - ''' - def geometry(self, key): - '''This function wraps access to the pos4d matrix. - - This is mostly a convenience method that gives access to vectors from the - ``pos4d`` matrix in familiar terms with string labels: - - - ``center``: The ``center`` is the origin of the local coordiante system - of the optical elemement. Typically, if will be the center of the - active plane, e.g. the center of the mirror surface. - - :math:`\hat v_{y,z}`: The box stretches in the y and z direction for - :math:`\pm \hat v_y` and :math:`\pm \hat v_z`. In optics, the relevant - interaction often happens on a surface, e.g. the surface of a mirror or - of a reflection grating. In the defaul configuration, this surface is in - the yz-plane. - - :math:`\hat v_x`: The thickness of the box is not as important in many - cases, - but is useful to e.g. render the geometry or to test if elements collide. - The box reaches from the yz-plane down to :math:`- \hat v_x`. Note, that - this definition means the size parallel to the y and z axis is - :math:`2 |\hat v_{y,z}|`, but only :math:`1 |\hat v_{x}|`. - - :math:`\hat e_{x,y,z}`: When an object is initialized, it automatically - adds unit vectors in the - direction of :math:`\hat v_{x,y,z}` called :math:`\hat e_{x,y,z}`. - - ``plane``: It also adds the homogeneous coordinates of the active plane - as ``plane``. - - Other labels get looked up in ``self._geometry`` and if the resulting - value is a 4-d vector, it gets transformed with ``pos4d``. - - Not all these labels make sense for every optical element (e.g. a curved - mirror does not really have a "plane"). - Access through this method is slower than direct indexing of ``self.pos4d``. - ''' + default_geometry = FinitePlane + '''If not geometry is passed in on initialization, an instance of this class will be used.''' - if key == 'center': - return self.pos4d[:, 3] - elif key in ['v_x', 'e_x']: - val = self.pos4d[:, 0] - elif key in ['v_y', 'e_y']: - val = self.pos4d[:, 1] - elif key in ['v_z', 'e_z']: - val = self.pos4d[:, 2] - elif key == 'plane': - return point_dir2plane(self.geometry('center'), self.geometry('e_x')) - else: - val = self._geometry[key] - if isinstance(val, np.ndarray) and (val.shape[-1] == 4): - val = np.dot(self.pos4d, val) - if key[:2] == 'e_': - return val / np.linalg.norm(val) - else: - return val + # SimulationSequenceElement has display none, but now we have a geometry + # so we don't need that default here. + display = {} def __init__(self, **kwargs): - self.pos4d = _parse_position_keywords(kwargs) - super(OpticalElement, self).__init__(**kwargs) - def intersect(self, dir, pos): - '''Calculate the intersection point between a ray and the element + geometry = kwargs.pop('geometry', self.default_geometry) + if isinstance(geometry, Geometry): + self.geometry = geometry + elif issubclass(geometry, Geometry): + self.geometry = geometry(kwargs) - Parameters - ---------- - dir : `numpy.ndarray` of shape (N, 4) - homogeneous coordinates of the direction of the ray - pos : `numpy.ndarray` of shape (N, 4) - homogeneous coordinates of a point on the ray + super(OpticalElement, self).__init__(**kwargs) - Returns - ------- - intersect : boolean array of length N - ``True`` if an intersection point is found. - interpos : `numpy.ndarray` of shape (N, 4) - homogeneous coordinates of the intersection point. Values are set - to ``np.nan`` if no intersection point is found. - interpos_local : `numpy.ndarray` of shape (N, 2) - y and z coordinates in the coordiante system of the active plane. - ''' - raise NotImplementedError + if not hasattr(self, "loc_coos_name"): + self.loc_coos_name = self.geometry.loc_coos_name + + @property + def pos4d(self): + return self.geometry.pos4d def process_photon(self, dir, pos, energy, polarization): '''Simulate interaction of optical element with a single photon. @@ -152,7 +98,7 @@ def process_photon(self, dir, pos, energy, polarization): return dir, pos, energy, polarization, probability, any, other, output, columns def __call__(self, photons): - intersect_out = self.intersect(photons['dir'].data, photons['pos'].data) + intersect_out = self.geometry.intersect(photons['dir'].data, photons['pos'].data) return self.process_photons(photons, *intersect_out) def process_photons(self, photons, intersect, interpos, intercoos): @@ -217,67 +163,7 @@ def process_photons(self, photons, intersect, interpos, intercoos): class FlatOpticalElement(OpticalElement): - '''Base class for geometrically flat optical elements. - - Compared with `OpticalElement` this adds methods to make the implementation of - flat elements easier. It adds a default `geometry`, a fast, vectorized method `intersect`, - and a template to `process_photons`. - - Derived classes have the option to implement their own `process_photons` or, alternatively, - they can implement a function called - ``specific_process_photons(self, photons, intersect, interpos, intercoos)`` that returns a dictionary - of the form ``{'column name': value, ...}`` where value is an array that holds one value for - each photon that intersects the optical element. In the special case of ``probability`` the - return value should only contain the probability assigned in **this** element. This value - will be multiplied with the previous probabilities of each photon automatically. - ''' - - display = {'shape': 'box'} - _geometry = {} - - loc_coos_name = ['y', 'z'] - '''name for output columns that contain the interaction point in local coordinates.''' - - def __init__(self, *args, **kwargs): - super(FlatOpticalElement, self).__init__(*args, **kwargs) - #copy class attribute to instance attribute - self._geometry = copy(self._geometry) - - def intersect(self, dir, pos): - '''Calculate the intersection point between a ray and the element - - Parameters - ---------- - dir : `numpy.ndarray` of shape (N, 4) - homogeneous coordinates of the direction of the ray - pos : `numpy.ndarray` of shape (N, 4) - homogeneous coordinates of a point on the ray - - Returns - ------- - intersect : boolean array of length N - ``True`` if an intersection point is found. - interpos : `numpy.ndarray` of shape (N, 4) - homogeneous coordinates of the intersection point. Values are set - to ``np.nan`` if no intersection point is found. - interpos_local : `numpy.ndarray` of shape (N, 2) - y and z coordinates in the coordiante system of the active plane - (not normalized to the dimensions of the element in question, but - in absolute units). - ''' - plucker = dir_point2line(h2e(dir), h2e(pos)) - interpos = intersect_line_plane(plucker, self.geometry('plane')) - vec_center_inter = - h2e(self.geometry('center')) + h2e(interpos) - ey = np.dot(vec_center_inter, h2e(self.geometry('e_y'))) - ez = np.dot(vec_center_inter, h2e(self.geometry('e_z'))) - intersect = ((np.abs(ey) <= np.linalg.norm(self.geometry('v_y'))) & - (np.abs(ez) <= np.linalg.norm(self.geometry('v_z')))) - if dir.ndim == 2: - interpos[~intersect, :3] = np.nan - # input of single photon. - elif (dir.ndim == 1) and not intersect: - interpos[:3] = np.nan - return intersect, interpos, np.vstack([ey, ez]).T + pass class FlatStack(FlatOpticalElement): diff --git a/marxs/optics/detector.py b/marxs/optics/detector.py index 5bfec13..148a2c1 100644 --- a/marxs/optics/detector.py +++ b/marxs/optics/detector.py @@ -3,12 +3,11 @@ import warnings import numpy as np -from transforms3d.affines import decompose44, compose +from transforms3d.affines import decompose44 from .base import FlatOpticalElement, OpticalElement -from ..math.utils import h2e from ..utils import SimulationSetupWarning - +from ..math.geometry import Cylinder class FlatDetector(FlatOpticalElement): '''Flat detector with square pixels @@ -31,6 +30,9 @@ class FlatDetector(FlatOpticalElement): pixsize : float size of pixels in mm + ignore_pixel_warning : bool + ignore warning if derived pixel number is not close to an integer + kwargs : see `args for optical elements` @@ -42,11 +44,13 @@ class FlatDetector(FlatOpticalElement): detpix_name = ['detpix_x', 'detpix_y'] '''name for output columns that contain this pixel number.''' + + display = {'color': (1.0, 1.0, 0.), 'shape': 'box', 'box-half': '+x'} - def __init__(self, pixsize=1, **kwargs): + def __init__(self, pixsize=1, ignore_pixel_warning=False, **kwargs): self.pixsize = pixsize super(FlatDetector, self).__init__(**kwargs) t, r, zoom, s = decompose44(self.pos4d) @@ -55,7 +59,7 @@ def __init__(self, pixsize=1, **kwargs): for i in (0, 1): z = zoom[i + 1] self.npix[i] = int(np.round(2. * z / self.pixsize)) - if np.abs(2. * z / self.pixsize - self.npix[i]) > 1e-2: + if (np.abs(2. * z / self.pixsize - self.npix[i]) > 1e-2) and not ignore_pixel_warning: warnings.warn('Detector size is not an integer multiple of pixel size in direction {0}. It will be rounded.'.format('xy'[i]), SimulationSetupWarning) self.centerpix[i] = (self.npix[i] - 1) / 2 @@ -94,169 +98,17 @@ class CircularDetector(OpticalElement): display = {'color': (1.0, 1.0, 0.), 'opacity': 0.7, - 'shape': 'surface', - 'coo1': np.linspace(0, 2 * np.pi, 50), - 'coo2': [-1, 1]} + } + + centerpix = [0, 0] + + default_geometry = Cylinder def __init__(self, pixsize=1, **kwargs): self.pixsize = pixsize - self.phi_offset = kwargs.pop('phi_offset', 0.) - self.inwards = kwargs.pop('inside', True) super(CircularDetector, self).__init__(**kwargs) - @classmethod - def from_rowland(cls, rowland, width): - '''Generate a `CircularDetector` from a `RowlandTorus`. - - Parameters - ---------- - rowland : `~marxs.design.RowlandTorus` - The circular detector is constructed to fit exactly into the - Rowland Circle defined by ``rowland``. - width : float - Half-width of the tube in the flat direction (z-axis) in mm - ''' - # Step 1: Get position and size from Rowland torus - pos4d_circ = compose([rowland.R, 0, 0], np.eye(3), [rowland.r, rowland.r, width]) - # Step 2: Transform to global coordinate system - pos4d_circ = np.dot(rowland.pos4d, pos4d_circ) - # Step 3: Make detector - return cls(pos4d=pos4d_circ, phi_offset=-np.pi) - - @property - def _inwardsoutwards(self): - 'Transform the self.inwards bool into [-1, +1]' - if self.inwards: - return 1. - else: - return -1. - - def intersect(self, dir, pos, transform=True): - '''Calculate the intersection point between a ray and the element - - Parameters - ---------- - dir : `numpy.ndarray` of shape (N, 4) - homogeneous coordinates of the direction of the ray - pos : `numpy.ndarray` of shape (N, 4) - homogeneous coordinates of a point on the ray - transform : bool - If ``True``, input is in global coordinates and needs to be transformed - here for the calculations; if ``False`` input is in local coordinates. - - Returns - ------- - intersect : boolean array of length N - ``True`` if an intersection point is found. - interpos : `numpy.ndarray` of shape (N, 4) - homogeneous coordinates of the intersection point. Values are set - to ``np.nan`` is no intersecton point is found. - interpos_local : `numpy.ndarray` of shape (N, 2) - phi, z coordiantes (in the local frame) for one of the intersection points. - If both intersection points are required, reset ``self.inner`` and call this - function again. - ''' - # This could be moved to a general function - if not np.all(dir[:, 3] == 0): - raise ValueError('First input must be direction vectors.') - # Could test pos, too... - if transform: - invpos4d = np.linalg.inv(self.pos4d) - dir = np.dot(invpos4d, dir.T).T - pos = np.dot(invpos4d, pos.T).T - - xyz = h2e(pos) - - # Solve quadratic equation in steps. a12 = (-xr +- sqrt(xr - r**2(x**2 - R**2))) - xy = xyz[:, :2] - r = dir[:, :2] - underroot = (np.einsum('ij,ij->i', xy, r))**2 - np.sum(r**2, axis=1) * (np.sum(xy**2, axis=1) - 1.) - intersect = (underroot >= 0) - i = intersect # just a shorthand because it's used so much below - - interpos_local = np.ones((pos.shape[0], 2)) - interpos_local[:] = np.nan - interpos = np.ones_like(pos) - interpos[:] = np.nan - - if intersect.sum() > 0: - b = np.sum(xy[i] * r[i], axis=1) - denom = np.sum(r[i]**2, axis=1) - a1 = (- b + np.sqrt(underroot[i])) / denom - a2 = (- b - np.sqrt(underroot[i])) / denom - x1 = xy[i, :] + a1[:, np.newaxis] * r[i, :] - apick = np.where(self._inwardsoutwards * np.sum(x1 * r[i, :], axis=1) >=0, a1, a2) - xy_p = xy[i, :] + apick[:, np.newaxis] * r[i, :] - phi = np.arctan2(xy_p[:, 1], xy_p[:, 0]) - # Shift phi by offset, then wrap to that it is in range [-pi, pi] - interpos_local[i, 0] = (phi - self.phi_offset + np.pi) % (2 * np.pi) - np.pi - # Those look like they hit in the xy plane. - # Still possible to miss if z axis is too large. - # Calculate z-coordiante at intersection - interpos_local[intersect, 1] = xyz[i, 2] + apick * dir[i, 2] - interpos[i, :2] = xy_p - interpos[i, 2] = interpos_local[i, 1] - interpos[i, 3] = 1 - # set those elements on intersect that miss in z to False - trans, rot, zoom, shear = decompose44(self.pos4d) - z_p = interpos[i, 2] - intersect[i.nonzero()[0][np.abs(z_p) > 1]] = False - # Now reset everything to nan that does not intersect - interpos_local[~i, :] = np.nan - # interpos_local in z direction is in local coordinates, i.e. - # the x coordiante is 0..1, but we want that in units of the - # global coordinate system. - interpos_local[:, 1] = interpos_local[:, 1] * zoom[2] - interpos[~i, :] = np.nan - - interpos = np.dot(self.pos4d, interpos.T).T - - return intersect, interpos, interpos_local - - def process_photons(self, photons, intersect, interpos, inter_local): - photons['pos'][intersect, :] = interpos[intersect, :] - self.add_output_cols(photons, self.loc_coos_name + self.detpix_name) - # Add ID number to ID col, if requested - if self.id_col is not None: - photons[self.id_col][intersect] = self.id_num - # Set position in different coordinate systems - photons['pos'][intersect] = interpos[intersect] - photons[self.loc_coos_name[0]][intersect] = inter_local[intersect, 0] - photons[self.loc_coos_name[1]][intersect] = inter_local[intersect, 1] - trans, rot, zoom, shear = decompose44(self.pos4d) - if np.isclose(zoom[0], zoom[1]): - photons[self.detpix_name[0]][intersect] = inter_local[intersect, 0] * zoom[0] / self.pixsize - else: - warnings.warn('Pixel coordinate for elliptical mirrors not implemented.', PixelSizeWarning) - photons[self.detpix_name[1]][intersect] = inter_local[intersect, 1] / self.pixsize - - return photons - - def parametric_surface(self, phi, z=np.array([-1, 1])): - '''Parametric description of the tube. - - This is just another way to obtain the shape of the tube, e.g. - for visualization. - - Parameters - ---------- - phi : np.array - ``phi`` is the angle around the tube profile. - z : np.array - The coordiantes along the radius coordinate. - - Returns - ------- - xyzw : np.array - Ring coordinates in global homogeneous coordinate system. - ''' - phi = np.asanyarray(phi) - z = np.asanyarray(z) - if (phi.ndim != 1) or (z.ndim != 1): - raise ValueError('input parameters have 1-dim shape.') - phi, z = np.meshgrid(phi, z) - x = np.cos(phi) - y = np.sin(phi) - w = np.ones_like(z) - coos = np.array([x, y, z, w]).T - return np.einsum('...ij,...j', self.pos4d, coos) + def specific_process_photons(self, photons, intersect, interpos, intercoos): + detx = intercoos[intersect, 0] * self.geometry['R'] / self.pixsize + self.centerpix[0] + dety = intercoos[intersect, 1] / self.pixsize + self.centerpix[1] + return {self.detpix_name[0]: detx, self.detpix_name[1]: dety} diff --git a/marxs/optics/filter.py b/marxs/optics/filter.py index c9b152f..79f92a8 100644 --- a/marxs/optics/filter.py +++ b/marxs/optics/filter.py @@ -4,9 +4,10 @@ import numpy as np from .base import OpticalElement, FlatOpticalElement +from ..base import MarxsElement +from ..math.geometry import NoGeometry - -class GlobalEnergyFilter(OpticalElement): +class GlobalEnergyFilter(MarxsElement): '''Energy dependent filter that globally affects all photons. This element is used on all photons in the list, there is no geometrical @@ -37,6 +38,7 @@ class GlobalEnergyFilter(OpticalElement): -------- marxs.optics.filter.EnergyFilter ''' + def __init__(self, **kwargs): self.filterfunc = kwargs.pop('filterfunc') super(GlobalEnergyFilter, self).__init__(**kwargs) diff --git a/marxs/optics/grating.py b/marxs/optics/grating.py index f8033dc..c092127 100644 --- a/marxs/optics/grating.py +++ b/marxs/optics/grating.py @@ -2,6 +2,7 @@ '''Gratings and efficiency files''' import math import numpy as np +from numpy.core.umath_tests import inner1d from ..math.utils import norm_vector, h2e, e2h from ..math.polarization import parallel_transport @@ -131,7 +132,7 @@ class FlatGrating(FlatOpticalElement): loc_coos_name = ['grat_y', 'grat_z'] '''name for output columns that contain the interaction point in local coordinates.''' - def order_sign_convention(self, p): + def order_sign_convention(self, p, e_perp_groove): '''Set sign convention for grating orders. This sets the following, somewhat arbitrary convention: @@ -152,7 +153,8 @@ def order_sign_convention(self, p): ---------- p : np.array Array of Eucleadian direction vectors - + e_perp_groove : np.array + Array of local groove directions at the positions where the photons hit ''' # -1 because n, l, d should be right-handed coordinate system # while n = e_x, l = e_x, and d = e_y would be left-handed. @@ -164,12 +166,34 @@ def __init__(self, **kwargs): if 'd' not in kwargs: raise ValueError('Input parameter "d" (Grating constant) is required.') self._d = kwargs.pop('d') - groove = kwargs.pop('groove_angle', 0.) + groove_angle = kwargs.pop('groove_angle', 0.) super(FlatGrating, self).__init__(**kwargs) - self._geometry['e_groove'] = np.array([0., math.sin(-groove), math.cos(-groove), 0.]) - self._geometry['e_perp_groove'] = np.array([0., math.cos(-groove), -math.sin(-groove), 0.]) + self.geometry._geometry['groove_angle'] = groove_angle + + def e_groove_coos(self, intercoos): + '''Get coordiante orthonormal coordinate system along groove direction. + + Parameters + ---------- + intercoos : `numpy.ndarray` of shape (N, 2) + coordinates in the coordiante system of the geometry (e.g. (x, y), or + (r, phi)). + + Returns + ------- + e_groove, e_perp_groove, n : `numpy.ndarray` of shape (N, 4) + Vectors pointing along the groove direction (parallel to the surface), + perpendicular to the groove direction (parallel to surface), + and normal to surface. + ''' + + groove = self.geometry['groove_angle'] + ex, ey, en = self.geometry.get_local_euklid_bases(intercoos) + e_groove = math.sin(-groove) * ex + math.cos(-groove) * ey + e_perp_groove = math.cos(-groove) * ex - math.sin(-groove) * ey + return e_groove, e_perp_groove, en def d(self, intercoos): '''Method that returns the grating constant at given positions. @@ -184,37 +208,47 @@ def d(self, intercoos): else: return self._d(intercoos) + def blaze_angle_modifier(self, intercoos): + '''Modify blaze angle + + In `diffract_photons` the blaze angle is calculated relative to the surface + of the grating. In cases where that number has to be modified (e.g. when the grating + bars are not perpendicualr to the surface) the blaze angle can be modified by + overriding this function. + ''' + return np.zeros(intercoos.shape[0]) + def diffract_photons(self, photons, intersect, interpos, intercoos): '''Vectorized implementation''' - p = norm_vector(h2e(photons['dir'].data[intersect])) - n = self.geometry('plane')[:3] - l = h2e(self.geometry('e_groove')) + p = norm_vector(photons['dir'].data[intersect]) + l, d, n = self.e_groove_coos(intercoos[intersect]) # Minus sign here because we want n, l, d to be a right-handed coordinate system - d = -h2e(self.geometry('e_perp_groove')) + d = - d wave = energy2wave / photons['energy'].data[intersect] # calculate angle between normal and (ray projected in plane perpendicular to groove) # -> this is the blaze angle - p_perp_to_grooves = norm_vector(p - np.dot(p, l)[:, np.newaxis] * l) + p_perp_to_grooves = norm_vector(p - inner1d(p, l)[:, None] * l) # Use abs here so that blaze angle is always in 0..pi/2 # independent of the relative orientation of p and n. - blazeangle = np.arccos(np.abs(np.dot(p_perp_to_grooves, n))) + blazeangle = np.arccos(np.abs(inner1d(p_perp_to_grooves, n))) + blazeangle += self.blaze_angle_modifier(intercoos[intersect, :]) m, prob = self.order_selector(photons['energy'].data[intersect], photons['polarization'].data[intersect], blazeangle) # The idea to calculate the components in the (d,l,n) system separately # is taken from MARX - sign = self.order_sign_convention(p) - p_d = np.dot(p, d) + sign * m * wave / self.d(intercoos[intersect, :]) - p_l = np.dot(p, l) + sign = self.order_sign_convention(p, d) + p_d = inner1d(p, d) + sign * m * wave / self.d(intercoos[intersect, :]) + p_l = inner1d(p, l) # The norm for p_n can be derived, but the direction needs to be chosen. p_n = np.sqrt(1. - p_d**2 - p_l**2) # Check if the photons have same direction compared to normal before - direction = np.sign(np.dot(p, n), dtype=np.float) + direction = np.sign(inner1d(p, n), dtype=np.float) if not self.transmission: direction *= -1 - dir = e2h(p_d[:, None] * d[None, :] + p_l[:, None] * l[None, :] + (direction * p_n)[:, None] * n[None, :], 0) + dir = p_d[:, None] * d + p_l[:, None] * l + (direction * p_n)[:, None] * n return dir, m, prob, blazeangle def specific_process_photons(self, photons, intersect, interpos, intercoos): @@ -236,15 +270,13 @@ class CATGrating(FlatGrating): ''' - def order_sign_convention(self, p): + def order_sign_convention(self, p, e_perp_groove): '''Convention to chose the sign for CAT grating orders Blazing happens on the side of the negative orders. Obviously, this convention is only meaningful if the photons do not arrive perpendicular to the grating. ''' - # Minus sign here because we want n, l, d to be a right-handed coordinate system - d = -h2e(self.geometry('e_perp_groove')) - dotproduct = np.dot(p, d) + dotproduct = inner1d(p, e_perp_groove) sign = np.sign(dotproduct) sign[sign == 0] = 1 return sign diff --git a/marxs/optics/mirror.py b/marxs/optics/mirror.py index 2b32794..7f740fd 100644 --- a/marxs/optics/mirror.py +++ b/marxs/optics/mirror.py @@ -33,7 +33,7 @@ def __init__(self, **kwargs): def specific_process_photons(self, photons, intersect, interpos, intercoos): # A ray through the center is not broken. # So, find out where a central ray would go. - focuspoints = h2e(self.geometry('center')) + self.focallength * norm_vector(h2e(photons['dir'][intersect])) + focuspoints = h2e(self.geometry['center']) + self.focallength * norm_vector(h2e(photons['dir'][intersect])) dir = e2h(focuspoints - h2e(interpos[intersect]), 0) pol = parallel_transport(photons['dir'].data[intersect, :], dir, photons['polarization'].data[intersect, :]) @@ -85,17 +85,17 @@ def __init__(self, **kwargs): super(ThinLens, self).__init__(**kwargs) def process_photon(self, dir, pos, energy, polarization): - intersect, h_intersect, loc_inter = self.intersect(dir, pos) + intersect, h_intersect, loc_inter = self.geometry.intersect(dir, pos) distance = distance_point_point(h_intersect, - self.geometry('center')[np.newaxis, :]) + self.geometry['center'][np.newaxis, :]) if distance == 0.: # No change of direction for rays through the origin. # Need to special case this, because rotation axis is not defined # in this case. new_ray_dir = h2e(dir) else: - delta_angle = distance / self.focallength - e_rotation_axis = np.cross(dir[:3], (h_intersect - self.geometry('center'))[:3]) + delta_angle = -distance / self.focallength + e_rotation_axis = np.cross(dir[:3], (h_intersect - self.geometry['center'])[:3]) # This is the first step that cannot be done on a stack of photons # Could have a for "i in photons", but might come up with better way rot = axangle2mat(e_rotation_axis, delta_angle) diff --git a/marxs/optics/multiLayerMirror.py b/marxs/optics/multiLayerMirror.py index 0b3e98f..89fd7a7 100644 --- a/marxs/optics/multiLayerMirror.py +++ b/marxs/optics/multiLayerMirror.py @@ -56,7 +56,7 @@ def specific_process_photons(self, photons, intersect, intersection, local): # - p is p polarization (in the plane of incidence) # First, make basis vectors. - v_s = np.cross(beam_dir, self.geometry('e_x')[0:3]) + v_s = np.cross(beam_dir, self.geometry['e_x'][0:3]) v_s /= np.linalg.norm(v_s, axis=1)[:, np.newaxis] v_p = np.cross(beam_dir, v_s) @@ -125,8 +125,8 @@ def interp_files(self, photons, local): tested_polarized_fraction = np.interp(photons['energy'], polarizedFile['Photon energy'] / 1000, polarizedFile['Polarization']) # find probability of being reflected due to position - local_x = local[:, 0] / np.linalg.norm(self.geometry('v_y')) - local_coords_in_file = reflectFile['X(mm)'] / np.linalg.norm(self.geometry('v_y')) - 1 + local_x = local[:, 0] / np.linalg.norm(self.geometry['v_y']) + local_coords_in_file = reflectFile['X(mm)'] / np.linalg.norm(self.geometry['v_y']) - 1 # interpolate 'Peak lambda', 'Peak' [reflectivity], and 'FWHM(nm)' to the actual photon positions peak_wavelength = np.interp(local_x, local_coords_in_file, reflectFile['Peak lambda']) max_refl = np.interp(local_x, local_coords_in_file, reflectFile['Peak']) / tested_polarized_fraction diff --git a/marxs/optics/tests/test_detector.py b/marxs/optics/tests/test_detector.py index 2d99fa7..1611827 100644 --- a/marxs/optics/tests/test_detector.py +++ b/marxs/optics/tests/test_detector.py @@ -7,6 +7,7 @@ from ...tests import closeornan from ...math.utils import h2e from ...design import RowlandTorus +from ...math.geometry import Cylinder def test_pixelnumbers(): pos = np.array([[0, 0., -0.25, 1.], @@ -35,123 +36,14 @@ def test_nonintegerwarning(recwarn): w = recwarn.pop() assert 'is not an integer multiple' in str(w.message) -def test_intersect_tube_miss(): - '''Rays passing at larger radius or through the center miss the tube.''' - circ = CircularDetector() - # Passing at larger radius - intersect, interpos, inter_local = circ.intersect(np.array([[1., 0., .1, 0], - [-1., -1., 0., 0.]]), - np.array([[0, -1.1, 0., 1.], - [2., 0., 0., 1.]])) - assert np.all(intersect == False) - assert np.all(np.isnan(interpos)) - assert np.all(np.isnan(inter_local)) - - # Through the center almost parallel to z-axis - intersect, interpos, inter_local = circ.intersect(np.array([[0., 0.1, 1., 0]]), - np.array([[0.5, 0.5, 0., 1.]])) - assert np.all(intersect == False) - assert np.all(np.isnan(interpos)) - assert np.all(np.isnan(inter_local)) - - # Repeat with a tube that's moved to make sure we did not mess up local and - # global coordinates in the implementation - circ = CircularDetector(position=[0.8, 0.8, 1.2]) - intersect, interpos, inter_local = circ.intersect(np.array([[1., 0., .0, 0], - [1., 0., 0., 0.]]), - np.array([[0., 0., 0., 1.], - [0., -0.3, 1., 1.]])) - assert np.all(intersect == False) - assert np.all(np.isnan(interpos)) - assert np.all(np.isnan(inter_local)) - -def test_intersect_tube_hitmiss_zoomed(): - '''Rays that hit a one tube, but miss a shorter one''' - dir = np.array([[-.1, 0., 0., 0], [1., 0, -0.2, 0.]]) - pos = np.array([[20., 0., 0.8, 1.], [2., .0, 0.4, 1.]]) - - circ = CircularDetector(zoom=[1, 1, 1]) - intersect, interpos, inter_local = circ.intersect(dir, pos) - assert np.all(intersect == True) - - circ = CircularDetector(zoom=[1, 1, .5]) - intersect, interpos, inter_local = circ.intersect(dir, pos) - assert np.all(intersect == False) - assert np.all(np.isnan(interpos)) - assert np.all(np.isnan(inter_local)) - -def test_intersect_tube_2points(): - dir = np.array([[-.1, 0., 0., 0], [-0.5, 0, 0., 0.], [-1., 0., 0., 0.]]) - pos = np.array([[50., 0., 0., 1.], [10., .5, 0., 1.], [2, 0, .3, 1.]]) - - circ = CircularDetector(phi_offset=-np.pi) - intersect, interpos, inter_local = circ.intersect(dir, pos) - assert np.all(intersect == True) - assert np.allclose(h2e(interpos), np.array([[-1., 0., 0.], - [-np.sqrt(0.75), 0.5, 0.], - [-1, 0., .3]])) - assert np.allclose(inter_local, np.array([[0., -np.arcsin(0.5), 0.], - [0., 0., 0.3]]).T) - -def test_intersect_tube_2points_zoomed(): - '''inter_local should be the same if we make the tube longer in z''' - dir = np.array([[-.1, 0., 0., 0], [-0.5, 0, 0., 0.], [-1., 0., 0., 0.]]) - pos = np.array([[50., 0., 0., 1.], [10., .5, 0., 1.], [2, 0, .3, 1.]]) - - circ = CircularDetector(phi_offset=-np.pi, zoom=[1, 1., 5.]) - intersect, interpos, inter_local = circ.intersect(dir, pos) - assert np.all(intersect == True) - assert np.allclose(h2e(interpos), np.array([[-1., 0., 0.], - [-np.sqrt(0.75), 0.5, 0.], - [-1, 0., .3]])) - assert np.allclose(inter_local, np.array([[0., -np.arcsin(0.5), 0.], - [0., 0., 0.3]]).T) - -def test_intersect_tube_2points_outside(): - '''Now hit the outside. Looks very similar, if we remove the phi_offset, too.''' - dir = np.array([[-.1, 0., 0., 0], [-0.5, 0, 0., 0.], [-1., 0., 0., 0.]]) - pos = np.array([[50., 0., 0., 1.], [10., .5, 0., 1.], [2, 0, .3, 1.]]) - - circ = CircularDetector(inside=False) - intersect, interpos, inter_local = circ.intersect(dir, pos) - assert np.all(intersect == True) - assert np.allclose(h2e(interpos), np.array([[1., 0., 0.], - [np.sqrt(0.75), 0.5, 0.], - [1, 0., .3]])) - assert np.allclose(inter_local, np.array([[0., np.arcsin(0.5), 0.], - [0., 0., 0.3]]).T) - -def test_intersect_tube_2points_translation(): - '''Repeat with a tube that's moved and zoomed to make sure we - did not mess up local and global coordinates in the implementation - ''' - circ = CircularDetector(position=[0.8, 0.8, 1.2], zoom=[1, 1, 2]) - intersect, interpos, inter_local = circ.intersect(np.array([[1., 0., .0, 0], - [1., 0., 0., 0.]]), - np.array([[0., 0., 0., 1.], - [0., 0., 1., 1.]])) - assert np.all(intersect == True) - assert np.allclose(h2e(interpos), np.array([[1.4, 0., 0.], - [1.4, 0., 1.]])) - -def test_tube_parametric(): - '''Generate points on surface using parametic. Then make rays for intersection - which should return those points as intersection points.''' - circ = CircularDetector(zoom=[2., 3., 4.]) - parametric = circ.parametric_surface(phi=[-.1, 0., 1.]) - parametric = parametric.reshape((6, 4)) - # select dir so that it's in the right direction to recover these points. - dir = np.tile([1, 0, 0, 0], (6, 1)) - intersect, interpos, inter_local = circ.intersect(dir, parametric) - assert np.allclose(parametric, interpos) def test_CircularDetector_from_Rowland(): '''If a circular detector is very narrow, then all points should be very close to the Rowland torus.''' rowland = RowlandTorus(R=6e4, r=5e4, position=[123., 345., -678.], orientation=transforms3d.euler.euler2mat(1, 2, 3, 'syxz')) - detcirc = CircularDetector.from_rowland(rowland, width=1e-6) + detcirc = CircularDetector(geometry=Cylinder.from_rowland(rowland, width=1e-6)) phi = np.mgrid[0:2.*np.pi:360j] - points = detcirc.parametric_surface(phi) + points = detcirc.geometry.parametric_surface(phi) # Quartic < 1e5 is very close for these large values of r and R. assert np.max(np.abs(rowland.quartic(h2e(points)))) < 1e5 diff --git a/marxs/optics/tests/test_grating.py b/marxs/optics/tests/test_grating.py index c1189d8..2240341 100644 --- a/marxs/optics/tests/test_grating.py +++ b/marxs/optics/tests/test_grating.py @@ -15,8 +15,9 @@ def test_zeros_order(): '''Photons diffracted into order 0 should just pass through''' - dir = random((10, 4)) * 10 - 5 - photons = Table({'pos': random((10, 4)) * 10 - 5, + dir = random((10, 4)) + dir[0, :] = dir[0, :] * 10 # make sure mostion is mostly along x so they hit grating + photons = Table({'pos': random((10, 4)) - 1, 'dir': dir, 'energy': random(10), 'polarization': polarization_vectors(dir, np.random.rand(10)), @@ -30,7 +31,7 @@ def test_zeros_order(): p = photons.copy() g0 = FlatGrating(d=1./500., order_selector=OrderSelector([0]), - zoom=np.array([1., 5., 5.])) + zoom=np.array([1., 50., 50.])) p = g0(p) # Direction unchanged d_in = h2e(photons['dir']) @@ -47,7 +48,7 @@ def test_zeros_order(): def test_translation_invariance(): '''For homogeneous gratings, the diffraction eqn does not care where the ray hits.''' photons = generate_test_photons(2) - photons['dir'] = np.tile(np.array([1., 2., 3., 0.]), (2, 1)) + photons['dir'] = np.tile(np.array([-1., 2., 3., 0.]), (2, 1)) photons['polarization'] = polarization_vectors(photons['dir'], [2., 4.]) pos = np.tile(np.array([1., 0., 0., 1.]), (2, 1)) delta_pos = np.array([0, .23, -.34, 0.]) @@ -93,7 +94,7 @@ def test_angle_dependence(): def test_order_dependence(): '''For small theta, the change in direction is m * dtheta''' photons = Table({'pos': np.ones((5,4)), - 'dir': np.tile([1., 0, 0, 0], (5,1)), + 'dir': np.tile([-1., 0, 0, 0], (5,1)), 'energy': np.ones(5), 'polarization': np.tile([0.,1.,0.,0.], (5,1)), 'probability': np.ones(5), @@ -112,7 +113,7 @@ def mock_order(x, y, z): def test_energy_dependence(): '''The grating angle should depend on the photon wavelength <-> energy.''' photons = Table({'pos': np.ones((5,4)), - 'dir': np.tile([1., 0, 0, 0], (5,1)), + 'dir': np.tile([-1., 0, 0, 0], (5,1)), 'energy': np.arange(1., 6), 'polarization': np.tile([0, 1., 0., 0.], (5, 1)), 'probability': np.ones(5), @@ -145,7 +146,8 @@ def test_groove_direction(): order1 = OrderSelector([1]) g = FlatGrating(d=1./500, order_selector=order1) - assert np.allclose(np.dot(g.geometry('e_groove'), g.geometry('e_perp_groove')), 0.) + e_groove, e_perp_groove, n = g.e_groove_coos(np.ones((2,1))) + assert np.isclose(np.dot(e_groove[0, :], e_perp_groove[0, :]), 0.) p = g(photons.copy()) g1 = FlatGrating(d=1./500, order_selector=order1, groove_angle=.3) @@ -212,7 +214,9 @@ def test_CAT_order_convention(): gm = CATGrating(d=1./5000, order_selector=OrderSelector([-5]), zoom=2) m5 = gm(photons.copy()) for g in [gm, gp]: - assert np.all(g.order_sign_convention(h2e(photons['dir'])) == np.array([1, -1, -1, 1, 1])) + e_groove, e_perp, n = g.e_groove_coos(np.zeros((5, 2))) + signs = g.order_sign_convention(photons['dir'], e_perp) + assert np.all(signs == np.array([1, 1, 1, -1, -1])) assert np.all(p5['dir'][1:3, 1] > 0) assert np.all(p5['dir'][3:, 1] < 0) assert np.all(m5['dir'][1:3, 1] < 0) @@ -320,7 +324,7 @@ def test_change_position_after_init(): # Make grating at some point g2 = FlatGrating(d=1./500, order_selector=OrderSelector([1]), position=trans) # then move it so that pos and groove direction match g1 - g2.pos4d = affines.compose(np.zeros(3), pos3d, 25.*np.ones(3)) + g2.geometry.pos4d = affines.compose(np.zeros(3), pos3d, 25.*np.ones(3)) p2 = g2(photons.copy()) assert np.allclose(p1['dir'], p2['dir']) @@ -335,4 +339,6 @@ def test_gratings_are_independent(): ''' g1 = FlatGrating(d=1./500, order_selector=OrderSelector([1]), groove_angle=.3) g2 = FlatGrating(d=1./500, order_selector=OrderSelector([1])) - assert not np.allclose(g1.geometry('e_groove'), g2.geometry('e_groove')) + groove1, perp1, n1 = g1.e_groove_coos(np.array([[-1, 2.3]])) + groove2, perp2, n2 = g2.e_groove_coos(np.array([[-1, 2.3]])) + assert not np.allclose(groove1, groove2) diff --git a/marxs/optics/tests/test_optics.py b/marxs/optics/tests/test_optics.py index 281b0c5..1b8fd30 100644 --- a/marxs/optics/tests/test_optics.py +++ b/marxs/optics/tests/test_optics.py @@ -34,18 +34,18 @@ def test_pos4d_no_transformation(): ''' oe = marxs.optics.aperture.RectangleAperture() assert np.all(oe.pos4d == np.eye(4)) - assert np.all(oe.geometry('center') == np.array([0, 0, 0, 1])) - assert np.all(oe.geometry('e_y') == np.array([0, 1, 0, 0])) - assert np.all(oe.geometry('e_z') == np.array([0, 0, 1, 0])) + assert np.all(oe.geometry['center'] == np.array([0, 0, 0, 1])) + assert np.all(oe.geometry['e_y'] == np.array([0, 1, 0, 0])) + assert np.all(oe.geometry['e_z'] == np.array([0, 0, 1, 0])) def test_pos4d_rotation(): rotation = axangle2aff(np.array([1, 1, 0]), np.deg2rad(90)) oe = marxs.optics.aperture.RectangleAperture(orientation=rotation[:3,:3]) assert np.all(oe.pos4d == rotation) - assert np.all(oe.geometry('center') == np.array([0, 0, 0, 1])) - assert np.allclose(oe.geometry('e_y'), np.array([.5, .5, 1. / np.sqrt(2), 0])) - assert np.allclose(oe.geometry('e_z'), np.array([1., -1, 0, 0]) / np.sqrt(2)) + assert np.all(oe.geometry['center'] == np.array([0, 0, 0, 1])) + assert np.allclose(oe.geometry['e_y'], np.array([.5, .5, 1. / np.sqrt(2), 0])) + assert np.allclose(oe.geometry['e_z'], np.array([1., -1, 0, 0]) / np.sqrt(2)) def test_pos4d_translation(): @@ -53,9 +53,9 @@ def test_pos4d_translation(): expectedpos4d = np.eye(4) expectedpos4d[:3, 3] = np.array([-2, 3, 4.3]) assert np.all(oe.pos4d == expectedpos4d) - assert np.all(oe.geometry('center') == np.array([-2, 3, 4.3, 1])) - assert np.all(oe.geometry('e_y') == np.array([0, 1, 0, 0])) - assert np.all(oe.geometry('e_z') == np.array([0, 0, 1, 0])) + assert np.all(oe.geometry['center'] == np.array([-2, 3, 4.3, 1])) + assert np.all(oe.geometry['e_y'] == np.array([0, 1, 0, 0])) + assert np.all(oe.geometry['e_z'] == np.array([0, 0, 1, 0])) @@ -123,9 +123,9 @@ def test_FlatStack(): elements=[marxs.optics.EnergyFilter, marxs.optics.FlatDetector], keywords=[{'filterfunc': lambda x: 0.5},{}]) # Check all layers are at the same position - assert np.allclose(fs.geometry('center'), np.array([0, 5, 0, 1])) - assert np.allclose(fs.elements[0].geometry('center'), np.array([0, 5, 0, 1])) - assert np.allclose(fs.elements[1].geometry('center'), np.array([0, 5, 0, 1])) + assert np.allclose(fs.geometry['center'], np.array([0, 5, 0, 1])) + assert np.allclose(fs.elements[0].geometry['center'], np.array([0, 5, 0, 1])) + assert np.allclose(fs.elements[1].geometry['center'], np.array([0, 5, 0, 1])) assert np.allclose(fs.pos4d, fs.elements[1].pos4d) p = generate_test_photons(5) diff --git a/marxs/simulator/__init__.py b/marxs/simulator/__init__.py index b3feb36..bea57e6 100644 --- a/marxs/simulator/__init__.py +++ b/marxs/simulator/__init__.py @@ -3,4 +3,6 @@ BaseContainer, Sequence, Parallel, ParallelCalculated, KeepCol, + propagate, + Propagator ) diff --git a/marxs/simulator/simulator.py b/marxs/simulator/simulator.py index 3bd9b50..b48c6c3 100644 --- a/marxs/simulator/simulator.py +++ b/marxs/simulator/simulator.py @@ -3,6 +3,7 @@ from transforms3d.affines import decompose44 from ..math.utils import translation2aff, zoom2aff, mat2aff, e2h, h2e +from ..math.geometry import Geometry from ..base import SimulationSequenceElement, _parse_position_keywords import transforms3d from transforms3d.utils import normalized_vector @@ -218,8 +219,18 @@ class Parallel(BaseContainer): of detectors on a detector wheel. The uncertainty is expressed as an affine transformation matrix. ''' + + @property + def pos4d(self): + return self.geometry.pos4d + def __init__(self, **kwargs): - self.pos4d = _parse_position_keywords(kwargs) + geometry = kwargs.pop('geometry', Geometry) + if isinstance(geometry, Geometry): + self.geometry = geometry + elif issubclass(geometry, Geometry): + self.geometry = geometry(kwargs) + self.id_num_offset = kwargs.pop('id_num_offset', 0) self.elem_class = kwargs.pop('elem_class') # Need to operate on a copy here, to avoid changing elem_args of outer level @@ -271,7 +282,7 @@ def move_center(self, change4d): change4d : np.array of shape (4,4) Homogeneous coordinate transformation matrix ''' - self.pos4d = np.dot(change4d, self.pos4d) + self.geometry.pos4d = np.dot(change4d, self.pos4d) inv = np.linalg.inv(change4d) self.elem_pos = [np.dot(inv, p) for p in self.elem_pos] self.generate_elements() @@ -517,3 +528,53 @@ def format_positions(self, atol=None): equal_nan=True): ind.append(i) return d[:, ind, :] + + +class Propagator(SimulationSequenceElement): + def __init__(self, **kwargs): + '''Move photons along the rays. + + This element moves all photons in the photon list by a specified amount + forwards or backwards along the ray. + Normally, photons only move forward in the simulation, but sometimes it + can be useful to reset the photons to an earlier position, e.g. to see + how the pattern on a detector changes for different detector positions. + + Parameters + ---------- + distance : float + Distance for the photons to move. Negative values move photons + backwards. + ''' + self.distance = kwargs.pop('distance') + super(Propagator, self).__init__(**kwargs) + + def __call__(self, photons): + photons['pos'] = photons['pos'] + self.distance * photons['dir'] + return photons + + +def propagate(photons, d): + '''Move photons along the rays. + + This function moves all photons in the photon list by a specified amount + forwards or backwards along the ray. + Normally, photons only move forward in the simulation, but sometimes it + can be useful to reset the photons to an earlier position, e.g. to see + how the pattern on a detector changes for different detector positions. + + Parameters + ---------- + photons : `astropy.table.Table` + Photon table with pos and dir entries + d : float + Distance for the photons to move. Negative values move photons + backwards. + + Returns + ------- + photons : `astropy.table.Table` + Photon table with pos and dir entries + ''' + photons['pos'] = photons['pos'] + d * photons['dir'] + return photons diff --git a/marxs/simulator/tests/test_parallel.py b/marxs/simulator/tests/test_parallel.py index d879f46..67a4541 100644 --- a/marxs/simulator/tests/test_parallel.py +++ b/marxs/simulator/tests/test_parallel.py @@ -79,8 +79,8 @@ def test_parallel_calculated_normals(): e1 = det.elements[0] e2 = det.elements[1] # These two edges should be close together: - edge1 = e1.geometry('center') + e1.geometry('v_y') - edge2 = e2.geometry('center') - e2.geometry('v_y') + edge1 = e1.geometry['center'] + e1.geometry['v_y'] + edge2 = e2.geometry['center'] - e2.geometry['v_y'] assert np.all(np.abs(edge1 - edge2) < 3.) def test_parallel_calculated_rotations(): diff --git a/marxs/source/source.py b/marxs/source/source.py index 4e8e1f9..687743f 100644 --- a/marxs/source/source.py +++ b/marxs/source/source.py @@ -401,7 +401,7 @@ def __init__(self, **kwargs): np.random.rand(n) * (self.func_par[0]**2 - self.func_par[1]**2)) super(DiskSource, self).__init__(**kwargs) -class GaussSource(RadialDistributionSource): +class GaussSource(AstroSource): '''Astrophysical source with a Gaussian brightness profile. This source uses a small angle approximation which is valid for radii less than @@ -413,11 +413,27 @@ class GaussSource(RadialDistributionSource): Gaussian sigma setting the width of the distribution ''' def __init__(self, **kwargs): - kwargs['func_par'] = kwargs.pop('sigma') - # Note: rand is in interavall [0..1[, so 1-rand is the same except for edges - kwargs['radial_distribution'] = lambda n: self.func_par * np.sqrt(np.log(np.random.rand(n))) + self.sigma = kwargs.pop('sigma') super(GaussSource, self).__init__(**kwargs) + def generate_photons(self, exposuretime): + '''Photon positions are generated in a frame that is centered on the + coordinates set in ``coords``, then they get transformed into the global sky + system. + ''' + photons = super(GaussSource, self).generate_photons(exposuretime) + + relative_frame = SkyOffsetFrame(origin=self.coords) + n = len(photons) + relative_coords = SkyCoord(np.random.normal(scale=self.sigma.value, size=n) * self.sigma.unit, + np.random.normal(scale=self.sigma.value, size=n) * self.sigma.unit, + frame=relative_frame) + origin_coord = relative_coords.transform_to(self.coords) + self.set_pos(photons, origin_coord) + + return photons + + class SymbolFSource(AstroSource): '''Source shaped like the letter F. diff --git a/marxs/source/tests/test_source.py b/marxs/source/tests/test_source.py index e7833ee..1cd0651 100644 --- a/marxs/source/tests/test_source.py +++ b/marxs/source/tests/test_source.py @@ -200,10 +200,10 @@ def test_disk_distribution(diskclass, diskpar): s = diskclass(coords=SkyCoord(213., -10., unit=u.deg), **diskpar) photons = s.generate_photons(1e5) - n = np.empty(10) + n = np.empty(20) for i in range(len(n)): circ = SkyCoord((213. + np.random.uniform(-0.1, .1)) * u.degree, - (- 10. + np.random.uniform(-0.1, .1))*u.degree) + (- 10. + np.random.uniform(-0.1, .1)) * u.degree) d = circ.separation(SkyCoord(photons['ra'], photons['dec'], unit='deg')) n[i] = (d < 5. * u.arcmin).sum() s, p = normaltest(n) diff --git a/marxs/tests/test_analysis.py b/marxs/tests/test_analysis.py index 882417f..18b04d0 100644 --- a/marxs/tests/test_analysis.py +++ b/marxs/tests/test_analysis.py @@ -20,7 +20,7 @@ def test_detector_position(): n = 1000 convergent_point = np.array([3., 5., 7.]) pos = np.random.rand(n, 3) * 100. + 10. - dir = pos - convergent_point[np.newaxis, :] + dir = - pos + convergent_point[np.newaxis, :] photons = Table({'pos': e2h(pos, 1), 'dir': e2h(dir, 0), 'energy': np.ones(n), 'polarization': np.ones(n), 'probability': np.ones(n)}) opt = find_best_detector_position(photons) diff --git a/marxs/tests/test_simulator.py b/marxs/tests/test_simulator.py index 8b8f3d2..88491e4 100644 --- a/marxs/tests/test_simulator.py +++ b/marxs/tests/test_simulator.py @@ -66,7 +66,7 @@ def test_list_elemargs(): for p in [p0, p1, p2]: assert len(p.elements) == 2 for e in p.elements: - assert np.linalg.norm(e.geometry('v_y')) == 3 + assert np.linalg.norm(e.geometry['v_y']) == 3 assert p.elements[0]._d == 0.001 assert p.elements[1]._d == 0.002 diff --git a/marxs/visualization/mayavi.py b/marxs/visualization/mayavi.py index 6d3ac25..9b01267 100644 --- a/marxs/visualization/mayavi.py +++ b/marxs/visualization/mayavi.py @@ -54,7 +54,7 @@ def triangulation(obj, display, viewer=None): '''Plot a plane with an inner hole such as an aperture.''' from mayavi import mlab - xyz, triangles = obj.triangulate() + xyz, triangles = obj.geometry.triangulate(display) t = mlab.triangular_mesh(xyz[:, 0], xyz[:, 1], xyz[:, 2], triangles, color=display['color']) return t @@ -65,12 +65,13 @@ def surface(surface, display, viewer=None): The parameter boundaries are taken from the ``coo1`` and ``coo2`` in the display dictionary. The plotting routine is generic. It calls the ``parametric_surface()`` method of the object that is plotted; see there - for a detailted description of parameters. + for a detailed description of parameters. ''' from mayavi import mlab - xyz = surface.parametric_surface(display.get('coo1', [-1, 1]), - display.get('coo2', [-1, 1])) + xyz = surface.geometry.parametric_surface(display.get('coo1', None), + display.get('coo2', None), + display) xyz = mutils.h2e(xyz) x = xyz[..., 0] y = xyz[..., 1] @@ -129,11 +130,11 @@ def cylinder(obj, display, viewer=None): ''' from mayavi import mlab - x0 = obj.geometry('center') - obj.geometry('v_x') - x1 = obj.geometry('center') + obj.geometry('v_x') + x0 = obj.geometry['center'] - obj.geometry['v_x'] + x1 = obj.geometry['center'] + obj.geometry['v_x'] c = mlab.plot3d([x0[0], x1[0]], [x0[1], x1[1]], [x0[2], x1[2]], color=display['color'], - tube_radius=np.linalg.norm(obj.geometry('v_y')), + tube_radius=np.linalg.norm(obj.geometry['v_y']), tube_sides=display.get('tube_sides', 20)) return c diff --git a/marxs/visualization/tests/test_utils.py b/marxs/visualization/tests/test_utils.py index 5cb2f36..887f2d4 100644 --- a/marxs/visualization/tests/test_utils.py +++ b/marxs/visualization/tests/test_utils.py @@ -4,7 +4,8 @@ from ..utils import (plane_with_hole, combine_disjoint_triangulations, get_color, color_tuple_to_hex, - MARXSVisualizationWarning) + MARXSVisualizationWarning, + DisplayDict) from ..mayavi import plot_object def test_hole_round(): @@ -131,3 +132,38 @@ class Display(): assert len(record) == 1 # check that the message matches assert '"shape" not set in display dict.' in record[0].message.args[0] + +def test_DiplayDict(): + '''Simple mock-up using a DisplayDict variable.''' + class A(object): + pass + + a = A() + b = A() + b.geometry = A() + b.geometry.value = 5 + + a.display = DisplayDict(a, r=7) + b.display = DisplayDict(b, value_1='w') + + assert a.display['r'] == 7 + assert b.display['value_1'] == 'w' + assert b.display['value'] == 5 + # value is set in display and in geometry + b.display['value'] = 6 + assert b.display['value'] == 6 + assert b.geometry.value == 5 + # values added to geometry later + b.geometry.value2 = 'q' + assert b.display['value2'] == 'q' + # Make sure there is the right error only when used on an object without geomtry + with pytest.raises(KeyError) as e: + temp = a.display['q'] + # Make sure there is the right error if something is not found + with pytest.raises(KeyError) as e: + temp = b.display['1'] + + # get works, too + assert b.display.get('value') == 6 + assert a.display.get('qwer') is None + assert a.display.get('asdf', 5) == 5 diff --git a/marxs/visualization/threejs.py b/marxs/visualization/threejs.py index 4125434..43dc549 100644 --- a/marxs/visualization/threejs.py +++ b/marxs/visualization/threejs.py @@ -213,7 +213,7 @@ def box(obs, display, outfile): @format_doc(doc_plot) def triangulation(obj, display, outfile): '''Output commands for a plane with an inner hole.''' - xyz, triangles = obj.triangulate() + xyz, triangles = obj.geometry.triangulate(display) materialspec = materialspec(display, 'MeshStandardMaterial') outfile.write('// {}\n'.format(obj.name)) outfile.write('var geometry = new THREE.BufferGeometry(); \n') diff --git a/marxs/visualization/threejsjson.py b/marxs/visualization/threejsjson.py index 58d4ff8..360a0ab 100644 --- a/marxs/visualization/threejsjson.py +++ b/marxs/visualization/threejsjson.py @@ -114,7 +114,7 @@ def box(obj, display): def triangulation(obj, display): '''Describe a plane with a hole, such as an aperture of baffle.''' - xyz, triangles = obj.triangulate() + xyz, triangles = obj.geometry.triangulate(display) out = {} out['n'] = 1 out['name'] = str(obj.name) diff --git a/marxs/visualization/utils.py b/marxs/visualization/utils.py index 78a9a51..e1b7ca6 100644 --- a/marxs/visualization/utils.py +++ b/marxs/visualization/utils.py @@ -25,8 +25,50 @@ def get_obj_name(obj): return str(obj) +class DisplayDict(dict): + '''A dictionary to store how an element is displayed in plotting. + + A dictionary of this type works just like a normal dictionary, except for an + additional look-up step for keys that are not found in the dictionary itself. + A ``DisplayDict`` is initialized with a reference to the object it describes + and any parameters accessed from ``DisplayDict`` that is not found in the + dictionary, will be searched in the objects geometry. This allows us to set + any and all display settings in the ``DisplayDict`` to custamize plotting in any + way, but for those values that are not set, fall back to the settings of the + geometry (e.g. the shape of an object is typically taken from the geometry, + while the color is not). + + Parameters + ---------- + parent : `marxs.base.MarxsElement` + Reference to the object that is described by this ``DisplayDict`` + args, kwargs: see `dict` + ''' + def __init__(self, parent, *args, **kwargs): + self.parent = parent + super(DisplayDict, self).__init__(*args, **kwargs) + + def __getitem__(self, key): + if (key not in self) and hasattr(self.parent, 'geometry'): + try: + return getattr(self.parent.geometry, key) + except AttributeError: + raise KeyError(key) + else: + return super(DisplayDict, self).__getitem__(key) + + def get(self, k, d=None): + '''D.get(k[,d]) -> D[k] if k in D, else d. d defaults to None.''' + try: + return self[k] + except: + KeyError + return d + + + def plot_object_general(plot_registry, obj, display=None, **kwargs): - '''Look up a plottig routine for an object and execute it. + '''Look up a plotting routine for an object and execute it. This function is not meant to be called directly by the user, instead, it is designed to simplify the implementation of new plotting backends. @@ -62,25 +104,25 @@ def plot_object_general(plot_registry, obj, display=None, **kwargs): MARXSVisualizationWarning) return None - out = None - if 'shape' in display: + try: shape = display['shape'] - shapes = [s.strip() for s in shape.split(';')] - for s in shapes: - if s == 'None': - break - elif s in plot_registry: - # turn into valid color tuple - display['color'] = get_color(display) - out = plot_registry[s](obj, display, **kwargs) - break - else: - warnings.warn('Skipping {0}: No function to plot {1}.'.format(get_obj_name(obj), shape), - MARXSVisualizationWarning) - else: + except KeyError: warnings.warn('Skipping {0}: "shape" not set in display dict.'.format(get_obj_name(obj)), MARXSVisualizationWarning) - return out + return None + + shapes = [s.strip() for s in shape.split(';')] + for s in shapes: + if s == 'None': + return None + elif s in plot_registry: + # turn into valid color tuple + display['color'] = get_color(display) + return plot_registry[s](obj, display, **kwargs) + else: + warnings.warn('Skipping {0}: No function to plot {1}.'.format(get_obj_name(obj), shape), + MARXSVisualizationWarning) + return None def get_color(d): diff --git a/setup.cfg b/setup.cfg index 86383e9..65d3b95 100644 --- a/setup.cfg +++ b/setup.cfg @@ -9,7 +9,7 @@ libdir = # Examples could be #srcdir = /melkor/d1/guenther/marx/dist/ -#libdir = /melkor/d1/guenther/marx/dev/lib/ +#libdir = /melkor/d1/guenther/marx/installed/devPIC/lib/ #srcdir = /Users/hamogu/code/dist/ #libdir = /Users/hamogu/marx/devPIC/lib/