First commit

main.py

@@ -0,0 +1,138 @@
+from typing import Any
+from numpy import array, diff, mean, median
+from numpy.linalg import norm
+from rich.console import Console
+from pathlib import Path
+from matplotlib.pyplot import figure, show
+
+from utils import diff_matrix, Directions, Axis, Vector
+
+from framel import frgetvect
+
+console = Console()
+
+
+def frgetvectN(
+    filename: str | list[str],
+    channels: list[str],
+    start: int = -1,
+    span: int = -1,
+    verbose: bool = False,
+) -> list[Any]:
+    return [frgetvect(filename, channel, start, span, verbose) for channel in channels]
+
+
+def main(
+    file: Path,
+    start: int = 15,
+    duration: int = 2 * 5 * 49,
+    factors: list = [1.00, 0.90, 1.95, 1.19],
+    step: int = 50,
+    _show: bool = False,
+):
+    factors = array(factors)
+    indexes = diff_matrix(3)
+    steps = step * factors
+
+    channels_name = ["FitPosY", "FitPosX"]
+    channels = [
+        "V1:Camera_Scatter_{}".format(channel_name) for channel_name in channels_name
+    ]
+
+    data = frgetvectN(str(file), channels)
+
+    if _show:
+        Figure = figure()
+        for i in range(len(data)):
+            Figure.gca().plot(
+                data[i][0][start : start + duration], label=channels_name[i]
+            )
+        Figure.gca().legend()
+        show()
+
+    medians = [
+        median(datum[0][start : start + duration].reshape((-1, 10)), axis=1)
+        for datum in data
+    ]
+    positions = [
+        [diff(median)[index] for index in indexes] for median in medians
+    ]  # liste de mouvements associé à un axe pour une direction
+
+    means: Any = list()
+
+    for i in range(len(positions[0])):
+        if i < 2:
+            means.append(mean(positions[0][i] / steps[i]))
+        else:
+            means.append(mean(positions[1][i] / steps[i]))
+
+    means = array(means)
+
+    temps = list()
+    for i in range(len(positions[0])):
+        temp = list()
+        for j in range(len(positions)):
+            temp.append(mean(positions[j][i]) / steps[i])
+        temps.append(Vector(x=temp[0], y=temp[1]))
+
+    directions = Directions(
+        x=Axis(pos=temps[0], neg=temps[1]),
+        y=Axis(pos=temps[2], neg=temps[3]),
+    )
+
+    i, k = 1, norm(directions.x.pos) / norm(directions.y.pos)
+    j = (
+        directions.y.neg.x * (-k * directions.y.pos.y - i * directions.x.pos.y)
+        + k * directions.y.neg.y * directions.y.pos.x
+        + i * directions.x.pos.x * directions.y.neg.y
+    ) / (
+        directions.x.neg.y * directions.y.neg.x
+        - directions.x.neg.x * directions.y.neg.y
+    )
+    l = -(
+        directions.x.neg.x * (-k * directions.y.pos.y - i * directions.x.pos.y)
+        + k * directions.x.neg.y * directions.y.pos.x
+        + i * directions.x.neg.y * directions.x.pos.x
+    ) / (
+        directions.x.neg.y * directions.y.neg.x
+        - directions.x.neg.x * directions.y.neg.y
+    )
+
+    if _show:
+        Figure = figure()
+        namess = ["X", "Y"]
+        names = ["X+", "X-", "Y+", "Y-"]
+        for i in range(len(positions)):
+            position = positions[i]
+            name = namess[i]
+            for j in range(len(position)):
+                name_ = names[j]
+                Figure.gca().plot(position[j], label="{} {}".format(name, name_))
+        Figure.gca().legend()
+        show()
+
+    new_factors = array([i, j, k, l])
+
+    console.print("movement for one step:")
+    console.print(directions)
+    console.print("factors are: {}".format(new_factors))
+
+    console.print("-------------------")
+    console.print(new_factors * means)
+    console.print(factors * means)
+
+
+main(
+    Path(
+        "/home/demagny/data/backscattermeter/camera/snail_N3_SXp50_old_1-1429443262-102.gwf"
+    ),
+)
+
+main(
+    Path(
+        "/home/demagny/data/backscattermeter/camera/snail_N3_SXp50_3-1429442495-124.gwf"
+    ),
+    96,
+    2 * 5 * 49,
+    [1.00, 0.90, 1.95, 1.404],
+)

mise.toml

@@ -0,0 +1,5 @@
+[tools]
+python = "3.12.9"
+
+[env]
+_.python.venv = { path = ".venv", create = true, uv_create_args = ["--seed"] }

pyproject.toml

@@ -0,0 +1,14 @@
+[project]
+name = "camera_analysis"
+version = "1.0.0"
+description = "Tools to check beam path on camera"
+authors = [{ name = "linarphy", email = "augustin.demagny@ligo.org" }]
+readme = "README.md"
+requires-python = ">3.12"
+dependencies = ["framel", "gwpy>=3.0.12", "pyqt6>=6.9.0", "rich>=14.0.0"]
+
+[dependency-groups]
+dev = ["basedpyright>=1.29.0", "ruff>=0.11.6"]
+
+[tool.uv.sources]
+framel = { path = "../../../Fr/custom_python" }

