메뉴
Astronomical data analysis with astropy and astroquery — FITS I/O, WCS transforms, catalog cross-matching, aperture photometry, and CMB power spectra.
기본적으로 소스를 먼저 확인하는 Prompt가 선택됩니다. 직접 명령으로 전환하거나 로컬 사본을 다운로드할 수도 있습니다.
설치 여부를 결정하기 전에 SKILL.md와 SkillsMP에 표시된 보조 파일을 읽어 보세요.
Codex 또는 Claude로 설치 이 Prompt를 복사해 Codex, Claude 또는 다른 어시스턴트에 붙여 넣으면 Skill 페이지를 검토하고 설치를 진행할 수 있습니다.
직접 명령은 검토 Prompt를 거치지 않습니다. 실행하기 전에 소스를 확인하세요.
npx skills add https://github.com/xjtulyc/awesome-rosetta-skills --skill astropy-astronomy명령은 한 줄로 유지됩니다. 복사하기 전에 가로로 스크롤해 전체 내용을 확인하세요.
로컬 사본을 원하시나요? SkillsMP에서 현재 제공할 수 있는 파일을 다운로드하세요.
SOC 직업 분류 기준
SKILL.md 표시 중
| name | astropy-astronomy |
| description | Astronomical data analysis with astropy and astroquery — FITS I/O, WCS transforms, catalog cross-matching, aperture photometry, and CMB power spectra. |
| tags | ["astronomy","astrophysics","fits","wcs","photometry","cmb"] |
| version | 1.0.0 |
| authors | [{"name":"Rosetta Skills Contributors","github":"@xjtulyc"}] |
| license | MIT |
| platforms | ["claude-code","codex","gemini-cli","cursor"] |
| dependencies | ["astropy>=5.3","astroquery>=0.4.6","healpy>=1.16","matplotlib>=3.7","numpy>=1.24"] |
| last_updated | 2026-03-17 |
| status | stable |
A comprehensive skill for astronomical data analysis using the Python astronomy ecosystem. Covers reading and writing FITS files, World Coordinate System (WCS) transforms, remote catalog queries, aperture photometry on imaging data, and CMB angular power spectrum analysis using HEALPix pixelisation.
Use this skill when you need to:
.fits filesThis skill is not appropriate for:
specutils skill)yt or pynbody skill)FITS is the standard file format in astronomy. A FITS file contains one or more Header/Data Units (HDUs). Each HDU has a text header of 80-character keyword cards and an optional data array or binary table. Common HDU types:
| HDU type | astropy class | Typical use |
|---|---|---|
| PrimaryHDU | fits.PrimaryHDU | 2-D image, metadata-only |
| ImageHDU | fits.ImageHDU | Additional image planes |
| BinTableHDU | fits.BinTableHDU | Source catalogs, spectra |
| CompImageHDU | fits.CompImageHDU | Tile-compressed image |
The FITS WCS standard maps pixel indices (x, y) to sky coordinates
(RA, Dec) via a projection (TAN, ZEA, CAR, …) plus a linear transform
encoded in the header keywords CRPIX, CRVAL, / .
reads these keywords and provides /
methods that return objects.
CDCDELTastropy.wcs.WCSpixel_to_worldworld_to_pixelSkyCoordastropy.coordinates.SkyCoord represents one or more positions on the sky in
any supported reference frame. Angle arithmetic, frame conversions, and
on-sky separations are all handled automatically. Key frames: ICRS (RA/Dec),
Galactic (l, b), FK5, AltAz (horizontal coordinates).
astroquery provides uniform Python interfaces to dozens of online services:
All queries return astropy.table.Table objects directly.
HEALPix (Hierarchical Equal Area isoLatitude Pixelisation) divides the sphere
into 12 * nside**2 equal-area pixels. healpy wraps the C++ HEALPix
library and provides map I/O (read_map, write_map), spherical harmonic
transforms (map2alm, alm2cl), and visualisation (mollview, gnomview).
The angular power spectrum C_l is the key observable for CMB cosmology.
photutils (companion to astropy) implements circular and elliptical apertures,
sky annuli for local background estimation, and source detection (DAOStarFinder,
IRAFStarFinder). The result is a flux in counts per second that can be converted
to magnitudes using a photometric zero-point.
pip install "astropy>=5.3" "astroquery>=0.4.6" "healpy>=1.16" \
"matplotlib>=3.7" "numpy>=1.24" photutils
On Apple Silicon / conda:
conda install -c conda-forge astropy astroquery healpy matplotlib numpy photutils
# Fetch a small SDSS r-band cutout for testing
python - <<'EOF'
from astroquery.skyview import SkyView
imgs = SkyView.get_images("M51", survey=["SDSSr"], pixels=512)
imgs[0].writeto("m51_r.fits", overwrite=True)
print("Saved m51_r.fits")
EOF
No API keys are required for the public Vizier / SIMBAD / SkyView services. If you use the ESO archive or proprietary data portals, store credentials as:
export ESO_USERNAME="<your-username>"
export ESO_PASSWORD=$(cat ~/.eso_passwd) # read from file, never hardcode
Access them in Python with:
import os
username = os.getenv("ESO_USERNAME", "")
password = os.getenv("ESO_PASSWORD", "")
from astropy.io import fits
import numpy as np
# Open without loading data into memory
with fits.open("m51_r.fits") as hdul:
hdul.info() # print HDU summary
header = hdul[0].header
data = hdul[0].data.astype(float) # pixel array (counts)
print(f"Image shape : {data.shape}")
print(f"Instrument : {header.get('INSTRUME', 'unknown')}")
print(f"Filter : {header.get('FILTER', 'unknown')}")
print(f"Exposure : {header.get('EXPTIME', 'unknown')} s")
print(f"Min / Max : {data.min():.1f} / {data.max():.1f} counts")
# Mask NaN / Inf pixels that can corrupt photometry
data = np.where(np.isfinite(data), data, 0.0)
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
import astropy.units as u
with fits.open("m51_r.fits") as hdul:
wcs = WCS(hdul[0].header)
data = hdul[0].data.astype(float)
ny, nx = data.shape
# Pixel centre of the image -> sky coordinates
cx, cy = nx / 2.0, ny / 2.0
sky_centre = wcs.pixel_to_world(cx, cy)
print(f"Image centre: RA={sky_centre.ra.deg:.4f} deg, "
f"Dec={sky_centre.dec.deg:.4f} deg")
# Convert to Galactic coordinates
gal = sky_centre.galactic
print(f"Galactic : l={gal.l.deg:.4f} deg, b={gal.b.deg:.4f} deg")
# Sky -> pixel (useful for placing apertures on known objects)
m51_coord = SkyCoord(ra=202.4696 * u.deg, dec=47.1952 * u.deg)
pix = wcs.world_to_pixel(m51_coord)
print(f"M51 nucleus at pixel ({pix[0]:.1f}, {pix[1]:.1f})")
# Compute pixel scale
pixel_scale = wcs.proj_plane_pixel_scales()
print(f"Pixel scale : {pixel_scale[0].to(u.arcsec):.3f} / pixel")
from astroquery.vizier import Vizier
from astroquery.simbad import Simbad
from astropy.coordinates import SkyCoord
import astropy.units as u
# -- Gaia DR3 point sources within 5 arcmin of M51 --
coord = SkyCoord(ra=202.4696 * u.deg, dec=47.1952 * u.deg, frame="icrs")
radius = 5 * u.arcmin
v = Vizier(columns=["Source", "RA_ICRS", "DE_ICRS", "Gmag", "BP-RP"],
row_limit=500)
result = v.query_region(coord, radius=radius, catalog="I/355/gaiadr3")
if result:
gaia = result[0]
print(f"Found {len(gaia)} Gaia DR3 sources")
print(gaia["Source", "RA_ICRS", "DE_ICRS", "Gmag"][:5])
else:
print("No Gaia sources returned (check network)")
# -- SIMBAD identification of M51 --
simbad = Simbad()
simbad.add_votable_fields("distance", "flux(V)", "sptype")
result_s = simbad.query_object("M51")
if result_s:
print("\nSIMBAD result for M51:")
print(result_s["MAIN_ID", "RA", "DEC", "FLUX_V"])
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
import astropy.units as u
from photutils.detection import DAOStarFinder
from photutils.aperture import (CircularAperture, CircularAnnulus,
aperture_photometry)
from photutils.background import MMMBackground
with fits.open("m51_r.fits") as hdul:
data = hdul[0].data.astype(float)
header = hdul[0].header
wcs = WCS(header)
# Estimate and subtract sky background
bkg_estimator = MMMBackground()
bkg = bkg_estimator(data)
data_sub = data - bkg
print(f"Background level: {bkg:.2f} counts")
# Detect sources with DAOStarFinder
from astropy.stats import sigma_clipped_stats
_, median, std = sigma_clipped_stats(data, sigma=3.0)
finder = DAOStarFinder(fwhm=3.5, threshold=5.0 * std)
sources = finder(data_sub)
print(f"Detected {len(sources)} sources")
# Perform circular aperture photometry (r=5 px, sky annulus 7-10 px)
positions = np.column_stack([sources["xcentroid"], sources["ycentroid"]])
aperture = CircularAperture(positions, r=5.0)
annulus = CircularAnnulus(positions, r_in=7.0, r_out=10.0)
phot_table = aperture_photometry(data_sub, [aperture, annulus])
# Local background correction
annulus_area = annulus.area
aperture_area = aperture.area
sky_per_pixel = phot_table["aperture_sum_1"] / annulus_area
sky_in_aper = sky_per_pixel * aperture_area
flux_corrected = phot_table["aperture_sum_0"] - sky_in_aper
# Instrumental magnitude (zero-point = 25.0 as example)
ZERO_POINT = 25.0
exptime = header.get("EXPTIME", 1.0)
flux_rate = flux_corrected / exptime
mag = ZERO_POINT - 2.5 * np.log10(np.abs(flux_rate))
phot_table["flux_corrected"] = flux_corrected
phot_table["mag_inst"] = mag
phot_table.sort("mag_inst")
print(phot_table["id", "xcenter", "ycenter", "flux_corrected", "mag_inst"][:10])
import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
# ---- Simulate a Gaussian CMB temperature map ----
nside = 512 # HEALPix resolution (npix = 12*nside^2)
lmax = 3 * nside - 1 # maximum multipole
npix = hp.nside2npix(nside)
# Planck 2018-like flat LCDM power spectrum (toy model)
ells = np.arange(lmax + 1)
Cl = np.zeros(lmax + 1)
# Simple Sachs-Wolfe plateau + first acoustic peak shape
Cl[2:] = 1e-10 * (1.0 / (ells[2:] * (ells[2:] + 1))) * (
1.0 + 6.0 * np.exp(-((ells[2:] - 200) / 40) ** 2)
)
# Generate a realisation
np.random.seed(42)
cmb_map = hp.synfast(Cl, nside=nside, lmax=lmax, verbose=False)
print(f"Map shape : {cmb_map.shape}, RMS = {cmb_map.std():.3e} K")
# Recover power spectrum from the map
alm = hp.map2alm(cmb_map, lmax=lmax)
Cl_rec = hp.alm2cl(alm)
ell_rec = np.arange(len(Cl_rec))
# Visualise
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
hp.mollview(cmb_map, title="Simulated CMB (K)", unit="K",
hold=False, fig=fig.number, sub=(1, 2, 1))
ax = axes[1]
ax.plot(ell_rec[2:], ell_rec[2:] * (ell_rec[2:] + 1) * Cl_rec[2:] / (2 * np.pi),
lw=1.5, label=r"Recovered $C_\ell$")
ax.set_xlabel(r"Multipole $\ell$")
ax.set_ylabel(r"$\ell(\ell+1)C_\ell / 2\pi$ [K$^2$]")
ax.set_title("Angular Power Spectrum")
ax.legend()
fig.tight_layout()
plt.savefig("cmb_power_spectrum.png", dpi=150)
print("Saved cmb_power_spectrum.png")
from astropy.coordinates import SkyCoord, match_coordinates_sky
import astropy.units as u
import numpy as np
# Suppose cat_a and cat_b are astropy Tables with RA/Dec columns
# Here we mock them with random positions for illustration
rng = np.random.default_rng(0)
n_a, n_b = 200, 5000
ra_a = rng.uniform(202.0, 203.0, n_a)
dec_a = rng.uniform(46.8, 47.5, n_a)
ra_b = rng.uniform(201.5, 203.5, n_b)
dec_b = rng.uniform(46.5, 47.8, n_b)
coords_a = SkyCoord(ra=ra_a * u.deg, dec=dec_a * u.deg, frame="icrs")
coords_b = SkyCoord(ra=ra_b * u.deg, dec=dec_b * u.deg, frame="icrs")
# For each source in A, find the nearest in B
idx, sep2d, _ = match_coordinates_sky(coords_a, coords_b)
# Keep only matches within 1 arcsec
match_mask = sep2d < 1.0 * u.arcsec
n_match = match_mask.sum()
print(f"Matched {n_match} / {n_a} sources within 1 arcsec")
matched_a_idx = np.where(match_mask)[0]
matched_b_idx = idx[match_mask]
separations = sep2d[match_mask].to(u.arcsec)
print(f"Median separation : {np.median(separations):.3f}")
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
ny, nx = 512, 512
wcs = WCS(naxis=2)
wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]
wcs.wcs.crpix = [nx / 2 + 0.5, ny / 2 + 0.5] # reference pixel (1-based)
wcs.wcs.crval = [202.4696, 47.1952] # RA, Dec of reference (deg)
wcs.wcs.cdelt = [-0.000277778, 0.000277778] # 1 arcsec/pixel
wcs.wcs.cunit = ["deg", "deg"]
# Create a synthetic Gaussian source
rng = np.random.default_rng(1)
image = rng.normal(100.0, 10.0, (ny, nx)) # sky background + noise
yy, xx = np.mgrid[:ny, :nx]
# Place a point source at pixel (256, 256) with peak 5000 counts
sigma = 2.5
image += 5000.0 * np.exp(-((xx - 256) ** 2 + (yy - 256) ** 2) / (2 * sigma ** 2))
# Write to FITS
header = wcs.to_header()
header["EXPTIME"] = (300.0, "Exposure time [s]")
header["FILTER"] = ("r", "Filter name")
hdu = fits.PrimaryHDU(data=image.astype(np.float32), header=header)
hdu.writeto("synthetic_field.fits", overwrite=True)
print("Saved synthetic_field.fits")
from astroquery.gaia import Gaia
# Download Gaia DR3 stars brighter than G=15 around M51
query = """
SELECT TOP 200
source_id, ra, dec, parallax, pmra, pmdec,
phot_g_mean_mag, phot_bp_mean_mag, phot_rp_mean_mag
FROM gaiadr3.gaia_source
WHERE CONTAINS(
POINT('ICRS', ra, dec),
CIRCLE('ICRS', 202.4696, 47.1952, 0.1)
) = 1
AND phot_g_mean_mag < 15
ORDER BY phot_g_mean_mag ASC
"""
job = Gaia.launch_job(query)
result = job.get_results()
print(f"Retrieved {len(result)} Gaia DR3 stars")
print(result["source_id", "ra", "dec", "phot_g_mean_mag"][:5])
from astropy.io import fits
from astropy.table import Table
import numpy as np
# Example: Chandra source catalogue (or any BinTable FITS)
# hdul = fits.open("chandra_sources.fits")
# For demonstration, create a mock BinTable
n = 100
rng = np.random.default_rng(2)
col_ra = fits.Column(name="RA", format="D", array=rng.uniform(0, 360, n))
col_dec = fits.Column(name="DEC", format="D", array=rng.uniform(-90, 90, n))
col_flux = fits.Column(name="FLUX", format="E",
array=rng.exponential(1e-13, n).astype(np.float32))
col_name = fits.Column(name="NAME", format="20A",
array=[f"CXO J{i:04d}" for i in range(n)])
hdu_table = fits.BinTableHDU.from_columns([col_ra, col_dec, col_flux, col_name])
hdu_table.name = "SOURCES"
primary = fits.PrimaryHDU()
hdul = fits.HDUList([primary, hdu_table])
hdul.writeto("mock_sources.fits", overwrite=True)
# Read back as an astropy Table (cleanest interface)
tbl = Table.read("mock_sources.fits", hdu="SOURCES")
print(tbl[:5])
print(f"Brightest source: {tbl[np.argmax(tbl['FLUX'])]['NAME']}")
astroquery times out or returns empty resultsVizier.TIMEOUT = 60 before the query.from astropy.table import Table
import os
CACHE = "gaia_cache.fits"
if os.path.exists(CACHE):
gaia = Table.read(CACHE)
else:
# ... run Vizier query ...
gaia.write(CACHE, overwrite=True)
WCS warnings about non-standard keywordsimport warnings
from astropy.wcs import FITSFixedWarning
warnings.filterwarnings("ignore", category=FITSFixedWarning)
photutils source detection finds too many / too few sourcesAdjust the detection threshold. A typical starting point:
from astropy.stats import sigma_clipped_stats
_, _, std = sigma_clipped_stats(data, sigma=3.0)
finder = DAOStarFinder(fwhm=3.0, threshold=5.0 * std) # raise threshold to reduce
finder = DAOStarFinder(fwhm=3.0, threshold=3.0 * std) # lower to detect fainter
Also ensure the PSF FWHM parameter matches the actual image seeing in pixels.
healpy crashes on import (macOS)On Apple Silicon, install via conda-forge rather than pip:
conda install -c conda-forge healpy
data = hdul[0].data
if data.dtype.byteorder not in ("=", "|"):
data = data.byteswap().newbyteorder()
"""
End-to-end aperture photometry pipeline:
1. Download SDSS r-band image of M51 via SkyView
2. Read FITS, solve WCS
3. Detect sources with DAOStarFinder
4. Perform aperture photometry
5. Cross-match detected sources against Gaia DR3
6. Compute instrumental zero-point from Gaia G magnitudes
7. Save a calibrated catalog to FITS
"""
import numpy as np
import astropy.units as u
from astropy.io import fits
from astropy.wcs import WCS
from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy.stats import sigma_clipped_stats
from astroquery.skyview import SkyView
from astroquery.vizier import Vizier
from photutils.detection import DAOStarFinder
from photutils.aperture import CircularAperture, CircularAnnulus, aperture_photometry
# --- 1. Download image ---
print("Downloading SDSS r-band image …")
imgs = SkyView.get_images("M51", survey=["SDSSr"], pixels=512)
hdul = imgs[0]
data = hdul[0].data.astype(float)
header = hdul[0].header
wcs = WCS(header)
exptime = float(header.get("EXPTIME", 53.91)) # SDSS typical 53.91 s
hdul.writeto("m51_sdssr.fits", overwrite=True)
# --- 2. Background subtraction ---
_, bkg_med, bkg_std = sigma_clipped_stats(data, sigma=3.0)
data_sub = data - bkg_med
# --- 3. Source detection ---
finder = DAOStarFinder(fwhm=3.5, threshold=5.0 * bkg_std)
sources = finder(data_sub)
print(f"Detected {len(sources)} sources")
# --- 4. Aperture photometry ---
positions = np.column_stack([sources["xcentroid"], sources["ycentroid"]])
aper = CircularAperture(positions, r=5.0)
ann = CircularAnnulus(positions, r_in=7.0, r_out=10.0)
phot = aperture_photometry(data_sub, [aper, ann])
sky_bkg = (phot["aperture_sum_1"] / ann.area) * aper.area
flux = phot["aperture_sum_0"] - sky_bkg
flux_rate = flux / exptime
# --- 5. WCS -> sky coordinates for each source ---
sky_pos = wcs.pixel_to_world(sources["xcentroid"], sources["ycentroid"])
sources["ra"] = sky_pos.ra.deg
sources["dec"] = sky_pos.dec.deg
# --- 6. Cross-match against Gaia DR3 ---
field_centre = wcs.pixel_to_world(data.shape[1] / 2, data.shape[0] / 2)
v = Vizier(columns=["Source", "RA_ICRS", "DE_ICRS", "Gmag"], row_limit=2000)
gaia_result = v.query_region(field_centre, radius=5 * u.arcmin,
catalog="I/355/gaiadr3")
gaia = gaia_result[0] if gaia_result else None
zp = 25.0 # fallback zero-point
if gaia is not None:
det_coords = SkyCoord(ra=sources["ra"] * u.deg, dec=sources["dec"] * u.deg)
gaia_coords = SkyCoord(ra=gaia["RA_ICRS"] * u.deg, dec=gaia["DE_ICRS"] * u.deg)
idx, sep, _ = det_coords.match_to_catalog_sky(gaia_coords)
good = sep < 1.0 * u.arcsec
if good.sum() > 5:
inst_mag = -2.5 * np.log10(np.abs(flux_rate[good]))
gaia_mag = gaia["Gmag"][idx[good]].data.astype(float)
zp = float(np.median(gaia_mag - inst_mag))
print(f"Photometric zero-point (Gaia G): {zp:.3f} mag")
# --- 7. Build calibrated catalog ---
mag_cal = zp - 2.5 * np.log10(np.abs(flux_rate))
catalog = Table()
catalog["id"] = sources["id"]
catalog["ra"] = sources["ra"]
catalog["dec"] = sources["dec"]
catalog["flux"] = flux_rate
catalog["mag"] = mag_cal
catalog.sort("mag")
catalog.write("m51_catalog.fits", overwrite=True)
print("Saved m51_catalog.fits with", len(catalog), "sources")
print(catalog["id", "ra", "dec", "mag"][:10])
"""
CMB power spectrum analysis:
1. Simulate a CMB temperature map with realistic Planck beam
2. Apply a Galactic mask (|b| < 20 deg)
3. Compute pseudo-Cl power spectrum
4. Correct for the beam and pixel window functions
5. Plot the final power spectrum
"""
import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
# --- Parameters ---
nside = 256
lmax = 2 * nside
fwhm_deg = 0.5 # Planck 143 GHz beam ~ 7 arcmin; use 30 arcmin here
fwhm_rad = np.radians(fwhm_deg)
npix = hp.nside2npix(nside)
ells = np.arange(lmax + 1)
# --- 1. Build input power spectrum (Planck 2018-like toy model) ---
Cl_in = np.zeros(lmax + 1)
Cl_in[2:] = 2.7e-9 / (ells[2:] * (ells[2:] + 1) / (2 * np.pi)) * (
1 + 5 * np.exp(-((ells[2:] - 200) / 50) ** 2)
+ 2 * np.exp(-((ells[2:] - 540) / 40) ** 2)
)
# --- 2. Simulate map with beam ---
np.random.seed(7)
beam_window = hp.gauss_beam(fwhm_rad, lmax=lmax)
alm_true = hp.synalm(Cl_in, lmax=lmax)
alm_beam = hp.almxfl(alm_true, beam_window)
cmb_map = hp.alm2map(alm_beam, nside=nside, verbose=False)
# --- 3. Build Galactic mask (cut |b| < 20 deg) ---
theta, phi = hp.pix2ang(nside, np.arange(npix))
b_rad = np.pi / 2 - theta # colatitude -> Galactic latitude proxy
mask = np.abs(b_rad) > np.radians(20)
fsky = mask.mean()
print(f"Sky fraction after masking: f_sky = {fsky:.3f}")
masked_map = cmb_map * mask
# --- 4. Pseudo-Cl estimation ---
alm_masked = hp.map2alm(masked_map, lmax=lmax)
Cl_pseudo = hp.alm2cl(alm_masked)
# --- 5. Beam and pixel window correction ---
pixel_win = hp.pixwin(nside)[:lmax + 1]
correction = (beam_window ** 2) * (pixel_win ** 2) * fsky
correction = np.where(correction > 1e-30, correction, 1e-30)
Cl_corr = Cl_pseudo / correction
# --- 6. Plot ---
fig, ax = plt.subplots(figsize=(10, 5))
norm = ells * (ells + 1) / (2 * np.pi)
ax.plot(ells[2:], norm[2:] * Cl_in[2:], lw=2, label="Input $C_\\ell$", color="k")
ax.plot(ells[2:], norm[2:] * Cl_corr[2:], lw=1.5, label="Recovered (beam-corrected)",
color="tab:blue", alpha=0.8)
ax.set_xlabel(r"Multipole $\ell$")
ax.set_ylabel(r"$\ell(\ell+1)C_\ell / 2\pi$ [$\mu$K$^2$]")
ax.set_xlim(2, lmax)
ax.set_yscale("log")
ax.legend()
ax.set_title(f"CMB Power Spectrum (nside={nside}, mask f_sky={fsky:.2f})")
fig.tight_layout()
plt.savefig("cmb_analysis.png", dpi=150)
print("Saved cmb_analysis.png")