astropy

What this skill does

# Astropy

## Overview

Astropy is the community standard Python library for astronomy, providing core functionality for astronomical data analysis and computation. This skill provides comprehensive guidance and tools for working with astropy's extensive capabilities across coordinate systems, file I/O, units and quantities, time systems, cosmology, modeling, and more.

## When to Use This Skill

This skill should be used when:
- Working with FITS files (reading, writing, inspecting, modifying)
- Performing coordinate transformations between astronomical reference frames
- Calculating cosmological distances, ages, or other quantities
- Handling astronomical time systems and conversions
- Working with physical units and dimensional analysis
- Processing astronomical data tables with specialized column types
- Fitting models to astronomical data
- Converting between pixel and world coordinates (WCS)
- Performing robust statistical analysis on astronomical data
- Visualizing astronomical images with proper scaling and stretching

## Core Capabilities

### 1. FITS File Operations

FITS (Flexible Image Transport System) is the standard file format in astronomy. Astropy provides comprehensive FITS support.

**Quick FITS Inspection**:
Use the included `scripts/fits_info.py` script for rapid file inspection:
```bash
python scripts/fits_info.py observation.fits
python scripts/fits_info.py observation.fits --detailed
python scripts/fits_info.py observation.fits --ext 1
```

**Common FITS workflows**:
```python
from astropy.io import fits

# Read FITS file
with fits.open('image.fits') as hdul:
    hdul.info()  # Display structure
    data = hdul[0].data
    header = hdul[0].header

# Write FITS file
fits.writeto('output.fits', data, header, overwrite=True)

# Quick access (less efficient for multiple operations)
data = fits.getdata('image.fits', ext=0)
header = fits.getheader('image.fits', ext=0)

# Update specific header keyword
fits.setval('image.fits', 'OBJECT', value='M31')
```

**Multi-extension FITS**:
```python
from astropy.io import fits

# Create multi-extension FITS
primary = fits.PrimaryHDU(primary_data)
image_ext = fits.ImageHDU(science_data, name='SCI')
error_ext = fits.ImageHDU(error_data, name='ERR')

hdul = fits.HDUList([primary, image_ext, error_ext])
hdul.writeto('multi_ext.fits', overwrite=True)
```

**Binary tables**:
```python
from astropy.io import fits

# Read binary table
with fits.open('catalog.fits') as hdul:
    table_data = hdul[1].data
    ra = table_data['RA']
    dec = table_data['DEC']

# Better: use astropy.table for table operations (see section 5)
```

### 2. Coordinate Systems and Transformations

Astropy supports ~25 coordinate frames with seamless transformations.

**Quick Coordinate Conversion**:
Use the included `scripts/coord_convert.py` script:
```bash
python scripts/coord_convert.py 10.68 41.27 --from icrs --to galactic
python scripts/coord_convert.py --file coords.txt --from icrs --to galactic --output sexagesimal
```

**Basic coordinate operations**:
```python
from astropy.coordinates import SkyCoord
import astropy.units as u

# Create coordinate (multiple input formats supported)
c = SkyCoord(ra=10.68*u.degree, dec=41.27*u.degree, frame='icrs')
c = SkyCoord('00:42:44.3 +41:16:09', unit=(u.hourangle, u.deg))
c = SkyCoord('00h42m44.3s +41d16m09s')

# Transform between frames
c_galactic = c.galactic
c_fk5 = c.fk5

print(f"Galactic: l={c_galactic.l.deg:.3f}, b={c_galactic.b.deg:.3f}")
```

**Working with coordinate arrays**:
```python
import numpy as np
from astropy.coordinates import SkyCoord
import astropy.units as u

# Arrays of coordinates
ra = np.array([10.1, 10.2, 10.3]) * u.degree
dec = np.array([40.1, 40.2, 40.3]) * u.degree
coords = SkyCoord(ra=ra, dec=dec, frame='icrs')

# Calculate separations
sep = coords[0].separation(coords[1])
print(f"Separation: {sep.to(u.arcmin)}")

# Position angle
pa = coords[0].position_angle(coords[1])
```

**Catalog matching**:
```python
from astropy.coordinates import SkyCoord
import astropy.units as u

catalog1 = SkyCoord(ra=[10, 11, 12]*u.degree, dec=[40, 41, 42]*u.degree)
catalog2 = SkyCoord(ra=[10.01, 11.02, 13]*u.degree, dec=[40.01, 41.01, 43]*u.degree)

# Find nearest neighbors
idx, sep2d, dist3d = catalog1.match_to_catalog_sky(catalog2)

# Filter by separation threshold
max_sep = 1 * u.arcsec
matched = sep2d < max_sep
```

