element_array.py

import numpy as np
import matplotlib.pyplot as plt
import csv
import os

# ---------------------------
# Parameters (change if needed)
# ---------------------------
Nx = 20                   # elements in x
Ny = 20                   # elements in y
d_mm = 3.85               # inter-element spacing in mm (from paper)
d = d_mm / 1000.0         # spacing in meters
freq = 27.5e9             # Hz (27.5 GHz)
c = 299792458.0           # speed of light in m/s
wavelength = c / freq
k = 2 * np.pi / wavelength

# Feed (horn) position: distance 0.25 m at incidence -30 deg in y-z plane
feed_distance_m = 0.25
incidence_deg = -30.0     # negative = tilted toward -y
inc_rad = np.deg2rad(incidence_deg)
r_f = np.array([
    0.0,                                  # x
    feed_distance_m * np.sin(-inc_rad),   # y (negative for -30deg)
    feed_distance_m * np.cos(-inc_rad)    # z
], dtype=float)

# Desired reflected beam direction (u0) - example uses +z (normal reflection)
u0 = np.array([0.0, 0.0, 1.0], dtype=float)

# Optional global constant reference phase (Delta phi)
delta_phi = 0.0  # radians

# ON/OFF phase values (unit-cell measured values).
# Paper aims for ~180 deg difference; default: on=0, off=pi.
# Replace these if you have HFSS-derived phi_on and phi_off (radians).
phi_on = -np.pi
phi_off = -2*np.pi

# ---------------------------
# Build element position vectors
# ---------------------------
# Put array center at (0,0,0). Indices i,j = 1..20 map to coordinates:
# x = (i - (Nx+1)/2) * d, y = (j - (Ny+1)/2) * d
ix = np.arange(1, Nx + 1)
jy = np.arange(1, Ny + 1)
# meshgrid with i->x, j->y
I, J = np.meshgrid(ix, jy, indexing='ij')

x_coords = (I - (Nx + 1) / 2.0) * d   # meters
y_coords = (J - (Ny + 1) / 2.0) * d   # meters
z_coords = np.zeros_like(x_coords)

# Flattened shapes for computations
positions = np.stack((x_coords, y_coords, z_coords), axis=-1)  # shape (Nx,Ny,3)

# ---------------------------
# Compute ideal phase phi_ij
# ---------------------------
# phi_ij = k * |r_e - r_f| - k*(u0 . r_e) + delta_phi
# We'll compute in radians and wrap to [-pi, pi] for clarity.
def wrap_phase_rad(phase_rad):
    """Wrap to [-pi, pi)."""
    return (phase_rad + np.pi) % (2 * np.pi) - np.pi

ideal_phase = np.zeros((Nx, Ny), dtype=float)
for i in range(Nx):
    for j in range(Ny):
        r_e = positions[i, j, :]
        dist = np.linalg.norm(r_e - r_f)
        term1 = k * dist
        term2 = k * np.dot(u0, r_e)
        phi = term1 - term2 + delta_phi
        ideal_phase[i, j] = wrap_phase_rad(phi)

# Convert to degrees for presentation
ideal_phase_deg = np.rad2deg(ideal_phase)

# ---------------------------
# Quantize to 1-bit using Eq (2)
# ---------------------------
# Choose ON if |phi_ij - phi_on| <= |phi_ij - phi_off| else OFF
# But account for circular nature of phase by wrapping differences to [-pi, pi].
def circular_diff(a, b):
    """Smallest signed angular difference (a - b) in radians wrapped to [-pi,pi)."""
    diff = a - b
    return wrap_phase_rad(diff)

quantized_phase = np.zeros_like(ideal_phase)
quantized_state = np.zeros((Nx, Ny), dtype=int)  # 1 => ON, 0 => OFF
for i in range(Nx):
    for j in range(Ny):
        phi = ideal_phase[i, j]
        d_on = abs(circular_diff(phi, phi_on))
        d_off = abs(circular_diff(phi, phi_off))
        if d_on <= d_off:
            quantized_phase[i, j] = phi_on
            quantized_state[i, j] = 1
        else:
            quantized_phase[i, j] = phi_off
            quantized_state[i, j] = 0

quantized_phase_deg = np.rad2deg(quantized_phase)

# ---------------------------
# Save CSVs
# ---------------------------
out_dir = "ris_phase_output"
os.makedirs(out_dir, exist_ok=True)

