compare_model_spectra_ASPCAP.ipynb Enhancement

[JosephKarpinski]

https://github.com/sdss/dr20_tutorials/blob/main/compare_model_spectra_ASPCAP.ipynb

Add enhancement code:

------------------------------------------------------------

NEW CELL: Auto-detect each element's diagnostic window + zoomed residual plot

------------------------------------------------------------

'''
Addresses the notebook's own stated limitation (Cell 16): the full-width
plots hide line-level detail. Rather than guessing exact line-center
wavelengths from memory, this finds where each element's model_flux
differs most from the full best-fit model -- those pixels are, by
construction, the ASPCAP fitting window for that element -- and zooms
into a region around them, with a residual panel underneath.
'''
full_model = aspcapStar["model_flux"][0]
window_half_width = 15 # Angstroms of context around the detected window

for elem in elems:
elem_model = aspcapStar[f"model_flux_{elem}_h"][0]
diff = np.abs(elem_model - full_model)

if not np.any(np.isfinite(diff)) or np.nanmax(diff) == 0:
    print(f"{elem}: no detectable difference from full model, skipping")
    continue

# locate the strongest single diagnostic feature for this element
peak_idx = np.nanargmax(diff)
center_wl = wavelength[peak_idx]
window = (wavelength > center_wl - window_half_width) & (wavelength < center_wl + window_half_width)
obs_window = (obs_wl > center_wl - window_half_width) & (obs_wl < center_wl + window_half_width) & mask

fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(9, 6), sharex=True,
                                gridspec_kw={'height_ratios': [3, 1]})
ax1.plot(obs_wl[obs_window], obs_flux[obs_window], color='k', label='observed', lw=1.2)
ax1.plot(wavelength[window], (full_model * continuum)[window], '--', label='full best-fit model')
ax1.plot(wavelength[window], (elem_model * continuum)[window], '--',
          label=f"{elem.title()} model")
ax1.axvline(center_wl, color='gray', lw=0.5, ls=':')
ax1.set_ylabel("flux")
ax1.set_title(fr"{elem.title()}: [{elem.title()}/H]={aspcapStar[f'{elem}_h'][0]:.3f} "
              fr"$\pm$ {aspcapStar[f'e_{elem}_h'][0]:.3f}, zoomed on strongest diagnostic pixel")
ax1.legend(fontsize=8)

# residual: element model vs full model, in normalized (continuum) units
ax2.plot(wavelength[window], (elem_model - full_model)[window], color='firebrick')
ax2.axhline(0, color='k', lw=0.7, ls='--')
ax2.set_ylabel("elem - full\nmodel")
ax2.set_xlabel("λ [Å]")
plt.tight_layout()
plt.show()

[JosephKarpinski]

Update:

ASPCAP employs a two-step process to extract abundances: first, the stellar atmospheric parameters are determined by fitting the entire APOGEE spectrum; and second, these parameters are then used to derive an individual element abundance by fitting over limited-wavelength regions of the spectrum (windows) that are dominated by spectral features associated with that element, and the DR16 spectroscopic stellar parameters (i.e., the FPARAM parameter array) are employed in the determination of the individual element abundances.
SDSS
This is exactly the mechanism that explains everything you've found:
model_flux (the "full" model) is the synthetic spectrum from the global fit — Teff, logg, [M/H], vturb, and the star's full multi-element abundance pattern all optimized together.
model_flux_mn_h, model_flux_na_h, etc. are each generated in the second-stage, per-element pass, where the stellar parameters (FPARAM) are held fixed from stage 1, but that element's abundance is re-optimized using only the pixels in its window — while every other element in that resynthesis is typically left at some baseline (often the overall [α/M]/[M/H] scaled-solar pattern), not at the star's actual individually-measured 20-element pattern.
That means model_flux_mn_h differs from model_flux everywhere the star's true abundance pattern departs from scaled-solar — every line of every element, not just Mn's — because only Mn was actually "unlocked" relative to the scaled-solar baseline in that resynthesis, while model_flux has everything unlocked. That's precisely why:
The diff is broadly spread across the whole spectrum rather than localized to each element's own lines.
Mn's and Na's model arrays disagree with each other almost as much as either disagrees with model_flux — because both differ from the "true" pattern in similar, correlated ways.
The single global nanargmax keeps landing on the same steep, abundance-sensitive feature (~15901 Å) regardless of which element column you feed it, since that's just wherever the scaled-solar vs. true-pattern mismatch happens to be largest overall.
Practical takeaway: the diff = elem_model - full_model approach can't be used to auto-locate each element's diagnostic window — it's not measuring "this element's line strength," it's measuring "how much this star's true abundances deviate from the second-stage baseline pattern," which is dominated by whatever's most metallicity-sensitive in general.

[johndonor3]

Unfortunately the time I can allocate for tutorials is very limited and enhancements are outside that scope until DR21, but feel free to send a pull request!