Separate out the geometry from the functional optics

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)

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()

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)

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

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.

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)

marxs/base/base.py

@@ -2,13 +2,15 @@
 from collections import OrderedDict
 import inspect
 import warnings
+from copy import deepcopy
 
 import numpy as np
 from transforms3d import affines
 
 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)

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,15 +68,15 @@ 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)
     r = np.sqrt(p['det_x'][ind]**2+p['det_y'][ind]**2.)
     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
 

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