import numpy as np
import panel as pn
import plotly.graph_objects as go
from scipy.special import wofz
from bokeh.models.formatters import BasicTickFormatter
# Ensure Panel extensions are loaded
pn.extension("plotly", "mathjax")
# Physical constants
c = 2.998e10 # Speed of light in cm/s
k_B = 1.38e-16 # Boltzmann constant in erg/K
h = 6.626e-27 # Planck's constant in erg·s
def voigt_profile(nu, nu_0, A21, T, m, v_shift):
"""Compute the Voigt profile given the input parameters."""
delta_nu_D = (nu_0 / c) * np.sqrt(2 * k_B * T / m) # Doppler broadening
gamma_L = A21 / (4 * np.pi) # Lorentzian broadening
x = (nu - nu_0 - v_shift * nu_0 / c) / delta_nu_D
a = gamma_L / delta_nu_D
phi_V = np.real(wofz(x + 1j * a)) / (delta_nu_D * np.sqrt(np.pi))
return phi_V
# Formatting for numerical display
fs = BasicTickFormatter(use_scientific=True, precision=4)
# Panel widgets
lambda_0_slider = pn.widgets.FloatInput(name="λ₀ (Å)", value=1215.67, step=1, format=fs)
lambda_delta_slider = pn.widgets.FloatInput(name="plot 𝛅λ = (λ-λ₀)/λ₀ ", value=1e-3, step=1, format=fs)
unit_selector = pn.widgets.RadioButtonGroup(name="X-Axis", options=["THz", "Angstrom"], value="Angstrom")
A21_slider = pn.widgets.FloatInput(name="A₂₁ (s⁻¹)", value=4.69e8, step=1, format=fs)
T_slider = pn.widgets.FloatInput(name="T (K)", value=10000, step=100, format=fs)
m_slider = pn.widgets.FloatInput(name="Mass (g)", value=1.67e-24, step=1e-16, format=fs)
v_shift_slider = pn.widgets.FloatInput(name="Velocity (km/s)", value=0, step=1)
N_slider = pn.widgets.FloatInput(name="Column Density N (cm⁻²)", value=1e18, step=1e12, format=fs)
# Create empty figures for interactive updates
fig_voigt = go.Figure()
fig_intensity = go.Figure()
# Create a Panel Markdown pane for text updates
text_pane_dynamic = pn.pane.Markdown("### Thermal line width:\n- **Loading...**", width=400)
# Panel Plotly panes
plot_voigt_pane = pn.pane.Plotly(fig_voigt, config={'responsive': True})
plot_intensity_pane = pn.pane.Plotly(fig_intensity, config={'responsive': True})
# Function to update text dynamically
def update_text(dnu, nu_0):
return fr"""
### Thermal Line Width:
- Doppler Width: **{dnu:.3e} Hz**
- Converted to Angstrom: **{(c*dnu/(nu_0*(nu_0+dnu))) * 1e8:.4g} Å**
- $$\delta \lambda/\lambda \sim $$ {(dnu/((nu_0+dnu))) * 1e8:.4g}
"""
def update_plot(event):
"""Update Plotly plots based on widget values."""
# Get current widget values
lambda_0 = lambda_0_slider.value
lambda_delta = lambda_delta_slider.value
unit = unit_selector.value
A21 = A21_slider.value
T = T_slider.value
m = m_slider.value
v_shift = v_shift_slider.value * 1e5 # Convert from km/s to cm/s
N = N_slider.value
# Convert wavelength (Å) to frequency (Hz)
nu_0 = c / (lambda_0 * 1e-8)
# Define frequency range
nu = np.linspace(nu_0 *(1 - lambda_delta), nu_0 *(1 + lambda_delta), 500)
# Compute Voigt profile
phi_V = voigt_profile(nu, nu_0, A21, T, m, v_shift)
# Compute absorption coefficient
g2og1 = 1 # Assume g2/g1 = 2 for simplicity
B12 = (c**2 / (2 * h * nu_0**3)) * (A21 * g2og1)
alpha_nu = (h * nu_0 / (4 * np.pi)) * B12 * phi_V
tau = N * alpha_nu # Optical depth
I = np.exp(-tau) # Intensity
# Convert x-axis based on selection
if unit == "Angstrom":
x_values = c / nu * 1e8 # Convert Hz to Å
xlabel = "Wavelength (Å)"
else:
x_values = nu / 1e12 # Convert Hz to THz
xlabel = "Frequency (THz)"
# Update Voigt profile figure
fig_voigt.data = []
fig_voigt.add_trace(go.Scatter(x=x_values, y=phi_V, mode='lines', name="Voigt Profile", line=dict(color='blue', shape='spline', smoothing=0.7)))
fig_voigt.update_layout(title="Absorption Line Profile", xaxis_title=xlabel, yaxis_title=r'line profile', template="plotly_white")
# Update Intensity figure
fig_intensity.data = []
fig_intensity.add_trace(go.Scatter(x=x_values, y=I, mode='lines', name="Intensity", line=dict(color='red', shape='spline', smoothing=0.7)))
fig_intensity.update_layout(title="Intensity Profile", xaxis_title=xlabel, yaxis_title="Intensity", template="plotly_white")
# Update Panel plot panes
plot_voigt_pane.object = fig_voigt
plot_intensity_pane.object = fig_intensity
# Add vertical line at λ₀
fig_voigt.add_shape(
type="line",
x0=lambda_0, x1=lambda_0,
y0=0, y1=1,
xref="x", yref="paper",
line=dict(color="black", width=2, dash="dot")
)
dnu = (nu_0 / c) * np.sqrt(2 * k_B * T / m) # Doppler broadening
print(dnu*c/nu_0**2 * 1e8)
# ✅ Update the text dynamically
text_pane_dynamic.object = update_text(dnu,nu_0)
# Attach widget callbacks
lambda_0_slider.param.watch(update_plot, 'value')
lambda_delta_slider.param.watch(update_plot, 'value')
unit_selector.param.watch(update_plot, 'value')
A21_slider.param.watch(update_plot, 'value')
T_slider.param.watch(update_plot, 'value')
m_slider.param.watch(update_plot, 'value')
v_shift_slider.param.watch(update_plot, 'value')
N_slider.param.watch(update_plot, 'value')
# Initialize the plot with default values
update_plot(None)
# Create Panel layout
dashboard = pn.Column(
pn.Row(lambda_0_slider, lambda_delta_slider, unit_selector),
pn.Row(A21_slider, T_slider),
pn.Row(m_slider, v_shift_slider),
pn.Row(N_slider,text_pane_dynamic),
pn.Row(plot_voigt_pane, plot_intensity_pane)
)
# Serve the dashboard
dashboard.servable()