Coverage for middlelayer/ngect/engine.py: 30.07%

1# -*- coding: utf-8 -*-

4# This file is part of ngVLA Exposure Calculator Tool (ngECT).

6# ngECT is free software: you can redistribute it and/or modify

7# it under the terms of the GNU General Public License as published by

8# the Free Software Foundation, either version 3 of the License, or

9# any later version.

10#

11# ngECT is distributed in the hope that it will be useful,

12# but WITHOUT ANY WARRANTY; without even the implied warranty of

13# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

14# GNU General Public License for more details.

15#

16# You should have received a copy of the GNU General Public License

17# along with ngECT. If not, see <https://www.gnu.org/licenses/>.

19"""Performance"""

20import logging

22import numpy as np

23import scipy.constants as const

24from astropy import units as u

25from scipy.interpolate import BSpline, CubicSpline, PPoly

27from .system import default_subarray_revision, read_receiver_data, read_subarray_data

28from .utils import get_si_prefix

30logger = logging.getLogger(__name__)

33class ArrayPerformance:

34 """Compute array performance metrics for ngVLA.

36 This class loads array configuration and receiver specifications and

37 provides helpers to map requested observing setups to supported

38 receiver bands and to compute sensitivity/observing time metrics.

40 Attributes

41 ----------

42 subarray : list[str]

43 Available subarray names loaded from configuration.

44 receiver : list[str]

45 Available receiver band names loaded from receiver data.

46 _ac_dt : dict

47 Subarray configuration dataset for the selected revision.

48 _rx_dt : dict

49 Receiver configuration dataset (latest by default).

50 _co_dt : dict

51 Correlator-related constants used in calculations (e.g., ``eta_c``, ``eta_Q``).

52 """

54 def __init__(self, ac_rev: str | None = None):

55 """Initialize the performance engine.

57 Parameters

58 ----------

59 ac_rev : str, optional

60 Subarray configuration revision to load. If not provided, the

61 repository's default revision is used (see

62 ``system.default_subarray_revision()``).

63 """

64 self._ac_rev = ac_rev if ac_rev else default_subarray_revision()

65 self._rx_dt = read_receiver_data()

66 self._ac_dt = read_subarray_data(self._ac_rev)

68 # eta_c = 0.98 # correlator efficiency

69 # eta_Q = 0.9625 # Digitizer quantization efficiency

70 self._co_dt = {"eta_c": 0.98, "eta_Q": 0.9625}

72 # ant_size (hardcoded)

73 for subarray in self._ac_dt:

74 if subarray == "sba":

75 self._ac_dt[subarray]["D"] = 6 * u.m

76 else:

77 self._ac_dt[subarray]["D"] = 18 * u.m

79 self.subarray = list(self._ac_dt.keys())

80 self.receiver = list(self._rx_dt.keys())

82 def freq_to_receiver(self, freq):

83 """Return receiver bands that support a requested frequency.

85 Parameters

86 ----------

87 freq : astropy.units.Quantity

88 Requested sky frequency (e.g., ``10 * u.GHz``). May also be a

89 plain float interpreted as GHz, but a Quantity is recommended.

91 Returns

92 -------

93 list[str]

94 List of receiver band names whose frequency coverage includes

95 the requested value. Multiple bands may include a given frequency

96 where they touch/overlap.

98 Raises

99 ------

100 ValueError

101 If the requested frequency is outside the range supported by all

102 receivers.

103 """

104

105 valid_rx = []

106 for rx in self.receiver:

107 band_freqs = self._rx_dt[rx]["freq"]

108 if freq <= np.max(band_freqs) * u.GHz and freq >= np.min(band_freqs) * u.GHz:

109 valid_rx.append(rx)

110

111 if not valid_rx:

112 error = "The input frequency {0} is outside of the range of the frequencies supported by the ngVLA".format(

113 freq

114 )

115 logger.error(error)

116 raise ValueError(error)

117

118 return valid_rx

119

120 def nearest_supported_frequency(self, freq: u.Quantity) -> dict:

