import numpy as np
            import matplotlib.pyplot as plt
            
            def z_scan_intensity(z, z_R, P_0, n_2):
                """Calculates the transmitted intensity of a Z-scan.
            
                Args:
                    z: Array of z-positions.
                    z_R: Rayleigh length.
                    P_0: Input power.
                    n_2: Nonlinear refractive index.
            
                Returns:
                    Array of transmitted intensities.
                """
            
                # Calculate the on-axis intensity
                I_0 = P_0 / (np.pi * z_R**2)
                I = I_0 * (1 + (z / z_R)**2)**(-1)
            
                # Calculate the nonlinear phase shift
                Δϕ = (2 * np.pi * n_2 * I_0 * L_eff) / λ
            
                # Calculate the transmitted intensity
                T = 1 - 2 * np.sin(Δϕ / 2)**2
            
                return T
            
            def plot_z_scan(z, T):
                """Plots the Z-scan curve.
            
                Args:
                    z: Array of z-positions.
                    T: Array of transmitted intensities.
                """
            
                plt.plot(z, T)
                plt.xlabel("z (mm)")
                plt.ylabel("Transmitted Intensity")
                plt.title("Z-Scan Curve")
                plt.grid(True)
                plt.show()
            
            # Example usage
            z_R = 10  # mm
            P_0 = 1000  # mW
            n_2 = 1e-16  # cm^2/W
            L_eff = 1  # cm
            λ = 532e-7  # cm
            
            z = np.linspace(-3 * z_R, 3 * z_R, 1000)
            T = z_scan_intensity(z, z_R, P_0, n_2)
            
            plot_z_scan(z, T)

import numpy as np

            def mueller_matrix(theta, phi, delta):
                """
                Calculates the Mueller matrix of a waveplate or polarizer.
            
                Args:
                    theta (float): The angle of the fast axis (or transmission axis) relative to the x-axis.
                    phi (float): The phase retardation (in radians) introduced by the waveplate.
                    delta (float): The extinction ratio of the polarizer (ratio of maximum to minimum transmittance).
            
                Returns:
                    np.ndarray: The 4x4 Mueller matrix.
                """
            
                cos_theta = np.cos(theta)
                sin_theta = np.sin(theta)
                cos_phi = np.cos(phi)
                sin_phi = np.sin(phi)
                cos_delta = np.cos(delta)
                sin_delta = np.sin(delta)
            
                # Calculate the Mueller matrix elements
                M11 = cos_theta**2 * cos_phi**2 + sin_theta**2 * sin_phi**2
                M12 = cos_theta**2 * sin_phi**2 + sin_theta**2 * cos_phi**2
                M13 = cos_theta * sin_theta * (cos_phi**2 - sin_phi**2)
                M14 = cos_theta * sin_theta * 2 * cos_phi * sin_phi
                M21 = M12
                M22 = M11
                M23 = -M13
                M24 = -M14
                M31 = M13
                M32 = -M13
                M33 = cos_theta**2 * sin_phi**2 + sin_theta**2 * cos_phi**2
                M34 = -cos_theta * sin_theta * 2 * cos_phi * sin_phi
                M41 = M14
                M42 = -M14
                M43 = M34
                M44 = cos_theta**2 * cos_phi**2 + sin_theta**2 * sin_phi**2
            
                # Construct the Mueller matrix
                M = np.array([[M11, M12, M13, M14],
                               [M21, M22, M23, M24],
                               [M31, M32, M33, M34],
                               [M41, M42, M43, M44]])
            
                return M
            
            # Example usage
            theta = np.pi / 4  # 45 degrees
            phi = np.pi / 2  # 90 degrees
            delta = 1e-6  # Small extinction ratio for a polarizer
            
            mueller_matrix = mueller_matrix(theta, phi, delta)
            print(mueller_matrix)