@@ -88,6 +88,12 @@ API Changes
 Bug Fixes
 ^^^^^^^^^
 
+- ``astropy.analytic_functions``
+
+  - Fixed the blackbody functions' handling of overflows on some platforms
+    (Windows with MSVC, older Linux versions) with a buggy ``expm1`` function.
+    [#4393]
+
 - ``astropy.config``
 
 - ``astropy.constants``
 
@@ -22,6 +22,14 @@
 FLAM = u.erg / (u.cm**2 * u.s * u.AA)
 
 
+# Some platform implementations of expm1() are buggy and Numpy uses
+# them anyways--the bug is that on certain large inputs it returns
+# NaN instead of INF like it should (it should only return NaN on a
+# NaN input
+# See https://github.com/astropy/astropy/issues/4171
+_has_buggy_expm1 = np.isnan(np.expm1(1000))
+
+
 def blackbody_nu(in_x, temperature):
     """Calculate blackbody flux per steradian, :math:`B_{\\nu}(T)`.
 
@@ -70,9 +78,24 @@ def blackbody_nu(in_x, temperature):
         warnings.warn('Input contains invalid wavelength/frequency value(s)',
                       AstropyUserWarning)
 
+    log_boltz = const.h * freq / (const.k_B * temp)
+    boltzm1 = np.expm1(log_boltz)
+
+    if _has_buggy_expm1:
+        # Replace incorrect nan results with infs--any result of 'nan' is
+        # incorrect unless the input (in log_boltz) happened to be nan to begin
+        # with.  (As noted in #4393 ideally this would be replaced by a version
+        # of expm1 that doesn't have this bug, rather than fixing incorrect
+        # results after the fact...)
+        boltzm1_nans = np.isnan(boltzm1)
+        if np.any(boltzm1_nans):
+            if boltzm1.isscalar and not np.isnan(log_boltz):
+                boltzm1 = np.inf
+            else:
+                boltzm1[np.where(~np.isnan(log_boltz) & boltzm1_nans)] = np.inf
+
     # Calculate blackbody flux
-    bb_nu = (2.0 * const.h * freq ** 3 /
-             (const.c ** 2 * np.expm1(const.h * freq / (const.k_B * temp))))
+    bb_nu = (2.0 * const.h * freq ** 3 / (const.c ** 2 * boltzm1))
     flux = bb_nu.to(FNU, u.spectral_density(freq))
 
     return flux / u.sr  # Add per steradian to output flux unit
 
@@ -3,7 +3,6 @@
 from __future__ import (absolute_import, division, print_function,
                         unicode_literals)
 
-import os
 import warnings
 # THIRD-PARTY
 import numpy as np
@@ -26,7 +25,6 @@
 __doctest_skip__ = ['*']
 
 
-@pytest.mark.skipif("os.environ.get('APPVEYOR')",  reason="fails on AppVeyor")
 @pytest.mark.skipif('not HAS_SCIPY')
 def test_blackbody_scipy():
     """Test Planck function.
@@ -47,7 +45,6 @@ def test_blackbody_scipy():
     np.testing.assert_allclose(intflux, ans.value, rtol=0.01)  # 1% accuracy
 
 
-@pytest.mark.skipif("os.environ.get('APPVEYOR')",  reason="fails on AppVeyor")
 def test_blackbody_overflow():
     """Test Planck function with overflow."""
     photlam = u.photon / (u.cm**2 * u.s * u.AA)
@@ -69,7 +66,6 @@ def test_blackbody_overflow():
     assert flux.value == 0
 
 
-@pytest.mark.skipif("os.environ.get('APPVEYOR')",  reason="fails on AppVeyor")
 def test_blackbody_synphot():
     """Test that it is consistent with IRAF SYNPHOT BBFUNC."""
     # Solid angle of solar radius at 1 kpc