121 """

122 Find the nearest supported receiver band edge to a requested frequency.

123

124 Parameters

125 ----------

126 freq : astropy.units.Quantity

127 Requested sky frequency. Must be a quantity with frequency units

128 (e.g., ``10 * u.GHz``).

129

130 Returns

131 -------

132 dict

133 A mapping with the following keys:

134

135 - ``unit``: always ``"GHz"``.

136 - ``band``: receiver band name (e.g., ``"BAND_1"``).

137 - ``freq``: nearest band-edge frequency as a float in GHz.

138 - ``range``: a mapping with keys ``low`` and ``high`` giving the band

139 limits in GHz.

140

141 Notes

142 -----

143 The search considers both low and high edges of every band and returns the

144 one closest to the requested frequency.

145

146 Examples

147 --------

148 >>> from ngect.engine import ArrayPerformance

149 >>> from astropy import units as u

150 >>> ap = ArrayPerformance()

151 >>> ap.nearest_supported_frequency(0.5 * u.GHz)["freq"]

152 1.2

153 """

154 f_req_ghz = freq.to_value(u.GHz)

155 nearest = None

156

157 for rx in self.receiver:

158 band_freqs = self._rx_dt[rx]["freq"] # array in GHz

159 f_low = float(np.min(band_freqs))

160 f_high = float(np.max(band_freqs))

161

162 # Check both band edges

163 for edge in (f_low, f_high):

164 distance = abs(f_req_ghz - edge)

165 if nearest is None or distance < nearest["dist"]:

166 nearest = {

167 "rx": rx,

168 "edge": edge,

169 "low": f_low,

170 "high": f_high,

171 "dist": distance,

172 }

173

174 # Return the nearest frequency found across all receivers

175 return {

176 "unit": "GHz",

177 "band": nearest["rx"],

178 "freq": nearest["edge"],

179 "range": {"low": nearest["low"], "high": nearest["high"]},

180 }

181

182 def get_receiver(self, rx):

183 """Get receiver properties.

184

185 Notes

186 -----

187 This method is a placeholder and may be implemented in a future

188 revision to return a structured view of receiver parameters.

189 """

190 return

191

192 def _validate_input(self, input):

193 """Validate input mapping for exposure calculations.

194

195 Parameters

196 ----------

197 input : dict

198 Input mapping expected to contain at least the following keys:

199 ``freq`` (Quantity), ``theta`` (Quantity), ``subarray`` (str),

200 ``specmode`` ("cont" or "line"), ``chanwidth`` (Quantity),

201 ``npol`` (int), ``tos`` (Quantity), ``elevation`` (Quantity),

202 ``ectmode`` (str), and ``pwv`` (str).

203

204 Notes

205 -----

206 This method currently acts as a placeholder and should be extended to

207 perform thorough validation and raise informative exceptions on

208 malformed inputs.

209 """

210

211 return

212

213 def ect(self, input):

214 """Compute exposure time/sensitivity metrics for a given setup.

215

216 Parameters

217 ----------

218 input : dict

219 Input parameters containing at least the following keys:

220 - ``freq`` (astropy.units.Quantity): requested observing frequency.

221 - ``theta`` (astropy.units.Quantity): target angular resolution.

222 - ``subarray`` (str): subarray name.

223 - ``specmode`` (str): "cont" or "line".

224 - ``chanwidth`` (astropy.units.Quantity): channel width (velocity or frequency).

225 - ``npol`` (int): number of polarizations.

226 - ``tos`` (astropy.units.Quantity): on-source integration time.

227 - ``elevation`` (astropy.units.Quantity): elevation angle.

228 - ``ectmode`` (str): calculator mode (e.g., "tos2rms").

229 - ``pwv`` (str): precipitable water vapor setting.

230

231 Returns

232 -------

233 dict

234 Mapping with computed quantities such as:

235 ``freq``, ``rx_name``, ``fov``, ``max_bw``, ``t_sys``,

236 ``chanwidth_freq``, ``chanwidth_velo``, ``tos``, ``sigma_fl``,

237 ``sigma_tb``, ``num_ant``, ``theta``, ``eta_A``, ``total_eff_A``,

238 and auxiliary efficiency factors.

239

240 Raises

241 ------

242 ValueError

243 If the requested frequency is unsupported by any receiver band.

244

245 Notes

246 -----

247 This method performs both instrument configuration lookups and several

248 numerical interpolations (e.g., cubic splines) over band-dependent

249 system temperature and aperture efficiency curves. The implementation is

250 adapted from the reference ngVLA sensitivity calculator.

251 """

252

253 self._validate_input(input)

254

255 freq = input["freq"]

256 chanwidth = input["chanwidth"]

257 theta = input["theta"]

258 tos = input["tos"]

259

260 rx = self.freq_to_receiver(freq)

261

262 rx_name = rx[0]

263

264 ar_specs = self._ac_dt[input["subarray"]]

265 rx_specs = self._rx_dt[rx_name]

266

267 ##########################################################

268 # !! need to refactor the following code block

