Compare galactic power simulations via galactic and local CS

import os
import sys
from radiocalibrationtoolkit import *

[INFO] LFmap: Import successful.

# some global plot settings
plt.rcParams["axes.labelweight"] = "bold"
plt.rcParams["font.weight"] = "bold"
plt.rcParams["font.size"] = 16
plt.rcParams["legend.fontsize"] = 14

plt.rcParams["xtick.major.width"] = 2
plt.rcParams["ytick.major.width"] = 2

plt.rcParams["xtick.major.size"] = 5
plt.rcParams["ytick.major.size"] = 5

plt.rcParams["xtick.labelsize"] = 14
plt.rcParams["ytick.labelsize"] = 14

layout_settings = dict(
    xaxis=dict(title="<b>LST</b>", tickprefix="<b>", ticksuffix="</b>", dtick=2),
    yaxis=dict(
        title="<b>frequency [MHz]</b>",
        tickprefix="<b>",
        ticksuffix="</b>",
        range=(30, 80),
        tick0=0,
        dtick=10,
        autorange=False,
    ),
    font=dict(
        size=20,
        color="black",
    ),
    coloraxis=dict(
        colorbar=dict(
            tickprefix="<b>",
            ticksuffix="</b>",
            title=dict(text="<b>Power [pW]</b>", side="right")),
        cmin=0,
        cmax=24,
    )
)

# This ensures Plotly output works in multiple places:
# plotly_mimetype: VS Code notebook UI
# notebook: "Jupyter: Export to HTML" command in VS Code
# See https://plotly.com/python/renderers/#multiple-renderers
import plotly.io as pio
pio.renderers.default = "plotly_mimetype+notebook"

To get the Power we need to just integrate the power spectrum density over the frequency.

\begin{equation*} P_{sky}(t,f) = \int_{f} \mathscr{P}_{sky}(t,f) df \end{equation*}

Thus, the full equation for the Sky Power Radio Emmision is:

\begin{equation*} P_{sky}(t,f) = \frac{k_{b}}{c^{2}} \int_{f} f^{2} \int_{\Omega} T_{sky}(t,f,\theta,\phi) \frac{|H(f,\theta,\phi)|^{2}Z_{0}}{R_{r}} df d\Omega \end{equation*}

where \(Z_0\) is the vacuum impedance, \(R_r\) is the antenna impedance, \(k_b\) the Boltzmann constant and \(c\) is the speed of light.

LATITUDE = -35.206667

antenna_inst = AntennaPattern(
    "./antenna_setup_files/SALLA_EW.xml",
)
galactic_map_inst = GlobalSkyModel2016(freq_unit="MHz")

lst_range = range(24)
freq_Mhz_range = range(30, 81, 1)

[INFO] Your keys are: ['EAHTheta_amp', 'EAHTheta_phase', 'EAHPhi_amp', 'EAHPhi_phase', 'absolute'] Use them as the 'quantity'

# Known impedance and HW response
hw_file_path = "./antenna_setup_files/HardwareProfileList_realistic.xml"
hw_dict = read_hw_file(hw_file_path, interp_args={'fill_value': 'extrapolate'})
lst_range = range(24)
freq_Mhz_range = range(30, 81, 1)
impedance_func = hw_dict["IImpedance"]["antenna_EW"]

<?xml version="1.0" encoding="iso-8859-1"?>
<Element HardwareProfileList at 0x7f3bf1d06e40>

# Calculate power via galactic CS

power_density_DF = calculate_power_spectral_density(
    antenna_inst=antenna_inst,
    galactic_map_inst=galactic_map_inst,
    lst_range=lst_range,
    freq_Mhz_range=freq_Mhz_range,
    latitude=LATITUDE,
    update_antenna_conventions={
        "shift_phi": -90,
        "flip_theta": True,
        "flip_phi": False,
        "in_degrees": True,
        "add_invisible_sky": True,
    },
    impedance_func=impedance_func,
)

power_via_galactic_cs_DF = integrate_spectral_density(
    power_density_DF, integrated_MHz_bands=power_density_DF.columns
)

100%|███████████████████████████████████████████| 51/51 [04:09<00:00,  4.89s/it]

fig = px.imshow(
    power_via_galactic_cs_DF.T * 1e12, width=600, aspect="cube", color_continuous_scale="jet"
)
fig.update_layout(**layout_settings)
fig.update_layout(
    title="<b>Simulated power dataset via galactic CS</b>",
    coloraxis=dict(
        cmin=0,  # Minimum color value
        cmax=0.1,  # Maximum color value
    )
)
fig.show()

# Calculate power via local CS

nside_transform = 128
power_density_DF = calculate_power_spectral_density(
    antenna_inst=antenna_inst,
    galactic_map_inst=galactic_map_inst,
    lst_range=lst_range,
    freq_Mhz_range=freq_Mhz_range,
    latitude=LATITUDE,
    update_antenna_conventions={
        "shift_phi": -90,
        "flip_theta": True,
        "flip_phi": False,
        "in_degrees": True,
        "add_invisible_sky": True,
    },
    impedance_func=impedance_func,
    # perform calculations in the local cs
    transform_hpmaps2_local_cs=True,
    nside_transform=nside_transform,
)

power_via_local_cs_DF = integrate_spectral_density(
    power_density_DF, integrated_MHz_bands=power_density_DF.columns
)

# power_via_local_cs_DF.to_csv('./power_via_local_CS_{}nside.csv'.format(nside_transform))

100%|███████████████████████████████████████████| 51/51 [04:46<00:00,  5.62s/it]

fig = px.imshow(
    power_via_local_cs_DF.T * 1e12, width=600, aspect="cube", color_continuous_scale="jet"
)
fig.update_layout(**layout_settings)
fig.update_layout(
    title="<b>Simulated power dataset via local CS</b>",
    coloraxis=dict(
        cmin=0,  # Minimum color value
        cmax=0.1,  # Maximum color value
    )
)
fig.show()

# Compare

fig = px.imshow(
    (power_via_galactic_cs_DF/power_via_local_cs_DF.values).T -1, width=650, aspect="cube", color_continuous_scale="jet"
)
fig.update_layout(**layout_settings)
fig.update_layout(
    title="<b>Power calculation comparison</b>",
    coloraxis=dict(
        colorbar=dict(title=dict(text="<b>(<sup>via galactic CS</sup>/<sub>via local CS</sub>)-1</b>", side="right",  font=dict(size=30))),
        cmin=None,
        cmax=None,
    )
)

fig.show()

# maximal difference
(np.abs(power_via_galactic_cs_DF/power_via_local_cs_DF.values-1)).values.max()

0.002198708764023616

# average difference
(power_via_galactic_cs_DF/power_via_local_cs_DF.values-1).values.mean()

-0.00019433273892603844