def save_csv(matrix_deg, filename):
    path = os.path.join(out_dir, filename)
    with open(path, "w", newline="") as f:
        writer = csv.writer(f)
        # header: j indices
        writer.writerow(["i\\j"] + [f"{j}" for j in range(1, Ny+1)])
        for i_idx in range(Nx):
            row = [f"{i_idx+1}"] + [f"{matrix_deg[i_idx, j_idx]:.3f}" for j_idx in range(Ny)]
            writer.writerow(row)
    print(f"Saved {path}")

save_csv(ideal_phase_deg, "ideal_phase_deg.csv")
save_csv(quantized_phase_deg, "quantized_phase_deg.csv")
save_csv(quantized_state, "quantized_state.csv")

# ---------------------------
# Print some example vectors & phases (like earlier message)
# ---------------------------
examples = [(1,1), (10,10), (10,11), (11,10), (20,20)]
print("\nFeed position (m):", r_f)
print(f"Wavelength = {wavelength*1e3:.6f} mm, k = {k:.6e} rad/m\n")

print("Example element positions (m) and phases:")
for (ii, jj) in examples:
    r_e = positions[ii-1, jj-1]
    print(f" Element (i={ii}, j={jj}): r_e = ({r_e[0]:.4f}, {r_e[1]:.4f}, {r_e[2]:.4f}) m")
    print(f"   ideal phase = {ideal_phase_deg[ii-1, jj-1]:.3f} deg")
    print(f"   quantized phase = {quantized_phase_deg[ii-1, jj-1]:.3f} deg ; state = {quantized_state[ii-1, jj-1]}")
    print("")

# ---------------------------
# Plot heatmaps
# ---------------------------
fig, axes = plt.subplots(1, 2, figsize=(12,5))
im1 = axes[0].imshow(ideal_phase_deg.T, origin='lower', cmap='twilight', 
                     extent=[1, Nx, 1, Ny], aspect='equal')
axes[0].set_title("Ideal phase (deg)")
axes[0].set_xlabel("i index")
axes[0].set_ylabel("j index")
plt.colorbar(im1, ax=axes[0], fraction=0.046, pad=0.04)

from matplotlib.colors import ListedColormap, BoundaryNorm

# two colors you previously chose
color_off = "#52c5b6"   # cyan/teal for OFF
color_on  = "#d73027"   # red for ON

# quantized_phase_deg should be in degrees (may contain negative values)
# normalize to [0, 360)
quant_plot = (quantized_phase_deg + 360) % 360   # shape (Nx, Ny)

# Convert phi_on/phi_off (radians) to degrees and normalize
phi_on_deg  = (np.rad2deg(phi_on)  + 360) % 360
phi_off_deg = (np.rad2deg(phi_off) + 360) % 360

# Compute circular midpoint from phi_on_deg -> phi_off_deg
# d is the angular distance going from phi_on to phi_off in [0,360)
d = (phi_off_deg - phi_on_deg + 360) % 360
mid_deg = (phi_on_deg + d/2.0) % 360   # midpoint on the circular interval

# We want mid_deg to sit at 180° so that the two bins are [0,180) and [180,360)
# Compute rotation amount (how much to subtract from data to rotate midpoint -> 180)
rot = (mid_deg - 180.0)  # degrees to subtract (can be negative)
# Apply rotation and renormalize into [0,360)
quant_rot = (quant_plot - rot + 360.0) % 360.0

# Build 2-color colormap and use simple boundaries [0,180,360]
cmap_2 = ListedColormap([color_off, color_on])
bounds = [0.0, 180.0, 360.0]
norm = BoundaryNorm(bounds, cmap_2.N)

# Plot the rotated/normalized quantized map
im2 = axes[1].imshow(quant_rot.T, origin='lower', cmap=cmap_2, norm=norm,
                     extent=[1, Nx, 1, Ny], aspect='equal')

axes[1].set_title("Quantized 1-bit phase (deg)")
axes[1].set_xlabel("i index")
axes[1].set_ylabel("j index")

# Colorbar: show labels corresponding to the *original* (unrotated) center values
tick_rot = np.array([90.0, 270.0])            # centers of the two bins after rotation
tick_orig = (tick_rot + rot) % 360.0          # map back to original degrees

# Replace 0 with 360 for display (so bar shows 360° instead of 0°)
tick_labels = []
for t in tick_orig:
    # treat very-small values as 0
    if abs(t) < 1e-6:
        tick_labels.append("360°")
    else:
        tick_labels.append(f"{int(round(t))}°")

cbar = plt.colorbar(im2, ax=axes[1], fraction=0.046, pad=0.04, ticks=tick_rot)
cbar.set_ticklabels(tick_labels)

plt.tight_layout()
plt.show()

# ---------------------------
# Done
# ---------------------------
print("Done. Outputs saved in directory:", out_dir)