Add cutoff blackbody model class - #20560

Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

• Do the proposed changes actually accomplish desired goals?
• Do the proposed changes follow the Astropy coding guidelines?
• Are tests added/updated as required? If so, do they follow the Astropy testing guidelines?
• Are docs added/updated as required? If so, do they follow the Astropy documentation guidelines?
• Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for rebase and squash.
• Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the bot.
• Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
• Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
• At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) before merge.

Hello and thanks.

scipy vs numpy

Is there any performance boost either way? I think that might be the main factor. Otherwise, if scipy is preferred, I think it would be acceptable to use but you will have to be more careful in handling "what if scipy is not installed" scenario.

I'd like to avoid hardcoding a discretized grid for this integration, which means that scipy vs. numpy implementations are going to be different and thus not directly comparable. Scipy makes it really easy with scipy.integrate.quad which just takes the function itself and internally chooses the evaluation points adaptively. With numpy I think the way to do it would be to use Legendre-Gauss quadrature via np.polynomial.legendre.leggauss for the nodes/weights. This does still require choosing a quadrature order (i.e., number of nodes), but we don't have to define the actual physical grid spacing or an arbitrary upper integration bound.

In practice, the two methods agree to within very small tolerances (max rel. error of ~1e-9 (1e-14 typ.) for N=64 nodes and ~1e-10 (1e-15 typ.) for N=128 nodes). The Gauss-Legendre (numpy) way is faster by a fair margin for both N=64 and N=128. Does that make it your preference, at the cost of probably slightly less obvious code/implementation?

I vaguely remember @karllark being passionate about our blackbody modeling, so maybe he has opinions. Also cc @astropy/modeling

I do not have an opinion on how the integral is computed. Thanks for thinking of me though.

Here's a scipy implementation for review. My thinking is that this is less obscure (thus easier to maintain), and wringing every ounce of performance out of this method isn't critical since it's not called at all during fitting or anything like that.

Sorry for the CI churn. I think this should basically be good to go now.