Crab Nebula example for low zenith range

from datetime import datetime, timezone
from pathlib import Path


from astroplan import Observer, FixedTarget
from astropy.time import Time
import astropy.units as u
from astropy.visualization import quantity_support
import matplotlib.pyplot as plt
from ipywidgets import interact, DatetimePicker

from iact_estimator import RESOURCES_PATH
from iact_estimator.io import read_yaml
from iact_estimator.core import (

    initialize_model,
    prepare_data,
    source_detection,
    calculate,
)
from iact_estimator.plots import plot_spectrum, plot_sed, plot_transit, plot_altitude_airmass

output_path=Path.cwd()
config = read_yaml(RESOURCES_PATH/"config.yml")
source_name = "Crab"

observer = Observer.at_site("Roque de los Muchachos")

crab = FixedTarget.from_name(source_name)

plotting_options = config["plotting_options"]
use_seaborn = config["use_seaborn"]
if use_seaborn:
    import seaborn as sns

    seaborn_options = config["seaborn_options"]
    sns.set_theme(**seaborn_options)


assumed_spectrum = initialize_model(config)

plot_energy_bounds = [
    u.Quantity(plotting_options["min_energy"]),
    u.Quantity(plotting_options["max_energy"]),
]

Source transit

target_source = FixedTarget.from_name(source_name)
observer = Observer.at_site("Roque de los Muchachos")

date_time = DatetimePicker(
    value=datetime.now(timezone.utc),
    description='Select a datetime',
    disabled=False
)

crab = FixedTarget.from_name("Crab")
plot_crab = True if (crab.coord == target_source.coord) else False

def interactive_plot_transit(date_time):
        with quantity_support():
                plot_transit(
                        config,
                        source_name,
                        target_source,
                        observer,
                        time = Time(date_time).utc,
                        merge_profiles=True,
                        plot_crab=False,
                        savefig=False,
                )

interact(interactive_plot_transit, date_time=date_time)
plt.show()

Altitude and airmass

date_time = DatetimePicker(
    value=datetime.now(timezone.utc),
    description='Select a datetime',
    disabled=False
)

def plot_alt(date_time):

    print(date_time)

    plot_altitude_airmass(
                config,
                source_name,
                target_source,
                observer,
                time=Time(date_time).utc,
                brightness_shading=True,
                airmass_yaxis=True,
                savefig=False,
            )

interact(plot_alt, date_time=date_time)
plt.show()

Spectrum

plot_spectrum(
            config,
        plot_energy_bounds,
        assumed_spectrum,
        source_name,
        plotting_options,
        savefig=False,
        )

Spectral energy distribution

energy_bins, gamma_rate, background_rate = prepare_data(config)

en, sed, dsed, sigmas, detected = calculate(
    energy_bins, gamma_rate, background_rate, config, assumed_spectrum
)

combined_significance = source_detection(
    sigmas, u.Quantity(config["observation_time"])
)

There are 14 energy bins
39.8 GeV-63.0 GeV GeV: exc. = 2455.3+-130.8 ev., SBR=22.34%, sigma = 19.4  DETECTION
63.0 GeV-100.0 GeV GeV: exc. = 9030.0+-159.0 ev., SBR=74.14%, sigma = 62.4  DETECTION
100.0 GeV-158.0 GeV GeV: exc. = 12870.0+-150.0 ev., SBR=178.01%, sigma = 101.7  DETECTION
158.0 GeV-251.0 GeV GeV: exc. = 10110.0+-110.8 ev., SBR=624.07%, sigma = 123.5  DETECTION
251.0 GeV-398.0 GeV GeV: exc. = 4080.0+-65.9 ev., SBR=2060.61%, sigma = 92.7  DETECTION
398.0 GeV-631.0 GeV GeV: exc. = 3660.0+-61.4 ev., SBR=4518.52%, sigma = 93.2  DETECTION
631.0 GeV-1000.0 GeV GeV: exc. = 2640.0+-51.9 ev., SBR=6616.54%, sigma = 80.7  DETECTION
1000.0 GeV-1585.0 GeV GeV: exc. = 1740.0+-42.0 ev., SBR=9830.51%, sigma = 66.5  DETECTION
1585.0 GeV-2512.0 GeV GeV: exc. = 900.0+-30.2 ev., SBR=11111.11%, sigma = 48.0  DETECTION
2512.0 GeV-3981.0 GeV GeV: exc. = 498.0+-22.5 ev., SBR=8300.00%, sigma = 35.4  DETECTION
3981.0 GeV-6310.0 GeV GeV: exc. = 279.0+-16.9 ev., SBR=6642.86%, sigma = 26.2  DETECTION
6310.0 GeV-10000.0 GeV GeV: exc. = 180.0+-14.1 ev., SBR=1304.35%, sigma = 18.5  DETECTION
10000.0 GeV-15848.9 GeV GeV: exc. = 53.6+-7.8 ev., SBR=1070.00%, sigma = 9.8  DETECTION
Combined significance (using the 13 data points shown in the SED) = 215.8
The source probably will be detected in 50.0 h.

annotation_options = {"rotation": 45,
        "xytext": (10, 10),
        "size": 15}

with quantity_support():
    plot_sed(
            config,
            sigmas,
            combined_significance,
            source_name,
            assumed_spectrum,
            en,
            sed,
            dsed,
            detected,
            savefig=False,
            annotation_options=annotation_options,
            )
    plt.ylim(1.e-12, 2.e-10)