Add a dispersion-compensator example and a modified michelson interfe…

…rometer

that includes arbitrary zernike polynomial distortions to create pretty ring
pattern.

examples/grating_dispersion_compensator_example.py

@@ -0,0 +1,93 @@
+import numpy as np
+
+from raypier.tracer import RayTraceModel
+from raypier.achromats import EdmundOptic45805
+from raypier.diffraction_gratings import RectangularGrating
+from raypier.mirrors import PECMirror, RectMirror
+from raypier.sources import BroadbandRaySource
+from raypier.beamstop import BeamStop
+from raypier.results import GroupVelocityDispersion
+from raypier.constraints import BaseConstraint
+from raypier.chirp_result import ChirpResult
+
+from traits.api import Range, on_trait_change
+from traitsui.api import View, Item
+
+
+direction = np.array([0.0,0.0,2.5]) - np.array([-29.685957,124.73850,1.5])
+
+source = BroadbandRaySource(origin=(-29.685957,124.73850,1.5),
+                            direction=tuple(direction),
+                            wavelength_start=0.76,
+                            wavelength_end = 0.80,
+                            number=260,
+                            uniform_deltaf=True,
+                            max_ray_len=300.0)
+
+
+grating = RectangularGrating(centre=(0,0,0),
+                             direction=(0,1,0),
+                             lines_per_mm=1400,
+                             order=-1)
+
+grating_init_rotation = 41.0
+grating.orientation = grating_init_rotation
+
+lens = EdmundOptic45805(centre=(0,75,0),
+                        direction=(0,-1,0))
+init_lens_rotation = lens.orientation
+
+mir = PECMirror(centre=(0,150,0),
+                direction=(0,1,0),
+                thickness=5)
+
+mir2 = RectMirror(centre=(-10.413843, 52.149718, -3.0),
+                  direction=(0.23150872699334185, -0.9727849212359605, -0.009654343160915735),
+                  width=6.0,
+                  length=10.0)
+
+bs = BeamStop(centre=(-35,150,15),
+              direction=(-0.2314983 ,  0.97131408,  0.05438228))
+
+gvd = GroupVelocityDispersion(source=source, target=bs.faces.faces[0])
+
+
+class GVDConstraint(BaseConstraint):
+    name = "Dispersion Compensator Adjustment"
+    focus_adjust = Range(70.0,80.0, value=73.5)
+    gdd_adjust = Range(-20.0,40.0, value=0.0)
+    lens_rotation = Range(-20.0, 20.0, value=0.0)
+    grating_angle = Range(-5.0,5.0, value=0.0)
+
+    traits_view = View(Item("focus_adjust"),
+                       Item("gdd_adjust"),
+                       Item("lens_rotation"),
+                       Item("grating_angle"),
+                       resizable=True)
+
+    @on_trait_change("focus_adjust, gdd_adjust")
+    def on_new_values(self):
+        lens.centre = (0.0, 75.0 + self.gdd_adjust, 0.0)
+        mir.centre = (0.0, 75.0 + self.gdd_adjust + self.focus_adjust, 0.0)
+        self.update = True
+
+    def _lens_rotation_changed(self):
+        lens.orientation = (init_lens_rotation + self.lens_rotation+180.0)%360.0 - 180.0
+        self.update = True
+
+    def _grating_angle_changed(self):
+        grating.orientation = grating_init_rotation + self.grating_angle
+        self.update = True
+
+
+gvd_cstr = GVDConstraint()
+
+
+chirp = ChirpResult(source=source, target=bs.faces.faces[0])
+
+model = RayTraceModel(optics=[grating, lens, mir, mir2, bs], 
+                      sources=[source,],
+                      results=[gvd, chirp],
+                      constraints=[gvd_cstr,])
+
+model.configure_traits(kind="live")

examples/zernike_demo_interferometer_example.py

