add f sweep

prsctrl/devices/lock_in/__init__.py

@@ -1,4 +1,4 @@
-from .base import Lock_In_Amp
+from .base import LockInAmp
 
 TYPENAME_DUMMY = "Dummy"
 TYPENAME_SR830 = "SR830"
@@ -16,7 +16,7 @@ def list_devices() -> dict[str,list[str]]:
         pass
     return devices
 
-def connect_device(type_name: str, device_name: str) -> Lock_In_Amp:
+def connect_device(type_name: str, device_name: str) -> LockInAmp:
     if type_name == TYPENAME_DUMMY:
         return DummyLockInAmp()
     elif type_name == TYPENAME_SR830:

prsctrl/devices/lock_in/base.py

@@ -1,7 +1,7 @@
 from abc import ABC, abstractmethod
 from typing import Callable
 
-class Lock_In_Amp(ABC):
+class LockInAmp(ABC):
     @abstractmethod
     def test_connection(self) -> None:
         """
@@ -19,16 +19,15 @@ class Lock_In_Amp(ABC):
         pass
     
     @abstractmethod
-    def read_value(self) -> tuple[float, float]:
+    def read_value(self, which:str) -> float:
         """
         Read a single value
 
-        Returns
+        :param which: X, Y, R, Theta
-        -------
-        [timestamp, voltage]
         """
         pass
 
+
     @abstractmethod
     def __str__(self):
         pass

prsctrl/devices/lock_in/impl/dummy.py

@@ -1,9 +1,9 @@
-from ..base import Lock_In_Amp
+from ..base import LockInAmp
 from typing import Callable
 from time import time as now
 import numpy as np
 
-class DummyLockInAmp(Lock_In_Amp):
+class DummyLockInAmp(LockInAmp):
     def __init__(self):
         super().__init__()
 
@@ -42,9 +42,8 @@ class DummyLockInAmp(Lock_In_Amp):
     def check_overloads(self) -> bool | str:
         return False
 
-    def read_value(self):
+    def read_value(self, which: str):
-        """Read the value of R"""
+        return -1.0
-        return float(self.query("OUTP? 3"))
 
     def reset(self):
         pass

prsctrl/devices/lock_in/impl/model7260.py

@@ -1,7 +1,7 @@
 import pyvisa
 # import pkg_resources
 
-from ..base import Lock_In_Amp
+from ..base import LockInAmp
 from prsctrl.util.visa import enumerate_devices
 
 import logging
@@ -9,7 +9,7 @@ log = logging.getLogger(__name__)
 
 import numpy as np
 
-class Model7260(Lock_In_Amp):
+class Model7260(LockInAmp):
     """
     Wrapper class for the Model 7260 DSP Lock-In controlled via pyvisa
     """

prsctrl/devices/lock_in/impl/sr830.py

@@ -4,13 +4,13 @@ import struct   # for converting bytes to float
 
 import numpy as np
 
-from ..base import Lock_In_Amp
+from ..base import LockInAmp
 from prsctrl.utility.visa import enumerate_devices
 
 import logging
 log = logging.getLogger(__name__)
 
-class SR830(Lock_In_Amp):
+class SR830(LockInAmp):
     """
     Wrapper class for the SR830 controlled via pyvisa
     """
@@ -152,10 +152,12 @@ class SR830(Lock_In_Amp):
             return "Output"
         return False
 
-    def measureTODO(self): pass
+    OUTP = ["X", "Y", "R", "theta"]
-    def read_value(self):
+    def read_value(self, which: str):
-        """Read the value of R"""
+        if which not in self.OUTP:
-        return float(self.query("OUTP? 3"))
+            raise ValueError(f"Invalid output: {which}. Must be one of {self.OUTP}")
+        outp = self.OUTP.index(which) + 1
+        return float(self.query(f"OUTP? {outp}"))
     
     def reset(self):
         self.instr.write("*RST")
@@ -230,14 +232,16 @@ class SR830(Lock_In_Amp):
         ofsl = int(self.query("OFSL?"))
         return SR830.OFSL[ofsl]
 
-    def get_wait_time_s(self):
+    def get_wait_time_s(self, time_const: float|None=None, filter_slope: int|None=None):
         """
         Get the wait time required to reach 99% of the final value.
+        :param time_const: If not passed, the value will be read from the device
+        :param filter_slop: If not passed, the value will be read from the device
         See Manual 3-21
         :return:
         """
-        time_const = self.get_time_constant_s()
+        if not time_const: time_const = self.get_time_constant_s()
-        filter_slope = self.get_filter_slope()
+        if not filter_slope: filter_slope = self.get_filter_slope()
         if filter_slope == 6: return 5 * time_const
         elif filter_slope == 12: return 7 * time_const
         elif filter_slope == 18: return 9 * time_const
@@ -327,22 +331,22 @@ class SR830(Lock_In_Amp):
         """
         Get the data from the buffer.
 
-        :return: np.ndarray
+        :return: np.ndarray | tuple[np.ndarray, np.ndarray]
-            Returns a numpy of shape (<1 if one channel, 2 if both channels>, n_points)
         """
         if self._buffer_length is None: raise RuntimeError(f"Buffer not set up, call buffer_setup() first.")
         self.run("PAUS")
         take_n_points = min(self.buffer_get_n_points(), self._buffer_length)  # there might be more points stored then was required
-        if CH1 and CH2:
+        # if CH1 and CH2:
-            data = np.empty((2, take_n_points), dtype=float)
+        #     data = (np.empty((2, take_n_points), dtype=float) for _ in range(2))
-        elif CH1 or CH2:
+        # elif CH1 or CH2:
-            data = np.empty((1, take_n_points), dtype=float)
+        #     data = np.empty((1, take_n_points), dtype=float)
-        else:
+        if (not CH1) and (not CH2):
             raise ValueError("Either CH1 or CH2 must be set True.")
+        data = []
         if CH1:
-            data[ 0, :] = self._buffer_get_data(1, 0, take_n_points)[:]
+            data.append(self._buffer_get_data(1, 0, take_n_points))
         if CH2:
-            data[-1, :] = self._buffer_get_data(2, 0, take_n_points)[:]
+            data.append(self._buffer_get_data(2, 0, take_n_points))
         return data
 
     def _buffer_get_data_slow(self, CH=1, start=0, n_points=None):
@@ -402,3 +406,20 @@ class SR830(Lock_In_Amp):
 
     def __str__(self):
         return "SR830"
+
+def set_measurement_params(lockin: SR830, p: dict={}, **kwargs):
+    params = p | kwargs
+    key_to_setter = {
+        "time_constant_s": lockin.set_time_constant_s,
+        "filter_slope": lockin.set_filter_slope,
+        "sync_filter": lockin.set_sync_filter,
+        "reserve": lockin.set_reserve,
+        "sensitivity_volt": lockin.set_sensitivity_volt,
+        "frequency_Hz": lockin.set_frequency_Hz,
+        "reference": lockin.set_reference,
+        "reference_trigger": lockin.set_reference_trigger,
+    }
+    for k, v in params.items():
+        if k not in key_to_setter.keys():
+            raise KeyError(f"Invalid parameter {k}")
+        key_to_setter[k](v)

prsctrl/devices/monochromator/impl/dummy.py

@@ -13,7 +13,7 @@ class DummyMonochromator(Monochromator):
         self.wavelength_nm = -1
 
     def set_wavelength_nm(self, wavelength_nm):
-        log.info("Dummy-Monochromator set to {wl} nm")
+        log.info(f"Dummy-Monochromator set to {wavelength_nm} nm")
         self.wavelength_nm = wavelength_nm
 
     def get_wavelength_nm(self):

prsctrl/measurement.py

@@ -4,133 +4,125 @@ Created on Fri Jan 24 15:18:31 2025
 
 @author: Matthias Quintern
 """
