Some work on the docs.

doc/source/gaussian_beamlet_propagation.rst

@@ -44,4 +44,114 @@ Gausslets have their own set of predefined source-objects, found in :py:mod:`ray
 Evaluating the E-field
 ======================
 
-The nice thing about Gausslet ray-tracing is that you can evaluate the E-field at any point in your model. 
+The nice thing about Gausslet ray-tracing is that you can evaluate the E-field at any point in your model. For script-based analysis,
+you can give any GaussletCollection object (obtained from a source-object after a tracing operation) to the 
+ :py:function:`raypier.core.fields.eval_Efield_from_gausslets` function.
+
+ .. py:module:: raypier.core.fields
+ 
+     .. py:function:: eval_Efield_from_gausslets(gc : GaussletCollection, points : ndarray[:,3], wavelengths=None, blending=1.0) -> Efield ndarray[:,3]
+        :canonical: raypier.core.fields.eval_Efield_from_gausslets
+
+        Calculates the vector E-field is each of the points given. The returned 
+        array of field-vectors will have the same length as `points` and 
+        has `numpy.complex128` dtype.
+
+        :param GaussletCollection gc: The set of Gausslets for which the field should be calculated
+        :param ndarray[N,3] points: An array of shape (N,3) giving the points at which the field will be evaluated.
+        :param ndarray[] wavelengths: A 1d array containing the wavelengths to be used for the field calculation, 
+                                        overriding the wavelengths data contained by the GaussletCollection object.
+        :param float blending: The 1/width of each Gaussian mode at the evaluation points. A value of unity (the default),
+                                means the parabasal rays are determined to be the 1/e point in the field amplitude.
+
+    .. py:class:: EFieldSummation(gc : GaussletCollection, points : ndarray[:,3], wavelengths=None, blending=1.0) -> EFieldSummation object
+
+        For situations where you wish to evaluate the E-field from a set of Gausslets with different sets of evaluation points,
+        this class provides a small optimisation by performing the maths to convert ray-intercepts to Gaussian mode parameters
+        up front.
+
+        .. py:method:: evaluate(points : ndarray[:,3]) -> Efield ndarray[:,3]
+        
+            Called to calculate the E-field for the given points.
+
+
+Beam Decomposition
+==================
+
+When a Beam-decomposition object intercepts a ray during the tracing operation, instead of immediately generating
+child rays as most other `Traceable` objects do, the decomposition objects simply store the intercepted ray.
+At the completion of tracing of the current ray generation (i.e. GaussletCollection), any decomposition-objects
+which have received one or more rays then perform their decomposition-algorithm to generate a new set of rays to
+be added to the other rays created in the last generation. The new rays created by the decomposition process
+will, in general, not join originate from the end-point of the input-rays.
+
+
+.. py:module:: raypier.gausslets
+
+High level `Optics` objects for beam decomposition are provided here.
+
+.. py:class:: PositionDecompositionPlane(BaseDecompositionPlane)
+    :canonical: raypier.gausslets.PositionDecompositionPlane
+
+    Defines a plane at which position-decomposition will be beformed.
+
+    .. py:attribute:: radius 
+        :type: float
+        Sets the radius used for capturing incoming rays. Rays outside of this will "miss"
+
+    .. py:attribute:: curvature
+        :type: float
+        An approximate radius-of-curvature to the beam focus. This is used to improve the 
+        phase-unwrapping of the wavefront. The default is zero, which means a plane-wave
+        is assumed. Negative values imply a focus behind the decomposition plane (i.e.
+        on the opposite side to the plane direction vector).
+
+    .. py:attribute:: resolution
+        :type: float
+        Sets the resampling density of the decomposition, in terms of the number 
+        of new rays per `radius` extent.
+
+    .. py:attribute:: blending
+        :type: float
+        Sets the blending values for the new rays. The new rays will have Gaussian
+        1/e**2 intensity widths equal to `spacing`/`blending`, where the `spacing`
+        value is `radius`/`resolution`.
+
+ .. py:class:: AngleDecomposition(BaseDecompositionPlane)
+ 
+    Defines a plane at which Gabor (angle)-decomposition is to be performed.
+
+    .. py:attribute:: sample_spacing
+        :type: float
+        Sets the sample-spacing at the decomposition plane, in microns.
+
+    .. py:attribute:: width
+        :type: int
+        A value in the range 1->512 to set the number of samples along the width of the sample-plane.
+
+    .. py:attribute:: height
+        :type: int
+        A value in the range 1->512 to set the number of samples along the height of the sample-plane.
+
+    .. py:attribute:: mask
+        :type: ndarray[:,:]
+        A 2d array with shape matching the (width, height) and dtype numpy.float64 . The array
+        values should be in the range 0.0 -> 1.0. This will be used to mask the input E-field.
+
+    .. py:attribute:: max_angle
+        :type: float
+        Limits the angular divergence of the outgoing rays.
+
+
+.. py:module:: raypier.core.gausslets
+
+The low-level beam-decomposition algorithms are found in this module. Two types of decomposition are available: position-decomposition
+and angle-decomposition. Use the former when the Gausslets are found to be too wide at a particular surface in the optical path 
+to re-sample the beam onto a set of more compact Gausslets. The later is used to simulate the effect of apertures much smaller than
+the Gausslet widths, such that each Gausslet can be treated like a plane-wave and the field-distribution found using a 2d Fourier
+transform.
+
+
+
+
+

doc/source/introduction.rst

@@ -6,16 +6,16 @@ design tools for modelling optical systems (cameras, imaging systems, telescopes
 
 The main features of ray-trace are:
  - Non-sequential tracing (no need to specify the order of optical components)
+ - Physical optics propagation with "Gausslet" tracing and beam decomposition
  - Nice visualisation of the traced result
  - Live update to the traced result as the user adjusts the model
  - Reasonable performance (tracing algorithms runs at C-speed using Cython)
  - STEP export of models for integration with CAD design (using PythonOCC)
  - Saving / Loading models in YAML format.
  - Trace rays with full polarisation and phase information
- - Physical optics propagation with "Gausslet" tracing and beam decomposition 
+ - Define Zernike Polynomial sequences and apply them as distortions to surfaces
  - Dielectric Materials with simple-coating supported, including dispersion
  - A basic library of materials (from RefractiveIndex.info)
- - 
  - Various analysis algorithms including E-field evaluation by sum-of-Gaussian-Modes, and
    dispersion calculations for ultra-fast optics applications.
 

raypier/core/fields.py

@@ -240,9 +240,12 @@ def evaluate(self, points):
         return E
 
 
-def eval_Efield_from_gausslets(gausslet_collection, points, wavelengths,
+def eval_Efield_from_gausslets(gausslet_collection, points, 
+                               wavelengths = None,
                                blending=1.0, **kwds):
     gc = gausslet_collection.copy_as_array() 
+    if wavelengths is None:
+        wavelengths = numpy.asarray(gc.wavelengths)
     rays, x, y, dx, dy = evaluate_neighbours_gc(gc)
     modes = evaluate_modes_c(x, y, dx, dy, blending=blending)
     _rays = RayCollection.from_array(rays)

raypier/core/gausslets.py

@@ -201,7 +201,7 @@ def decompose_angle(origin, direction, axis1, E_field, input_spacing, max_angle,
 
     if E_max is not None and pos_max is not None:
         print("E_max:", E_max, "pos_max:", pos_max)
-        e_test = eval_Efield_from_gausslets(rays, pos_max[None,:], wl, blending=1.0)
+        e_test = eval_Efield_from_gausslets(rays, pos_max[None,:], blending=1.0)
         print("e_test:", e_test)
         scaling = (E_max/e_test).mean()
         #ray_data['E1_amp'] *= scaling
@@ -344,10 +344,11 @@ def decompose_position(input_rays, origin, direction, axis1, radius, resolution,
     ray_data['E2_amp'] = E2_amp
 
     gausslets = GaussletCollection.from_rays(ray_data)
+    gausslets.wavelengths = wavelengths
     gausslets.config_parabasal_rays(wavelengths, spacing/blending, 0.0)
     apply_mode_curvature(gausslets, -A, -B, -C)
 
-    E_test = eval_Efield_from_gausslets(gausslets, origins, wavelengths)
+    E_test = eval_Efield_from_gausslets(gausslets, origins)
 
     power_scaling = (E_in.real**2 + E_in.imag**2).sum() / (E_test.real**2 + E_test.imag**2).sum()
 

raypier/fields.py

@@ -171,7 +171,7 @@ def evaluate(self, src_list):
 
             if isinstance(rays, GaussletCollection):
                 n_list.append(rays.base_rays.refractive_index.real)
-                E = eval_Efield_from_gausslets(rays, points2, wavelengths, 
+                E = eval_Efield_from_gausslets(rays, points2,
                                                blending=self.blending)
             else:
                 n_list.append(rays.refractive_index.real)

raypier/gausslets.py

@@ -162,7 +162,7 @@ def evaluate_decomposed_rays(self, input_rays):
         points += origin[None,None,:]
 
         flat_points = points.reshape(-1,3)
-        E_field = eval_Efield_from_gausslets(input_rays, flat_points, wavelengths).reshape(*points.shape)
+        E_field = eval_Efield_from_gausslets(input_rays, flat_points).reshape(*points.shape)
         self._E_field = E_field
         mask = self._mask
         if mask is not None: