astropy

What this skill does

# Astropy - Astronomy & Astrophysics

Astropy is more than a library; it's a community standard. It ensures that astronomical calculations are reproducible and physically consistent by linking numerical arrays with physical units and celestial reference frames.

## When to Use

- Handling physical units and constants in calculations.
- Working with celestial coordinates (Right Ascension, Declination) and frame transformations.
- Reading/writing FITS (Flexible Image Transport System), ASDF, and VO (Virtual Observatory) files.
- Managing astronomical time scales (UTC, TAI, TDB, Julian Dates).
- World Coordinate System (WCS) mapping (pixels to sky coordinates).
- Cosmological calculations (distances, ages, Hubble parameters).
- Tabular data with attached units and metadata (astropy.table).

## Reference Documentation

**Official docs**: https://docs.astropy.org/  
**Learn Astropy**: https://learn.astropy.org/  
**Search patterns**: `astropy.units.Quantity`, `astropy.coordinates.SkyCoord`, `astropy.io.fits`, `astropy.time.Time`

## Core Principles

### Units and Quantities
Everything in Astropy should be a Quantity — a number or array paired with a Unit. This prevents errors like adding meters to feet.

### Celestial Coordinates
Coordinates are represented by SkyCoord objects, which handle the complex math of spherical trigonometry and frame precession automatically.

### Tables with Units
While Pandas is great for dataframes, astropy.table.QTable is superior for physics because it preserves units during operations.

## Quick Reference

### Installation

```bash
pip install astropy
```

### Standard Imports

```python
import numpy as np
from astropy import units as u
from astropy import constants as const
from astropy.coordinates import SkyCoord
from astropy.time import Time
from astropy.table import QTable
from astropy.io import fits
from astropy.wcs import WCS
```

### Basic Pattern - Physical Calculation

```python
# Calculate Schwarzschild radius of the Sun
mass_sun = 1.0 * u.M_sun
rs = (2 * const.G * mass_sun) / (const.c**2)

print(f"Radius: {rs.to(u.km):.2f}")
# Output: Radius: 2.95 km
```

## Critical Rules

### ✅ DO

- **Use astropy.units** - Always attach units to raw numbers from the start.
- **Prefer QTable** - Use QTable instead of Table to keep columns as Quantity objects.
- **Use SkyCoord for logic** - Never manually calculate angular separations; use `c1.separation(c2)`.
- **Vectorize coordinates** - SkyCoord can hold arrays of positions, which is much faster than loops.
- **Specify Time Scales** - Always define the scale (utc, tdb) when creating Time objects.
- **Use with for FITS** - Always use context managers to ensure file handles are closed.

### ❌ DON'T

- **Strip units early** - Avoid using `.value` until you absolutely need to pass data to a non-astropy function.
- **Hardcode constants** - Don't use 3e8 for speed of light; use `const.c`.
- **Assume J2000** - Be explicit about the frame (ICRS, FK5, Galactic).
- **Ignore WCS** - Don't assume pixels are linear; use the WCS object for mapping.

## Anti-Patterns (NEVER)

```python
# ❌ BAD: Manual unit conversion (Error prone)
distance_pc = 10.0
distance_m = distance_pc * 3.086e16 

# ✅ GOOD: Automatic conversion
dist = 10.0 * u.pc
dist_m = dist.to(u.m)

# ❌ BAD: Stripping units to use with NumPy functions
# (Astropy Quantities are NumPy-compatible, no need to strip!)
import numpy as np
vals = [1, 2, 3] * u.m
res = np.mean(vals.value) # ❌ Unit info lost

# ✅ GOOD: Keep units
res = np.mean(vals) # ✅ Returns <Quantity 2. m>

# ❌ BAD: Manual RA/Dec math
# sep = np.sqrt((ra1-ra2)**2 + (dec1-dec2)**2) # ❌ Incorrect for spherical geometry

# ✅ GOOD: Proper separation
c1 = SkyCoord(ra=10*u.deg, dec=20*u.deg)
c2 = SkyCoord(ra=11*u.deg, dec=21*u.deg)
sep = c1.separation(c2)
```

## Units and Constants (astropy.units)

### Complex Conversions

```python
# Equivalencies (e.g., converting wavelength to frequency)
wav = 500 * u.nm
freq = wav.to(u.THz, equivalencies=u.spectral())

# Flux conversions
flux_density = 15 * u.Jy
obs_freq = 1.4 * u.GHz
brightness_temp = flux_density.to(u.K, equivalencies=u.brightness_temperature(obs_freq, beam_area=1*u.sr))

# Solar units
lum = 3.8e26 * u.W
print(lum.to(u.L_sun))
```