@@ -0,0 +1,73 @@
+"""
+A modification of the original Michelson interferometer example where the spherical face has been
+replaced by a plane-flat with a Zernike polynomial distortion applied.
+
+Any number of Zernike poly terms can be applied. The ANSI standard J-indexing scheme is used.
+Use the "Faces" tab of the Propterties window of the Distorted Mirror. Zernike coefficient
+amplitudes can be added, updated or removed "live" in the GUI.
+
+"""
+
+from raypier.api import RayTraceModel, UnpolarisingBeamsplitterCube, CollimatedGaussletSource, CircleShape, GeneralLens,\
+        PlanarFace, GaussletCapturePlane, EFieldPlane, IntensitySurface
+from raypier.intensity_image import IntensityImageView
+from raypier.distortions import ZernikeSeries
+from raypier.faces import DistortionFace
+
+
+src=CollimatedGaussletSource(origin=(-30,0,0),
+                             direction=(1,0,0),
+                             radius=5.0,
+                             beam_waist=10.0,
+                             resolution=10.0,
+                             E_vector=(0,1,0),
+                             wavelength=1.0,
+                             display="wires",
+                             opacity=0.05,
+                             max_ray_len=50.0)
+
+
+bs = UnpolarisingBeamsplitterCube(centre=(0,0,0),
+                                  size=10.0)
+
+shape = CircleShape(radius=10.0)
+
+f1 = PlanarFace(mirror=True)
+
+dist = ZernikeSeries(unit_radius=5.0, coefficients=[(3,0.0005),(7,0.0001)])
+f2 = PlanarFace(mirror=True) #SphericalFace(curvature=2000.0, mirror=True)
+f2 = DistortionFace(base_face=f2, mirror=True, distortion=dist)
+
+m1 = GeneralLens(name="Flat Mirror",
+                 shape=shape,
+                 surfaces=[f1],
+                 centre=(0,20,0),
+                 direction=(0,-1,0))
+
+m2 = GeneralLens(name="Distorted Mirror",
+                 shape=shape,
+                 surfaces=[f2],
+                 centre=(20,0,0),
+                 direction=(-1,0,0))
+
+cap = GaussletCapturePlane(centre=(0,-20,0),
+                           direction=(0,1,0))
+
+field = EFieldPlane(centre=(0,-20,0),
+                    direction=(0,1,0),
+                    detector=cap,
+                    align_detector=True,
+                    width=10.0,
+                    height=10.0,
+                    size=100)
+
+image = IntensityImageView(field_probe=field)
+surf = IntensitySurface(field_probe=field)
+
+
+model = RayTraceModel(optics=[bs, m1, m2],
+                      sources=[src],
+                      probes=[field, cap],
+                      results=[image, surf])
+
+model.configure_traits()

raypier/results.py

@@ -107,9 +107,8 @@ def evaluate_phase(all_wavelengths, traced_rays, target_face,
     returns - (freq, phase) #freq in THz
     """
     c = 2.99792458e8 * 1e-9 #convert to mm/ps
-    all_rays = [r.copy_as_array() for r in reversed(traced_rays)]
     idx = target_face.idx
-    last = all_rays[0]
+    last = traced_rays[-1].copy_as_array() #all_rays[0]
     selected_idx = numpy.argwhere(last['end_face_idx']==idx).ravel()
     wavelengths = all_wavelengths[last['wavelength_idx'][selected_idx]]
     sort_idx = numpy.argsort(wavelengths)[::-1]
@@ -119,13 +118,9 @@ def evaluate_phase(all_wavelengths, traced_rays, target_face,
     idx = len(selected_idx)//2
     phase = last['phase'][selected_idx].copy()
     phase -= phase.mean()
-    total = numpy.zeros(len(selected_idx), 'd')
-    for ray in all_rays:
-        selected = ray[selected_idx]
-        total += selected['length'] * selected['refractive_index'].real
-        selected_idx = selected['parent_idx']
-
-    #print "Phase:", phase
+
+    total = last['accumulated_path'][selected_idx]
+
     if len(total) < 6:
         raise ValueError("Not enough rays to evaluate 2nd and 3rd derivatives")
     ave_path = total.mean()
@@ -137,7 +132,6 @@ def evaluate_phase(all_wavelengths, traced_rays, target_face,
     fs_total -= fs_total.mean()
     total += fs_total
 
-    #phase += 2*numpy.pi*total/(0.001*wavelengths) #total*omega/c
     phase += total*f*((2*numpy.pi)/c)
     return (f, phase)