Adapted ideal_n_photon_source.py to the Fock backedn

quasi/_math/fock/ops.py

@@ -548,3 +548,27 @@ def mix(state, n):
     out_str = [indices[: 2 * n]]
     einstr = "".join(left_str + [","] + right_str + ["->"] + out_str)
     return np.einsum(einstr, state, state.conj())
+
+
+@njit
+def matrix_power(matrix, n):
+    """
+    Multiplies a square matrix by itself n times using Numba, ensuring both matrices have the same data type.
+
+    Parameters:
+    matrix (np.ndarray): A square NumPy array representing the matrix, can be float64 or complex128.
+    n (int): The number of times the matrix is to be multiplied by itself.
+
+    Returns:
+    np.ndarray: The resulting matrix after n multiplications, in complex128 if input is complex or float64 otherwise.
+    """
+    # Ensure both matrices are of the same complex data type if matrix contains complex numbers
+    if np.iscomplexobj(matrix):
+        matrix = matrix.astype(np.complex128)
+        result = np.eye(matrix.shape[0], dtype=np.complex128)  # Initialize result as a complex identity matrix
+    else:
+        result = np.eye(matrix.shape[0], dtype=matrix.dtype)  # Use the same dtype as the input matrix
+
+    for _ in range(n):
+        result = np.dot(result, matrix)
+    return result

quasi/devices/sources/ideal_n_photon_source.py

@@ -17,7 +17,7 @@
 from quasi.simulation import Simulation, SimulationType, ModeManager
 
 from quasi.experiment import Experiment
-from quasi._math.fock.ops import adagger
+from quasi._math.fock.ops import adagger, matrix_power
 
 
 
@@ -75,15 +75,19 @@ def simulate_fock(self):
         """
         Fock Simulation
         """
+        print(self.ports["photon_num"].signal.contents)
+        photon_num=4
         # Get the Experiment object reference
         exp = Experiment.get_instance()
         # Get the mode manager
         mm = ModeManager()
         mode = mm.create_new_mode()
 
+
         # Generate the creation operator
-        ad = adagger(exp.cutoff)
-        exp.add_operation(ad, [mm.get_mode_index(mode)])
+        operator = matrix_power(adagger(exp.cutoff), photon_num)
+
+        exp.add_operation(operator, [mm.get_mode_index(mode)])
         self.ports["output"].signal.set_contents(
             timestamp=0,
             mode_id=mode)