Examples

Basic Analysis of a CO Cube

This example shows how to perform basic analysis on a CO data cube:

from casa_cube import Cube
import matplotlib.pyplot as plt

# Load the data
cube = Cube('co_cube.fits')

# Create a figure with multiple panels
fig, axes = plt.subplots(2, 2, figsize=(12, 12))

# Plot integrated intensity (M0)
cube.plot(moment=0, ax=axes[0,0], title='Integrated Intensity')

# Plot velocity field (M1)
cube.plot(moment=1, ax=axes[0,1], title='Velocity Field')

# Plot velocity dispersion (M2)
cube.plot(moment=2, ax=axes[1,0], title='Velocity Dispersion')

# Plot peak intensity (M8)
cube.plot(moment=8, ax=axes[1,1], title='Peak Intensity')

plt.tight_layout()
plt.show()

Channel Map Analysis

This example demonstrates how to create and analyze channel maps:

# Create channel maps
cube.plot_channels(n=20, ncols=5, vmin=-10, vmax=10)

# Get line profile
cube.plot_line(x_axis="velocity", threshold=0.1)

# Calculate and plot moment maps with threshold
m0 = cube.get_moment_map(moment=0, threshold=0.1)
m1 = cube.get_moment_map(moment=1, v0=0, threshold=0.1)

# Plot with custom parameters
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

cube.plot(moment=0, ax=ax1, threshold=0.1,
          cmap='viridis', colorbar=True)
ax1.set_title('Integrated Intensity')

cube.plot(moment=1, ax=ax2, v0=0, threshold=0.1,
          cmap='RdBu_r', colorbar=True)
ax2.set_title('Velocity Field')

plt.tight_layout()
plt.show()

High-Pass Filtering

This example shows how to apply high-pass filtering to remove large-scale structures:

# Apply high-pass filter
filtered_map = cube.get_high_pass_filter_map(
    moment=0,
    w0=5.0,  # 5 arcsec scale
    gamma=0.5  # radial stretch
)

# Plot original and filtered maps
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

cube.plot(moment=0, ax=ax1, cmap='viridis')
ax1.set_title('Original')

ax2.imshow(filtered_map, cmap='viridis', origin='lower')
ax2.set_title('High-Pass Filtered')

plt.tight_layout()
plt.show()

Advanced Analysis

This example demonstrates some advanced features:

# Convert to brightness temperature
Tb = cube._Jybeam_to_Tb(cube.get_moment_map(moment=8))

# Calculate turbulent velocity
vturb = cube.get_vturb(mol_weight=28)  # for CO

# Create a cut through the cube
x, y, z = cube.make_cut(100, 100, 200, 200, num=100)

# Plot results
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

ax1.imshow(Tb, cmap='viridis', origin='lower')
ax1.set_title('Brightness Temperature')

ax2.imshow(vturb, cmap='viridis', origin='lower')
ax2.set_title('Turbulent Velocity')

plt.tight_layout()
plt.show()

# Plot the cut
plt.figure(figsize=(8, 6))
plt.plot(x, z)
plt.xlabel('Position')
plt.ylabel('Intensity')
plt.title('Cut Through Cube')
plt.show()