269 ##########################################################

270

271 max_bw = rx_specs["max_bw"] * u.GHz

272 band_freqs = rx_specs["freq"]

273 fC = rx_specs["freq_center"]

274 all_T_s = rx_specs["tSys"]

275 all_eta_A = rx_specs["antEff"]

276

277 fL, fH = np.min(band_freqs), np.max(band_freqs)

278 cs_T_s = CubicSpline(band_freqs, all_T_s)

279 cs_eta_A = CubicSpline(band_freqs, all_eta_A)

280

281 b_max = ar_specs["b_max"]

282 b_min = ar_specs["b_min"]

283 N_ant = ar_specs["N_ant"]

284

285 # nu is the observed frequency in GHz, which could be different from the requested (input) frequency

286 # see the potential modification below in the 'continuum' mode.

287 nu = freq.to_value(u.GHz)

288

289 if input["specmode"] == "line":

290 T_s = cs_T_s(nu)

291 eta_A = cs_eta_A(nu)

292 else:

293 if rx_name != "BAND_6":

294 x_freqs = np.linspace(fL, fH, 100)

295 nu = np.mean([fL, fH])

296 logger.debug(

297 "using entire continuum bandwidth ({1}-{2} GHz) at center frequency {0} GHz".format(nu, fL, fH)

298 )

299 else:

300 if nu - 10 < fL:

301 logger.debug(

302 "continuum bandwidth extends beyond receiver edges, shifting center frequency from {0} to {1}".format(

303 nu, fL + 10

304 )

305 )

306 nu = fL + 10

307 logger.debug(

308 "using continuum bandwidth ({1}-{2} GHz) at center frequency {0} GHz".format(

309 nu, nu - 10, nu + 10

310 )

311 )

312 elif nu + 10 > fH:

313 logger.debug(

314 "continuum bandwidth extends beyond receiver edges, shifting center frequency from {0} to {1}".format(

315 nu, fH - 10

316 )

317 )

318 nu = fH - 10

319 logger.debug(

320 "using entire continuum bandwidth ({1}-{2} GHz) at center frequency {0} GHz".format(

321 nu, nu - 10, nu + 10

322 )

323 )

324 x_freqs = np.linspace(nu - 10, nu + 10, 100)

325

326 T_s = np.mean(cs_T_s(x_freqs))

327 eta_A = np.mean(cs_eta_A(x_freqs))

328

329 logger.debug(

330 "using interpolated Tsys: {0} K at frequency {2} GHz (band average: {1})".format(T_s, rx_specs["Tsys"], nu)

331 )

332 logger.debug(

333 "using interpolated eta_A: {0} at frequency {2} GHz (band average: {1})".format(

334 eta_A, rx_specs["eta_A"], nu

335 )

336 )

337

338 A = np.pi * (ar_specs["D"] / 2.0) ** 2 # [m**2]

339 Field_view = (

340 1.02 * (206265 / 60.0) * const.c / (ar_specs["D"].to_value(u.m) * nu * 1e9)

341 ) # arcmin Uniform illumination

342 total_eff_A = eta_A * N_ant * A

343 Res_max_base = (206265.0) * const.c * 1e3 / (nu * 1e9 * b_max) * u.mas

344 LAS = (206265.0) * const.c / (nu * 1e9 * b_min) * u.arcsec

345

346 SEFD = (2 * const.k * T_s / (self._co_dt["eta_Q"] * eta_A * A.to_value(u.m * u.m))) / 1e-26 * u.Jy

347

348 avg_SEFD = SEFD.to_value(u.Jy) * (rx_specs["Tsys"] / T_s) * (rx_specs["eta_A"] / eta_A)

349

350 logger.debug(

351 "calculated SEFD: {0} Jy at frequency {1} GHz (band average: {2})".format(SEFD.to_value(u.Jy), nu, avg_SEFD)

352 )

353

354 if chanwidth.unit.physical_type == "speed":

355 delta_nu = chanwidth.to_value(u.m / u.s) / const.c * nu * 1e9 * u.Hz # line width in Hz

356 delta_v = chanwidth

357 else:

358 delta_nu = chanwidth

359 delta_v = chanwidth.to_value(u.Hz) / nu * 1e9 * const.c * u.m / u.s

360

361 # delta_nu: in Hz

362 if input["specmode"] == "cont":

363 delta_nu = max_bw

364 elif delta_nu > max_bw:

365 # ensure that we constrain to the maximum bandwidth

366 delta_nu = max_bw

