Chapter 8: Matplotlib and Visual Evidence

import numpy as np

# Question 1: What shape will this array have?
x = np.linspace(0, 10, 100)
y = np.sin(x)
# Shape of y: _______

# Question 2: What will this boolean mask select?
data = np.array([1, 5, 3, 8, 2, 9])
mask = data > 5
selected = data[mask]
# selected contains: _______

# Question 3: What's wrong with this code?
# matrix = np.array([[1, 2], [3, 4]])
# result = matrix + [10, 20, 30]  # Will this work?
# Your answer: _______

# Question 4: What will this create?
X, Y = np.meshgrid(np.arange(5), np.arange(3))
# Shape of X and Y: _______

# Create a figure with explicit control over components
fig = plt.figure(figsize=(10, 6))

# Add axes manually to see the structure
# [left, bottom, width, height] in figure coordinates (0-1)
ax_main = fig.add_axes([0.1, 0.3, 0.7, 0.6])   # Main plot
ax_zoom = fig.add_axes([0.65, 0.6, 0.2, 0.2])  # Inset zoom
ax_residual = fig.add_axes([0.1, 0.05, 0.7, 0.2])  # Residual panel

# Generate data
x = np.linspace(0, 10, 100)
y = np.sin(x) * np.exp(-x/10)
y_model = np.sin(x) * np.exp(-x/10.5)  # Slightly different model

# Main plot
ax_main.plot(x, y, 'ko', markersize=3, alpha=0.5, label='Data')
ax_main.plot(x, y_model, 'r-', linewidth=2, label='Model')
ax_main.set_ylabel('Signal', fontsize=12)
ax_main.legend(loc='upper right')
ax_main.grid(True, alpha=0.3)

# Zoom inset
zoom_mask = (x > 3) & (x < 5)
ax_zoom.plot(x[zoom_mask], y[zoom_mask], 'ko', markersize=2)
ax_zoom.plot(x[zoom_mask], y_model[zoom_mask], 'r-', linewidth=1)
ax_zoom.set_xlim(3, 5)
ax_zoom.grid(True, alpha=0.3)
ax_zoom.set_title('Zoom', fontsize=9)

# Residuals
residuals = y - y_model
ax_residual.scatter(x, residuals, s=5, alpha=0.5, color='blue')
ax_residual.axhline(y=0, color='red', linestyle='--', linewidth=1)
ax_residual.set_xlabel('Time', fontsize=12)
ax_residual.set_ylabel('Residuals', fontsize=10)
ax_residual.grid(True, alpha=0.3)

fig.suptitle('Explicit Control Over Figure Components', fontsize=14, fontweight='bold')
plt.show()

from matplotlib.gridspec import GridSpec

# Create figure with shared x-axes using GridSpec
fig = plt.figure(figsize=(10, 8))
gs = GridSpec(3, 1, hspace=0.05)

# Create subplots with shared x-axis
ax1 = fig.add_subplot(gs[0])
ax2 = fig.add_subplot(gs[1], sharex=ax1)
ax3 = fig.add_subplot(gs[2], sharex=ax1)

# Generate multi-wavelength data
time = np.linspace(0, 10, 200)
rng = np.random.default_rng(42)
optical = np.sin(2 * np.pi * time / 3) + rng.normal(0, 0.1, 200)
xray = 2 * np.sin(2 * np.pi * time / 3 + 0.5) + rng.normal(0, 0.2, 200)
radio = 0.5 * np.sin(2 * np.pi * time / 3 - 0.3) + rng.normal(0, 0.05, 200)

# Plot with shared axes
ax1.plot(time, optical, 'b-', alpha=0.7)
ax1.set_ylabel('Optical')
ax1.tick_params(labelbottom=False)  # Hide x labels except bottom

ax2.plot(time, xray, 'r-', alpha=0.7)
ax2.set_ylabel('X-ray')
ax2.tick_params(labelbottom=False)

ax3.plot(time, radio, 'g-', alpha=0.7)
ax3.set_ylabel('Radio')
ax3.set_xlabel('Time (days)')

fig.suptitle('Multi-wavelength Observations with Shared Time Axis')
plt.show()

# Generate different types of relationships
x = np.logspace(-1, 3, 100)  # 0.1 to 1000

# Different mathematical relationships
linear = 2 * x + 5
quadratic = 0.5 * x**2
exponential = np.exp(x/100)
power_law = 10 * x**(-1.5)
logarithmic = 50 * np.log10(x) + 10