## Tables (astropy.table)

### Managing Astronomical Catalogs

```python
# Create a QTable (Quantity Table)
tbl = QTable()
tbl['source_id'] = ['Star A', 'Star B', 'Star C']
tbl['flux'] = [1.2, 2.5, 0.8] * u.mJy
tbl['distance'] = [10, 25, 100] * u.pc

# Filtering
nearby = tbl[tbl['distance'] < 50 * u.pc]

# New column based on existing ones
tbl['luminosity'] = 4 * np.pi * tbl['distance']**2 * tbl['flux']

# Exporting
tbl.write('catalog.fits', overwrite=True)
tbl.write('catalog.csv', format='ascii.csv')
```

## Coordinates (astropy.coordinates)

### Frame Transformations

```python
# Create coordinate from string
c = SkyCoord("05h35m17.3s", "-05d23m28s", frame='icrs')

# Convert to Galactic coordinates
gal = c.galactic
print(f"L: {gal.l:.2f}, B: {gal.b:.2f}")

# Transform to AltAz (Horizontal) for a specific location and time
from astropy.coordinates import EarthLocation, AltAz
loc = EarthLocation(lat=45*u.deg, lon=-110*u.deg, height=2000*u.m)
time = Time.now()

altaz = c.transform_to(AltAz(obstime=time, location=loc))
print(f"Altitude: {altaz.alt:.2f}")
```

## Time (astropy.time)

### Precise Timekeeping

```python
# Julian Date vs ISO
t = Time("2023-12-01T12:00:00", scale='utc')
print(t.jd)
print(t.mjd) # Modified Julian Date

# Time deltas
dt = 24 * u.hour
future_t = t + dt

# Siderial Time
print(t.sidereal_time('mean', longitude=-110*u.deg))
```

## FITS and WCS (astropy.io.fits, astropy.wcs)

### Image Data Handling

```python
# Reading FITS
with fits.open("galaxy.fits") as hdul:
    data = hdul[0].data
    header = hdul[0].header
    
    # Initialize WCS from header
    w = WCS(header)

# Convert Pixel to World (Sky)
pix_coords = [[100, 100], [200, 200]]
world_coords = w.pixel_to_world(pix_coords[0][0], pix_coords[0][1])

# Convert Sky to Pixel
target = SkyCoord(ra=150.0*u.deg, dec=2.0*u.deg, frame='icrs')
px, py = w.world_to_pixel(target)
```

## Practical Workflows

### 1. Distance Calculation from Parallax

```python
def get_absolute_magnitude(mag, parallax_mas):
    """Calculate absolute magnitude given apparent mag and parallax in mas."""
    dist = (parallax_mas * u.mas).to(u.pc, equivalencies=u.parallax())
    # Absolute mag M = m - 5*log10(d/10pc)
    abs_mag = mag - 5 * np.log10(dist / (10 * u.pc))
    return abs_mag

# Star with mag 15 and 2 mas parallax
M = get_absolute_magnitude(15, 2)
```

### 2. Aperture Photometry Units Correction

```python
def flux_to_mag(flux, zero_point=25.0):
    """Convert raw flux to magnitudes with unit checking."""
    if not isinstance(flux, u.Quantity):
        raise ValueError("Flux must have units (e.g., u.electron / u.s)")
    
    # Assuming zero_point is defined for 1 electron/s
    mag = -2.5 * np.log10(flux / (1 * u.electron / u.s)) + zero_point
    return mag
```

### 3. Cosmological Distance Evolution

```python
from astropy.cosmology import FlatLambdaCDM

# Define cosmology (Planck 2018 defaults)
cosmo = FlatLambdaCDM(H0=67.4, Om0=0.315)

z = np.linspace(0, 5, 100)
dist = cosmo.luminosity_distance(z) # Returns distance in Mpc

print(f"Lookback time at z=1: {cosmo.lookback_time(1).to(u.Gyr):.2f}")
```

## Performance Optimization

### Using astropy.utils.console for Long Tasks

```python
from astropy.utils.console import ProgressBar

with ProgressBar(len(large_catalog)) as bar:
    for source in large_catalog:
        # Complex calculation
        bar.update()
```

### Fast I/O with C Engine

```python
# When reading huge CSV catalogs
tbl = QTable.read('massive_data.csv', format='ascii.fast_csv')
```

## Common Pitfalls and Solutions

### The "Unit Stripping" bug

```python
# ❌ Problem: Using a library that doesn't support Quantities
# solution = some_external_lib(data * u.m) # Crashes!

# ✅ Solution: Wrap the function or convert ex