-from .measurement_device.base import VoltageMeasurementDevice
+from .devices.lock_in.base import LockInAmp
-from .led_control_device.base import LedControlDevice
+from .devices.shutter.base import Shutter
-from .led_script import LedScript
+from .devices.monochromator import Monochromator
-from .utility.prsdata import DataCollector
+from .utility.prsdata import PrsData
 
 import time
 import datetime
 from queue import Queue
+import numpy as np
 
 import logging
 log = logging.getLogger(__name__)
 
-def measure(
+def get_wavelengths_values(wl_range: tuple | list):
-    vm_dev: VoltageMeasurementDevice, 
+    """
-    led_dev: LedControlDevice, 
+    :param wl_range:
-    led_script: LedScript,
+        if tuple, return list(range(*wl_range))
-    data: DataCollector,
+        if list, return copy of wl_range
-    delta_t: float=0.1, 
+    """
-    flush_after:int|None=None,
+    if isinstance(wl_range, tuple):
-    use_buffer=False,
+        wavelengths = list(range(*wl_range))
-    max_measurements: int=None,
+    elif isinstance(wl_range, list) or isinstance(wl_range, np.ndarray):
-    stop_on_script_end: bool=False,
+        wavelengths = wl_range.copy()
-    verbose: bool=False,
+    else:
-    command_queue: None|Queue=None,
+        raise ValueError(f"Invalid type for 'wavelengths_nm': {type(wl_range)}")
-    data_queue: None|Queue=None,
+    return wavelengths
-    add_measurement_info_to_metadata=True
+
+
+def measure_spectrum(
+        monochromator: Monochromator,
+        lockin: LockInAmp,
+        shutter: Shutter,
+        data: PrsData,
+        measurement_params:dict,
+        aux_DC="Aux In 4",
+        command_queue: None | Queue = None,
+        data_queue: None | Queue = None,
+        add_measurement_info_to_metadata=True
 ):
-    """
-    Perform a measurement
 
-    Parameters
+    import pyvisa
-    ----------
+    def run_lockin_cmd(cmd, n_try=2):
-    vm_dev : VoltageMeasurementDevice
+        com_success = n_try
-        DESCRIPTION.
+        e = None
-    led_dev : LedControlDevice
+        while com_success > 0:
-        DESCRIPTION.
+            try:
-    led_script : LedScript
+                return cmd()
-        DESCRIPTION.
+            except pyvisa.VisaIOError as e:
-    data : DataCollector
+                # TODO: retry if status bit is set
-        DESCRIPTION.
+                lockin.try_recover_from_communication_error(e)
-    delta_t : float, optional
+            com_success -= 1
-        Target interval between measurements and led updates. The default is 0.1.
+        raise e
-    flush_after : int|None, optional
+
-        If int, flush values to disk after <flush_after>. The default is None.
+    default_measurement_params = {
-    use_buffer : TYPE, optional
+        "measurement_time_s": 30,
-        If True, use the buffer measurement mode. The default is False.
+        "sample_rate_Hz": 512,
-    max_measurements : int, optional
+        "wait_time_s": 1,
-        Number of measurements to perform before returning.
+        "wavelengths_nm": (390, 720, 1),
-        Note: If use_buffer=True, a few more than max_measurements might be performed
+    }
-        The default is None.
+    measurement_params = default_measurement_params | measurement_params
-    stop_on_script_end : bool, optional
+    wait_time_s = measurement_params["wait_time_s"]
-        Stop when the script end is reached.
+    sample_rate_Hz = measurement_params["sample_rate_Hz"]
-    verbose : bool, optional
+    measurement_time_s = measurement_params["measurement_time_s"]
-        If True, print some messages. The default is False.
+    n_bins = sample_rate_Hz * measurement_time_s
-    command_queue : None|Connection, optional
+    wavelengths = get_wavelengths_values(measurement_params["wavelengths_nm"])
-        A queue to receive to commands from. 
+    print(wavelengths)
-        Commands may be:
+
-            "stop" -> stops the measurement
+    timeout_s = 3 * measurement_time_s
-            ("led_script", <LedScript object>) a new led script to use
+    timeout_interval = 0.5
-        The default is None.
-    data_queue : None|Queue, optional
-        A queue to put data in. The default is None.
-    add_measurement_info_to_metadata : bool, optional
-        If True, add measurement info to the metadata:
-        time, measurement_interval, measurement_use_buffer, measurement_voltage_device, measurement_led_device
-        The default is True.
-    Returns
-    -------
-    None.
 
-    """
     get_time = lambda: datetime.datetime.now().strftime("%Y-%m-%d_%H:%M:%S")
     if add_measurement_info_to_metadata:
-        data.metadata["measurement_interval"] = str(delta_t) + " s"
+        data.metadata["device_lock-in"] = str(lockin)
-        data.metadata["measurement_use_buffer"] = str(use_buffer)
+        data.metadata["device_monochromator"] = str(monochromator)
-        data.metadata["measurement_voltage_measurement_device"] = str(vm_dev)
-        data.metadata["measurement_led_control_device"] = str(led_dev)
-        led_name = led_dev.get_led_name()
-        if led_name: data.metadata["measurement_led_lamp"] = led_name
         data.metadata["measurement_time_start"] = get_time()
     # write metadata to disk
     data.write_metadata()
-    vm_dev.reset(True)
+    print(wait_time_s)
-    if use_buffer:
+    data.metadata["messages"] = []
-        vm_dev.buffer_measure(delta_t, verbose=True)
+
-    # allow 0 instead of None
-    if max_measurements == 0: max_measurements = None
-    if flush_after == 0: flush_after = None
     try:
-        i = 0
+        shutter.open()
-        led_val = led_script.start()
+        for i_wl, wl in enumerate(wavelengths):
-        try:
+            log.info(f"Measuring at lambda={wl} nm")
-            led_dev.set_level(led_val)
+            run_lockin_cmd(lambda: lockin.buffer_setup(CH1="R", CH2=aux_DC, length=n_bins, sample_rate=sample_rate_Hz))
-        except Exception as e:
+
-            log.error(f"Error setting led to {led_val:03}%: {e}")
+            data_queue.put(("set_wavelength", wl))
-            raise e
+            monochromator.set_wavelength_nm(wl)
-        t_iter_start = time.time()
+            data_queue.put(("wait_stable", ))
-        while True:
+            # wait the wait time
-            # using while True and if, to be able to log the stop reason
+            time.sleep(wait_time_s)
-            if max_measurements is not None and i >= max_measurements:
+            overload = run_lockin_cmd(lambda: lockin.check_overloads())
-                log.info(f"Reached maximum number of measurements ({i}{max_measurements}), stopping measurement")
+            if overload:
-                break
+                msg = f"Overload of {overload} at {wl} nm"
-            # 1) read value(s)
+                log.warning(msg)
-            if use_buffer:
+                data.metadata["messages"].append(msg)
-                try:
+            theta = []
-                    values = vm_dev.buffer_read_new_values()
+            measure_phase = lambda: theta.append(run_lockin_cmd(lambda: lockin.read_value("theta")))
-                except ValueError as e:
+            data_queue.put(("measuring", ))
-                    # print(f"Error in buffer measurement {i}:", e)
+            measure_phase()
-                    values = []
+            run_lockin_cmd(lambda: lockin.buffer_start_fill())
-            else:
+            # check if its done
-                values = [vm_dev.read_value()]
+            t = timeout_s
-            # print(values)
+            while t > 0:
-            # 2) process value(s)
+                t -= timeout_interval
-            for (tval, vval) in values:
+                time.sleep(timeout_interval)
-                if i == 0:
+                measure_phase()
-                    t0 = tval
+                if run_lockin_cmd(lambda: lockin.buffer_is_done()):
-                tval -= t0
+                    break
-                current_data = (i, tval, vval, led_val)
+            if t < 0: raise RuntimeError("Timed out waiting for buffer measurement to finish")
-                data.add_data(*current_data)
+            # done
-                # 3) write data
+            dR_raw, R_raw = run_lockin_cmd(lambda: lockin.buffer_get_data(CH1=True, CH2=True))
-                if verbose: print(f"n = {i:6d}, t = {tval: .2f} s, U = {vval: .5f} V, LED = {led_val:03}%" + " "*10, end='\r')
+            data[wl] = {}
-                if flush_after is not None and (i+1) % flush_after == 0:
+            data[wl]["dR_raw"] = dR_raw * 2 * np.sqrt(2)  # convert RMS to Peak-Peak
-                    data.flush(verbose=verbose)
+            data[wl]["R_raw"] = R_raw
-                # if a queue was given, put the data
+            data[wl]["theta_raw"] = np.array(theta)
-                if data_queue is not None:
+            spec_data = data.get_spectrum_data([wl])
-                    data_queue.put(current_data)
+            data.write_partial_file(wl)
-                i += 1
+            data_queue.put(("data", spec_data))
 
             # if a pipe was given, check for messages
             if command_queue is not None and command_queue.qsize() > 0:
@@ -138,49 +130,56 @@ def measure(
                 if recv == "stop":
                     log.info(f"Received 'stop', stopping measurement")
                     break
-                elif type(recv) == tuple and recv[0] == "led_script":
-                    log.info(f"Received 'led_script', replacing script")
-                    led_script = recv[1]
                 elif type(recv) == tuple and recv[0] == "metadata":
                     log.info(f"Received 'metadata', updating metadata")
                     data.metadata |= recv[1]
                     data.write_metadata()
                 else:
                     log.error(f"Received invalid message: '{recv}'")
-
-            # 4) sleep
-            # subtract the execution time from the sleep time for a more
-            # accurate frequency
-            dt_sleep = delta_t - (time.time() - t_iter_start)
-            if dt_sleep > 0:
-                # print(f"Sleeping for {dt_sleep}")
-                time.sleep(dt_sleep)
-            t_iter_start = time.time()
-            # 5) update LED
-            if stop_on_script_end and led_script.is_done(t_iter_start):
-                log.info("Reached led script end, stopping measurement")
-                break
-            new_led_val = led_script.get_state(t_iter_start)
-            if new_led_val != led_val:
-                try:
-                    led_dev.set_level(new_led_val)
-                    led_val = new_led_val
-                except Exception as e:
-                    log.error(f"Error setting led to {new_led_val:03}%: {e}")
-                    raise e
-
     except KeyboardInterrupt:
         log.info("Keyboard interrupt, stopping measurement")
     except Exception as e:
         log.critical(f"Unexpected error, stopping measurement. Error: {e}")
         if command_queue is not None:
             command_queue.put(("exception", e))
+        raise e
+
     if add_measurement_info_to_metadata:
         data.metadata["measurement_time_stop"] = get_time()
     # Write again after having updated the stop time
     data.write_metadata()
-    data.flush()
+    data.write_full_file()
-    led_dev.off()
 
-    
+
-    
+def set_offsets_laser_only(
+        lockin: LockInAmp,
+        shutter: Shutter,
+        wait_time_s,
+        R=True,
+        phase=True,
+        data_queue: None | Queue = None,
+    ):
+    """
+    Set the R offset from the signal when only the laser is on.
+    This signal should be stray laser light and laser induced PL
+    :param phase: If True, use the Auto-Phase function to offset the phase
+    :param R: If True, use the Auto-Offset function to offset R
+    :return: Offset as percentage of the full scale R, Phase offset in degrees
+    """
+    log.info("Setting offset when the lamp is off.")
+    shutter.close()
+    if data_queue:
+        data_queue.put(("wait_stable", ))
+    time.sleep(wait_time_s + 10)
+    # TODO: generalize for other lock-ins
+    if data_queue:
+        data_queue.put(("set_offsets", ))
+    lockin.run("AOFF 3")  # auto offset R
+    # R must come before phase, because after auto-phase the signal needs to stabilize again
+    if R:
+        R_offset_fs = float(lockin.query("OEXP? 3").split(",")[0])  # returns R offset and expand
+    if phase:
+        lockin.run("APHS")
+    phase_offset_deg = float(lockin.query("PHAS? 3"))  # returns R offset and expand
+    if data_queue: data_queue.put(("offsets", R_offset_fs, phase_offset_deg))
+    return R_offset_fs, phase_offset_deg

prsctrl/prsctrl_interactive.py

@@ -2,10 +2,13 @@
 run this before using this library:
 ipython -i prctrl_interactive.py
 """
+from toolz import frequencies
+
 version = "0.1"
 
 import numpy as np
 import matplotlib.pyplot as plt
+import time
 
 from datetime import datetime as dtime
 from os import path, makedirs
@@ -30,13 +33,14 @@ from .devices import lock_in as mod_lock_in
 from .devices import lamp as mod_lamp
 from .devices import monochromator as mod_monochromator
 # import base classes
-from .devices.lock_in import Lock_In_Amp
+from .devices.lock_in import LockInAmp
 from .devices.shutter import Shutter
 from .devices.lamp import Lamp
 from .devices.monochromator import Monochromator
 
-# from .measurement import measure as _measure
+from .devices.lock_in.impl.sr830 import set_measurement_params
-from .utility.data_collector import PrsDataCollector
+from .measurement import measure_spectrum as _measure_spectrum, set_offsets_laser_only, get_wavelengths_values
+from .utility.prsdata import PrsData, plot_spectrum
 from .utility.config_file import ConfigFile
 from .utility.device_select import select_device_interactive, connect_device_from_config_or_interactive
 from .update_funcs import Monitor
@@ -64,12 +68,11 @@ test = False
 
 # DEVICES
 # global variable for the instrument/client returned by pyvisa/bleak
-lockin: Lock_In_Amp|None = None
+lockin: LockInAmp | None = None
 shutter: Shutter|None = None
 lamp: Lamp|None = None
 mcm: Monochromator|None = None
 
-data_collector = PrsDataCollector(data_path=settings["datadir"], data_name="interactive", dirname="interactive_test", add_number_if_dir_exists=True)
 t0 = 0
 data = None
 md = None
@@ -78,97 +81,179 @@ from .test_measurement import  _measure_both_sim
 def measure_both_sim(**kwargs):
     return _measure_both_sim(mcm, lockin, shutter, **kwargs)
 
-def monitor(script: str|int=0, interval: float|None=None, metadata:dict={}, flush_after: int|None=None, use_buffer: bool|None=None, max_measurements=None, stop_on_script_end: bool=False, max_points_shown=None):
+default_lockin_params = {
-    """
+    "time_constant_s": 10,
-    Monitor the voltage with matplotlib.
+    # "time_constant_s": 300e-3,
-    - Opens a matplotlib window and takes measurements depending on settings["interval"]
+    "sensitivity_volt": 500e-6,
-    - Waits for the user to press a key
+    "filter_slope": 12,
-    
+    "sync_filter": 1,
-    If use_buffer=False, uses python's time.sleep() for waiting the interval, which is not very precise. 
+    "reserve": "Normal",
-    With use_buffer=True, the timing of the voltage data readings will be very precise, however,
+    "reference": "Internal",
-    the led updates may deviate up to <interval>.
+    "reference_trigger": "Falling Edge",
-    
+    "frequency_Hz": 173,
-    The data is automatically saved to "<date>_<time>_<name>" in the data directory.
+}
-    
+default_measurement_params = {
-    Parameters
+    "measurement_time_s": 30,
-    ----------
+    "sample_rate_Hz": 512,
-    script : str|int
+    "wait_time_s": 0,
-        Path to a led script file, or a constant value between 0 and 100 for the LED.
+    "wavelengths_nm": (550, 555, 1),
-    interval : float|None
+}
-        Time between measurements. 
+
-        If None, the value is taken from the settings.
+
-    metadata : dict
+def get_time_estimate(lockin_params: dict, measurement_params: dict, offset_with_laser_only=True, extra_wait_time_s=10) -> float:
-        Metadata to append to the data header.
+    global lockin
-        The set interval as well as the setting for 'name' and 'led' are automatically added.
+
-    flush_after : int|None
+    lockin_params = default_lockin_params | lockin_params
-        Flush the data to disk after <flush_after> readings
+    measurement_params = default_measurement_params | {"wait_time_s":  lockin.get_wait_time_s(lockin_params["time_constant_s"], lockin_params["filter_slope"]) + extra_wait_time_s} | measurement_params
-        If None, the value is taken from the settings.
+    wls = get_wavelengths_values(measurement_params["wavelengths_nm"])
-    use_buffer : bool
+
-        If True, use the voltage measurement device's internal buffer for readings, which leads to more accurate timings.
+    wait_time_s = measurement_params["wait_time_s"]
-        If None, the value is taken from the settings.
+    measurement_time_s = measurement_params["measurement_time_s"]
-    max_points_shown : int|None
+
-        how many points should be shown at once. None means infinite
+    t = 0
-    max_measurements : int|None
+    if offset_with_laser_only:
-        maximum number of measurements. None means infinite
+        t += 2 * (wait_time_s + 10)   # 2 * (at beginning and end), 10 for auto functions
-    stop_on_script_end : bool, optional
+    t += len(wls) * (5 + wait_time_s + measurement_time_s + 10)     # 5 for setting wavelength, 10 for data transfers
-        Stop measurement when the script end is reached
+    return t
-    """
+
-    global _runtime_vars, data_collector, dev, led
+def measure_spectrum(metadata:dict={},
+        lockin_params={},
+        measurement_params={},
+        aux_DC="Aux In 4",
+        offset_with_laser_only=True,
+        extra_wait_time_s=10,
+        name:str|None=None,
+        dirname=None,
+        save_spectrum=True
+    ):
+    global _runtime_vars
     global data, md
+    global lockin, shutter, lamp, mcm
     _runtime_vars["last_measurement"] = dtime.now().isoformat()
-    if interval is None: interval = settings["interval"]
+
-    if flush_after is None: flush_after = settings["flush_after"]
+    lockin_params = default_lockin_params | lockin_params
-    if use_buffer is None: use_buffer = settings["use_buffer"]
+    set_measurement_params(lockin, lockin_params)
+
+    # must come after lockin settings, since its using get_wait_time_s
+    measurement_params = default_measurement_params | {"wait_time_s":  lockin.get_wait_time_s() + extra_wait_time_s} | measurement_params
+
+    # get the frequency, if not given
+    refkey = "frequency_Hz"
+    if not refkey in measurement_params:
+        measurement_params[refkey] = lockin.get_frequency_Hz()
+
+    # trigger on the falling edge, since the light comes through when the ref signal is low
+    # could of course also trigger on rising and apply 180° shift
+    lockin.run("RSLP 2")
+
     # set metadata
-    metadata["interval"] = str(interval)
+    metadata["lock-in_settings"] = lockin_params
-    metadata["name"] = settings["name"]
+    metadata["measurement_parameters"] = measurement_params
-    metadata["led"] = settings["led"]
+    if name is None:
-    metadata["led_script"] = str(script)
+        name = settings["name"]
-    print(f"Starting measurement with:\n\tinterval = {interval}s\nUse <C-c> to stop. Save the data using 'data.save_csv()' afterwards.")
+    metadata["name"] = name
+    print(f"Starting measurement with: Use <C-c> to stop. Save the data using 'data.save_csv()' afterwards.")
     plt.ion()
-    plt_monitor = Monitor(max_points_shown=max_points_shown)
+    plt_monitor = Monitor(r"$\lambda$ [nm]", [
-    led_script = LedScript(script=script, auto_update=True, verbose=True)
+        dict(ax=0, ylabel=r"$\Delta R$", color="green"),
-    data_collector = PrsDataCollector(metadata=metadata, data_path=settings["datadir"], data_name=settings["name"])
+        dict(ax=1, ylabel=r"$\sigma_{\Delta R}$", color="green"),
+        dict(ax=2, ylabel=r"$R$", color="blue"),
+        dict(ax=3, ylabel=r"$\sigma_R$", color="blue"),
+        dict(ax=4, ylabel=r"$\Delta R/R$", color="red"),
+        dict(ax=5, ylabel=r"$\sigma_{\Delta R/R}$", color="red"),
+        dict(ax=6, ylabel=r"$\theta$", color="aqua"),
+        dict(ax=7, ylabel=r"$\sigma_{\theta}$", color="aqua"),
+    ])
+    plt_monitor.set_fig_title(f"Turn on laser and plug detector into A and {aux_DC} ")
+    data = PrsData(data={}, metadata=metadata, write_data_path=settings["datadir"], write_data_name=settings["name"], write_dirname=dirname)
+
+    # measure/set offset
+    full_scale_voltage = lockin_params["sensitivity_volt"]
+    def set_offsets(name):
+        shutter.close()
+        plt_monitor.set_fig_title(f"Measuring baseline with lamp off")
+        R_offset_fs, phase_offset_deg = set_offsets_laser_only(lockin, shutter, measurement_params["wait_time_s"])
+        R_offset_volt = R_offset_fs * full_scale_voltage
+        data.metadata[f"R_offset_volt_{name}"] = R_offset_volt
+        data.metadata[f"phase_offset_deg_{name}"] = phase_offset_deg
+        print(f"R_offset_volt_{name} {R_offset_volt}")
+        print(f"phase_offset_deg_{name}: {phase_offset_deg}")
+    if offset_with_laser_only: set_offsets("before")
+
     # data_collector.clear()
     data_queue = mp.Queue()
     command_queue = mp.Queue()
     # Argument order must match the definition
-    proc_measure = mt.Thread(target=_measure, args=(dev,
+    proc_measure = mt.Thread(target=_measure_spectrum, args=(
-        led, 
+        mcm,
-        led_script, 
+        lockin,
-        data_collector,
+        shutter,
-        interval,
+        data,
-        flush_after, 
+        measurement_params,
-        use_buffer,
+        aux_DC,
-        max_measurements,
-        stop_on_script_end,
-        False,  # verbose
         command_queue,
-        data_queue
+        data_queue,
+        True,  # add metadata
     ))
     proc_measure.start()
     try:
         while proc_measure.is_alive():
             while not data_queue.empty():
                 # print(data_queue.qsize(), "\n\n")
-                current_data = data_queue.get(block=False)
+                msg = data_queue.get(block=False)
-                i, tval, vval, led_val = current_data
+                if msg[0] == "data":
-                plt_monitor.update(i, tval, vval, led_val)
+                    # data is as returned by PrsData.get_spectrum_data()
+                    plt_monitor.update(*msg[1][0,:])
+                elif msg[0] == "set_wavelength":
+                    plt_monitor.set_ax_title("Setting wavelength...")
+                    plt_monitor.set_fig_title(f"$\\lambda = {msg[1]}$ nm")
+                elif msg[0] == "wait_stable":
+                    plt_monitor.set_ax_title("Waiting until signal is stable...")
+                elif msg[0] == "measuring":
+                    plt_monitor.set_ax_title("Measuring...")
+                else:
+                    log.error(f"Invalid tuple received from measurement thread: {msg[0]}")
+            time.sleep(0.01)
 
     except KeyboardInterrupt:
         pass
     command_queue.put("stop")
     proc_measure.join()
+    plt_monitor.set_fig_title("")
+    plt_monitor.set_ax_title("")
+    try:
+        if offset_with_laser_only: set_offsets("before")
+        data["reference_freq_Hz_after"] = lockin.get_frequency_Hz()
+    except Exception as e:
+        print(e)
     print("Measurement stopped" + " "*50)
-    led_script.stop_updating()  # stop watching for file updates (if enabled)
+    if save_spectrum:
-    data_collector.save_csv(verbose=True)
+        plt.ioff()
-    data, metadata = data_collector.get_data()
+        fig = plot_spectrum(data, title=name, what=["dR_R", "theta"])
-    fig = data_plot(data, CPD=True, LED=True)
+        fig_path = path.join(data.path, data.dirname + ".pdf")
-    plt.ioff()
+        fig.savefig(fig_path)
-    fig_path = path.join(data_collector.path, data_collector.dirname + ".pdf")
+        return fig
-    fig.savefig(fig_path)
 
 
+def sweep_ref():
+    wavelenghts = [500, 550, 650, 660, 670, 680]
+    frequencies = list(range(27, 500, 5))
+    frequencies = [111, 444]
+    lockin_params = {
+        "time_constant_s": 10,
+    }
+    measurement_params = {
+        "wavelengths_nm": wavelenghts,
+    }
+    time_est = len(frequencies) * get_time_estimate(lockin_params=lockin_params, measurement_params=measurement_params, offset_with_laser_only=True, extra_wait_time_s=10)
+    print(f"Estimated time: {time_est}")
+    return
+    for f in frequencies:
+        dirname = f"2025-05-07_f-scan_f={f}_Hz"
+        lockin_params["frequency"] = f
+        measure_spectrum(lockin_params=lockin_params, measurement_params=measurement_params, dirname=dirname, name="Frequency scan $f = {f}$ Hz")
+        plt.close('all')
+
 # DATA
 def data_load(dirname:str) -> tuple[np.ndarray, dict]:
     """
@@ -185,7 +270,7 @@ def data_load(dirname:str) -> tuple[np.ndarray, dict]:
         dirpath = dirname
     else:
         dirpath = path.join(settings["datadir"], dirname)
-    data, md = PrsDataCollector.load_data_from_dir(dirpath, verbose=True)
+    data, md = PrsData.load_data_from_dir(dirpath, verbose=True)
     
 # SETTINGS
 def set(setting, value):
@@ -305,7 +390,7 @@ Enter 'help()' for a list of commands""")
     log_path = path.expanduser(config_file.get_or("path_log", "~/.cache/prsctrl-interactive.log"))
     makedirs(path.dirname(log_path), exist_ok=True)
     logging.basicConfig(
-        level=logging.WARN,
+        level=logging.INFO,
         format="%(asctime)s [%(levelname)s] [%(name)s] %(message)s",
         handlers=[
             logging.FileHandler(log_path),

prsctrl/test_measurement.py

@@ -11,7 +11,7 @@ from prsctrl.devices.lamp import Lamp
 from prsctrl.devices.shutter import Shutter
 from prsctrl.devices.monochromator import Monochromator
 from .update_funcs import Monitor
-from prsctrl.devices.lock_in import Lock_In_Amp
+from prsctrl.devices.lock_in import LockInAmp
 from prsctrl.devices.lock_in.impl.sr830 import SR830
 
 def set_measurement_params(lockin: SR830, p: dict={}, **kwargs):
@@ -69,6 +69,7 @@ def _measure_both_sim(monochromator: Monochromator, lockin: SR830, shutter: Shut
     measurement_params = {
         "measurement_time_s": 30,
         "sample_rate_Hz": 512,
+        "extra_wait_time_s": 10,
     }
     if laser_power_mW:
         measurement_params["laser_power_mW"] = laser_power_mW
@@ -98,9 +99,10 @@ def _measure_both_sim(monochromator: Monochromator, lockin: SR830, shutter: Shut
             com_success -= 1
         raise e
 
+    extra_wait_time_s = measurement_params["extra_wait_time_s"]
     # 5s for setting buffer,
     # 5s for get values and plot
-    print(f"Time estimate {(measurement_time_s + wait_time_s + 10 + 5 + 5)/60 * ((wl_range[1]-wl_range[0])/wl_range[2])} minutes")
+    print(f"Time estimate {(measurement_time_s + wait_time_s + extra_wait_time_s + 5 + 5)/60 * ((wl_range[1]-wl_range[0])/wl_range[2])} minutes")
     input("Make sure the laser is turned on and press enter > ")
     mon = monitor if monitor is not None else Monitor(r"$\lambda$ [nm]", [
         dict(ax=0, ylabel=r"$\Delta R$", color="green"),
@@ -139,7 +141,7 @@ def _measure_both_sim(monochromator: Monochromator, lockin: SR830, shutter: Shut
         monochromator.set_wavelength_nm(wl)
         mon.set_fig_title(f"Waiting for signal to stabilize")
         # wait the wait time
-        sleep(wait_time_s + 10)
+        sleep(wait_time_s + extra_wait_time_s)
         overload = run_lockin_cmd(lambda: lockin.check_overloads())
         if overload:
             msg = f"Overload of {overload} at {wl} nm"

prsctrl/tests/sweep_frequency.py

@@ -0,0 +1,65 @@
+from prsctrl.utility.prsdata import PrsData
+import os
+import re
+import numpy as np
+import matplotlib.pyplot as plt
+
+def process_results(data_dir, dir_regex=r"202.-..-.._f-scan_f=(\d+)_Hz", out_dir=None):
+    data_dir = os.path.expanduser(data_dir)
+    if out_dir is None:
+        out_dir = data_dir
+    paths = os.listdir(data_dir)
+    data_dirs = []
+    for p in paths:
+        full_path = os.path.join(data_dir, p)
+        if not os.path.isdir(full_path): continue
+        m = re.fullmatch(dir_regex, p)
+        if m:
+            data_dirs.append(full_path)
+        else:
+            print(f"Unmatched directory {p}")
+    assert len(data_dirs) > 0
+    data_dirs.sort()
+
+    frequencies = []
+    data = {}
+    shape = None
+    wls = None
+    for d in data_dirs:
+        print(f"Getting data from {d}")
+        pd = PrsData(load_data_path=d)
+        f = pd.metadata["lock-in_settings"]["frequency_Hz"]
+        # print(d, f)
+        sdata = pd.get_spectrum_data()
+        print(pd.wavelengths)
+        print(pd.data.keys())
+        if wls is None: wls = sdata[:,0]
+        if shape is None: shape = sdata.shape
+        else:
+            if shape != sdata.shape:
+                print(f"ERROR Shape mismatch for f={f}: {shape} != {sdata.shape}")
+                continue
+                # raise ValueError(f"Shape mismatch for {d}: {shape} != {sdata.shape}")
+        frequencies.append(f)
+        data[f] = sdata
+    data_per_wl_and_f = np.empty((shape[0], len(frequencies), shape[1]))
+    frequencies.sort()
+    for i in range(shape[0]):
+        for j, f in enumerate(frequencies):
+            data_per_wl_and_f[i, j, :] = data[f][i,:]
+    print(f"Found wavelengths: {wls}")
+    n_cols = 2
+    for qty in ["theta", "stheta", "dR_R", "sdR_R"]:
+        fig, axs = plt.subplots(wls.shape[0]//n_cols, n_cols, sharex=True, figsize=(8, 8))
+        axs = axs.flatten()
+        qty_idx = PrsData.default_spectrum_columns.index(qty)
+        fig.suptitle(f"Frequency scan: {PrsData.key_names[qty]}")
+        axs[-1].set_xlabel("Modulation Frequency $f$ [Hz]")
+        for i, wl in enumerate(wls):
+            ax = axs[i]
+            ax.set_ylabel(PrsData.labels[qty])
+            ax.plot(frequencies, data_per_wl_and_f[i, :, qty_idx])
+            ax.set_title(f"$\\lambda = {wl}$ nm")
+        fig.tight_layout()
+        fig.savefig(out_dir + f"result_{qty}.pdf")
+    print(frequencies)

prsctrl/utility/data_collector.py

@@ -9,7 +9,7 @@ from abc import abstractmethod
 log = logging.getLogger(__name__)
 
 from ..utility.file_io import get_next_filename, sanitize_filename
-from ..utility.prsdata import PrsData, FLUSH_TYPE, FLUSH_PREFIX, METADATA_FILENAME
+from ..utility.prsdata import PrsData, FLUSH_TYPE, PARTIAL_PREFIX, METADATA_FILENAME
 
 """
 Wollen:
@@ -76,7 +76,7 @@ class DataCollector:
            The full data and the metadata
         """
         if self.full_data is None:
-            self.full_data = PrsData(path=self.dirpath, metadata=self.metadata)
+            self.full_data = PrsData(load_data_path=self.dirpath, metadata=self.metadata)
         return self.full_data
 
     def save_csv_in_dir(self, sep=",", verbose=False):

prsctrl/utility/prsdata.py

@@ -3,59 +3,191 @@ import numpy as np
 import os
 import matplotlib.pyplot as plt
 import pickle
+import datetime
 import logging
 log = logging.getLogger(__name__)
 
 from ..utility.file_io import get_next_filename, sanitize_filename
 
 FLUSH_TYPE = "pickle-ndarray"
-FLUSH_PREFIX = "PART_"
+PARTIAL_PREFIX = "PART_"
-METADATA_FILENAME = FLUSH_PREFIX + "measurement_metadata.pkl"
+METADATA_FILENAME = PARTIAL_PREFIX + "measurement_metadata.pkl"
 
 
 class PrsData:
     """
     Class managing data and metadata.
-    Can be initialized from data directly, or a file or directory path.
+    Can be initialized from data directly, or from a file or directory path.
+
+    Keys:
+        - dR: delta R, the change in reflectivity. (This is the signal amplitude "R" from the lock-in)
+        - R: the baseline reflectivity. (DC signal measured using Aux In of the lock-in)
+        - theta: phase
+        - <qty>_raw: The raw measurement data (all individual samples)
+
+    data is a dictionary with:
+        - key: wavelength as int
+        - value: dictionary with:
+            - key: quantity as string: ["dR_raw", "R_raw", "theta_raw", "dR", ...]
+            - value: quantity value, or array if "_raw"
     """
-    def __init__(self, path:str|None=None, data:tuple|None=None, metadata:dict|None=None, verbose=False):
+    def __init__(self, data: dict | None=None, metadata: dict | None=None,
-        self.data = data
+                 load_data_path: str|None=None,
+                 write_data_path: str|None=None,
+                 write_data_name: str = "PRS",
+                 write_dirname: str | None = None,
+                 add_number_if_dir_exists=True,
+                 ):
         if type(metadata) == dict:
             self.metadata = metadata
         else:
             self.metadata = {}
-        if data is None and path is None:
+        self.data = data
+        if data is None and load_data_path is None:
             raise ValueError("Either path or data must be defined.")
-        if data is not None and path is not None:
+        if data is not None and load_data_path is not None:
             raise ValueError("Either path or data must be defined, but not both.")
-        if path is not None:  # load from file
+        if load_data_path is not None:  # load from file
-            if os.path.isdir(path):
+            if os.path.isdir(load_data_path):
-                self.data, md = PrsData.load_data_from_dir(path, verbose=verbose)
+                self.data, md = PrsData.load_data_from_dir(load_data_path)
                 self.metadata |= md
-            elif os.path.isfile(path):
+            elif os.path.isfile(load_data_path):
-                if path.endswith(".csv"):
+                if load_data_path.endswith(".csv"):
-                    self.data, md = PrsData.load_data_from_csv(path)
+                    self.data, md = PrsData.load_data_from_csv(load_data_path)
                     self.metadata |= md
-                elif path.endswith(".pkl"):
+                elif load_data_path.endswith(".pkl"):
-                    self.data, md = PrsData.load_data_from_pkl(path)
+                    self.data, md = PrsData.load_data_from_pkl(load_data_path)
                     self.metadata |= md
                 else:
                     raise NotImplementedError(f"Only .csv and .pkl files are supported")
             else:
-                raise FileNotFoundError(f"Path '{path}' is neither a file nor a directory.")
+                raise FileNotFoundError(f"Path '{load_data_path}' is neither a file nor a directory.")
-        else:
+        self.wavelengths = []
-            self.data = data
+        keys = list(self.data.keys())
+        for key in keys:
+            # for some reason, the wavelengths can end up as string
+            try:
+                wl = int(key)
+                self.wavelengths.append(wl)
+                self.data[wl] = self.data[key]
+                del self.data[key]
+            except ValueError:
+                pass
+        self.wavelengths.sort()
 
-    # Convert data
+        # INIT WRITE MODE
-    def to_dataframe(self):
+        self.dirname = None
-        df = pd.DataFrame(self.data, columns=self.columns)
+        self.path = None
-        df.meta = str(self.metadata)
+        self.name = write_data_name
-        return df
+        if write_data_path:
+            self.path = os.path.abspath(os.path.expanduser(write_data_path))
+            if write_dirname is None:
+                self.dirname = sanitize_filename(datetime.datetime.now().strftime("%Y-%m-%d_%H-%M") + "_" + self.name)
+            else:
+                self.dirname = sanitize_filename(write_dirname)
+            self.dirpath = os.path.join(self.path, self.dirname)
+
+            if os.path.exists(self.dirpath):
+                if not add_number_if_dir_exists:
+                    raise Exception(f"Directory '{self.dirname}' already exists. Provide a different directory or pass `add_number_if_dir_exists=True` to ignore this")
+                else:
+                    i = 1
+                    dirpath = f"{self.dirpath}-{i}"
+                    while os.path.exists(dirpath):
+                        i += 1
+                        dirpath = f"{self.dirpath}-{i}"
+                        print(f"Directory '{self.dirname}' already exists. Trying '{dirpath}' instead")
+                    self.dirpath = dirpath
+            self._assert_directory_exists()
+
+    def __setitem__(self, key, value):
+        self.data[key] = value
+        try:
+            wl = int(key)
+            self.wavelengths.append(wl)
+            self.wavelengths.sort()
+        except ValueError:
+            pass
+
+    def __getitem__(self, key):
+        return self.data[key]
+
+
+    def _check_has_wavelength(self, wl):
+        if wl not in self.data: raise KeyError(f"No data for wavelength '{wl}'")
+
+    def get_for_wl(self, wl, key) -> float:
+        self._check_has_wavelength(wl)
+        if key in self.data[wl]:
+            return self.data[wl][key]
+        elif key == "dR_R":
+            return self.calculate_dR_R_for_wl(wl)[0]
+        elif key == "sdR_R":
+            return self.calculate_dR_R_for_wl(wl)[1]
+        elif f"{key}_raw" in self.data[wl]:
+            vals = self.data[wl][f"{key}_raw"]
+            mean = np.mean(vals)
+            self.data[wl][key] = mean
+            return mean
+        elif key.startswith("s") and f"{key[1:]}_raw" in self.data[wl]:
+            vals = self.data[wl][f"{key[1:]}_raw"]
+            err = np.std(vals)
+            self.data[wl][key] = err
+            return err
+        raise KeyError(f"No '{key}' data for wavelength '{wl}'")
+
+    def calculate_dR_R_for_wl(self, wl) -> tuple[float, float]:
+        dR, sdR = self.get_for_wl(wl, "dR"), self.get_for_wl(wl, "sdR")
+        R, sR = self.get_for_wl(wl, "R"), self.get_for_wl(wl, "sR")
+        dR_R = dR / R
+        sdR_R = np.sqrt((sdR / R)**2 + (dR * sR/R**2)**2)
+        self.data[wl]["dR_R"] = dR_R
+        self.data[wl]["sdR_R"] = sdR_R
+        return dR_R, sdR_R
+
+    key_names = {
+        "wl": "Wavelength [nm]",
+        "dR": "dR [V]",
+        "R": "R [V]",
+        "sR": "sigma(R) [V]",
+        "sdR": "sigma(dR) [V]",
+        "dR_R": "dR/R",
+        "sdR_R": "sigma(dR/R)",
+        "theta": "theta [°]",
+        "stheta": "sigma(theta) [°]"
+    }
+    default_spectrum_columns=["wl", "dR", "sdR", "R", "sR", "dR_R", "sdR_R", "theta", "stheta"]
+    labels = {
+        "dR_R": r"$\Delta R/R$",
+        "dR": r"$\Delta R$ [V]",
+        "R": r"$R$ [V]",
+        "theta": r"$\theta$ [°]",
+    }
+    labels |= {f"s{k}": f"$\sigma_{{{v[1:]}}}" for k, v in labels.items()}
+    def get_spectrum_data(self, wavelengths=None, keys=None) -> np.ndarray:
+        """
+        Return the spectral data for the specified keys and wavelengths.
+        :param wavelengths: List of wavelengths, or None to use all wavelengths.
+        :param keys: List of keys, or None to use dR, R, dR_R, Theta
+        :return: numpy.ndarray where the first index is (wavelength=0, <keys>...=1...) and the second is the wavelengths.
+        """
+        if keys is None: keys = self.default_spectrum_columns
+        if wavelengths is None:
+            wavelengths = self.wavelengths
+        data = np.empty((len(wavelengths), len(keys)), dtype=float)
+        # this might be slow but it shouldnt be called often
+        for j, wl in enumerate(wavelengths):
+            for i, key in enumerate(keys):
+                if key == "wl":
+                    data[j][i] = wl
+                else:
+                    val = self.get_for_wl(wl, key)
+                    data[j][i] = val
+        return data
 
     def to_csv(self, sep=","):
         # self.to_dataframe().to_csv(os.path.join(self.path, self.name + ".csv"), index=False, metadata=True)
-        return PrsData.get_csv(self.data, self.metadata, sep=sep)
+        return PrsData.get_csv(self.get_spectrum_data(), self.metadata, sep=sep)
-
 
     def save_csv_at(self, filepath, sep=",", verbose=False):
         if verbose: print(f"Writing csv to {filepath}")
@@ -68,20 +200,77 @@ class PrsData:
         filepath = os.path.join(self.path, self.dirname + ".csv")
         self.save_csv_at(filepath, sep, verbose)
 
+    # FILE IO
+    def _check_write_mode(self):
+        if self.dirpath is None:
+            raise RuntimeError(f"Can not write data because {__class__.__name__} is not in write mode.")
+
+    def _assert_directory_exists(self):
+        if not os.path.isdir(self.dirpath):
+            os.makedirs(self.dirpath)
+
+
+    def write_partial_file(self, key):
+        self._check_write_mode()
+        if key not in self.data:
+            raise KeyError(f"Invalid key '{key}'")
+        filename = sanitize_filename(PARTIAL_PREFIX + str(key)) + ".pkl"
+        self._assert_directory_exists()
+        filepath = os.path.join(self.dirpath, filename)
+        log.info(f"Writing data '{key}' to {filepath}")
+        with open(filepath, "wb") as file:
+            pickle.dump(self.data[key], file)
+
+    def write_full_file(self):
+        filename = sanitize_filename("full-data") + ".pkl"
+        self._assert_directory_exists()
+        filepath = os.path.join(self.dirpath, filename)
+        log.info(f"Writing data to {filepath}")
+        with open(filepath, "wb") as file:
+            pickle.dump((self.data, self.metadata), file)
+
+    def write_metadata(self):
+        f"""
+        Write the metadata to the disk as '{METADATA_FILENAME}'
+        """
+        filepath = os.path.join(self.dirpath, METADATA_FILENAME)
+        log.debug(f"Writing metadata to {filepath}")
+        with open(filepath, "wb") as file:
+            pickle.dump(self.metadata, file)
+
     # STATIC CONVERTER
     @staticmethod
-    def get_csv(data, metadata, sep=","):
+    def get_csv(data: np.ndarray, metadata: dict, columns: list, sep=","):
         csv = ""
-        for k, v in metadata.items():
+        # metadata
-            csv += f"# {k}: {v}\n"
+        md_keys = list(metadata.keys())
-        csv += "".join(f"{colname}{sep}" for colname in PrsData.columns).strip(sep) + "\n"
+        md_keys.sort()
+        for k in md_keys:
+            v = metadata[k]
+            if type(v) == dict:
+                csv += f"# {k}:\n"
+                for kk, vv in v.items():
+                    csv += f"#     {kk}:{vv}\n"
+            if type(v) == list:
+                csv += f"# {k}:\n"
+                for vv in v:
+                    csv += f"#     - {vv}\n"
+            else:
+                csv += f"# {k}: {v}\n"
+        # data header
+        csv += "".join(f"{colname}{sep}" for colname in columns).strip(sep) + "\n"
+        # data
         for i in range(data.shape[0]):
-            csv += f"{i}{sep}{data[i,1]}{sep}{data[i,2]}{sep}{data[i,3]}\n"
+            csv += f"{data[i, 0]}"
+            for j in range(1, data.shape[1]):
+                csv += f"{sep}{data[i,j]}"
+            csv += "\n"
         return csv.strip("\n")
 
     # STATIC LOADERS
+    # TODO
     @staticmethod
-    def load_data_from_csv(filepath:str, sep: str=",") -> tuple[np.ndarray, dict]:
+    def load_data_from_csv(filepath:str, sep: str=",") -> tuple[dict, dict]:
         """
         Loads data from a single csv file.
         Lines with this format are interpreted as metadata:
@@ -119,25 +308,18 @@ class PrsData:
         return data, metadata
 
     @classmethod
-    def load_data_from_pkl(cls, filepath:str) -> tuple[np.ndarray, dict]:
+    def load_data_from_pkl(cls, filepath:str) -> tuple[dict, dict]:
         """
-        Loads data from a single csv file.
+        Loads data from a single pkl file.
-        Lines with this format are interpreted as metadata:
-            # key: value
-        Lines with this format are interpreted as data:
-            index, timestamp [s], CPD [V], LED [%]
         Parameters
         ----------
-        filepath
+        :param filepath Path to the file.
-           Path to the csv file.
+        :return
-        Returns
+            data
-        -------
+                2D numpy array with shape (n, 4) where n is the number of data points.
-        data
+            metadata
-            2D numpy array with shape (n, 4) where n is the number of data points.
+                Dictionary with metadata.
-        metadata
-            Dictionary with metadata.
         """
-        data = None
         metadata = {}
         with open(filepath, "rb") as f:
             obj = pickle.load(f)
@@ -146,72 +328,63 @@ class PrsData:
                 raise ValueError(f"Pickle file is a tuple with length {len(obj)}, however it must be 2: (data, metadata)")
             data = obj[0]
             metadata = obj[1]
-            if not isinstance(data, np.ndarray):
+            if not isinstance(data, dict):
-                raise ValueError(f"First object in tuple is not a numpy.ndarray")
+                raise ValueError(f"First object in tuple is not a dictionary but {type(data)}")
-        elif isinstance(obj, np.ndarray):
+        elif isinstance(obj, dict):
             data = obj
         else:
-            raise ValueError(f"Pickled object must be either numpy.ndarray or (numpy.ndarray, dict), but is of type {type(obj)}")
+            raise ValueError(f"Pickled object must be either dict=data or (dict=data, dict=metadata), but is of type {type(obj)}")
         # must be loaded by now
-        if not len(data.shape) == 2 and data.shape[1] == len(cls.columns):
-            raise ValueError(f"numpy.ndarray has invalid shape: {data.shape}, however the shape must be (N, {len(cls.columns)})")
         if not isinstance(metadata, dict):
             raise ValueError(f"Metadata is not a of type dict")
         return data, metadata
 
     @staticmethod
-    def load_data_from_dir(dirpath:str, verbose:bool=False) -> tuple[np.ndarray, dict]:
+    def load_data_from_dir(dirpath:str) -> tuple[dict, dict]:
         """
-        Combines all data files with the FLUSH_PREFIX from a directory into a numpy array
+        Combines all data files with the PARTIAL_PREFIX from a directory into a dictionary
 
-        Parameters
+        :param dirpath Path to the data directory
-        ----------
+        :return data, metadata
-        dirpath : str
-            Path to the data directory
-        verbose : bool, optional
-            If True, print a message for every file that is opened. The default is False.
-
-        Raises
-        ------
-        NotImplementedError
-            DESCRIPTION.
-
-        Returns
-        -------
-        data : ndarray
-            First index: Measurement
-            Second index: (index, timestamp [s], CPD [V], LED [%])
         """
         files = os.listdir(dirpath)
         files.sort()
-        data = np.empty((0, 4))
+        data = {}
         metadata = {}
         for filename in files:
             filepath = os.path.join(dirpath, filename)
-            if filename.startswith(FLUSH_PREFIX):
+            if filename.startswith(PARTIAL_PREFIX):
-                if filename.endswith(".csv"):
+                log.debug(f"Loading {filename}")
-                    if verbose: print(f"Opening {filepath} as csv")
+                # must be first
-                    df = pd.read_csv(filepath)
+                if filename == METADATA_FILENAME:  # Metadata filename must also start with FLUSH_PREFIX
-                    arr = df.to_numpy()
-                    data = np.concatenate((data, arr))
-                elif filename.endswith(".ndarray.pkl"):
-                    with open(filepath, "rb") as file:
-                        arr = pickle.load(file)
-                        if len(arr.shape) != 2 or arr.shape[1] != 4:
-                            print(f"Skipping file '{filepath}' with invalid array shape: {arr.shape}")
-                            continue
-                        data = np.concatenate((data, arr))
-                elif filename == METADATA_FILENAME:  # Metadata filename must also start with FLUSH_PREFIX
                     with open(filepath, "rb") as file:
                         metadata = pickle.load(file)
+                elif filename.endswith(".csv"):
+                    raise NotADirectoryError(f"Partial .csv files are not supported: '{filename}'")
+                elif filename.endswith(".pkl"):
+                    key = filename.strip(PARTIAL_PREFIX).strip(".pkl")
+                    with open(filepath, "rb") as file:
+                        val = pickle.load(file)
+                    data[key] = val
                 else:
                     raise NotImplementedError(f"Unknown file extension for file '{filepath}'")
             else:
                 log.info(f"Skipping unknown file: '{filepath}'")
         return data, metadata
 
+    def plot_raw_for_wl(self, wl, what=["dR", "R", "theta"]):
+        self._check_has_wavelength(wl)
+        fig, ax = plt.subplots(len(what))  # no sharex since theta might have less points
+        ax[-1].set_xlabel("Index")
+        fig.suptitle(f"Raw data for $\\lambda = {wl}$ nm")
+        for i, qty in enumerate(what):
+            ax[i].set_ylabel(PrsData.labels[qty])
+            ax[i].plot(self.data[wl][f"{qty}_raw"])
+        fig.tight_layout()
+        return fig
 
-def plot_cpd_data(data: str or pd.DataFrame or np.ndarray, t: str="seconds", title: str="", CPD:bool=True, LED:bool=True):
+
+def plot_spectrum(data: str or pd.DataFrame or np.ndarray, title: str="", errors=False, what=["dR_R"]):
     """
     Plot recorded data 
 
@@ -219,56 +392,41 @@ def plot_cpd_data(data: str or pd.DataFrame or np.ndarray, t: str="seconds", tit
     ----------
     data : str or np.ndarray
         Path to the data directory or 
-        numpy array with columns (idx, t [s], V [V], LED [%])
+        numpy array with columns PrsData.default_spectrum_columns
-    t : str, optional
+    :param title : str, optional
-        Which timescale to use for the x axis:
-        Must be one of "seconds", "mintutes", "hours".
-        The default is "seconds".
-    title : str, optional
         Title for the plot. The default is "".
-    CPD : bool, optional
+    :return fig Matplotlib figure object.
-        Wether to plot the voltage (CPD) line. The default is True.
-    LED : bool, optional
-        Wether to plot the LED state line. The default is False.
-
-    Returns
-    -------
-    fig : TYPE
-        Matplotlib figure object.
     """
     if type(data) == str:
-        _data, _ = PrsData.load_data_from_dir(data)
+        prsdata = PrsData(data)
+        _data = prsdata.get_spectrum_data()
+    elif isinstance(data, PrsData):
+        _data = data.get_spectrum_data()
     else:
         _data = data
-    fig, ax = plt.subplots()
+    n_plots = len(what)
-    xdata = _data[:,1].copy()
+    if errors: n_plots *= 2
-    xlabel = "t [s]"
+
-    if t == "minutes":
+    fig, ax = plt.subplots(n_plots, 1, sharex=True)
-        xdata /= 60
+    ax[-1].set_xlabel("Wavelength [nm]")
-        xlabel = "t [minutes]"
+    i = 0
-    elif t == "hours":
+    colors = {
-        xdata /= 3600
+        "dR_R": "red",
-        xlabel = "t [hours]"
+        "dR": "green",
-    ax.set_xlabel(xlabel)
+        "R": "blue",
-    ax_cpd = ax
+        "theta": "magenta",
-    ax_led = ax
+    }
-    if CPD and LED:
+
-        ax_led = ax.twinx()
+    wl_idx =  PrsData.default_spectrum_columns.index("wl")
-    if CPD:
+    for key in what:
-        ax_cpd = ax
+        key_and_err = [key, f"s{key}"] if errors else [key]
-        ax_cpd.set_ylabel("CPD [V]")
+        for k in key_and_err:
-        ax_cpd.plot(xdata, _data[:,2], color="blue", label="CPD")
+            data_idx = PrsData.default_spectrum_columns.index(k)
-    if LED:
+            ax[i].plot(_data[:,wl_idx], _data[:,data_idx], color=colors[key])
-        ax_led.set_ylabel("LED [%]")
+            ax[i].set_ylabel(PrsData.labels[k])
-        ax_led.plot(xdata, _data[:,3], color="orange", label="LED")
+            i += 1
-        ax_led.set_ylim(-2, 102)
-        ax_led.set_yticks([0, 20, 40, 60, 80, 100])
-    if CPD and LED:
-        # ax_led.legend()
-        # ax_cpd.legend()
-        pass
     if title:
-        ax.set_title(title)
+        fig.suptitle(title)
     fig.tight_layout()
     return fig
 

utility/testing.py

@@ -1,17 +0,0 @@
-import numpy as np
-
-def testcurve(x, frequency=10, peak_width=2, amplitude=20, bias=0):
-    # 0 = pk - width
-    # 2pi = pk + width
-    # want peak at n*time == frequency
-    nearest_peak = np.round(x / frequency, 0)
-    # if not peak at 0 and within peak_width
-    if nearest_peak > 0 and np.abs((x - nearest_peak * frequency)) < peak_width:
-        # return sin that does one period within 2*peak_width
-        return amplitude * np.sin(2*np.pi * (x - nearest_peak * frequency - peak_width) / (2*peak_width)) + bias
-    else:
-        return bias
-
-def get_testcurve(frequency=10, peak_width=2, amplitude=20, bias=0):
-    return np.vectorize(lambda x: testcurve(x, frequency=frequency, peak_width=peak_width, amplitude=amplitude, bias=bias))
-