Repository: ehpor/hcipy
Issue: 'gaussian_mode' is not defined

Looks like the mode basis 'gaussian_mode' is missing. Maybe it was removed? It is called by default in the SingleModeFiber Detector.

Changed files:
hcipy/optics/__init__.py, hcipy/optics/fiber.py, tests/test_optics.py

diff --git a/hcipy/optics/__init__.py b/hcipy/optics/__init__.py
index f5e27ba..57c9384 100644
--- a/hcipy/optics/__init__.py
+++ b/hcipy/optics/__init__.py
@@ -17,8 +17,10 @@ __all__ = [
 	'Detector',
 	'NoiselessDetector',
 	'NoisyDetector',
+	'make_gaussian_fiber_mode',
 	'StepIndexFiber',
 	'SingleModeFiber',
+	'SingleModeFiberInjection',
 	'SingleModeFiberArray',
 	'GaussianBeam',
 	'make_sellmeier_glass',
diff --git a/hcipy/optics/fiber.py b/hcipy/optics/fiber.py
index 915f143..87e22e4 100644
--- a/hcipy/optics/fiber.py
+++ b/hcipy/optics/fiber.py
@@ -3,7 +3,36 @@ from .detector import Detector
 from .optical_element import OpticalElement, AgnosticOpticalElement, make_agnostic_forward, make_agnostic_backward
 from .wavefront import Wavefront
 from ..mode_basis import ModeBasis, make_lp_modes
-from ..field import Field
+from ..field import Field, CartesianGrid, RegularCoords
+from ..dev import deprecated
+
+def make_gaussian_fiber_mode(mode_field_diameter):
+	'''The Gaussian approximation of a fiber mode.
+
+	Parameters
+	----------
+	mode_field_diameter : scalar
+		The mode field diameter of the fiber.
+
+	Returns
+	-------
+	Field generator
+		The Gaussian fiber mode.
+	'''
+	def func(grid):
+		if grid.is_('cartesian'):
+			r2 = grid.x**2 + grid.y**2
+		else:
+			r2 = grid.as_('polar').r**2
+
+		# Compute and normalize the Gaussian function.
+		mode_field_radius = 0.5 * mode_field_diameter
+		res = np.exp(-r2 / mode_field_radius**2)
+		res *= np.sqrt(2 / np.pi) / mode_field_radius
+
+		return Field(res, grid)
+
+	return func
 
 class StepIndexFiber(AgnosticOpticalElement):
 	'''A step-index fiber.
@@ -18,23 +47,26 @@ class StepIndexFiber(AgnosticOpticalElement):
 		The numerical aperture of the fiber.
 	fiber_length : scalar
 		The length of the optical fiber.
+	position : array_like
+		The 2D position of the fiber in the optical plane.
 	'''
 	def __init__(self, core_radius, NA, fiber_length, position=None):
+		super().__init__(False, True)
+
 		self._core_radius = core_radius
 		self._NA = NA
 		self.fiber_length = fiber_length
 
 		self._position = np.zeros((2,)) if position is None else np.array(position)
 
-		AgnosticOpticalElement.__init__(self, False, True)
-
 	def make_instance(self, instance_data, input_grid, output_grid, wavelength):
 		monochromatic_V = self.V(wavelength)
+
 		instance_data.NA = self.evaluate_parameter(self._NA, input_grid, output_grid, wavelength)
 		instance_data.fiber_modes, instance_data.beta = make_lp_modes(input_grid.shifted(-self._position), monochromatic_V, self.core_radius, return_betas=True)
 
 	def num_modes(self, wavelength):
-		'''The approximate amount of modes of the fiber.
+		'''The approximate number of modes of the fiber.
 		'''
 		V = self.V(wavelength)
 		return V**2 / 2
@@ -68,7 +100,7 @@ class StepIndexFiber(AgnosticOpticalElement):
 
 	@property
 	def position(self):
-		'''The position of the fiber.
+		'''The 2D position of the fiber.
 		'''
 		return self._position
 
@@ -129,7 +161,6 @@ class StepIndexFiber(AgnosticOpticalElement):
 		Wavefront
 			The wavefront that exits the fiber.
 		'''
-
 		M = instance_data.fiber_modes.transformation_matrix
 		mode_coefficients = np.einsum('...i, ij->...j', wavefront.electric_field * wavefront.grid.weights, M.conj())
 		output_electric_field = np.einsum('...i, ij->...j', mode_coefficients * np.exp(1j * instance_data.beta * self.fiber_length), M.T.conj())
@@ -158,33 +189,113 @@ class StepIndexFiber(AgnosticOpticalElement):
 
 		return Wavefront(output_electric_field, wavefront.wavelength)
 
+@deprecated('This class operated counter to HCIPy best practices and will be removed in HCIPy 0.7 or later. Please use a combination of the SingleModeFiberInjection and Detector classes instead.')
 class SingleModeFiber(Detector):
 	def __init__(self, input_grid, mode_field_diameter, mode=None):
 		self.input_grid = input_grid
 		self.mode_field_diameter = mode_field_diameter
 
 		if mode is None:
-			mode = gaussian_mode
+			mode = make_gaussian_fiber_mode(mode_field_diameter)
 
-		self.mode = mode(self.input_grid, mode_field_diameter)
-		self.mode /= np.sum(np.abs(self.mode)**2 * self.input_grid.weights)
+		try:
+			self.mode = mode(self.input_grid, mode_field_diameter)
+		except TypeError:
+			self.mode = mode(self.input_grid)
+
+		self.mode /= np.sqrt(np.sum(np.abs(self.mode)**2 * self.input_grid.weights))
 		self.intensity = 0
 
 	def integrate(self, wavefront, dt, weight=1):
-		self.intensity += weight * dt * (np.dot(wavefront.electric_field * wavefront.electric_field.grid.weights, self.mode))**2
+		self.intensity += weight * dt * np.abs(np.dot(wavefront.electric_field.conj() * wavefront.electric_field.grid.weights, self.mode))**2
 
 	def read_out(self):
 		intensity = self.intensity
 		self.intensity = 0
 		return intensity
 
-# This implementation assumes orthogonality of the different fibers.
-# Forward() should be changed if this is added in the future.
-# Also, the modes are independent of wavelength.
+class SingleModeFiberInjection(OpticalElement):
+	'''Injection into a single mode fiber.
+
+	Parameters
+	----------
+	input_grid : Grid
+		The grid of the incoming wavefront.
+	mode : Field generator
+		The mode of the fiber.
+	position : array_like
+		The 2D position of the fiber. If this is None, the fiber will be located at (0, 0).
+	'''
+	def __init__(self, input_grid, mode, position=None):
+		self.input_grid = input_grid
+
+		if position is None:
+			position = np.zeros(2)
+		else:
+			position = np.array(position)
+
+		self.output_grid = CartesianGrid(RegularCoords([1, 1], [1, 1], position))
+		self.output_grid.weights = 1
+
+		# Compute the fiber mode and normalize.
+		self.mode = mode(self.input_grid.shifted(-position))
+		self.mode.grid = self.input_grid
+		self.mode /= np.sqrt(np.sum(np.abs(self.mode)**2 * self.input_grid.weights))
+
+	def forward(self, wavefront):
+		'''Forward propagate the light through the fiber.
+
+		Parameters
+		----------
+		wavefront : Wavefront
+			The incoming wavefront.
+
+		Returns
+		-------
+		Wavefront
+			The complex amplitude of the single-mode fiber.
+		'''
+		output = np.dot(wavefront.electric_field.conj() * self.input_grid.weights, self.mode)
+
+		if not np.isscalar(output):
+			output = np.array(output)
+		else:
+			# The incoming wavefront was polarized, so we need to add a dimension to
+			# make it a valid Field.
+			output = output[..., np.newaxis]
+
+		output = Field(output, self.output_grid)
+
+		return Wavefront(output, wavefront.wavelength, wavefront.input_stokes_vector)
+
+	def backward(self, wavefront):
+		'''Backwards propagate the light through the fiber.
+
+		Parameters
+		----------
+		wavefront : Wavefront
+			The complex amplitude for the single-mode fiber.
+
+		Returns
+		-------
+		Wavefront
+			The outgoing wavefront.
+		'''
+		output = wavefront.electric_field * self.mode
+		output = Field(output, self.input_grid)
+
+		return Wavefront(output, wavefront.wavelength, wavefront.input_stokes_vector)
+
 class SingleModeFiberArray(OpticalElement):
 	def __init__(self, input_grid, fiber_grid, mode, *args, **kwargs):
 		'''An array of single-mode fibers.
 
+		This implementation assumes orthogonality of the different fibers. It
+		does not attempt to perform any orthogonalization of the fibers. This
+		may lead to creation of light power.
+
+		Additionally, all fiber modes are assumed to be independent of wavelength.
+
 		Parameters
 		----------
 		input_grid : Grid
@@ -205,15 +316,33 @@ class SingleModeFiberArray(OpticalElement):
 		self.projection_matrix = self.fiber_modes.transformation_matrix
 
 	def forward(self, wavefront):
-		res = self.projection_matrix.T.dot(wavefront.electric_field * self.input_grid.weights)
+		'''Forward propagate the light through the fiber array.
+
+		Parameters
+		----------
+		wavefront : Wavefront
+			The incoming wavefront.
+
+		Returns
+		-------
+		Wavefront
+			The complex amplitudes of each fiber.
+		'''
+		res = self.projection_matrix.T.dot(wavefront.electric_field.conj() * self.input_grid.weights)
 		return Wavefront(Field(res, self.fiber_grid), wavefront.wavelength)
 
 	def backward(self, wavefront):
-		res = self.projection_matrix.dot(wavefront.electric_field)
-		return Wavefront(Field(res, self.input_grid), wavefront.wavelength)
+		'''Backwards propagate the light through the fiber array.
 
-	def get_transformation_matrix_forward(self, wavelength=1):
-		return self.projection_matrix.T
+		Parameters
+		----------
+		wavefront : Wavefront
+			The complex amplitudes for each of the single-mode fibers.
 
-	def get_transformation_matrix_backward(self, wavelength=1):
-		return self.projection_matrix
+		Returns
+		-------
+		Wavefront
+			The outgoing wavefront.
+		'''
+		res = self.projection_matrix.dot(wavefront.electric_field)
+		return Wavefront(Field(res, self.input_grid), wavefront.wavelength)
diff --git a/tests/test_optics.py b/tests/test_optics.py
index a295a6b..f2747e9 100644
--- a/tests/test_optics.py
+++ b/tests/test_optics.py
@@ -636,6 +636,56 @@ def test_step_index_fiber():
 		assert forward_wf.total_power <= img.total_power * (1 + 1e-4)
 		assert backward_wf.total_power <= forward_wf.total_power * (1 + 1e-4)
 
+def test_gaussian_fiber_mode_power():
+	for q in [16, 32, 64]:
+		for mfd in [0.25, 0.5, 1]:
+			grid = make_focal_grid(q, 4)
+			mode = make_gaussian_fiber_mode(mfd)(grid)
+
+			wf = Wavefront(mode)
+
+			assert np.allclose(wf.total_power, 1)
+
+def test_gaussian_fiber_mode_mfd():
+	for mfd in [0.25, 0.5, 1]:
+		grid = CartesianGrid(UnstructuredCoords(([0, mfd / 2], [0, 0])))
+		mode = make_gaussian_fiber_mode(mfd)(grid)
+
+		assert np.allclose(mode[1] / mode[0], np.exp(-1))
+
+def test_single_mode_fiber_injection():
+	grid = make_focal_grid(16, 4)
+	mode_1 = make_gaussian_fiber_mode(1)
+
+	def mode_2(grid):
+		return mode_1(grid) * 1j
+
+	for mode in [mode_1, mode_2]:
+		smf = SingleModeFiberInjection(grid, mode)
+
+		wf = Wavefront(mode(grid))
+		wf.total_power = 1
+
+		post_fiber_wf = smf.forward(wf)
+
+		# This wavefront should fully couple into the fiber.
+		assert np.allclose(post_fiber_wf.electric_field, 1)
+		assert np.allclose(post_fiber_wf.total_power, 1)
+
+		pre_fiber_wf = smf.backward(post_fiber_wf)
+
+		# Power should be conserved when injecting a perfect mode.
+		assert np.allclose(pre_fiber_wf.electric_field, wf.electric_field)
+		assert np.allclose(pre_fiber_wf.total_power, 1)
+
+		wf = Wavefront(mode(grid) * grid.x)
+		wf.total_power = 1
+
+		post_fiber_wf = smf.forward(wf)
+
+		# This wavefront should not couple at all into the fiber.
+		assert np.allclose(post_fiber_wf.total_power, 0)
+
 def check_thin_lens(wf, lens, num_steps, focal_length):
 	dz = focal_length / (num_steps / 2)
 	prop = FresnelPropagator(wf.electric_field.grid, dz)