1  #!/usr/bin/env python3
  2  import numpy as np
  3  import matplotlib.pyplot as plt
  4  
  5  c0 = 299_792_458.0
  6  f0 = 10500000000.0
  7  lam0 = c0/f0
  8  k0 = 2*np.pi/lam0
  9  
 10  M = 16
 11  N = 32
 12  dy = 14.275831333333334/1e3
 13  dz = 16.915/1e3
 14  
 15  theta0_deg = 0.0
 16  phi0_deg   = 0.0
 17  theta0 = np.deg2rad(theta0_deg)
 18  phi0   = np.deg2rad(phi0_deg)
 19  
 20  beta = 1.65
 21  wy = np.ones(16, float)
 22  wz = np.kaiser(32, beta)
 23  wz /= wz.max()
 24  
 25  m_idx = np.arange(M) - (M-1)/2
 26  n_idx = np.arange(N) - (N-1)/2
 27  y_positions = m_idx * dy
 28  z_positions = n_idx * dz
 29  
 30  def element_factor(theta_rad, phi_rad):
 31      return np.abs(np.cos(theta_rad))
 32  
 33  def array_factor(theta_rad, phi_rad, y_positions, z_positions, wy, wz, theta0_rad, phi0_rad):
 34      k0 = 2*np.pi/(299_792_458.0/10500000000.0)
 35      ky = k0*np.sin(theta_rad)*np.sin(phi_rad)
 36      kz = k0*np.sin(theta_rad)*np.cos(phi_rad)
 37      ky0 = k0*np.sin(theta0_rad)*np.sin(phi0_rad)
 38      kz0 = k0*np.sin(theta0_rad)*np.cos(phi0_rad)
 39      Ay = np.sum(wy[:,None] * np.exp(1j * y_positions[:,None]*(ky[None,:]-ky0)), axis=0)
 40      Az = np.sum(wz[:,None] * np.exp(1j * z_positions[:,None]*(kz[None,:]-kz0)), axis=0)
 41      return Ay * Az
 42  
 43  def cut_curve(phi_deg, num_pts=721):
 44      th_deg = np.linspace(0, 90, num_pts)
 45      th = np.deg2rad(th_deg)
 46      ph = np.deg2rad(phi_deg) * np.ones_like(th)
 47      AF = array_factor(th, ph, y_positions, z_positions, wy, wz, theta0, phi0)
 48      PAT = np.abs(AF) * element_factor(th, ph)
 49      PAT /= PAT.max()
 50      return th_deg, 20*np.log10(PAT + 1e-15)
 51  
 52  def hpbw_deg(theta_deg, pat_db):
 53      p = np.argmax(pat_db)
 54      peak = pat_db[p]
 55      mask = pat_db >= (peak - 3.0)
 56      idx = np.where(mask)[0]
 57      if len(idx) < 2:
 58          return np.nan
 59      return float(theta_deg[idx[-1]] - theta_deg[idx[0]])
 60  
 61  thE_deg, patE_db = cut_curve(0.0)
 62  bwE = hpbw_deg(thE_deg, patE_db)
 63  plt.figure(figsize=(7,5), dpi=130)
 64  plt.plot(thE_deg, patE_db, linewidth=1.5)
 65  plt.grid(True, linestyle='--', alpha=0.5)
 66  plt.xlabel('Theta (deg)')
 67  plt.ylabel('Normalized Gain (dB)')
 68  plt.title(f'E-plane (phi=0°)  |  -3 dB BW ≈ {bwE:.2f}°')
 69  plt.tight_layout()
 70  plt.savefig('E_plane_Kaiser25dB_like.png', bbox_inches='tight')
 71  plt.show()
 72  
 73  thH_deg, patH_db = cut_curve(90.0)
 74  bwH = hpbw_deg(thH_deg, patH_db)
 75  plt.figure(figsize=(7,5), dpi=130)
 76  plt.plot(thH_deg, patH_db, linewidth=1.5)
 77  plt.grid(True, linestyle='--', alpha=0.5)
 78  plt.xlabel('Theta (deg)')
 79  plt.ylabel('Normalized Gain (dB)')
 80  plt.title(f'H-plane (phi=90°) |  -3 dB BW ≈ {bwH:.2f}°')
 81  plt.tight_layout()
 82  plt.savefig('H_plane_Kaiser25dB_like.png', bbox_inches='tight')
 83  plt.show()
 84  
 85  theta_deg = np.linspace(0.0, 90.0, 121)
 86  phi_deg   = np.linspace(-90.0, 90.0, 121)
 87  TH, PH = np.meshgrid(theta_deg, phi_deg, indexing='xy')
 88  PAT_db = np.empty_like(TH, dtype=float)
 89  for i in range(TH.shape[0]):
 90      th = np.deg2rad(TH[i, :])
 91      ph = np.deg2rad(PH[i, :])
 92      AF = array_factor(th, ph, y_positions, z_positions, wy, wz, theta0, phi0)
 93      EF = element_factor(th, ph)
 94      pat = np.abs(AF)*EF
 95      pat /= pat.max()
 96      PAT_db[i, :] = 20*np.log10(pat + 1e-15)
 97  
 98  plt.figure(figsize=(7,5), dpi=130)
 99  extent = [theta_deg[0], theta_deg[-1], phi_deg[0], phi_deg[-1]]
100  plt.imshow(PAT_db, origin='lower', extent=extent, aspect='auto')
101  plt.colorbar(label='Normalized Gain (dB)')
102  plt.xlabel('Theta (deg)')
103  plt.ylabel('Phi (deg)')
104  plt.title('Array Pattern Heatmap (|AF·EF|, dB) — Kaiser ~-25 dB')
105  plt.tight_layout()
106  plt.savefig('Heatmap_Kaiser25dB_like.png', bbox_inches='tight')
107  plt.show()
108  
109  print('Saved: E_plane_Kaiser25dB_like.png, H_plane_Kaiser25dB_like.png, Heatmap_Kaiser25dB_like.png')