From 0f87cfe114fdd67d16fc505454c86c5293fa65d4 Mon Sep 17 00:00:00 2001
From: Viljar Femoen <viljar.femoen@hotmail.no>
Date: Mon, 2 Dec 2024 15:15:16 +0100
Subject: [PATCH 1/4] Refactor simulation, cleanup initialization

---
 src/instamatic/controller.py                  | 32 +++++++++++--------
 src/instamatic/microscope/microscope.py       |  2 +-
 src/instamatic/simulation/__init__.py         | 29 +++++++++++++++++
 src/instamatic/simulation/camera.py           | 27 +++-------------
 .../microscope.py}                            |  2 +-
 5 files changed, 53 insertions(+), 39 deletions(-)
 rename src/instamatic/{microscope/interface/simu_microscope.py => simulation/microscope.py} (99%)

diff --git a/src/instamatic/controller.py b/src/instamatic/controller.py
index 66a022d3..b00c5033 100644
--- a/src/instamatic/controller.py
+++ b/src/instamatic/controller.py
@@ -16,6 +16,7 @@
 from instamatic.microscope import components
 from instamatic.microscope.base import MicroscopeBase
 from instamatic.microscope.microscope import get_microscope
+from instamatic.simulation import get_simulation_instances
 
 _ctrl = None  # store reference of ctrl so it can be accessed without re-initializing
 