367 delta_v = chanwidth.to_value(u.Hz) / delta_nu * 1e9 * const.c * u.m / u.s

368

369 # emit a message about it

370 logger.debug(

371 "resetting delta_nu (currently %s) to max_bw %s since it exceeds the maximum", delta_nu, max_bw

372 )

373

374 # We assume t_int=1000s as a fiducial value to generate sigma_rms/sigma_T

375 # and the reported values would be scaled from these values.

376 t_int_fiducial = 1000.0 # sec

377 sigma_ps = (

378 SEFD.to_value(u.Jy)

379 / (

380 self._co_dt["eta_c"]

381 * np.sqrt(input["npol"] * delta_nu.to_value(u.Hz) * t_int_fiducial * N_ant * (N_ant - 1))

382 )

383 ) / 1e-6 # uJy

384

385 if theta < 0:

386 eta_w = 1.0

387 theta = (ar_specs["theta_nat_30_GHz"] * u.mas) * 30.0 / nu

388 logger.debug("using native resolution {0} mas at frequency {1} GHz".format(theta.to_value(u.mas), nu))

389 else:

390 if ar_specs["spline_type"] == "cubic":

391 cs = PPoly(*ar_specs["spline_params"])

392 elif ar_specs["spline_type"] == "univariate":

393 cs = BSpline(*ar_specs["spline_params"])

394 eta_w = float(cs(np.log10(theta.to_value(u.mas) * nu / 30.0)))

395 logger.debug(

396 "using inefficiency factor eta_w: {0} at frequency {1} GHz and resolution {2} mas".format(

397 eta_w, nu, theta.to_value(u.mas)

398 )

399 )

400

401 sigma_rms = eta_w * sigma_ps # uJy

402 sigma_T = 1.216 * (sigma_rms / (nu**2 * theta.to_value(u.arcsec) ** 2)) # K when sigma_ps in uJy

403

404 if input["ectmode"] == "tos2rms":

405 tos_scale = tos.to_value(u.s) / t_int_fiducial

406 rms_scale = 1 / tos_scale**0.5

407 if input["ectmode"] == "rms2tos":

408 if input["sigma"].unit.physical_type == "temperature":

409 rms_scale = input["sigma"].to_value(u.K) / sigma_T

410 if (input["sigma"].unit * u.beam).physical_type == "spectral flux density":

411 rms_scale = input["sigma"].to_value(u.uJy / u.beam) / sigma_rms

412 tos_scale = 1 / rms_scale**2.0

413

414 sigma_fl_value = sigma_rms * rms_scale / 1e6 # value in units of Jy/beam

415 sigma_fl_prefix, sigma_fl_scale = get_si_prefix(sigma_fl_value, select="num ", lztol=1)

416

417 sigma_tb_value = sigma_T * rms_scale # value in units of K

418 sigma_tb_prefix, sigma_tb_scale = get_si_prefix(sigma_tb_value, select="m ", lztol=1)

419

420 theta_value = theta.to_value(u.arcsec)

421 theta_prefix, theta_scale = get_si_prefix(theta.to_value(u.arcsec), select="m ", lztol=1)

422 theta_unit = u.Unit(theta_prefix + "arcsec")

423 if theta_scale == "m":

424 theta_unit = u.Unit("mas")

425 else:

426 theta_unit = u.Unit(theta_prefix + "arcsec")

427

428 output = {

429 "freq": nu * u.GHz,

430 "sigma_fl": sigma_fl_value / sigma_fl_scale * u.Unit(sigma_fl_prefix + "Jy/beam"),

431 "sigma_tb": sigma_tb_value / sigma_tb_scale * u.Unit(sigma_tb_prefix + "K"),

432 "tos": t_int_fiducial * u.s * tos_scale,

433 "rx_name": rx_name.replace("_", " "),

434 "rx_freq_min": fL * u.GHz,

435 "rx_freq_max": fH * u.GHz,

436 "max_bw": max_bw,

437 "fov": Field_view * u.arcmin,

438 "eta_A": eta_A,

439 "total_eff_A": total_eff_A,

440 "t_sys": T_s * u.K,

441 "sefd": SEFD,

442 "theta_max": Res_max_base,

443 "theta": theta_value / theta_scale * theta_unit,

444 "chanwidth_freq": delta_nu,

445 "chanwidth_velo": delta_v,

446 "eta_W": eta_w,

447 "LAS": LAS,

448 "num_ant": N_ant,

449 }

450

451 return output

Coverage for middlelayer / ngect / engine.py: 30.07%

143 statements