# Create a comprehensive comparison
fig, axes = plt.subplots(5, 4, figsize=(14, 16))

datasets = [
    (linear, 'Linear: y = 2x + 5'),
    (quadratic, 'Quadratic: y = 0.5x^2'),
    (exponential, 'Exponential: y = e^(x/100)'),
    (power_law, 'Power Law: y = 10x^(-1.5)'),
    (logarithmic, 'Logarithmic: y = 50log(x) + 10')
]

for i, (data, title) in enumerate(datasets):
    # Linear-linear
    axes[i, 0].plot(x, data, 'b-', linewidth=2)
    axes[i, 0].set_title(f'{title}\nLinear-Linear')
    axes[i, 0].grid(True, alpha=0.3)

    # Log-log: both axes must be positive. Never take abs(data) to force this.
    if np.all(data > 0):
        axes[i, 1].loglog(x, data, 'r-', linewidth=2)
        axes[i, 1].set_title('Log-Log')
        axes[i, 1].grid(True, alpha=0.3, which='both')
    else:
        axes[i, 1].text(0.5, 0.5, 'Invalid: y contains\nzero or negative values',
                        ha='center', va='center', transform=axes[i, 1].transAxes)
        axes[i, 1].set_title('Log-Log: not valid')

    # Semilog-x
    axes[i, 2].semilogx(x, data, 'g-', linewidth=2)
    axes[i, 2].set_title('Semilog-X')
    axes[i, 2].grid(True, alpha=0.3)

    # Semilog-y has the same positive-y requirement.
    if np.all(data > 0):
        axes[i, 3].semilogy(x, data, 'm-', linewidth=2)
        axes[i, 3].set_title('Semilog-Y')
        axes[i, 3].grid(True, alpha=0.3)
    else:
        axes[i, 3].text(0.5, 0.5, 'Invalid: y contains\nzero or negative values',
                        ha='center', va='center', transform=axes[i, 3].transAxes)
        axes[i, 3].set_title('Semilog-Y: not valid')

# Highlight which scale reveals linearity
axes[0, 0].set_facecolor('#E8F4F8')  # Linear shows linear
axes[2, 3].set_facecolor('#F8E8E8')  # Semilogy shows exponential
axes[3, 1].set_facecolor('#F8F8E8')  # Loglog shows power law
axes[4, 2].set_facecolor('#E8F8E8')  # Semilogx shows logarithmic

fig.suptitle('Choose the Scale that Reveals Your Data\'s Nature',
             fontsize=16, fontweight='bold')
plt.tight_layout()
plt.show()

print("Key insight: The 'right' scale makes relationships linear!")
print("- Power laws -> log-log")
print("- Exponential growth -> semilog-y")
print("- Logarithmic growth -> semilog-x")

# Synthetic, shape-only Schechter-style luminosity-function example.
# It illustrates scale choice; it is not a calibrated number-density estimate.
rng = np.random.default_rng(42)

# Schechter function parameters (Schechter 1976)
L_star = 1e10  # Characteristic luminosity
alpha = -1.25  # Faint-end slope

# Generate a relative luminosity-function shape and multiplicative scatter.
L = np.logspace(7, 12, 1000)
relative_phi = (L/L_star)**alpha * np.exp(-L/L_star)

# This is synthetic scatter, not observational uncertainty.
observed_relative_phi = relative_phi * rng.lognormal(0, 0.2, len(relative_phi))

fig, axes = plt.subplots(2, 2, figsize=(10, 10))

# Linear scale - useless for this data
axes[0, 0].plot(L, observed_relative_phi, '.', markersize=1, alpha=0.5)
axes[0, 0].set_xlabel(r'Luminosity [$L_\odot$]')
axes[0, 0].set_ylabel('Relative number density (arbitrary normalization)')
axes[0, 0].set_title('Linear Scale: Cannot See Structure')

# Log-log - reveals power law at faint end
axes[0, 1].loglog(L, observed_relative_phi, '.', markersize=1, alpha=0.5,
                  label='Synthetic sample')
axes[0, 1].loglog(L, relative_phi, 'r-', linewidth=2,
                  label='Schechter-style shape')
axes[0, 1].set_xlabel(r'Luminosity [$L_\odot$]')
axes[0, 1].set_ylabel('Relative number density (arbitrary normalization)')
axes[0, 1].set_title('Log-Log: Reveals Power Law + Exponential Cutoff')
axes[0, 1].legend()