@@ -46,22 +47,25 @@ def initialize(
     ctrl : `TEMController`
         Return TEM control object
     """
-    print(f"Microscope: {tem_name}{' (server)' if use_tem_server else ''}")
-    tem = get_microscope(tem_name, use_server=use_tem_server)
-
-    if cam_name:
-        if use_cam_server:
-            cam_tag = ' (server)'
-        elif stream:
-            cam_tag = ' (stream)'
-        else:
-            cam_tag = ''
-
-        print(f'Camera    : {cam_name}{cam_tag}')
 
-        cam = Camera(cam_name, as_stream=stream, use_server=use_cam_server)
+    if tem_name == 'simulate' or cam_name == 'simulate':
+        cam, tem = get_simulation_instances()
     else:
-        cam = None
+        print(f"Microscope: {tem_name}{' (server)' if use_tem_server else ''}")
+        tem = get_microscope(tem_name, use_server=use_tem_server)
+        if cam_name:
+            if use_cam_server:
+                cam_tag = ' (server)'
+            elif stream:
+                cam_tag = ' (stream)'
+            else:
+                cam_tag = ''
+
+            print(f'Camera    : {cam_name}{cam_tag}')
+
+            cam = Camera(cam_name, as_stream=stream, use_server=use_cam_server)
+        else:
+            cam = None
 
     global _ctrl
     ctrl = _ctrl = TEMController(tem=tem, cam=cam)
diff --git a/src/instamatic/microscope/microscope.py b/src/instamatic/microscope/microscope.py
index c7a0d8ab..dde3aa9d 100644
--- a/src/instamatic/microscope/microscope.py
+++ b/src/instamatic/microscope/microscope.py
@@ -19,7 +19,7 @@ def get_microscope_class(interface: str) -> 'type[MicroscopeBase]':
             raise PermissionError('Access to the TEM interface requires admin rights.')
 
     if simulate or interface == 'simulate':
-        from .interface.simu_microscope import SimuMicroscope as cls
+        from ..simulation.microscope import MicroscopeSimulation as cls
     elif interface == 'jeol':
         from .interface.jeol_microscope import JeolMicroscope as cls
     elif interface == 'fei':
diff --git a/src/instamatic/simulation/__init__.py b/src/instamatic/simulation/__init__.py
index e69de29b..75cc4c27 100644
--- a/src/instamatic/simulation/__init__.py
+++ b/src/instamatic/simulation/__init__.py
@@ -0,0 +1,29 @@
+from __future__ import annotations
+
+from typing import Tuple
+
+from instamatic import config
+from instamatic.camera.camera_client import CamClient
+from instamatic.camera.videostream import VideoStream
+from instamatic.microscope.client import MicroscopeClient
+from instamatic.simulation.camera import CameraSimulation
+from instamatic.simulation.microscope import MicroscopeSimulation
+
+
+def get_simulation_instances() -> Tuple[CameraSimulation, MicroscopeSimulation]:
+    """Initialize simulated camera and microscope.
+
+    Returns
+    -------
+    Tuple[CameraSimulation, MicroscopeSimulation]
+    """
+    if config.settings.use_tem_server:
+        tem = MicroscopeClient(interface='simulate')
+    else:
+        tem = MicroscopeSimulation()
+    if config.settings.use_cam_server:
+        cam = CamClient(name='simulate', interface=config.camera.interface)
+    else:
+        camsim = CameraSimulation(tem=tem)
+        cam = VideoStream(camsim)
+    return cam, tem
diff --git a/src/instamatic/simulation/camera.py b/src/instamatic/simulation/camera.py
index 29ccd2f8..b74e8d48 100644
--- a/src/instamatic/simulation/camera.py
+++ b/src/instamatic/simulation/camera.py
@@ -5,16 +5,17 @@
 from numpy import ndarray
 
 from instamatic.camera.camera_base import CameraBase
+from instamatic.microscope.base import MicroscopeBase
 from instamatic.simulation.stage import Stage
 
 
 class CameraSimulation(CameraBase):
     streamable = True
 
-    def __init__(self, name: str = 'simulate'):
+    def __init__(self, tem: MicroscopeBase, name: str = 'simulate'):
         super().__init__(name)
 
-        self.ready = False
+        self.tem = tem
 
         # TODO put parameters into config
         self.stage = Stage()
@@ -23,30 +24,10 @@ def __init__(self, name: str = 'simulate'):
     def establish_connection(self):
         pass
 
-    def actually_establish_connection(self):
-        if self.ready:
-            return
-        import time
-
-        time.sleep(2)
-        from instamatic.controller import get_instance
-
-        ctrl = get_instance()
-        self.tem = ctrl.tem
-
-        ctrl.stage.set(z=0, a=0, b=0)
-        print(self.tem.getStagePosition())
-        print(self.stage.samples[0].x, self.stage.samples[0].y)
-
-        self.ready = True
-
     def release_connection(self):
-        self.tem = None
-        self.ready = False
+        pass
 
     def get_image(self, exposure: float = None, binsize: int = None, **kwargs) -> ndarray:
-        self.actually_establish_connection()
-
         if exposure is None:
             exposure = self.default_exposure
         if binsize is None:
diff --git a/src/instamatic/microscope/interface/simu_microscope.py b/src/instamatic/simulation/microscope.py
similarity index 99%
rename from src/instamatic/microscope/interface/simu_microscope.py
rename to src/instamatic/simulation/microscope.py
index be37c70e..b3471da5 100644
--- a/src/instamatic/microscope/interface/simu_microscope.py
+++ b/src/instamatic/simulation/microscope.py
@@ -33,7 +33,7 @@
 MIN = 0
 
 
-class SimuMicroscope(MicroscopeBase):
+class MicroscopeSimulation(MicroscopeBase):
     """Simulates a microscope connection.
 
     Has the same variables as the real JEOL/FEI equivalents, but does

From d93c299c4c57808835eb3170ca2e456b6be93251 Mon Sep 17 00:00:00 2001
From: Viljar Femoen <viljar.femoen@hotmail.no>
Date: Wed, 15 Jan 2025 15:52:10 +0100
Subject: [PATCH 2/4] Properly handle rotation and eucentric height

---
 src/instamatic/simulation/stage.py | 23 +++++++++++++++++------
 1 file changed, 17 insertions(+), 6 deletions(-)

diff --git a/src/instamatic/simulation/stage.py b/src/instamatic/simulation/stage.py
index a5b57036..b1cd1f9f 100644
--- a/src/instamatic/simulation/stage.py
+++ b/src/instamatic/simulation/stage.py
@@ -38,9 +38,8 @@ def __init__(
         self.z = 0
         self.alpha_tilt = 0
         self.beta_tilt = 0
-        self.in_plane_rotation = 0  # TODO change this with focus/magnification
+        self.in_plane_rotation = 10  # TODO change this with focus/magnification
         self.rotation_matrix = np.eye(3)
-        self.origin = np.array([0, 0, 0])
 
         # TODO parameters
         self.grid = Grid()
@@ -65,6 +64,10 @@ def __init__(
             for _ in range(num_crystals)
         ]
 
+    @property
+    def origin(self) -> np.ndarray:
+        return np.array([self.x, self.y, self.z])
+
     def set_position(
         self,
         x: float = None,
@@ -135,21 +138,21 @@ def image_extent_to_sample_coordinates(
                 'Tilting is not fully implemented yet', NotImplementedWarning, stacklevel=2
             )
         # https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection
-        n = self.rotation_matrix @ np.array([0, 0, 1])
-        p0 = self.origin
         l = np.array([0, 0, 1])  # noqa: E741
+        n = self.rotation_matrix @ l
+        p0 = self.origin
         l0 = np.array(
             [
                 p.flatten()
                 for p in np.meshgrid(
                     np.linspace(x_min, x_max, shape[1]),
                     np.linspace(y_min, y_max, shape[0]),
-                    [0],
+                    [self.z],
                 )
             ]
         )
 
-        p = l0 + np.array([0, 0, 1])[:, np.newaxis] * np.dot(-l0.T + p0, n) / np.dot(l, n)
+        p = l0 + l[:, np.newaxis] * np.dot(-l0.T + p0, n) / np.dot(l, n)
 
         x, y, z = self.rotation_matrix.T @ p
         x = x.reshape(shape)
@@ -188,6 +191,10 @@ def get_image(
         x, y = self.image_extent_to_sample_coordinates(
             shape=shape, x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max
         )
+        x_min = x.min()
+        x_max = x.max()
+        y_min = y.min()
+        y_max = y.max()
 
         grid_mask = self.grid.array_from_coords(x, y)
 
@@ -239,6 +246,10 @@ def get_diffraction_pattern(
         x, y = self.image_extent_to_sample_coordinates(
             shape=shape, x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max
         )
+        x_min = x.min()
+        x_max = x.max()
+        y_min = y.min()
+        y_max = y.max()
         d_min = 1.0
 
         grid_mask = self.grid.array_from_coords(x, y)

From 5db04e6557782e77786416415319f425e2c4db0f Mon Sep 17 00:00:00 2001
From: Viljar Femoen <viljar.femoen@hotmail.no>
Date: Wed, 15 Jan 2025 16:50:24 +0100
Subject: [PATCH 3/4] Use diffsims for simulations

---
 src/instamatic/simulation/crystal.py | 114 ++++++++++++++++++---------
 src/instamatic/simulation/stage.py   |  27 ++-----
 2 files changed, 80 insertions(+), 61 deletions(-)

diff --git a/src/instamatic/simulation/crystal.py b/src/instamatic/simulation/crystal.py
index 13875be3..f98b3501 100644
--- a/src/instamatic/simulation/crystal.py
+++ b/src/instamatic/simulation/crystal.py
@@ -4,13 +4,28 @@
 
 import numpy as np
 from diffpy import structure as diffpy
+from diffsims.crystallography._diffracting_vector import DiffractingVector
+from diffsims.generators.simulation_generator import (
+    Simulation2D,
+    SimulationGenerator,
+    get_kinematical_intensities,
+)
+from orix.crystal_map import Phase
+from orix.quaternion import Rotation
 
 Crystal_T = TypeVar('Crystal_T', bound='Crystal')
 
 
 class Crystal:
     def __init__(
-        self, a: float, b: float, c: float, alpha: float, beta: float, gamma: float
+        self,
+        a: float,
+        b: float,
+        c: float,
+        alpha: float,
+        beta: float,
+        gamma: float,
+        space_group: int = 1,
     ) -> None:
         """Simulate a primitive crystal given the unit cell. No additional
         symmetry is imposed.
@@ -31,6 +46,8 @@ def __init__(
             Angle between a and c, in degrees
         gamma : float
             Angle between a and b, in degrees
+        space_group: int
+            Space group number. Defaults to 1 (P1)
         """
         self.a = a
         self.b = b
@@ -41,9 +58,10 @@ def __init__(
 
         self.lattice = diffpy.Lattice(self.a, self.b, self.c, self.alpha, self.beta, self.gamma)
         self.structure = diffpy.Structure(
-            atoms=[diffpy.Atom(xyz=[0, 0, 0])],
+            atoms=[diffpy.Atom('Au', xyz=[0, 0, 0])],
             lattice=self.lattice,
         )
+        self.phase = Phase(space_group=space_group, structure=self.structure)
 
     @property
     def a_vec(self) -> np.ndarray:
@@ -133,11 +151,11 @@ def diffraction_pattern_mask(
         shape: tuple[int, int],
         d_min: float,
         rotation_matrix: np.ndarray,
-        wavelength: float,
+        acceleration_voltage: float,
         excitation_error: float,
     ) -> np.ndarray:
         """Get a diffraction pattern with a given shape, up to a given
-        resolution, in a given orientation and wavelength.
+        resolution, in a given orientation and acceleration voltage.
 
         Parameters
         ----------
@@ -147,8 +165,8 @@ def diffraction_pattern_mask(
             Minimum d-spacing, in Å
         rotation_matrix : np.ndarray
             Orientation
-        wavelength : float
-            Wavelength of incident beam, in Å
+        acceleration_voltage : float
+            acceleration_voltage of incident beam, in kV
         excitation_error : float
             Excitation error used for intensity calculation, in reciprocal Å
 
@@ -157,39 +175,57 @@ def diffraction_pattern_mask(
         np.ndarray
             Diffraction pattern
         """
-        # TODO calibration
-        out = np.zeros(shape, dtype=bool)
-
-        # TODO this depends on convergence angle
-        spot_radius = 3  # pixels
-
-        vecs = self.reciprocal_space_lattice(d_min)
-        d = np.sum(vecs**2, axis=1)
-        vecs = vecs[d < d_min**2]
-
-        k = 2 * np.pi / wavelength
-        k_vec = rotation_matrix @ np.array([0, 0, -k])
-
-        # Find intersect with Ewald's sphere
-        q_squared = np.sum((vecs - k_vec) ** 2, axis=1)
-        vecs = vecs[
-            (q_squared > (k - excitation_error) ** 2)
-            & (q_squared < (k + excitation_error) ** 2)
-        ]
-
-        # Project onto screen
-        vecs_xy = (rotation_matrix.T @ vecs.T).T[:, :-1]  # ignoring curvature
-
-        # Make image
-        for vec in vecs_xy:
-            x = int(vec[0] * d_min * shape[1] / 2) + shape[1] // 2
-            y = int(vec[1] * d_min * shape[0] / 2) + shape[0] // 2
-            min_x = max(0, x - spot_radius)
-            max_x = min(shape[1], x + spot_radius)
-            min_y = max(0, y - spot_radius)
-            max_y = min(shape[0], y + spot_radius)
-            out[min_y:max_y, min_x:max_x] = 1
-        return out
+        gen = SimulationGenerator(accelerating_voltage=acceleration_voltage)
+        wavelength = gen.wavelength
+
+        # Rotate using all the rotations in the list
+        recip = DiffractingVector.from_min_dspacing(
+            self.phase,
+            min_dspacing=d_min,
+            include_zero_vector=False,
+        )
+        rotation = Rotation.from_matrix(rotation_matrix)
+        # Calculate the reciprocal lattice vectors that intersect the Ewald sphere.
+        (
+            intersected_vectors,
+            hkl,
+            shape_factor,
+        ) = gen.get_intersecting_reflections(
+            recip,
+            rotation,
+            wavelength,
+            max_excitation_error=excitation_error,
+            with_direct_beam=False,
+        )
+
+        # Calculate diffracted intensities based on a kinematic model.
+        intensities = get_kinematical_intensities(
+            self.structure,
+            hkl,
+            intersected_vectors.gspacing,
+            prefactor=shape_factor,
+            scattering_params=gen.scattering_params,
+        )
+
+        # Threshold peaks included in simulation as factor of zero beam intensity.
+        peak_mask = intensities > np.max(intensities) * gen.minimum_intensity
+        intensities = intensities[peak_mask]
+        intersected_vectors = intersected_vectors[peak_mask]
+        intersected_vectors.intensity = intensities
+
+        # Create a simulation object
+        sim = Simulation2D(
+            phases=self.phase,
+            coordinates=intersected_vectors,
+            rotations=rotation,
+            simulation_generator=gen,
+            reciprocal_radius=1 / d_min,
+        )
+
+        # Simulate diffraction pattern
+        return sim.get_diffraction_pattern(
+            shape, sigma=1, calibration=1 / d_min / (shape[0] / 2)
+        )
 
     def __str__(self) -> str:
         return f'{self.__class__.__name__}(a = {self.a}, b = {self.b}, c = {self.c}, alpha = {self.alpha}, beta = {self.beta}, gamma = {self.gamma})'
diff --git a/src/instamatic/simulation/stage.py b/src/instamatic/simulation/stage.py
index b1cd1f9f..28f59b0d 100644
--- a/src/instamatic/simulation/stage.py
+++ b/src/instamatic/simulation/stage.py
@@ -258,12 +258,7 @@ def get_diffraction_pattern(
             # no transmission
             return np.zeros(shape, dtype=int)
 
-        reflections = np.zeros(shape, dtype=bool)
-
-        # Direct beam
-        reflections[
-            shape[0] // 2 - 4 : shape[0] // 2 + 4, shape[1] // 2 - 4 : shape[1] // 2 + 4
-        ] = 1
+        reflections = np.zeros(shape, dtype=float)
 
         for sample in self.samples:
             if not sample.range_might_contain_crystal(
@@ -275,30 +270,18 @@ def get_diffraction_pattern(
                 # Crystal is completely on the grid
                 continue
 
-            reflections |= self.crystal.diffraction_pattern_mask(
+            reflections += self.crystal.diffraction_pattern_mask(
                 shape,
                 d_min=d_min,
                 rotation_matrix=self.rotation_matrix @ sample.rotation_matrix,
-                wavelength=0.02,
+                acceleration_voltage=200,
                 excitation_error=0.01,
             )
         # TODO diffraction shift
 
         # TODO noise
 
-        # Simple scaling
-        # TODO improve, proper form factors maybe
-        # TODO camera length
-        kx, ky = np.meshgrid(
-            np.linspace(-1 / d_min, 1 / d_min, shape[1]),
-            np.linspace(-1 / d_min, 1 / d_min, shape[0]),
-        )
-        k_squared = kx**2 + ky**2
-        scale = 1 / (3 * k_squared + 1)
-
-        scale[~reflections] = 0
-
         # Convert to int array
-        scale = (scale * 0x8000).astype(int)
+        reflections = (reflections * 0x8000).astype(int)
 
-        return scale
+        return reflections

From 9f5e5ea4f90696f6cbe77f816609d3d717128a54 Mon Sep 17 00:00:00 2001
From: Viljar Femoen <viljar.femoen@hotmail.no>
Date: Thu, 4 Sep 2025 11:16:06 +0200
Subject: [PATCH 4/4] Speed up kinematic simulations

---
 pyproject.toml                                |  3 +
 .../calibrate/calibrate_beamshift.py          |  2 +-
 src/instamatic/camera/fakevideostream.py      |  2 +-
 src/instamatic/controller.py                  | 41 ++++++++--
 .../experiments/autocred/experiment.py        | 19 +++--
 src/instamatic/experiments/red/experiment.py  |  6 +-
 src/instamatic/gui/autocred_frame.py          |  2 -
 src/instamatic/gui/ctrl_frame.py              | 18 +++++
 src/instamatic/processing/find_crystals.py    |  4 +-
 src/instamatic/simulation/__init__.py         |  5 +-
 src/instamatic/simulation/camera.py           | 13 +++-
 src/instamatic/simulation/crystal.py          | 74 +++++++++++--------
 src/instamatic/simulation/microscope.py       |  2 +
 src/instamatic/simulation/stage.py            | 39 +++++-----
 14 files changed, 155 insertions(+), 75 deletions(-)

diff --git a/pyproject.toml b/pyproject.toml
index 7c856aee..93086a1b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -88,6 +88,9 @@ publishing = [
     "wheel",
     "build",
 ]
+simulation = [
+    "diffsims",
+]
 
 [project.scripts]
 "instamatic" = "instamatic.main:main"
diff --git a/src/instamatic/calibrate/calibrate_beamshift.py b/src/instamatic/calibrate/calibrate_beamshift.py
index 1fa69725..4182070a 100644
--- a/src/instamatic/calibrate/calibrate_beamshift.py
+++ b/src/instamatic/calibrate/calibrate_beamshift.py
@@ -108,7 +108,7 @@ def plot(self, to_file=None, outdir=''):
 
     def center(self, ctrl):
         """Return beamshift values to center the beam in the frame."""
-        pixel_center = [val / 2.0 for val in ctrl.cam.get_image_dimensions()]
+        pixel_center = [val / 2.0 for val in ctrl.cam.get_camera_dimensions()]
 
         beamshift = self.pixelcoord_to_beamshift(pixel_center)
         if ctrl:
diff --git a/src/instamatic/camera/fakevideostream.py b/src/instamatic/camera/fakevideostream.py
index a253caa7..605bec71 100644
--- a/src/instamatic/camera/fakevideostream.py
+++ b/src/instamatic/camera/fakevideostream.py
@@ -51,7 +51,7 @@ def __getattr__(self, attrname):
     def get_image(self, exposure=None, binsize=None):
         frame = self.cam.get_image(exposure=exposure, binsize=binsize)
 
-        self.frame, scale = autoscale(frame, maxdim=self.display_dim)
+        # self.frame, scale = autoscale(frame, maxdim=self.display_dim)
 
         return frame
 
diff --git a/src/instamatic/controller.py b/src/instamatic/controller.py
index b00c5033..3e757d1a 100644
--- a/src/instamatic/controller.py
+++ b/src/instamatic/controller.py
@@ -388,11 +388,11 @@ def find_eucentric_height(
 
         def one_cycle(tilt: float = 5, sign=1) -> list:
             angle1 = -tilt * sign
-            self.stage.a = angle1
+            self.stage.set(a=angle1)
             img1 = self.get_rotated_image()
 
             angle2 = +tilt * sign
-            self.stage.a = angle2
+            self.stage.set(a=angle2)
             img2 = self.get_rotated_image()
 
             if sign < 1:
@@ -402,7 +402,7 @@ def one_cycle(tilt: float = 5, sign=1) -> list:
 
             return shift
 
-        self.stage.a = 0
+        self.stage.set(a=0)
         # self.stage.z = 0 # for testing
 
         zc = self.stage.z
@@ -414,7 +414,7 @@ def one_cycle(tilt: float = 5, sign=1) -> list:
         sign = 1
 
         for i, z in enumerate(zs):
-            self.stage.z = z
+            self.stage.set(z=z)
             if verbose:
                 print(f'z = {z:.1f} nm')
 
@@ -423,9 +423,7 @@ def one_cycle(tilt: float = 5, sign=1) -> list:
 
             sign *= -1
 
-        mean_shift = shifts[-1] + shifts[0]
-        mean_shift = mean_shift / np.linalg.norm(mean_shift)
-        ds = np.dot(shifts, mean_shift)
+        ds = np.sum(shifts, axis=1)
 
         p = np.polyfit(zs, ds, 1)  # linear fit
         alpha, beta = p
@@ -596,6 +594,35 @@ def get_rotated_image(self, exposure: float = None, binsize: int = None) -> np.n
 
         return arr
 
+    def get_rotated_raw_image(self, exposure: float = None, binsize: int = None) -> np.ndarray:
+        """Simplified function equivalent to `get_image` that returns the
+        rotated image array.
+
+        Parameters
+        ----------
+        exposure: float
+            Exposure time in seconds
+        binsize: int
+            Binning to use for the image, must be 1, 2, or 4, etc
+        mode : str
+            Magnification mode
+        mag : int
+            Magnification value
+
+        Returns
+        -------
+        arr : np.array
+            Image as 2D numpy array.
+        """
+
+        mag = self.magnification.value
+        mode = self.mode.get()
+
+        arr = self.get_raw_image(exposure=exposure, binsize=binsize)
+        arr = rotate_image(arr, mode=mode, mag=mag)
+
+        return arr
+
     def get_image(
         self,
         exposure: float = None,
diff --git a/src/instamatic/experiments/autocred/experiment.py b/src/instamatic/experiments/autocred/experiment.py
index 5089aef3..8860e7c1 100644
--- a/src/instamatic/experiments/autocred/experiment.py
+++ b/src/instamatic/experiments/autocred/experiment.py
@@ -252,6 +252,8 @@ def __init__(
             except BaseException:
                 print('Is DIALS server running? Connection failed.')
                 self.s_c = 0
+        else:
+            self.s_c = 0
 
         if use_vm:
             self.s2 = socket.socket()
@@ -264,6 +266,12 @@ def __init__(
             except BaseException:
                 print('Is VM server running? Connection failed.')
                 self.s2_c = 0
+        else:
+            self.s2_c = 0
+
+    @property
+    def magnification(self) -> float:
+        return self.ctrl.magnification.value
 
     def image_cropper(self, img, window_size=0):
         crystal_pos, r = find_defocused_image_center(
@@ -1152,7 +1160,7 @@ def raster_scan(self):
             config.calibration['mag1']['pixelsize'][self.magnification] / 1000
         )  # nm -> um
         xdim, ydim = self.ctrl.cam.get_camera_dimensions()
-        box_x, box_y = self.pixelsize_mag1 * xdim, self.pixelsize_mag1 * ydim
+        box_x, box_y = pixelsize_mag1 * xdim, pixelsize_mag1 * ydim
 
         # Make negative to reflect config change 2019-07-03 to make omega more in line with other software
         rot_axis = -config.camera.camera_rotation_vs_stage_xy
@@ -1185,7 +1193,6 @@ def raster_scan(self):
                 try:
                     img, h = self.ctrl.get_image(exposure=self.expt, header_keys=None)
                     if img.mean() > 10:
-                        self.magnification = self.ctrl.magnification.value
                         crystal_positions = find_crystals_timepix(
                             img, self.magnification, spread=self.spread, offset=self.offset
                         )
@@ -1380,8 +1387,6 @@ def start_collection_point(self):
                 y=self.calib_beamshift.reference_shift[1],
             )
 
-            self.magnification = self.ctrl.magnification.value
-
             self.rotation_direction = self.eliminate_backlash_in_tiltx()
             img, h = self.ctrl.get_image(exposure=self.expt, header_keys=header_keys)
             # img, h = Experiment.apply_corrections(img, h)
@@ -1690,9 +1695,9 @@ def start_collection(self):
                 input(
                     'No z-height adjustment found. Please find an area with particles! Press Enter to continue auto adjustment of z height>>>'
                 )
-                x_zheight, y_zheight = center_z_height_HYMethod(
-                    self.ctrl, spread=self.spread, offset=self.offset
-                )
+                # x_zheight, y_zheight = center_z_height_HYMethod(
+                #     self.ctrl, spread=self.spread, offset=self.offset
+                # )
                 xpoint, ypoint, zpoint, aaa, bbb = self.ctrl.stage.get()
                 self.logger.info(
                     f'Stage position: x = {xpoint}, y = {ypoint}. Z height adjusted to {zpoint}. Tilt angle x {aaa} deg, Tilt angle y {bbb} deg'
diff --git a/src/instamatic/experiments/red/experiment.py b/src/instamatic/experiments/red/experiment.py
index 3b01c1a2..ca3c10d5 100644
--- a/src/instamatic/experiments/red/experiment.py
+++ b/src/instamatic/experiments/red/experiment.py
@@ -98,12 +98,12 @@ def start_collection(self, exposure_time: float, tilt_range: float, stepsize: fl
             ctrl.cam.block()
 
         # for i, a in enumerate(tilt_positions):
-        for i, angle in enumerate(tqdm(tilt_positions)):
-            ctrl.stage.a = angle
+        for i, angle in enumerate((tilt_positions)):
+            ctrl.stage.set(a=angle)
 
             j = i + self.offset
 
-            img, h = self.ctrl.get_image(exposure_time)
+            img, h = ctrl.get_image(exposure_time, header_keys='')
 
             self.buffer.append((j, img, h))
 
diff --git a/src/instamatic/gui/autocred_frame.py b/src/instamatic/gui/autocred_frame.py
index 857ee48f..daa028af 100644
--- a/src/instamatic/gui/autocred_frame.py
+++ b/src/instamatic/gui/autocred_frame.py
@@ -385,8 +385,6 @@ def acquire_data_autocRED(controller, **kwargs):
     )
     cexp.start_collection()
 
-    stop_event.clear()
-    stop_event_experiment.clear()
     controller.log.info('Finish autocRED experiment')
 
 
diff --git a/src/instamatic/gui/ctrl_frame.py b/src/instamatic/gui/ctrl_frame.py
index 1d1c8e4d..21199c9c 100644
--- a/src/instamatic/gui/ctrl_frame.py
+++ b/src/instamatic/gui/ctrl_frame.py
@@ -50,6 +50,8 @@ def __init__(self, parent):
         frame.pack(side='top', fill='x', padx=10, pady=10)
 
         frame = Frame(self)
+        stage_reset_btn = Button(frame, text='Reset stage', command=self.reset_stage)
+        stage_reset_btn.grid(row=0, column=0, sticky='W')
 
         Label(frame, text='Angle (-)', width=20).grid(row=1, column=0, sticky='W')
         Label(frame, text='Angle (0)', width=20).grid(row=2, column=0, sticky='W')
@@ -161,6 +163,13 @@ def __init__(self, parent):
         )
         slider.grid(row=12, column=0, columnspan=3, sticky='EW')
 
+        # Magnification
+        Label(frame, text='Magnification', width=20).grid(row=14, column=0, sticky='W')
+        mag_inc_btn = Button(frame, text='+', command=self.increase_mag)
+        mag_inc_btn.grid(row=14, column=1)
+        mag_dec_btn = Button(frame, text='-', command=self.decrease_mag)
+        mag_dec_btn.grid(row=14, column=2)
+
         frame.pack(side='top', fill='x', padx=10, pady=10)
 
         frame = Frame(self)
@@ -228,6 +237,15 @@ def set_brightness(self, event=None):
     def get_brightness(self, event=None):
         self.var_brightness.set(self.ctrl.brightness.get())
 
+    def increase_mag(self):
+        self.ctrl.magnification.increase()
+
+    def decrease_mag(self):
+        self.ctrl.magnification.decrease()
+
+    def reset_stage(self):
+        self.ctrl.stage.set(0, 0, 0, 0, 0)
+
     def set_difffocus(self, event=None):
         self.var_difffocus.set(self.var_difffocus.get())
         self.q.put(('ctrl', {'task': 'difffocus.set', 'value': self.var_difffocus.get()}))
diff --git a/src/instamatic/processing/find_crystals.py b/src/instamatic/processing/find_crystals.py
index c58bbd84..4c8bda7d 100644
--- a/src/instamatic/processing/find_crystals.py
+++ b/src/instamatic/processing/find_crystals.py
@@ -82,7 +82,7 @@ def segment_crystals(img, r=101, offset=5, footprint=5, remove_carbon_lacing=Tru
     # remove carbon lines
     if remove_carbon_lacing:
         arr = morphology.remove_small_objects(arr, min_size=8 * 8, connectivity=0)
-        arr = morphology.remove_small_holes(arr, min_size=32 * 32, connectivity=0)
+        arr = morphology.remove_small_holes(arr, area_threshold=32 * 32, connectivity=0)
     arr = morphology.binary_dilation(arr, morphology.disk(footprint))  # dilation
 
     # get background pixels
@@ -227,7 +227,7 @@ def main_entry():
 
     for fn in args:
         img, h = read_image(fn)
-
+        print(h)
         crystals = find_crystals_timepix(img, h['exp_magnification'], plot=True)
 
         for crystal in crystals:
diff --git a/src/instamatic/simulation/__init__.py b/src/instamatic/simulation/__init__.py
index 75cc4c27..6138efbd 100644
--- a/src/instamatic/simulation/__init__.py
+++ b/src/instamatic/simulation/__init__.py
@@ -4,7 +4,7 @@
 
 from instamatic import config
 from instamatic.camera.camera_client import CamClient
-from instamatic.camera.videostream import VideoStream
+from instamatic.camera.fakevideostream import VideoStream
 from instamatic.microscope.client import MicroscopeClient
 from instamatic.simulation.camera import CameraSimulation
 from instamatic.simulation.microscope import MicroscopeSimulation
@@ -17,6 +17,9 @@ def get_simulation_instances() -> Tuple[CameraSimulation, MicroscopeSimulation]:
     -------
     Tuple[CameraSimulation, MicroscopeSimulation]
     """
+    tem = MicroscopeSimulation()
+    cam = CameraSimulation(tem=tem)
+    return cam, tem
     if config.settings.use_tem_server:
         tem = MicroscopeClient(interface='simulate')
     else:
diff --git a/src/instamatic/simulation/camera.py b/src/instamatic/simulation/camera.py
index b74e8d48..da96442c 100644
--- a/src/instamatic/simulation/camera.py
+++ b/src/instamatic/simulation/camera.py
@@ -2,7 +2,7 @@
 
 from typing import Tuple
 
-from numpy import ndarray
+from numpy import ndarray, random
 
 from instamatic.camera.camera_base import CameraBase
 from instamatic.microscope.base import MicroscopeBase
@@ -27,6 +27,12 @@ def establish_connection(self):
     def release_connection(self):
         pass
 
+    def block(self):
+        print('Blocking')
+
+    def unblock(self):
+        print('Unblocking')
+
     def get_image(self, exposure: float = None, binsize: int = None, **kwargs) -> ndarray:
         if exposure is None:
             exposure = self.default_exposure
@@ -76,9 +82,12 @@ def get_image(self, exposure: float = None, binsize: int = None, **kwargs) -> nd
                 shape=shape, x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max
             )
         else:
-            return self.stage.get_image(
+            img = self.stage.get_image(
                 shape=shape, x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max
             )
+            # Add some noise
+            img += random.random_integers(1, 10, shape).astype(img.dtype)
+            return img
 
     def _mag_to_ranges(self, mag: float) -> Tuple[float, float, float, float]:
         # assume 50x = 2mm full size
diff --git a/src/instamatic/simulation/crystal.py b/src/instamatic/simulation/crystal.py
index f98b3501..1abef961 100644
--- a/src/instamatic/simulation/crystal.py
+++ b/src/instamatic/simulation/crystal.py
@@ -58,10 +58,17 @@ def __init__(
 
         self.lattice = diffpy.Lattice(self.a, self.b, self.c, self.alpha, self.beta, self.gamma)
         self.structure = diffpy.Structure(
-            atoms=[diffpy.Atom('Au', xyz=[0, 0, 0])],
+            atoms=[diffpy.Atom('C', xyz=[0, 0, 0])],
             lattice=self.lattice,
         )
         self.phase = Phase(space_group=space_group, structure=self.structure)
+        self.recip = DiffractingVector.from_min_dspacing(
+            self.phase,
+            min_dspacing=1,
+            include_zero_vector=False,
+        )
+        self.recip.sanitise_phase()
+        self.recip.calculate_structure_factor()
 
     @property
     def a_vec(self) -> np.ndarray:
@@ -153,6 +160,7 @@ def diffraction_pattern_mask(
         rotation_matrix: np.ndarray,
         acceleration_voltage: float,
         excitation_error: float,
+        intensity_scale: float = 1,
     ) -> np.ndarray:
         """Get a diffraction pattern with a given shape, up to a given
         resolution, in a given orientation and acceleration voltage.
@@ -169,62 +177,68 @@ def diffraction_pattern_mask(
             acceleration_voltage of incident beam, in kV
         excitation_error : float
             Excitation error used for intensity calculation, in reciprocal Å
+        intensity_scale : float
+            Multiply final pattern by this value, defaults to 1.0
 
         Returns
         -------
         np.ndarray
             Diffraction pattern
         """
-        gen = SimulationGenerator(accelerating_voltage=acceleration_voltage)
+        gen = SimulationGenerator(
+            accelerating_voltage=acceleration_voltage, shape_factor_model='sin2c'
+        )
         wavelength = gen.wavelength
 
-        # Rotate using all the rotations in the list
-        recip = DiffractingVector.from_min_dspacing(
-            self.phase,
-            min_dspacing=d_min,
-            include_zero_vector=False,
-        )
+        max_excitation_error = excitation_error
+
         rotation = Rotation.from_matrix(rotation_matrix)
+        from diffsims.generators.simulation_generator import (
+            Vector3d,
+            get_intersection_with_ewalds_sphere,
+        )
+
+        optical_axis = rotation * Vector3d.zvector()
+
         # Calculate the reciprocal lattice vectors that intersect the Ewald sphere.
-        (
-            intersected_vectors,
-            hkl,
-            shape_factor,
-        ) = gen.get_intersecting_reflections(
-            recip,
-            rotation,
+        intersection, excitation_error = get_intersection_with_ewalds_sphere(
+            self.recip,
+            optical_axis,
             wavelength,
-            max_excitation_error=excitation_error,
-            with_direct_beam=False,
+            max_excitation_error,
+            gen.precession_angle,
         )
+        # Select intersected reflections
+        intersected_vectors = self.recip[intersection].rotate_with_basis(rotation)
+        excitation_error = excitation_error[intersection]
+        r_spot = intersected_vectors.norm
 
-        # Calculate diffracted intensities based on a kinematic model.
-        intensities = get_kinematical_intensities(
-            self.structure,
-            hkl,
-            intersected_vectors.gspacing,
-            prefactor=shape_factor,
-            scattering_params=gen.scattering_params,
-        )
+        # Calculate shape factor
+        shape_factor = gen.get_shape_factor(excitation_error, max_excitation_error, r_spot)
+        # Calculate intensity
+        f_hkls = self.recip.structure_factor[intersection]
+        intensities = (f_hkls * f_hkls.conjugate()).real * shape_factor
 
         # Threshold peaks included in simulation as factor of zero beam intensity.
         peak_mask = intensities > np.max(intensities) * gen.minimum_intensity
         intensities = intensities[peak_mask]
         intersected_vectors = intersected_vectors[peak_mask]
-        intersected_vectors.intensity = intensities
+        intersected_vectors.intensity = intensities * intensity_scale
 
-        # Create a simulation object
         sim = Simulation2D(
             phases=self.phase,
             coordinates=intersected_vectors,
             rotations=rotation,
-            simulation_generator=gen,
+            simulation_generator=self,
             reciprocal_radius=1 / d_min,
         )
-
         # Simulate diffraction pattern
         return sim.get_diffraction_pattern(
-            shape, sigma=1, calibration=1 / d_min / (shape[0] / 2)
+            shape,
+            sigma=1,
+            calibration=1 / d_min / (shape[0] / 2),
+            fast=False,
+            normalize=False,
         )
 
     def __str__(self) -> str:
diff --git a/src/instamatic/simulation/microscope.py b/src/instamatic/simulation/microscope.py
index b3471da5..ca1b4cc4 100644
--- a/src/instamatic/simulation/microscope.py
+++ b/src/instamatic/simulation/microscope.py
@@ -56,6 +56,8 @@ def __init__(self, name: str = 'simulate'):
 
         self.BeamShift_x = random.randint(MIN, MAX)
         self.BeamShift_y = random.randint(MIN, MAX)
+        self.BeamShift_x = 0
+        self.BeamShift_y = 0
 
         self.BeamTilt_x = random.randint(MIN, MAX)
         self.BeamTilt_y = random.randint(MIN, MAX)
diff --git a/src/instamatic/simulation/stage.py b/src/instamatic/simulation/stage.py
index 28f59b0d..03a7294c 100644
--- a/src/instamatic/simulation/stage.py
+++ b/src/instamatic/simulation/stage.py
@@ -16,7 +16,7 @@ def __init__(
         self,
         num_crystals: int = 100_000,
         min_crystal_size: float = 100,
-        max_crystal_size: float = 1000,
+        max_crystal_size: float = 3_000,
         random_seed: int = 100,
     ) -> None:
         """Handle many samples on a grid.
@@ -32,14 +32,9 @@ def __init__(
         random_seed : int, optional
             Seed for random number generation, by default 100
         """
-        # TODO make this settable
-        self.x = 0
-        self.y = 0
-        self.z = 0
-        self.alpha_tilt = 0
-        self.beta_tilt = 0
+
         self.in_plane_rotation = 10  # TODO change this with focus/magnification
-        self.rotation_matrix = np.eye(3)
+        self.set_position(0, 0, 0, 0, 0)
 
         # TODO parameters
         self.grid = Grid()
@@ -47,9 +42,11 @@ def __init__(
         self.rng = np.random.Generator(np.random.PCG64(random_seed))
 
         # TODO parameters
-        # TODO multiple phases
         # TODO amorphous phase
-        self.crystal = Crystal(*self.rng.uniform(5, 25, 3), *self.rng.uniform(80, 110, 3))
+        self.crystals = [
+            Crystal(20, 20, 20, 90, 90, 90, 221),
+            Crystal(*self.rng.uniform(5, 25, 3), *self.rng.uniform(80, 110, 3)),
+        ]
 
         self.samples = [
             Sample(
@@ -104,7 +101,7 @@ def set_position(
             degrees=True,
         ).as_matrix()
 
-    def image_extent_to_sample_coordinates(
+    def image_extent_to_stage_coordinates(
         self,
         shape: tuple[int, int],
         x_min: float,
@@ -151,10 +148,11 @@ def image_extent_to_sample_coordinates(
                 )
             ]
         )
-
         p = l0 + l[:, np.newaxis] * np.dot(-l0.T + p0, n) / np.dot(l, n)
 
-        x, y, z = self.rotation_matrix.T @ p
+        # Rotate around sample-holder-constant axes, not stage
+        p0[2] = 0
+        x, y, z = p0[:, np.newaxis] + self.rotation_matrix.T @ (p - p0[:, np.newaxis])
         x = x.reshape(shape)
         y = y.reshape(shape)
         return x, y
@@ -188,7 +186,7 @@ def get_image(
         np.ndarray
             Image
         """
-        x, y = self.image_extent_to_sample_coordinates(
+        x, y = self.image_extent_to_stage_coordinates(
             shape=shape, x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max
         )
         x_min = x.min()
@@ -198,14 +196,16 @@ def get_image(
 
         grid_mask = self.grid.array_from_coords(x, y)
 
-        sample_data = np.ones(shape, dtype=int) * 1000
+        sample_data = np.full(shape, fill_value=0xF000, dtype=np.uint32)
         for ind, sample in enumerate(self.samples):
             if not sample.range_might_contain_crystal(
                 x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max
             ):
                 continue
             # TODO better logic here
-            sample_data[sample.pixel_contains_crystal(x, y)] = 1000 * (1 - sample.thickness)
+            sample_data[sample.pixel_contains_crystal(x, y)] = np.round(
+                0xF000 * (1 - sample.thickness)
+            )
 
         sample_data[grid_mask] = 0
 
@@ -243,7 +243,7 @@ def get_diffraction_pattern(
         np.ndarray
             diffraction pattern
         """
-        x, y = self.image_extent_to_sample_coordinates(
+        x, y = self.image_extent_to_stage_coordinates(
             shape=shape, x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max
         )
         x_min = x.min()
@@ -270,18 +270,19 @@ def get_diffraction_pattern(
                 # Crystal is completely on the grid
                 continue
 
-            reflections += self.crystal.diffraction_pattern_mask(
+            reflections += self.crystals[sample.crystal_index].diffraction_pattern_mask(
                 shape,
                 d_min=d_min,
                 rotation_matrix=self.rotation_matrix @ sample.rotation_matrix,
                 acceleration_voltage=200,
                 excitation_error=0.01,
+                intensity_scale=0xFFFF,
             )
         # TODO diffraction shift
 
         # TODO noise
 
         # Convert to int array
-        reflections = (reflections * 0x8000).astype(int)
+        reflections = (reflections).astype(np.uint32)
 
         return reflections