utils.py

@@ -0,0 +1,164 @@
+from numpy import zeros
+from numpy.typing import NDArray
+from typing import NamedTuple, Any
+
+
+def compute_matrix(number: int):
+    """
+    Compute the index order (from 1 because it will be used by MATLAB) of measurmenet taken with a "snail" shape.
+    Used for the backscattermeter analysis MATLAB script
+
+    Parameters
+    ----------
+    number: int
+        Number of block from center to the extremity of the map
+    """
+    matrix = zeros((number * 2 + 1, number * 2 + 1), dtype=int)
+    cursor = [number, number]
+    counter = 1
+    matrix[*cursor] = counter
+    for i in range(number * 2):
+        for _ in range(i + 1):
+            if i % 2 == 0:
+                cursor[1] += 1
+            else:
+                cursor[1] -= 1
+            counter += 1
+            matrix[*cursor] = counter
+        for _ in range(i + 1):
+            if i % 2 == 0:
+                cursor[0] -= 1
+            else:
+                cursor[0] += 1
+            counter += 1
+            matrix[*cursor] = counter
+    for i in range(number):
+        cursor[1] += 1
+        counter += 1
+        matrix[*cursor] = counter
+    for i in range(number):
+        cursor[1] += 1
+        counter += 1
+        matrix[*cursor] = counter
+    return matrix
+
+
+def diff_matrix(number: int) -> list[list[int]]:
+    posX, negX, posY, negY = [], [], [], []
+    index = 0
+    for i in range(number * 2):
+        for _ in range(i + 1):
+            if i % 2 == 0:
+                posX.append(index)
+            else:
+                negX.append(index)
+            index += 1
+        for _ in range(i + 1):
+            if i % 2 == 0:
+                posY.append(index)
+            else:
+                negY.append(index)
+            index += 1
+    for _ in range(number):
+        posX.append(index)
+        index += 1
+    for _ in range(number):
+        posY.append(index)
+        index += 1
+    return [posX, negX, posY, negY]
+
+
+def python2matlab(array: NDArray[Any]):
+    """
+    Return the MATLAB representation of a 2d python numpy array
+
+    Parameters
+    ----------
+    array: NDArray[Any]
+        2d array
+    """
+    if len(array.shape) != 2:
+        raise ValueError("Can only show the representation of 2D numpy array")
+    result = "matrix = [ "
+    for column_index in range(array.shape[1]):
+        for row_index in range(array.shape[0]):
+            result += "{} ".format(array[row_index, column_index])
+        if column_index != array.shape[1] - 1:
+            result += "; "
+    result += "]';"
+    return result
+
+
+class Vector(NamedTuple):
+    x: float
+    y: float
+
+    def __truediv__(self, other: "Vector | float | int") -> "Vector":
+        if type(other) is type(self):
+            return Vector(x=self.x / other.x, y=self.y / other.y)
+        return Vector(x=self.x / other, y=self.y / other)
+
+    def __mul__(self, other: "Vector | float | int") -> "Vector":
+        if type(other) is type(self):
+            return Vector(x=self.x * other.x, y=self.y * other.y)
+        return Vector(x=self.x * other, y=self.y * other)
+
+    def add(self, other: "Point | Vector | float | int") -> "Vector":
+        if type(other) == type(self) or type(other) is Point:
+            return Vector(x=self.x + other.x, y=self.y + other.y)
+        return Vector(x=self.x + other, y=self.y + other)
+
+
+class Point(NamedTuple):
+    x: float
+    y: float
+
+    def add(self, other: "Point | Vector | float | int") -> "Point":
+        if type(other) == type(self) or type(other) is Vector:
+            return Point(x=self.x + other.x, y=self.y + other.y)
+        return Point(x=self.x + other, y=self.y + other)
+
+
+class Axis(NamedTuple):
+    pos: Vector
+    neg: Vector
+
+
+class Directions(NamedTuple):
+    x: Axis
+    y: Axis
+
+
+def convert_to_base(base: tuple[Vector, Vector], vector: Vector) -> Vector:
+    y = (vector.y - (base[0].y * vector.x / base[0].x)) / (
+        base[1].y - (base[0].y * base[1].x / base[0].x)
+    )
+    x = vector.x / base[0].x - y * (base[1].x / base[0].x)
+    return Vector(x=x, y=y)
+
+
+def move_snail(
+    n_grid: int, point: Point, directions: Directions
+) -> tuple[list[Vector], Point]:
+    movement: list[Vector] = []
+    for i in range(2 * n_grid):
+        for _ in range(i + 1):
+            if (i % 2) == 0:
+                movement.append(directions.x.pos)
+            else:
+                movement.append(directions.x.neg)
+            point = point.add(movement[-1])
+        for _ in range(i + 1):
+            if (i % 2) == 0:
+                movement.append(directions.y.pos)
+            else:
+                movement.append(directions.y.neg)
+            point = point.add(movement[-1])
+    # go back
+    for i in range(n_grid):
+        movement.append(directions.x.pos)
+        point = point.add(movement[-1])
+    for i in range(n_grid):
+        movement.append(directions.y.pos)
+        point = point.add(movement[-1])
+    return movement, point