# Semilog-y - emphasizes exponential cutoff
axes[1, 0].semilogx(L, observed_relative_phi, '.', markersize=1, alpha=0.5)
axes[1, 0].semilogx(L, relative_phi, 'r-', linewidth=2)
axes[1, 0].set_xlabel(r'Luminosity [$L_\odot$]')
axes[1, 0].set_ylabel('Relative number density (arbitrary normalization)')
axes[1, 0].set_title('Semilog-X: Emphasizes Bright End Cutoff')

# Custom: log-log with ratio to model
ratio = observed_relative_phi / relative_phi
axes[1, 1].semilogx(L, ratio, '.', markersize=1, alpha=0.5)
axes[1, 1].axhline(y=1, color='red', linestyle='--', linewidth=2)
axes[1, 1].set_xlabel(r'Luminosity [$L_\odot$]')
axes[1, 1].set_ylabel('Observed / Model')
axes[1, 1].set_title('Residual Plot: Shows Systematic Deviations')
axes[1, 1].set_ylim(0.1, 10)

fig.suptitle('Galaxy Luminosity Function: Scale Choice Matters', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# Example: A reusable light curve plotting function
def plot_light_curve(time, flux, flux_err=None, period=None,
                     title=None, figsize=(12, 8)):
    """
    Create a publication-quality light curve plot with optional phase folding.

    Parameters
    ----------
    time : array-like
        Time values (days)
    flux : array-like
        Flux or magnitude values
    flux_err : array-like, optional
        Flux uncertainties
    period : float, optional
        Period for phase folding (days)
    title : str, optional
        Plot title
    figsize : tuple, optional
        Figure size (width, height)

    Returns
    -------
    fig, axes : matplotlib figure and axes objects
    """
    if period is not None:
        fig, axes = plt.subplots(2, 2, figsize=figsize)
        axes = axes.flatten()
    else:
        fig, axes = plt.subplots(1, 1, figsize=(figsize[0], figsize[1]/2))
        axes = [axes]  # Make it iterable

    # Main light curve
    ax = axes[0]
    if flux_err is not None:
        ax.errorbar(time, flux, yerr=flux_err, fmt='k.', markersize=2,
                   alpha=0.5, elinewidth=0.5, capsize=0)
    else:
        ax.scatter(time, flux, s=1, alpha=0.5, color='black')

    ax.set_xlabel('Time (days)', fontsize=11)
    ax.set_ylabel('Relative Flux', fontsize=11)
    ax.set_title('Light Curve' if title is None else title, fontsize=12)
    ax.grid(True, alpha=0.3)

    # If period provided, add phase-folded plot
    if period is not None:
        phase = (time % period) / period

        # Phase folded
        ax = axes[1]
        if flux_err is not None:
            ax.errorbar(phase, flux, yerr=flux_err, fmt='b.', markersize=2,
                       alpha=0.3, elinewidth=0.5, capsize=0)
        else:
            ax.scatter(phase, flux, s=1, alpha=0.3, color='blue')
        ax.set_xlabel('Phase', fontsize=11)
        ax.set_ylabel('Relative Flux', fontsize=11)
        ax.set_title(f'Phase Folded (P = {period:.3f} days)', fontsize=12)
        ax.set_xlim(0, 1)
        ax.grid(True, alpha=0.3)

        # Double-plotted phase folded
        ax = axes[2]
        phase_double = np.concatenate([phase, phase + 1])
        flux_double = np.concatenate([flux, flux])
        if flux_err is not None:
            err_double = np.concatenate([flux_err, flux_err])
            ax.errorbar(phase_double, flux_double, yerr=err_double,
                       fmt='r.', markersize=2, alpha=0.3,
                       elinewidth=0.5, capsize=0)
        else:
            ax.scatter(phase_double, flux_double, s=1, alpha=0.3, color='red')
        ax.set_xlabel('Phase', fontsize=11)
        ax.set_ylabel('Relative Flux', fontsize=11)
        ax.set_title('Double Phase Plot', fontsize=12)
        ax.set_xlim(0, 2)
        ax.grid(True, alpha=0.3)

        # Binned phase curve - more efficient version
        ax = axes[3]
        n_bins = 50
        binned_flux, bin_edges = np.histogram(phase, bins=n_bins, weights=flux)
        bin_counts, _ = np.histogram(phase, bins=n_bins)
        bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

        # Estimate each occupied bin's sample mean and standard error.
        # A one-point bin has no sample scatter estimate, so do not draw it
        # as though it had a measured standard error.
        valid_bins = bin_counts > 1
        binned_mean = np.full(n_bins, np.nan)
        binned_mean[valid_bins] = binned_flux[valid_bins] / bin_counts[valid_bins]

        # The sample variance needs the n - 1 correction before the SEM.
        flux_squared, _ = np.histogram(phase, bins=n_bins, weights=flux**2)
        sample_variance = np.full(n_bins, np.nan)
        sample_variance[valid_bins] = (
            flux_squared[valid_bins]
            - bin_counts[valid_bins] * binned_mean[valid_bins]**2
        ) / (bin_counts[valid_bins] - 1)
        binned_err = np.full(n_bins, np.nan)
        binned_err[valid_bins] = np.sqrt(
            np.maximum(sample_variance[valid_bins] / bin_counts[valid_bins], 0)
        )

        # Only plot bins with a stated uncertainty estimate.
        ax.errorbar(bin_centers[valid_bins], binned_mean[valid_bins],
                   yerr=binned_err[valid_bins],
                   fmt='go-', markersize=4, linewidth=1, capsize=3)
        ax.set_xlabel('Phase', fontsize=11)
        ax.set_ylabel('Relative Flux', fontsize=11)
        ax.set_title(f'Binned ({n_bins} bins)', fontsize=12)
        ax.set_xlim(0, 1)
        ax.grid(True, alpha=0.3)

    plt.tight_layout()
    return fig, axes

# Test the function with synthetic data
rng = np.random.default_rng(42)
time = np.linspace(0, 30, 500)
period = 2.3456  # days
phase = 2 * np.pi * time / period
flux = 1.0 - 0.01 * np.sin(phase)**2  # Transit-like signal
flux += rng.normal(0, 0.002, len(time))  # Add noise
flux_err = np.ones_like(flux) * 0.002

fig, axes = plot_light_curve(time, flux, flux_err, period=period,
                            title='Exoplanet Transit Light Curve')
plt.show()

print("This reusable function encodes best practices:")
print("- Automatic phase folding when period is provided")
print("- Double phase plot to see continuity")
print("- Binned version to see average behavior")
print("- Consistent styling throughout")

# my_plots.py - Your personal plotting library
import numpy as np
import matplotlib.pyplot as plt

def plot_spectrum(wavelength, flux, flux_err=None, lines=None,
                  title=None, figsize=(12, 5)):
    """Plot a stellar spectrum with optional line identification."""
    fig, ax = plt.subplots(figsize=figsize)

    # Plot spectrum
    ax.plot(wavelength, flux, 'k-', linewidth=0.8, label='Spectrum')

    if flux_err is not None:
        ax.fill_between(wavelength, flux - flux_err, flux + flux_err,
                        alpha=0.3, color='gray', label='Uncertainty')

    # Mark spectral lines
    if lines is not None:
        for wave, name in lines:
            ax.axvline(x=wave, color='red', linestyle=':', alpha=0.5)
            ax.text(wave, ax.get_ylim()[1]*0.95, name,
                   rotation=90, va='top', ha='right', fontsize=8)

    ax.set_xlabel('Wavelength (Å)', fontsize=11)
    ax.set_ylabel('Normalized Flux', fontsize=11)
    ax.set_title('Spectrum' if title is None else title, fontsize=12)
    ax.grid(True, alpha=0.3)
    ax.legend(loc='best')

    return fig, ax

def plot_cmd(color, magnitude, title=None, figsize=(8, 10)):
    """Create a color-magnitude diagram."""
    fig, ax = plt.subplots(figsize=figsize)

    # 2D histogram for density
    h = ax.hist2d(color, magnitude, bins=50, cmap='YlOrBr',
                  density=True)

    ax.set_xlabel('Color (B-V)', fontsize=11)
    ax.set_ylabel('Absolute Magnitude', fontsize=11)
    ax.invert_yaxis()  # Astronomical convention
    ax.set_title('Color-Magnitude Diagram' if title is None else title,
                fontsize=12)

    plt.colorbar(h[3], ax=ax, label='Density')

    return fig, ax

def plot_finding_chart(ra, dec, image=None, sources=None,
                       figsize=(10, 10)):
    """Create a finding chart with marked sources."""
    fig, ax = plt.subplots(figsize=figsize,
                           subplot_kw={'projection': 'rectilinear'})

    if image is not None:
        im = ax.imshow(image, cmap='gray_r', origin='lower',
                      extent=[ra.min(), ra.max(), dec.min(), dec.max()])
        plt.colorbar(im, ax=ax, label='Intensity')

    if sources is not None:
        for source in sources:
            circle = plt.Circle((source['ra'], source['dec']),
                               source.get('radius', 0.1),
                               fill=False, color='red', linewidth=2)
            ax.add_patch(circle)
            ax.text(source['ra'], source['dec'] + 0.15,
                   source.get('name', ''),
                   ha='center', color='red')

    ax.set_xlabel('RA (degrees)', fontsize=11)
    ax.set_ylabel('Dec (degrees)', fontsize=11)
    ax.set_title('Finding Chart', fontsize=12)
    ax.invert_xaxis()  # RA increases to the left

    return fig, ax

# Example usage
print("Your personal plotting library is ready!")
print("Available functions:")
print("  - plot_light_curve(): For time series photometry")
print("  - plot_spectrum(): For spectroscopic data")
print("  - plot_cmd(): For color-magnitude diagrams")
print("  - plot_finding_chart(): For sky position plots")

# Synthetic stellar-spectrum teaching model, not a calibrated observation.
rng = np.random.default_rng(42)
wavelength = np.linspace(4000, 7000, 1000)  # Angstroms
continuum = 1 - 0.0001 * (wavelength - 5500)**2 / 1e6  # Continuum shape

# Add toy absorption features at H-alpha, H-beta, and H-gamma.
lines = [(6563, 50, 0.3), (4861, 30, 0.25), (4340, 25, 0.2)]  # center, width, depth
spectrum = continuum.copy()
for center, width, depth in lines:
    spectrum *= (1 - depth * np.exp(-0.5 * ((wavelength - center) / width)**2))

# Add noise
spectrum += rng.normal(0, 0.01, len(wavelength))

# Create publication-quality spectrum plot
fig, ax = plt.subplots(figsize=(12, 5))

ax.plot(wavelength, spectrum, 'k-', linewidth=0.8, label='Synthetic flux')
ax.axhline(y=1, color='gray', linestyle='--', alpha=0.5, label='Continuum')

# Mark important lines
for center, width, depth in lines:
    ax.axvline(x=center, color='red', linestyle=':', alpha=0.5)
    ax.text(center, 0.65, f'{center}Å', rotation=90,
            va='bottom', ha='right', fontsize=9)

ax.set_xlabel('Wavelength (Å)', fontsize=12)
ax.set_ylabel('Normalized Flux', fontsize=12)
ax.set_title('Synthetic Stellar Spectrum with Balmer Lines', fontsize=14)
ax.set_xlim(4000, 7000)
ax.set_ylim(0.6, 1.1)
ax.grid(True, alpha=0.2)
ax.legend(loc='lower right')

plt.tight_layout()
plt.show()

# Demonstrating proper error bar visualization
rng = np.random.default_rng(42)

# Generate data with stated one-standard-deviation measurement uncertainties.
x = np.linspace(1, 10, 20)
y_err = 0.5 + 0.1 * x + rng.uniform(0, 0.5, 20)
y_true = 2 * x + 3
y = y_true + rng.normal(0, y_err)  # Simulated noise matches the displayed y_err

fig, axes = plt.subplots(1, 3, figsize=(14, 5))

# Standard error bars
axes[0].errorbar(x, y, yerr=y_err, fmt='o', capsize=3, capthick=1,
                 elinewidth=1, markersize=5, label='Data')
axes[0].set_title('Standard Error Bars')
axes[0].set_xlabel('X')
axes[0].set_ylabel('Y')
axes[0].grid(True, alpha=0.3)

# Error bars with a fit-derived uncertainty band
axes[1].errorbar(x, y, yerr=y_err, fmt='o', capsize=0,
                 elinewidth=0.5, markersize=5, alpha=0.7)
# `polyfit` expects inverse-standard-deviation weights, not inverse variance.
coeffs, parameter_covariance = np.polyfit(x, y, 1, w=1 / y_err, cov='unscaled')
model_x = np.linspace(0, 11, 100)
model_y = np.polyval(coeffs, model_x)
design_matrix = np.column_stack([model_x, np.ones_like(model_x)])
model_sigma = np.sqrt(
    np.einsum('ij,jk,ik->i', design_matrix, parameter_covariance, design_matrix)
)
axes[1].plot(model_x, model_y, 'r-', label='Weighted linear fit')
axes[1].fill_between(model_x, model_y - model_sigma, model_y + model_sigma,
                     alpha=0.3, color='red',
                     label='1σ parameter uncertainty (not prediction error)')
axes[1].set_title('Fit with Parameter-Uncertainty Band')
axes[1].set_xlabel('X')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# Weighted fit visualization
weights = 1 / y_err
# Fit weighted linear model with the same stated convention.
coeffs = np.polyfit(x, y, 1, w=weights)
fit_y = np.polyval(coeffs, x)

axes[2].errorbar(x, y, yerr=y_err, fmt='none', capsize=3,
                 elinewidth=1, alpha=0.6)
axes[2].scatter(x, y, s=20 + 40 * weights / weights.max(), alpha=0.6,
                label='Data (marker area ∝ 1 / σ)')
axes[2].plot(x, fit_y, 'g-', linewidth=2, label='Weighted fit')
axes[2].set_title('Weighted Fitting')
axes[2].set_xlabel('X')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

fig.suptitle('Error Bar Visualization Best Practices', fontsize=14)
plt.tight_layout()
plt.show()

# Create synthetic galaxy image
rng = np.random.default_rng(42)
x = np.linspace(-5, 5, 200)
y = np.linspace(-5, 5, 200)
X, Y = np.meshgrid(x, y)

# Exponential disk profile
R = np.sqrt(X**2 + Y**2)
disk = np.exp(-R / 1.5)

# Add spiral arms (simplified)
theta = np.arctan2(Y, X)
spiral = 1 + 0.3 * np.sin(2 * theta - R)
galaxy = disk * spiral

# Add noise and background
galaxy += rng.normal(0, 0.02, galaxy.shape)
galaxy += 0.1  # Sky background

# Display with different scaling
fig, axes = plt.subplots(2, 3, figsize=(12, 8))

# Linear scale
im1 = axes[0, 0].imshow(galaxy, cmap='gray', origin='lower')
axes[0, 0].set_title('Linear Scale')
plt.colorbar(im1, ax=axes[0, 0], fraction=0.046)

# Log scale
from matplotlib.colors import LogNorm
galaxy_positive = galaxy - galaxy.min() + 1e-3  # Ensure positive
im2 = axes[0, 1].imshow(galaxy_positive, cmap='gray',
                         norm=LogNorm(), origin='lower')
axes[0, 1].set_title('Log Scale')
plt.colorbar(im2, ax=axes[0, 1], fraction=0.046)

# Histogram equalization (adaptive)
from matplotlib.colors import PowerNorm
im3 = axes[0, 2].imshow(galaxy, cmap='gray',
                         norm=PowerNorm(gamma=0.5), origin='lower')
axes[0, 2].set_title('Power Scale (gamma=0.5)')
plt.colorbar(im3, ax=axes[0, 2], fraction=0.046)

# Different colormaps
im4 = axes[1, 0].imshow(galaxy, cmap='viridis', origin='lower')
axes[1, 0].set_title('Viridis Colormap')
plt.colorbar(im4, ax=axes[1, 0], fraction=0.046)

im5 = axes[1, 1].imshow(galaxy, cmap='hot', origin='lower')
axes[1, 1].set_title('Hot Colormap')
plt.colorbar(im5, ax=axes[1, 1], fraction=0.046)

im6 = axes[1, 2].imshow(galaxy, cmap='cividis', origin='lower')
axes[1, 2].set_title('Cividis Colormap')
plt.colorbar(im6, ax=axes[1, 2], fraction=0.046)

for ax in axes.flat:
    ax.set_xlabel('X (pixels)')
    ax.set_ylabel('Y (pixels)')

fig.suptitle('Galaxy Image with Different Displays', fontsize=14)
plt.tight_layout()
plt.show()

# Simulated coordinate-display example: no FITS WCS is present here.
fig = plt.figure(figsize=(10, 8))

# Create synthetic data
ra_range = np.linspace(150, 151, 100)  # RA in degrees
dec_range = np.linspace(2, 3, 100)     # Dec in degrees
RA, DEC = np.meshgrid(ra_range, dec_range)

# Synthetic astronomical image
image = np.exp(-((RA - 150.5)**2 + (DEC - 2.5)**2) / 0.01)

# Standard matplotlib plot (pixel coordinates)
ax1 = fig.add_subplot(121)
im1 = ax1.imshow(image, origin='lower', cmap='viridis')
ax1.set_xlabel('X (pixels)')
ax1.set_ylabel('Y (pixels)')
ax1.set_title('Pixel Coordinates')
plt.colorbar(im1, ax=ax1)

# Approximate RA/Dec display; `extent` is not a WCS transform.
ax2 = fig.add_subplot(122)
im2 = ax2.imshow(image, origin='lower', cmap='viridis',
                 extent=[ra_range[0], ra_range[-1],
                        dec_range[0], dec_range[-1]])
ax2.set_xlabel('RA (degrees)')
ax2.set_ylabel('Dec (degrees)')
ax2.set_title('Approximate RA/Dec Display (not WCS)')
ax2.invert_xaxis()  # RA increases to the left
plt.colorbar(im2, ax=ax2)

plt.suptitle('Pixel Coordinates vs an Approximate Coordinate Display')
plt.tight_layout()
plt.show()

from astropy.io import fits
from astropy.wcs import WCS

with fits.open('science_image.fits') as hdul:
    image = hdul[0].data
    wcs = WCS(hdul[0].header)

fig, ax = plt.subplots(subplot_kw={'projection': wcs})
ax.imshow(image, origin='lower', cmap='viridis')
ax.coords.grid(color='white', ls='solid', alpha=0.4)
ax.set(xlabel='Right Ascension', ylabel='Declination')

# Example of journal-specific styling
# Note: SciencePlots package provides pre-configured styles

# Manual compact serif style; adapt it to the figure's final placement.
def apply_compact_serif_style():
    """Apply a compact serif style for a draft figure."""
    plt.rcParams.update({
        'font.size': 10,
        'font.family': 'serif',
        'font.serif': ['Computer Modern Roman'],
        'text.usetex': False,  # Set True if LaTeX is installed
        'figure.figsize': (3.5, 3.5),  # Single column
        'figure.dpi': 300,
        'savefig.dpi': 300,
        'savefig.bbox': 'tight',
        'axes.linewidth': 0.5,
        'lines.linewidth': 1.0,
        'lines.markersize': 4,
        'xtick.major.width': 0.5,
        'ytick.major.width': 0.5,
    })

# Example plot with journal style
apply_compact_serif_style()
fig, ax = plt.subplots()
x = np.linspace(0, 2*np.pi, 100)
ax.plot(x, np.sin(x), label='sin(x)')
ax.plot(x, np.cos(x), label='cos(x)')
ax.set_xlabel('Phase')
ax.set_ylabel('Amplitude')
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()

# Reset to defaults
plt.rcParams.update(plt.rcParamsDefault)

# Conceptual example of yt usage (requires yt installation)
"""
# yt installation: pip install yt
import yt

# Load simulation data (supports many formats)
ds = yt.load("galaxy_simulation/output_0050")

# Create a slice plot through the simulation volume
slc = yt.SlicePlot(ds, 'z', 'density')
slc.set_cmap('density', 'viridis')
slc.annotate_timestamp()
slc.save()

# Volume rendering
sc = yt.create_scene(ds, 'density')
sc.save()

# Phase plots
plot = yt.PhasePlot(ds.all_data(),
                    'density', 'temperature', 'cell_mass')
plot.save()
"""

# Demonstrate the concept with synthetic data
from matplotlib.colors import LogNorm

fig, axes = plt.subplots(2, 3, figsize=(14, 9))

# Simulate different yt-style visualizations
# 1. Density slice
x = np.linspace(-10, 10, 100)
y = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(x, y)
density = np.exp(-(X**2 + Y**2)/10) + 0.5*np.exp(-((X-3)**2 + Y**2)/5)

im1 = axes[0, 0].imshow(density, cmap='magma', origin='lower')
axes[0, 0].set_title('Density Slice (z=0)')
plt.colorbar(im1, ax=axes[0, 0], label='Synthetic density (arbitrary units)')

# 2. Temperature slice
rng = np.random.default_rng(42)
temperature = 1e4 * density**0.5 + rng.normal(0, 100, density.shape)
im2 = axes[0, 1].imshow(temperature, cmap='hot', origin='lower')
axes[0, 1].set_title('Temperature Slice')
plt.colorbar(im2, ax=axes[0, 1], label='T [K]')

# 3. Velocity field
u = -Y / np.sqrt(X**2 + Y**2 + 0.1)
v = X / np.sqrt(X**2 + Y**2 + 0.1)
axes[0, 2].streamplot(x, y, u, v, density=1, color='blue', alpha=0.7)
axes[0, 2].set_title('Velocity Streamlines')
axes[0, 2].set_xlim(-10, 10)
axes[0, 2].set_ylim(-10, 10)

# 4. Phase plot
n_points = 1000
dens_sample = rng.lognormal(0, 1, n_points)
temp_sample = 1e4 * dens_sample**0.5 + rng.normal(0, 500, n_points)
h = axes[1, 0].hist2d(np.log10(dens_sample), np.log10(temp_sample),
                       bins=30, cmap='YlOrBr')
axes[1, 0].set_xlabel('log(Density)')
axes[1, 0].set_ylabel('log(Temperature)')
axes[1, 0].set_title('Phase Diagram')
plt.colorbar(h[3], ax=axes[1, 0])

# 5. Projection (column density)
column_density = np.sum(np.stack([density]*20), axis=0)  # Fake projection
im5 = axes[1, 1].imshow(column_density, cmap='viridis',
                         origin='lower', norm=LogNorm())
axes[1, 1].set_title('Column Density Projection')
plt.colorbar(im5, ax=axes[1, 1], label='Synthetic column sum (arbitrary units)')

# 6. Multi-field composite
axes[1, 2].imshow(density, cmap='Blues', alpha=0.7, origin='lower')
contour = axes[1, 2].contour(temperature, levels=5, colors='red',
                              linewidths=2, alpha=0.8)
axes[1, 2].set_title('Multi-field Visualization')
axes[1, 2].clabel(contour, inline=True, fontsize=8)

fig.suptitle('yt-style Visualizations for Simulation Data',
             fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

print("\nyt advantages for simulation data:")
print("✓ Supports many simulation-data frontends; verify your output format")
print("✓ Handles AMR and particle data seamlessly")
print("✓ Unit-aware (automatic unit conversions)")
print("✓ Parallel processing for large datasets")
print("✓ Volume rendering capabilities")
print("✓ On-the-fly derived fields")
print("\nLearn more at: https://yt-project.org/")

# Plotly example (requires installation: pip install plotly)
import plotly.express as px
import plotly.graph_objects as go
import numpy as np

# Create interactive 3D scatter plot
rng = np.random.default_rng(42)
n = 500
x = rng.normal(0, 1, n)
y = rng.normal(0, 1, n)
z = x**2 + y**2 + rng.normal(0, 0.2, n)
colors = np.sqrt(x**2 + y**2)

fig = go.Figure(data=[go.Scatter3d(
    x=x, y=y, z=z,
    mode='markers',
    marker=dict(size=3, color=colors, colorscale='Viridis', opacity=0.8)
)])

fig.update_layout(
    title='Interactive 3D Scatter (drag to rotate!)',
    scene=dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z')
)
fig.show()

print("Plotly excels at:")
print("✓ Interactive exploration (zoom, pan, hover)")
print("✓ 3D visualization")
print("✓ Web deployment and dashboards")
print("✓ Animations")

# Astropy visualization (requires: pip install astropy)
from astropy.visualization import astropy_mpl_style, quantity_support
from astropy.coordinates import SkyCoord
from astropy import units as u
import matplotlib.pyplot as plt

# Enable Astropy's matplotlib style
plt.style.use(astropy_mpl_style)
quantity_support()

# Create sky coordinates
ra = [10.5, 15.2, 20.8, 25.1, 30.0] * u.degree
dec = [-25.3, -30.1, -22.7, -28.5, -24.2] * u.degree
coords = SkyCoord(ra=ra, dec=dec, frame='icrs')

# Plot with proper coordinate handling
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='aitoff')
ax.scatter(coords.ra.wrap_at(180*u.deg).radian, coords.dec.radian,
           s=50, c='red', marker='*')
ax.set_title('All-Sky Projection (Aitoff)')
ax.grid(True)
plt.show()

print("Astropy visualization excels at:")
print("✓ WCS (World Coordinate System) support")
print("✓ All-sky projections")
print("✓ Unit-aware plotting")
print("✓ FITS image display")