**Horizontal coordinates (Alt/Az)**:
```python
from astropy.coordinates import SkyCoord, EarthLocation, AltAz
from astropy.time import Time
import astropy.units as u

location = EarthLocation(lat=40*u.deg, lon=-70*u.deg, height=300*u.m)
obstime = Time('2023-01-01 03:00:00')
target = SkyCoord(ra=10*u.degree, dec=40*u.degree, frame='icrs')

altaz_frame = AltAz(obstime=obstime, location=location)
target_altaz = target.transform_to(altaz_frame)

print(f"Alt: {target_altaz.alt.deg:.2f}°, Az: {target_altaz.az.deg:.2f}°")
```

**Available coordinate frames**:
- `icrs` - International Celestial Reference System (default, preferred)
- `fk5`, `fk4` - Fifth/Fourth Fundamental Katalog
- `galactic` - Galactic coordinates
- `supergalactic` - Supergalactic coordinates
- `altaz` - Horizontal (altitude-azimuth) coordinates
- `gcrs`, `cirs`, `itrs` - Earth-based systems
- Ecliptic frames: `BarycentricMeanEcliptic`, `HeliocentricMeanEcliptic`, `GeocentricMeanEcliptic`

### 3. Units and Quantities

Physical units are fundamental to astronomical calculations. Astropy's units system provides dimensional analysis and automatic conversions.

**Basic unit operations**:
```python
import astropy.units as u

# Create quantities
distance = 5.2 * u.parsec
velocity = 300 * u.km / u.s
time = 10 * u.year

# Convert units
distance_ly = distance.to(u.lightyear)
velocity_mps = velocity.to(u.m / u.s)

# Arithmetic with units
wavelength = 500 * u.nm
frequency = wavelength.to(u.Hz, equivalencies=u.spectral())
```

**Working with arrays**:
```python
import numpy as np
import astropy.units as u

wavelengths = np.array([400, 500, 600]) * u.nm
frequencies = wavelengths.to(u.THz, equivalencies=u.spectral())

fluxes = np.array([1.2, 2.3, 1.8]) * u.Jy
luminosities = 4 * np.pi * (10*u.pc)**2 * fluxes
```

**Important equivalencies**:
- `u.spectral()` - Convert wavelength ↔ frequency ↔ energy
- `u.doppler_optical(rest)` - Optical Doppler velocity
- `u.doppler_radio(rest)` - Radio Doppler velocity
- `u.doppler_relativistic(rest)` - Relativistic Doppler
- `u.temperature()` - Temperature unit conversions
- `u.brightness_temperature(freq)` - Brightness temperature

**Physical constants**:
```python
from astropy import constants as const

print(const.c)      # Speed of light
print(const.G)      # Gravitational constant
print(const.M_sun)  # Solar mass
print(const.R_sun)  # Solar radius
print(const.L_sun)  # Solar luminosity
```

**Performance tip**: Use the `<<` operator for fast unit assignment to arrays:
```python
# Fast
result = large_array << u.m

# Slower
result = large_array * u.m
```

### 4. Time Systems

Astronomical time systems require high precision and multiple time scales.

**Creating time objects**:
```python
from astropy.time import Time
import astropy.units as u

# Various input formats
t1 = Time('2023-01-01T00:00:00', format='isot', scale='utc')
t2 = Time(2459945.5, format='jd', scale='utc')
t3 = Time(['2023-01-01', '2023-06-01'], format='iso')

# Convert formats
print(t1.jd)    # Julian Date
print(t1.mjd)   # Modified Julian Date
print(t1.unix)  # Unix timestamp
print(t1.iso)   # ISO format

# Convert time scales
print(t1.tai)   # International Atomic Time
print(t1.tt)    # Terrestrial Time
print(t1.tdb)   # Barycentric Dynamical Time
```

**Time arithmetic**:
```python
from astropy.time import Time, TimeDelta
import astropy.units as u

t1 = Time('2023-01-01T00:00:00')
dt = TimeDelta(1*u.day)

t2 = t1 + dt
diff = t2 - t1
print(diff.to(u.hour))

# Array of times
times = t1 + np.arange(10) * u.day
```

**Astronomical time calculations**:
```python
from astropy.time import Time
from astrop