diff --git a/docs/sphinx/source/whatsnew/v0.6.0.rst b/docs/sphinx/source/whatsnew/v0.6.0.rst
index dd5ff4208b..9dc9a13313 100644
--- a/docs/sphinx/source/whatsnew/v0.6.0.rst
+++ b/docs/sphinx/source/whatsnew/v0.6.0.rst
@@ -125,6 +125,7 @@ Bug fixes
   by lonlat bounding box (:issue:`521`)
 * Fix issue with unbounded clearness index calculation in disc. (:issue:`540`)
 * Limit pvwatts_ac results to be greater than or equal to 0. (:issue:`541`)
+* Fix bug in get_relative_airmass(model='youngirvine1967'). (:issue:`545`)
 
 
 Documentation
@@ -155,6 +156,7 @@ Testing
 * Use pytest-mock to ensure that ModelChain DC model is set up correctly.
 * Add Python 3.7 to build matrix
 * Make test_forecast.py more robust. (:issue:`293`)
+* Improve test_atmosphere.py. (:issue:`158`)
 
 
 Contributors
diff --git a/pvlib/atmosphere.py b/pvlib/atmosphere.py
index 66a1aed11e..fb8775a190 100644
--- a/pvlib/atmosphere.py
+++ b/pvlib/atmosphere.py
@@ -229,8 +229,8 @@ def get_relative_airmass(zenith, model='kastenyoung1989'):
         am = (1.0 / (np.sin(np.radians(90 - z +
               244.0 / (165 + 47.0 * (90 - z) ** 1.1)))))
     elif 'youngirvine1967' == model:
-        am = ((1.0 / np.cos(zenith_rad)) *
-              (1 - 0.0012*((1.0 / np.cos(zenith_rad)) ** 2) - 1))
+        sec_zen = 1.0 / np.cos(zenith_rad)
+        am = sec_zen * (1 - 0.0012 * (sec_zen * sec_zen - 1))
     elif 'young1994' == model:
         am = ((1.002432*((np.cos(zenith_rad)) ** 2) +
               0.148386*(np.cos(zenith_rad)) + 0.0096467) /
diff --git a/pvlib/test/test_atmosphere.py b/pvlib/test/test_atmosphere.py
index 17cf227f94..b77d565c9a 100644
--- a/pvlib/test/test_atmosphere.py
+++ b/pvlib/test/test_atmosphere.py
@@ -1,69 +1,77 @@
 import itertools
 
 import numpy as np
+from numpy import nan
 from numpy.testing import assert_allclose
 import pandas as pd
+from pandas.util.testing import assert_series_equal
 import pytest
 
 from pvlib import atmosphere
-from pvlib import solarposition
 from pvlib._deprecation import pvlibDeprecationWarning
 
 from conftest import fail_on_pvlib_version
 
-latitude, longitude, tz, altitude = 32.2, -111, 'US/Arizona', 700
-
-times = pd.date_range(start='20140626', end='20140626', freq='6h', tz=tz)
-
-ephem_data = solarposition.get_solarposition(times, latitude, longitude)
-
-
-# need to add physical tests instead of just functional tests
 
 def test_pres2alt():
-    atmosphere.pres2alt(100000)
-
-
-def test_alt2press():
-    atmosphere.pres2alt(1000)
-
-
-@pytest.mark.parametrize("model",
-    ['simple', 'kasten1966', 'youngirvine1967', 'kastenyoung1989',
-     'gueymard1993', 'young1994', 'pickering2002'])
-def test_airmass(model):
-    out = atmosphere.get_relative_airmass(ephem_data['zenith'], model)
-    assert isinstance(out, pd.Series)
-    out = atmosphere.get_relative_airmass(ephem_data['zenith'].values, model)
-    assert isinstance(out, np.ndarray)
+    out = atmosphere.pres2alt(np.array([10000, 90000, 101325]))
+    expected = np.array([15797.638, 988.637, 0.124])
+    assert_allclose(out, expected, atol=0.001)
+
+
+def test_alt2pres():
+    out = atmosphere.alt2pres(np.array([-100, 0, 1000, 8000]))
+    expected = np.array([102532.073, 101324.999,  89874.750,  35600.496])
+    assert_allclose(out, expected, atol=0.001)
+
+
+@pytest.fixture
+def zeniths():
+    return np.array([100, 89.9, 80, 0])
+
+
+@pytest.mark.parametrize("model,expected",
+                         [['simple', [nan, 572.958,   5.759,   1.000]],
+                          ['kasten1966', [nan, 35.365,  5.580,  0.999]],
+                          ['youngirvine1967', [
+                                 nan, -2.251358367165932e+05, 5.5365, 1.0000]],
+                          ['kastenyoung1989', [nan, 36.467,  5.586,  1.000]],
+                          ['gueymard1993', [nan, 36.431,  5.581,  1.000]],
+                          ['young1994', [nan, 30.733,  5.541,  1.000]],
+                          ['pickering2002', [nan, 37.064,  5.581,  1.000]]])
+def test_airmass(model, expected, zeniths):
+    out = atmosphere.get_relative_airmass(zeniths, model)
+    expected = np.array(expected)
+    assert_allclose(out, expected, equal_nan=True, atol=0.001)
+    # test series in/out. index does not matter
+    # hits the isinstance() block in get_relative_airmass
+    times = pd.DatetimeIndex(start='20180101', periods=len(zeniths), freq='1s')
+    zeniths = pd.Series(zeniths, index=times)
+    expected = pd.Series(expected, index=times)
+    out = atmosphere.get_relative_airmass(zeniths, model)
+    assert_series_equal(out, expected, check_less_precise=True)
 
 
 def test_airmass_scalar():
     assert not np.isnan(atmosphere.get_relative_airmass(10))
 
 
-def test_airmass_scalar_nan():
-    assert np.isnan(atmosphere.get_relative_airmass(100))
-
-
 def test_airmass_invalid():
     with pytest.raises(ValueError):
-        atmosphere.get_relative_airmass(ephem_data['zenith'], 'invalid')
+        atmosphere.get_relative_airmass(0, 'invalid')
 
 
 def test_get_absolute_airmass():
-    relative_am = atmosphere.get_relative_airmass(ephem_data['zenith'],
-                                                  'simple')
-    atmosphere.get_absolute_airmass(relative_am)
-    atmosphere.get_absolute_airmass(relative_am, pressure=100000)
-
-
-def test_get_absolute_airmass_numeric():
-    atmosphere.get_absolute_airmass(2)
-
-
-def test_get_absolute_airmass_nan():
-    np.testing.assert_equal(np.nan, atmosphere.get_absolute_airmass(np.nan))
+    # input am
+    relative_am = np.array([nan, 40, 2, .999])
+    # call without pressure kwarg
+    out = atmosphere.get_absolute_airmass(relative_am)
+    expected = np.array([nan, 40., 2., 0.999])
+    assert_allclose(out, expected, equal_nan=True, atol=0.001)
+    # call with pressure kwarg
+    out = atmosphere.get_absolute_airmass(relative_am, pressure=90000)
+    expected = np.array([nan, 35.529, 1.776, 0.887])
+    assert_allclose(out, expected, equal_nan=True, atol=0.001)
 
 
 @fail_on_pvlib_version('0.7')