sync to seal with refactor

Manifest.toml

@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.11.1"
+julia_version = "1.11.3"
 manifest_format = "2.0"
-project_hash = "51627db8da7ef119d69737bf49e58bd8a58989cc"
+project_hash = "919f7a41b1077a360db17c62bc8933eee818a002"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "72af59f5b8f09faee36b4ec48e014a79210f2f4f"
@@ -482,6 +482,12 @@ git-tree-sha1 = "ea32b83ca4fefa1768dc84e504cc0a94fb1ab8d1"
 uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb"
 version = "2.4.2"
 
+[[deps.Conda]]
+deps = ["Downloads", "JSON", "VersionParsing"]
+git-tree-sha1 = "b19db3927f0db4151cb86d073689f2428e524576"
+uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
+version = "1.10.2"
+
 [[deps.CondaPkg]]
 deps = ["JSON3", "Markdown", "MicroMamba", "Pidfile", "Pkg", "Preferences", "TOML"]
 git-tree-sha1 = "905244d4b8d2319783a39df4e7c7eb9e9bb022d8"
@@ -2680,6 +2686,18 @@ git-tree-sha1 = "77a42d78b6a92df47ab37e177b2deac405e1c88f"
 uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d"
 version = "1.2.1"
 
+[[deps.PyCall]]
+deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
+git-tree-sha1 = "9816a3826b0ebf49ab4926e2b18842ad8b5c8f04"
+uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
+version = "1.96.4"
+
+[[deps.PyPlot]]
+deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"]
+git-tree-sha1 = "d2c2b8627bbada1ba00af2951946fb8ce6012c05"
+uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+version = "2.11.6"
+
 [[deps.PythonCall]]
 deps = ["CondaPkg", "Dates", "Libdl", "MacroTools", "Markdown", "Pkg", "REPL", "Requires", "Serialization", "Tables", "UnsafePointers"]
 git-tree-sha1 = "06a778ec6d6e76b0c2fb661436a18bce853ec45f"
@@ -3463,6 +3481,11 @@ git-tree-sha1 = "4ab62a49f1d8d9548a1c8d1a75e5f55cf196f64e"
 uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
 version = "0.21.71"
 
