def uniaxial_reflectivity(q: NDArray, layers: NDArray, tensor: NDArray, energy: float):
    """
    EMpy implementation of the uniaxial 4x4 matrix formalism.

    for calculating reflectivity from a stratified
    medium.

    Uses the implementation developed by FSRStools -
    https://github.com/ddietze/FSRStools - written by Daniel Dietze

    Parameters
    ----------
    q: np.ndarray
        q values for which to calculate the reflectivity.
        Units: 1/Angstroms

    layers: np.ndarray
        coefficients required for the calculation, has shape (2 + N, 4),
        where N is the number of layers
        - layers[0, 1] - SLD of fronting (/1e-6 Angstrom**-2)
        - layers[0, 2] - iSLD of fronting (/1e-6 Angstrom**-2)
        - layers[N, 0] - thickness of layer N
        - layers[N, 1] - SLD of layer N (/1e-6 Angstrom**-2)
        - layers[N, 2] - iSLD of layer N (/1e-6 Angstrom**-2)
        - layers[N, 3] - roughness between layer N-1/N
        - layers[-1, 1] - SLD of backing (/1e-6 Angstrom**-2)
        - layers[-1, 2] - iSLD of backing (/1e-6 Angstrom**-2)
        - layers[-1, 3] - roughness between backing and last layer

    tensor: np.ndarray
        contains the 1x3x3 dimensions
        First dimension may change in teh fiture to account for multi-energy
        currently it will just cycle
    scale: float
        Multiply all reflectivities by this value.
    bkg: float
        Linear background to be added to all reflectivities

    Returns
    -------
    Reflectivity: np.ndarray
        Calculated reflectivity values for each q value.
    """
    # Plane of incidence - required to define polarization vectors
    OpticAxis = np.array([0.0, 0.0, 1.0])
    phi = 0  # This is a uniaxial calculation

    ##Organize qvals into proper order
    qvals = np.asarray(q)
    flatq = qvals.ravel()
    numpnts = flatq.size  # Number of q-points

    ##Grab the number of layers
    nlayers = layers.shape[0]

    ##hc has been converted into eV*Angstroms
    wl = (
        hc / energy
    )  ##calculate the wavelength array in Aangstroms for layer calculations
    k0 = 2 * np.pi / (wl)

    # Convert optical constants into dielectric tensor
    tensor = np.conj(np.eye(3) - 2 * tensor[:, :, :])  # From tensor[:,0,:,:]

    # freq = 2*np.pi * c/wls #Angular frequency
    theta_exp = np.zeros(numpnts, dtype=float)
    theta_exp = np.pi / 2 - np.arcsin((flatq[:] / (2 * k0)).clip(-1, 1))

    ##Generate arrays of data for calculating transfer matrix
    ##Scalar values ~~
    ## Special cases!
    ##kx is constant for each wavelength but changes with angle
    ## Dimensionality ##
    ## (angle)
    # kx = np.zeros(numpnts, dtype=complex)
    # ky = np.zeros(numpnts, dtype=complex) #Used to keep the k vectors three
    # components later on for cross / dot products
    kx = k0 * np.sin(theta_exp) * np.cos(phi)
    ky = k0 * np.sin(theta_exp) * np.sin(phi)

    ## Calculate the eigenvalues corresponding to kz ~~ Each one has 4 solutions
    ## Dimensionality ##
    ## (angle, #layer, solution)
    kz = np.zeros((numpnts, nlayers, 4), dtype=complex)

    ## Calculate the eignvectors corresponding to each kz ~~ polarization of D and H
    ## Dimensionality ##
    ## (angle, #layers, solution, vector)
    Dpol = np.zeros(
        (numpnts, nlayers, 4, 3), dtype=complex
    )  ##The polarization of the displacement field
    Hpol = np.zeros(
        (numpnts, nlayers, 4, 3), dtype=complex
    )  ##The polarization of the magnetic field

    # Cycle through the layers and calculate kz
    for j, epsilon in enumerate(
        tensor
    ):  # Each layer will have a different epsilon and subsequent kz
        kz[:, j, :] = calculate_kz_uni(epsilon, kx, ky, k0, opticaxis=OpticAxis)
        Dpol[:, j, :, :], Hpol[:, j, :, :] = calculate_Dpol_uni(
            epsilon,
            kx,
            ky,
            kz[:, j, :],
            k0,
            opticaxis=OpticAxis,  # type: ignore
        )

    ##Make matrices for the transfer matrix calculation
    ##Dimensionality ##
    ##(angles, #layers, Matrix (4,4)

    ## Propogation Matrix
    P = calculate_P(
        numpnts, nlayers, kz[:, :, :], layers[:, 0]
    )  ##layers[k,0] is the thicknes of layer k
    ##Nevot-Croche roughness matrix
    W = calculate_W(numpnts, nlayers, kz[:, :, :], kz[:, :, :], layers[:, 3])
    ##Dynamic Matrix and inverse
    D, Di = calculate_D(numpnts, nlayers, Dpol[:, :, :, :], Hpol[:, :, :, :])

    ##Calculate the full system transfer matrix
    ##Dimensionality ##
    ##(angles, Matrix (4,4))
    M = np.ones((numpnts, 4, 4), dtype=complex)
    # Make a (numpnts x 4x4) identity matrix for the TMM -
    M = np.einsum("...ij,ij->...ij", M, np.identity(4))  # type: ignore
    M = calculate_TMM(numpnts, nlayers, M, D, Di, P, W)
    ##Calculate the final outputs and organize into the appropriate waves for later
    refl, tran = calculate_output(numpnts, M)

    return refl, tran, kx, ky, kz, Dpol, Hpol, D, Di, P, W, M