Galactic power simulation¶

[1]:

import os
import sys
from radiocalibrationtoolkit import *

[INFO] LFmap: Import successful.

[2]:

# 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,
    )
)

[3]:

# 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.

[4]:

LATITUDE = -35.206667

[5]:

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

lst_range = np.asarray(list(range(24))) + 0.5
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'

[6]:

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,
    },
)

[INFO] No impedance function, assuming it is equal to '1'

100%|███████████████████████████████████████████| 51/51 [03:12<00:00,  3.77s/it]

[7]:

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

[8]:

fig = px.imshow(
    power_DF.T * 1e12, width=600, aspect="cube", color_continuous_scale="jet"
)
fig.update_layout(**layout_settings)
fig.update_layout(
    title="<b>Simulated power dataset at ant. level (Z=1)</b>",
)

fig.show()

[9]:

# Known impedance and HW response

[10]:

hw_file_path = "./antenna_setup_files/HardwareProfileList_realistic.xml"
hw_dict = read_hw_file(hw_file_path, interp_args={'fill_value':'extrapolate'})

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

[11]:

hw_dict

[11]:

defaultdict(<function radiocalibrationtoolkit.common.helpfunctions.read_hw_file.<locals>.<lambda>()>,
            {'RResponse': defaultdict(<function radiocalibrationtoolkit.common.helpfunctions.read_hw_file.<locals>.<lambda>()>,
                         {'LNA': <scipy.interpolate._interpolate.interp1d at 0x7f95c21b2160>,
                          'cable_fromLNA2digitizer': <scipy.interpolate._interpolate.interp1d at 0x7f95c0d4e840>,
                          'digitizer': <scipy.interpolate._interpolate.interp1d at 0x7f95c0d4d080>,
                          'impedance_matching_EW': <scipy.interpolate._interpolate.interp1d at 0x7f95c0d4c400>,
                          'impedance_matching_NS': <scipy.interpolate._interpolate.interp1d at 0x7f95c0d4d030>}),
             'IImpedance': defaultdict(<function radiocalibrationtoolkit.common.helpfunctions.read_hw_file.<locals>.<lambda>()>,
                         {'antenna_EW': <scipy.interpolate._interpolate.interp1d at 0x7f95c0d4cea0>,
                          'antenna_NS': <scipy.interpolate._interpolate.interp1d at 0x7f95c0d4d6c0>})})

[12]:

lst_range = np.asarray(list(range(24))) + 0.5
freq_Mhz_range = range(30, 81, 1)
impedance_func = hw_dict["IImpedance"][
    "antenna_EW"
]

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,
)

100%|███████████████████████████████████████████| 51/51 [03:07<00:00,  3.68s/it]

[13]:

power_DF = integrate_spectral_density(power_density_DF, integrated_MHz_bands=power_density_DF.columns.values)

[14]:

fig = px.imshow(
    power_DF.T * 1e12, width=600, aspect="cube", color_continuous_scale="jet"
)
fig.update_layout(**layout_settings)
fig.update_layout(
    title="<b>Simulated power dataset at antenna level</b>",
    coloraxis=dict(
        colorbar=dict(
            title=dict(
            text="\u03C1",
            side="top",
            )
        ),
        cmin=0,  # Minimum color value
        cmax=0.1,  # Maximum color value
    )
)
fig.show()

[15]:

hw_reponse_1 = dB2PowerAmp(
    hw_dict["RResponse"]["LNA"](power_DF.columns)
)
hw_reponse_2 = dB2PowerAmp(
    hw_dict["RResponse"]["digitizer"](power_DF.columns)
)
hw_reponse_3 = dB2PowerAmp(
    hw_dict["RResponse"]["cable_fromLNA2digitizer"](power_DF.columns)
)
hw_reponse_4 = dB2PowerAmp(
    hw_dict["RResponse"][
        "impedance_matching_EW"
    ](power_DF.columns)
)

power_in_HW_DF = power_DF.multiply(
    hw_reponse_1 * hw_reponse_2 * hw_reponse_3 * hw_reponse_4
)

[16]:

# power_in_HW_DF.to_csv('power_GSM16.csv')

[17]:

fig = px.imshow(
    power_in_HW_DF.T * 1e12, width=600, aspect="cube", color_continuous_scale="jet"
)
fig.update_layout(**layout_settings)
fig.update_layout(
    title="<b>Simulated power dataset at digitzer level</b>",
)
fig.show()

[18]:

# for just the voltage square spectral density
# do not use impedance and do not integrate!

lst_range = range(24)
freq_Mhz_range = range(10, 125, 1)

voltage2_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,
    },
)
voltage2_density_DF.index.name = "LST"
voltage2_density_DF.columns.name = "freq_MHz"

[INFO] No impedance function, assuming it is equal to '1'

100%|█████████████████████████████████████████| 115/115 [06:50<00:00,  3.57s/it]

[19]:

hw_reponse_1 = dB2PowerAmp(hw_dict["RResponse"]["LNA"](voltage2_density_DF.columns))
hw_reponse_2 = dB2PowerAmp(
    hw_dict["RResponse"]["digitizer"](voltage2_density_DF.columns)
)
hw_reponse_3 = dB2PowerAmp(
    hw_dict["RResponse"]["cable_fromLNA2digitizer"](voltage2_density_DF.columns)
)
hw_reponse_4 = dB2PowerAmp(
    hw_dict["RResponse"]["impedance_matching_EW"](voltage2_density_DF.columns)
)

voltage2_density_in_HW_DF = voltage2_density_DF.multiply(
    hw_reponse_1 * hw_reponse_2 * hw_reponse_3 * hw_reponse_4
)

[20]:

# voltage2_density_DF.to_csv('./voltage2_density/voltage2_density_gsm16.csv')
# voltage2_density_in_HW_DF.to_csv('./voltage2_density/voltage2_density_in_HW_gsm16.csv')