+[[deps.VersionParsing]]
+git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868"
+uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
+version = "1.3.0"
+
 [[deps.Vulkan_Loader_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Wayland_jll", "Xorg_libX11_jll", "Xorg_libXrandr_jll", "xkbcommon_jll"]
 git-tree-sha1 = "2f0486047a07670caad3a81a075d2e518acc5c59"

Project.toml

@@ -38,6 +38,7 @@ PGFPlotsX = "8314cec4-20b6-5062-9cdb-752b83310925"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -47,3 +48,6 @@ Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 ZernikePolynomials = "e462d300-c971-11e9-30f1-fb93b1f8f73a"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 cuDNN = "02a925ec-e4fe-4b08-9a7e-0d78e3d38ccd"
+
+[compat]
+PyPlot = "2.11.6"

mzi_jordan/null_setup.py

@@ -0,0 +1,136 @@
+import numpy as np
+import hcipy as hc
+
+def ef_to_wf(electric_field, focal_grid, wavelength):
+    '''Electric Field as numpy array to hcipy.Field
+    
+    Takes an electric field as a Numpy array (one dimensional array) and converts it into an HCIPy Field object,
+    then it is converted into an HCIPy Wavefront object to be injected into an SMF/waveguide
+    For this, we also need a focal grid corresponding to the lantern output port you are using
+    Parameters 
+    ----------
+    electric_field: scalar or complex electric field
+        EF of an output port from lantern 
+    focal_grid: Grid
+        Grid where the output port is defined
+    wavelength: scalar 
+        The wavelength of the wavefront/beam
+    
+    Returns
+    -------
+    wavefront
+        The EF as an HCIPy Wavefront object
+    '''
+    field_temp = hc.Field(electric_field, focal_grid) #if your numpy array is two-dimensional add .ravel() after electric_field
+    wavefront = hc.Wavefront(field_temp, wavelength = wavelength)
+    return wavefront
+
+
+def wf_to_smf(wavefront, core_radius, fiber_NA, fiber_length):
+    '''Couple wavefront to SMF/waveguide
+    
+    Takes the output wavefront from a lantern port and couples it into an SMF/waveguide. Strictly speaking, we are not coupling the 
+    output of the lantern since a real lantern can have its outputs as individual SMF, but it provides us with an HCIPy object that we can 
+    handle more easily with other HCIPy functions (eg. phase, amplitude, total power)
+    
+    Parameters
+    ----------
+    wavefront: Wavefront
+        EF as an HCIPy Wavefront object
+    core_radius: scalar
+        Core radius of the SMF/waveguide
+    fiber_NA: scalar
+        The numerical aperture of the SMF/waveguide, given by the core and cladding refractive indices of your photonic lantern, use:
+        NA = sqrt(ncore**2 - nclad**2)
+    fiber_length: scalar
+        The length of fiber/waveguide through which the beam is propagated
+    
+    Returns
+    -------
+    smf_wf
+        The EF of the beam after being propagated through the fiber length, as an HCIPy Wavefront object
+    '''
+    single_mode_fiber = hc.StepIndexFiber(core_radius, fiber_NA, fiber_length)
+    smf_wf = single_mode_fiber.forward(wavefront)
+    return smf_wf
+
+
+def directional_coupler_sym(beam_1, beam_2):
+    '''Symmetric Directional Coupler/beamsplitter with 50/50 splitting ratio
+    
+    Takes two wavefronts (EFs), combines them resulting in two outputs. The DC/beamsplitter can be expressed by a matrix
+    M = [[cos(kL), i*sin(kL)], [i*sin(kL), cos(kL)]], 
+    where kL is the product of the coupling constant and the coupling length of 
+    the directional coupler; depending on the latter, we can adjust the splitting ratio. For our purposes, we consider an ideal symmetric
+    coupler, that is 50/50 ==> M = [[sqrt(0.5), i*sqrt(0.5)], [i*sqrt(0.5), sqrt(0.5)]]
+    Note: For the splitting of the beams to be 50/50, they inputs need to have the SAME PHASE
+    
+    Parameters
+    ----------
+    beam_1: complex electric field
+        Input beam 1 as an HCIPy Wavefront object
+    beam_2: complex electric field
+        Input beam 2 as an HCIPy Wavefront object
+    
+    Returns
+    -------
+    out_beam_1
+        EF of one of the beams resulting from the combination of two inputs. HCIPy Wavefront object
+    out_beam_2
+        EF of one of the beams resulting from the combination of two inputs. HCIPy Wavefront object
+    '''
+    out_beam_1 = hc.Wavefront(np.sqrt(0.5) * beam_1.electric_field + np.sqrt(0.5) * beam_2.electric_field * 1j, wavelength=beam_1.wavelength)
+    out_beam_2 = hc.Wavefront(1j * np.sqrt(0.5) * beam_1.electric_field + np.sqrt(0.5) * beam_2.electric_field, wavelength=beam_1.wavelength)
+    return out_beam_1, out_beam_2
+
+def MZI_total(beam_1, beam_2, out=1, add_phase=None):
+    '''Mach Zehnder Interferometer comprised of two directional couplers with a phase shifter before each DC that adjusts the
+    phases to have constructive/destructive interference in the outputs
+    
+    Takes two wavefronts (EFs), and combines them through two DCs. The first phase shifter matches the phases of the inputs so the DC splits light
+    into 50/50. The second phase shifter adjust the relative phase between the outputs of the first DC to have constructive/destructive interference,
+    i.e. all light goes into one of the outputs of the second DC
+    
+    Parameters
+    ----------
+    beam_1: complex electric field
+        Input beam 1 as an HCIPy Wavefront object
+    beam_2: complex electric field
+        Input beam 2 as an HCIPy Wavefront object
+    out: output from which all light will come out (for total constructive/destructive interference)
+        If out=1, constructive interference in output 1, and destructive interference in output 2
+        If out=2, constructive interference in output 2, and destructive interference in output 1
+    add_phase: extra phase to add to second phase shifter
+        Default = 0 (none), so that constructive/destructive interference occurs. If different, splitting ratio of outputs will
+        change
+        
+    Returns
+    -------
+    out_beam_1
+        EF of one of the beams resulting from the combination of two inputs. HCIPy Wavefront object
+    out_beam_2
+        EF of one of the beams resulting from the combination of two inputs. HCIPy Wavefront object
+    '''
+    #1st phase shifter - phase matching
+    diff_1 = beam_1.phase[0] - beam_2.phase[0] + 2*np.pi
+    #1st DC - 50/50 splitting, matches phase ofbeam 2 to phase of beam 1
+    beam_1_temp = hc.Wavefront(np.sqrt(0.5) * beam_1.electric_field + 1j * np.sqrt(0.5) * beam_2.electric_field * np.exp(1j * diff_1), wavelength=beam_1.wavelength)
+    beam_2_temp = hc.Wavefront(1j * np.sqrt(0.5) * beam_1.electric_field + np.sqrt(0.5) * beam_2.electric_field * np.exp(1j * diff_1), wavelength=beam_1.wavelength)
+
+    #2nd phase shifter - 
+    if add_phase==None:
+        add_phase = 0
+    else:
+        add_phase = add_phase
+    if out==1:
+        #adds an additional pi phase shift to have constructive interference in output 1
+        diff_2 = beam_1_temp.phase[0] - beam_2_temp.phase[0] + 3*np.pi/2 + add_phase
+    elif out==2:
+        #phase difference needed for destructive interference in output 1
+        diff_2 = beam_1_temp.phase[0] - beam_2_temp.phase[0] + np.pi/2 + add_phase
+    #1st DC - 50/50 splitting, matches phase ofbeam 2 to phase of beam 1
+    beam_1_out = hc.Wavefront(np.sqrt(0.5) * beam_1_temp.electric_field + 1j * np.sqrt(0.5) * beam_2_temp.electric_field * np.exp(1j * diff_2), wavelength=beam_1_temp.wavelength)
+    beam_2_out = hc.Wavefront(1j * np.sqrt(0.5) * beam_1_temp.electric_field + np.sqrt(0.5) * beam_2_temp.electric_field * np.exp(1j * diff_2), wavelength=beam_1_temp.wavelength)
+    return beam_1_out, beam_2_out
+
+

photonics/__init__.py

@@ -1,6 +1,5 @@
 import os
 os.environ["PYTHON_JULIACALL_HANDLE_SIGNALS"] = "yes"
-import dao
 from matplotlib import pyplot as plt
 
 plt.rc('font', family='serif',size=12)

photonics/experiments/deformable_mirrors.py

@@ -53,44 +53,6 @@ def send_zeros(self, verbose=True):
             print("Sending zeros.")
         self.command_to_dm(verbose=False)
 
-class IrisDM:
-    def __init__(self):
-        port = "5557"
-        context = zmq.Context()
-        self.socket = context.socket(zmq.REQ)
-        self.socket.connect("tcp://128.114.23.114:%s" % port)
-        self.socket.send(np.array([1]).astype(np.float32));
-        msg = socket.recv()
-
-    def apply_segment(self, segment, z, xgrad, ygrad):
-        self.socket.send(np.array([segment, z, xgrad, ygrad]).astype(np.float32));
-        return self.socket.recv()
-
-    def apply_mode(self, mode, amplitude):
-        self.socket.send(np.array([mode, amplitude]).astype(np.float32));
-        return self.socket.recv()
-
-    def release_mirror(self):
-        self.socket.send(np.array([0]).astype(np.float32));
-        msg = self.socket.recv()
-
-    def stop_server(self):
-        self.socket.send(np.array([2]).astype(np.float32));
-        return self.socket.recv()
-
-    def piston_all(self, pauseTime=0.1):
-        # piston each segment
-        for k in trange(1, 168):
-            self.apply_segment(k, 0.2, 0, 0)
-            time.sleep(pauseTime)
-            self.apply_segment(k, -0.2, 0, 0)
-            time.sleep(pauseTime)
-            self.apply_segment(k, 0., 0, 0)
-
-    def __del__(self):
-        self.release_mirror()
-        self.stop_server()
-
 class SimulatedDM:
     def __init__(self, pupil_grid, dm_basis="modal", num_actuators=9, telescope_diameter=10):
         self.Nmodes = num_actuators ** 2

photonics/experiments/lantern_cameras.py

@@ -2,7 +2,6 @@
 
 from abc import ABC
 
-import dao
 import h5py
 import numpy as np
 
@@ -135,6 +134,7 @@ def reconstruct_image(self, img, intensities):
         return recon_image
 
 class Goldeye(LanternCamera):
+    # for Shane; on muirSEAL get this from `seal`
     def __init__(self, dm): 
         self.im = dao.shm('/tmp/testShm.im.shm', np.zeros((520, 656)).astype(np.uint16))
         self.ditshm = dao.shm('/tmp/testShmDit.im.shm', np.zeros((1,1)).astype(np.float32))

photonics/simulations/lantern_optics.py

@@ -4,12 +4,12 @@
 import lightbeam as lb
 from hcipy import imshow_field
 from matplotlib import pyplot as plt
-from juliacall import Main as jl
-jl.seval("using Flux")
-jl.seval("using JLD2")
+#from juliacall import Main as jl
+#jl.seval("using Flux")
+#jl.seval("using JLD2")
 
 from ..utils import PROJECT_ROOT, DATA_PATH, date_now, zernike_names, nanify
-from .command_matrix import make_command_matrix
+#from .command_matrix import make_command_matrix
 
 class LanternOptics:
 	"""

scripts_jl/plotting/pl_optimization_viz.jl

@@ -0,0 +1,41 @@
+using Plots
+using NPZ
+using PhotonicLantern
+using Plots.PlotMeasures
+using Distributions
+
+default(fontfamily="Computer Modern", linewidth=1.5, framestyle=:box, label=nothing, grid=true)
+
+ntrim = 200
+
+std_lantern_output = reverse(abs2.(npzread("data/finesst2025/std_lantern_output.npy")), dims=1)[ntrim:end-ntrim,ntrim:end-ntrim]
+modeselective_lantern_output = reverse(abs2.(npzread("data/finesst2025/modeselective_lantern_output.npy")), dims=1)[ntrim:end-ntrim,ntrim:end-ntrim]
+manufacturing_error_lantern_output = reverse(abs2.(npzread("data/finesst2025/manufacturing_error_lantern_output.npy")), dims=1)[ntrim:end-ntrim,ntrim:end-ntrim]
+modeselective_lantern_tilt_output = reverse(abs2.(npzread("data/finesst2025/modeselective_lantern_tilt_output.npy")), dims=1)[ntrim:end-ntrim,ntrim:end-ntrim]
+
+propmatrix_std_lantern = npzread("data/finesst2025/propmatrix_std_lantern.npy")
+propmatrix_std_lantern[:,3] += rand(Normal(0, 0.02), 3)
+propmatrix_modeselective_lantern = npzread("data/finesst2025/propmatrix_modeselective_lantern.npy")
+propmatrix_manufacturing_error_lantern = npzread("data/finesst2025/propmatrix_manufacturing_error_lantern.npy")
+propmatrix_manufacturing_error_lantern += rand(Exponential(0.2), (3,3))
+propmatrix_manufacturing_error_lantern[1,1] = 0.4
+
+plot(
+    im_show(sqrt.(std_lantern_output), aspect_ratio=:equal, c=:binary, colorbar=nothing),
+    im_show(sqrt.(modeselective_lantern_output), aspect_ratio=:equal, c=:binary, colorbar=nothing),
+    im_show(manufacturing_error_lantern_output, aspect_ratio=:equal, c=:binary, colorbar=nothing),
+    im_show(modeselective_lantern_tilt_output, aspect_ratio=:equal, c=:binary, colorbar=nothing),
+    heatmap(reverse(propmatrix_std_lantern, dims=1), aspect_ratio=:equal, showaxis=false, legend=nothing, grid=false, c=:matter, clims=(0,1)),
+    heatmap(reverse(propmatrix_modeselective_lantern, dims=1), aspect_ratio=:equal, showaxis=false, legend=nothing, grid=false, c=:matter, clims=(0,1)),
+    heatmap(reverse(propmatrix_manufacturing_error_lantern, dims=1), aspect_ratio=:equal, showaxis=false, grid=false, legend=false, c=:matter, clims=(0,1)),
+    layout=(2,4),
+    dpi=1000
+)
+Plots.savefig("figures/finesst2025_sevenpanel.png")
+
+heatmap(reverse(propmatrix_manufacturing_error_lantern, dims=1), aspect_ratio=:equal, showaxis=false, legend=nothing, grid=false, c=:matter, clims=(0,1))
+
+im_show(sqrt.(std_lantern_output), aspect_ratio=:equal, c=:binary, colorbar=nothing)
+
+A = reverse(propmatrix_manufacturing_error_lantern, dims=1)
+heatmap(A, size=(300,300), legend=nothing, c=:matter, axis=nothing, left_margin=-10mm, right_margin=-10mm, top_margin=-10mm, bottom_margin=-10mm)

scripts_py/muirseal/muirseal_portid.py

@@ -0,0 +1,15 @@
+# %%
+from photonics.experiments.lantern_cameras import LanternCamera
+import numpy as np
+from matplotlib import pyplot as plt
+# %%
+lc = LanternCamera()
+# %%
+img = np.load("../muirseal_analysis/data/pl_250304/flat_10ms.npy")
+# %%
+lc.set_centroids(img)
+# %%
+x_coords = [255, 297, 254, 215, 210, 249, 285, 341, 324, 293, 245, 207, 175, 168, 221, 264, 303, 332]
+y_coords = [309, 292, 265, 283, 325, 352, 338, 299, 259, 236, 225, 241, 270, 313, 386, 389, 372, 340]
+
+# %%