Repository: astropy/astropy

Issue:
Make units be consistent with versioned constants
This motivation for this issue is fixing a subtle but important inconsistency that can now appear because constants are versioned (new feature in v2.0).  The basic problem is that some units are effectively tied to constants.  For example `units.solMass` exists, but so does `constants.M_sun`.  Right now with the default units they are consistent, but *if* someone uses the constant versioning machinery to update e.g. M_sun, they become inconsistent.

There are a couple possible fixes here:
* Tie the units to the constant values so that if a particular constant is set as the "current" version, the unit is assumed to have that value (I think that is more "right" but probably requires subtle and/or major changes to `units`).
* Make unit conversion factors versionable, and provide helper functions or just instructions for how to sync them up (possibly easier, but means manual intervention is required).

Discussion:
mhvk: Interesting! I certainly like that it is more obvious where the default is stored, and that by setting the constants right up front, the units are automatically OK.
stargaser: @mhvk I think with the commit I just pushed modifying `constants/si.py`, it works. I need to put in several tests like this (after I modified my `astropy.cfg` to point to `codata2010` and `iau2012`):

```
In [1]: import astropy.units as u

In [2]: import astropy.constants as const

In [3]: (u.M_earth / const.M_earth).to(u.dimensionless_unscaled)
Out[3]: <Quantity 1.>
```

It will be critical to make sure the `constants.config` is modified before the first import of `units`.
mhvk: That makes sense. It will be equally critical to document this clearly - while I do think it will be a feature that is rarely enough used that we can expect users to jump through some hoops, the hoops should be clear! And it is wonderful that this way it is possible to have consistent constants and units.
mhvk: p.s.  Except...  I guess if we introduce new units based on new constants, we have a bit of a problem. Though it mostly means that some of the units code will have to check the constants configuration.
codecov[bot]: # [Codecov](https://codecov.io/gh/astropy/astropy/pull/8517?src=pr&el=h1) Report
> Merging [#8517](https://codecov.io/gh/astropy/astropy/pull/8517?src=pr&el=desc) into [master](https://codecov.io/gh/astropy/astropy/commit/1508303ccd3c1e93dd13269e91555405a68780c6?src=pr&el=desc) will **decrease** coverage by `<.01%`.
> The diff coverage is `76.92%`.

[![Impacted file tree graph](https://codecov.io/gh/astropy/astropy/pull/8517/graphs/tree.svg?width=650&token=eEhbJxwU0j&height=150&src=pr)](https://codecov.io/gh/astropy/astropy/pull/8517?src=pr&el=tree)

```diff
@@            Coverage Diff            @@
##           master   #8517      +/-   ##
=========================================
- Coverage   86.81%   86.8%   -0.01%     
=========================================
  Files         386     386              
  Lines       58164   58179      +15     
  Branches     1060    1060              
=========================================
+ Hits        50493   50502       +9     
- Misses       7056    7062       +6     
  Partials      615     615
```


| [Impacted Files](https://codecov.io/gh/astropy/astropy/pull/8517?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [astropy/constants/si.py](https://codecov.io/gh/astropy/astropy/pull/8517/diff?src=pr&el=tree#diff-YXN0cm9weS9jb25zdGFudHMvc2kucHk=) | `100% <100%> (ø)` | :arrow_up: |
| [astropy/constants/\_\_init\_\_.py](https://codecov.io/gh/astropy/astropy/pull/8517/diff?src=pr&el=tree#diff-YXN0cm9weS9jb25zdGFudHMvX19pbml0X18ucHk=) | `85.18% <75%> (-9.69%)` | :arrow_down: |
| [astropy/modeling/blackbody.py](https://codecov.io/gh/astropy/astropy/pull/8517/diff?src=pr&el=tree#diff-YXN0cm9weS9tb2RlbGluZy9ibGFja2JvZHkucHk=) | `89.83% <0%> (ø)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/astropy/astropy/pull/8517?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/astropy/astropy/pull/8517?src=pr&el=footer). Last update [1508303...9a5d890](https://codecov.io/gh/astropy/astropy/pull/8517?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

mhvk: @astrofrog - as noted by @stargaser, using a specific set of constants is easy only if it is the first thing you select, before importing all the units that rely on those constants. For this, a configuration item seems the easiest solution. Or can you see how to do it on-the-fly?

Note that this likely will be quite rare, so I think it is OK for there to be some hoops to jump through; in the end, for reproducibility, one would arguably be better of using an environment with fixed versions of all the software used.


stargaser: The code coverage check has passed! A documentation section has been added too. I've removed the WIP annotations -- I think it is ready for @mhvk to take a closer look.
stargaser: For the larger question: I could add a config item "constants_version" which a user could set to 'astropyconst13' or 'astropyconst20'.  I think an exception should be raised if a user sets "constants_version" and "physical_constants" or "astronomical_constants".
mhvk: Alternatively, it could be allowed to set both `physical_constants` and `astronomical_constants` to "astropy20" etc.
stargaser: It gets tricky to implement "constants_version". When a user specifies only `constants_version='astropyconst13'`, the value of `conf.physical_constants` is still `'codata2014'` and the value of `conf.astronomical_constants` is still `'iau2015'`. We don't know if the user put either of those in the config file, or simply did not specify them.

mhvk: So, how about my suggestion of allowing `physical_constants = 'astropy20'`? Then at least the answer should be unique...
stargaser: Oh, got it! Don't add the extra config item, but allow the two items to take on values of `'astropyconst13'` and `'astropyconst20'` (to keep the naming consistent with the context manager), in addition to setting to `'codata2010'` and `'iau2012'`. Yes, that is a good idea. I will go implement that now.
mhvk: @stargaser - the failure on circle-ci looks related to what you have ...
stargaser: Same failure on the doctest for the context manager. There's a very similar test in the module which passes. Unfortunately, I'm having trouble running doctests in my development environment -- `python setup.py test` stops with "path too long" error on my Macbook. I'm at a loss right now for how to fix this.
pllim: There is still one failure on 32-bit
```
______________________ [doctest] docs/constants/index.rst ______________________
[gw0] linux -- Python 3.7.3 /opt/python/cp37-cp37m/bin/python
123 
124     >>> from astropy import constants as const
125     >>> with const.set_enabled_constants('astropyconst13'):
126     ...     print(const.h)
127       Name   = Planck constant
128       Value  = 6.62606957e-34
129       Uncertainty  = 2.9e-41
130       Unit  = J s
131       Reference = CODATA 2010
132     >>> print(const.h)
Differences (unified diff with -expected +actual):
    @@ -1,5 +1,5 @@
       Name   = Planck constant
    -  Value  = 6.62607004e-34
    -  Uncertainty  = 8.1e-42
    +  Value  = 6.62606957e-34
    +  Uncertainty  = 2.9e-41
       Unit  = J s
    -  Reference = CODATA 2014
    +  Reference = CODATA 2010
```
stargaser: Working on it...Looks like I have to debug this in a 32-bit container. I set up a new development environment on a 64-bit Linux machine and all the tests passed with Python 3.7.3.
stargaser: I am running the same container as CircleCI 32-bit and I have installed everything the same. All the tests pass for me. 

I am seeing a problem on all platforms in which  `import astropy.constants as const; print(const.h)` raises ValueErrors with "J is not a valid unit". If I first import astropy.units then everything is fine. I'll try a few more things but I am despairing of getting it all sorted today. 
mhvk: @stargaser - I'll try to have a look, that error is very strange
mhvk: Adding another "trial" import of units right before importing `si` solves it, but somehow I cannot just remove the one on top...
mhvk: @stargaser - I think I have a solution: do the configuration in a separate submodule and import from there in `si` and `cgs`; I'll make a PR to your PR momentarily (or, if you are OK, can just push to your branch).
stargaser: @mhvk -- that's great news! It's fine for you to just push to my branch. I'll add you.
mhvk: @stargaser - well, I may have spoken too quickly. All the configuration tests now seem to fail...
mhvk: OK, now things are OK again...
mhvk: Now pushed to your branch. It solves the `print()` issue - hopefully also the 32 bit one!
mhvk: @stargaser - as long as you allow pushes to your branch (which is the default for astropy), no need to invite me...
stargaser: oh no! I saw the 32-bit build fail and I was trying to restart it -- and it restarted all. Apologies!
bsipocz: > oh no! I saw the 32-bit build fail and I was trying to restart it -- and it restarted all. Apologies!

no worries. To restart only one job you need to select it from the dropdown.
bsipocz: @mhvk - could you use the `merge-when-ci-passes` label when this is ready to go? I'll be online until quite late tonight so will be able to merge stuff.
astrofrog: I’m still concerned that we shouldn’t have settings in the config file that affect the scientific results of a script (and that was the motivation for ScienceState). Is there no way around that in this case?
mhvk: @astrofrog - The trouble is that we need to do something before we load units. In principle, I guess `ScienceState` could be used, if loaded before anything else, and we somehow check that if units is already loaded, it just fails.

p.s. I do not feel this is so hugely pressing that we could not postpone to 4.0; better to have a system we all are happy with! Though I think independently of anything else, we should figure out why the 32-bit build is not working!
stargaser: 32-bit failed again, with what looks like a problem with astropy-helpers. I tested the branch in a 32-bit container on my laptop with the same commands and it worked.

I've started a rerun of the job with ssh enabled -- not sure if I'll be online when it finishes but I thought I'd get it in the queue.
stargaser: I'm okay with pushing this off to 4.0 -- maybe we can figure out better ways to do the imports to make `ScienceState` work.
astrofrog: I think the 32-bit failure can be ignored.

What if the ScienceState lived at the root of the package, in astropy/__init__.py?
stargaser: @astrofrog -- great suggestion! I think putting `ScienceState` in the Astropy root has a good chance of giving us what we want. I think we need something outside both constants and units, because once we fully import one or the other, it is too late to make units consistent with non-default constants.  The configuration system does its thing before those imports, and that is why we were able to get the configuration file to work.

My new goals for this PR are to have a solution using `ScienceState` ready for review close to World Metrology Day on 2019 May 20th, when the new SI takes effect and I am expecting renewed interest in the constants package.
eteq: @stargaser - just wanted to check here: we've been waiting on the release to get the new values, but any progress here?  I.e., is there any chance of "fast-tracking" this one for 3.2, or should we continue to aim this for 4.0?
stargaser: I'm afraid I haven't found time to push this forward lately. I've been assuming this feature is too big to put in the late stages of the 3.2 release cycle.
stargaser: Starting work on putting a ScienceState class for physical constants into `astropy/__init__.py`...
stargaser: I'm having trouble writing tests, because `units` and `constants` have already been imported when pytest gets to my test file. Would @mhvk or @astrofrog have suggestions for testing?
stargaser: The functionality is there. I can do
```
>>> import astropy
>>> astropy.physical_constants.set('codata2010')
<ScienceState physical_constants: 'codata2010'>
>>> astropy.astronomical_constants.set('iau2012')
<ScienceState astronomical_constants: 'iau2012'>
>>> import astropy.units as u
>>> import astropy.constants as const
>>> (const.M_sun / u.M_sun).to(u.dimensionless_unscaled) - 1
<Quantity 0.>
>>> const.M_sun
<<class 'astropy.constants.iau2012.IAU2012'> name='Solar mass' value=1.9891e+30 uncertainty=5e+25 unit='kg' reference="Allen's Astrophysical Quantities 4th Ed.">
```

I can't seem to write that in a working test, because `units` and `constants` are already imported when pytest gets to my test file. Accordingly, the tests are simple.
pllim: CircleCI 32-bit is weird. It appears to be running under your own CircleCI instance, not Astropy's. I think @nden ran into this before and the problem went away when she removed her own CircleCI instance (or something like that).
nden: Yes, the only way I could "correct" this is to close astropy testing on my Circleci account.
mhvk: @pllim is right about the documentation: Probably should have something in `units` as well.
stargaser: Thanks @pllim and @mhvk for all the review comments! I believe I have addressed the code comments. The documentation is still to do. I've fixed the CircleCI problem. Tests are all passing now!

After the docs, I really need to squash commits. I'll save a copy of this branch first.
bsipocz: I'm puzzled:
```
======================================================= FAILURES =======================================================
_____________________________________ [doctest] docs/units/constants_versions.rst ______________________________________
018 >>> import astropy
019 >>> astropy.physical_constants.set('codata2010')
020 <ScienceState physical_constants: 'codata2010'>
021 >>> astropy.astronomical_constants.set('iau2012')
022 <ScienceState astronomical_constants: 'iau2012'>
023 >>> import astropy.units as u
024 >>> import astropy.constants as const
025 >>> (const.M_sun / u.M_sun).to(u.dimensionless_unscaled) - 1   # doctest: +FLOAT_CMP
026 <Quantity 0.>
027 >>> print(const.M_sun)
Differences (unified diff with -expected +actual):
    @@ -1,5 +1,5 @@
       Name   = Solar mass
    -  Value  = 1.9891e+30
    -  Uncertainty  = 5e+25
    +  Value  = 1.9884754153381438e+30
    +  Uncertainty  = 9.236140093538353e+25
       Unit  = kg
    -  Reference = Allen's Astrophysical Quantities 4th Ed.
    +  Reference = IAU 2015 Resolution B 3 + CODATA 2014

/private/var/folders/dc/hsm7tqpx2d57n7vb3k1l81xw0000gq/T/astropy-test-78uso7fm/lib/python3.7/site-packages/docs/units/constants_versions.rst:27: DocTestFailure
```

while if I copy paste the code snippet into ipython I get the same output (from Allen's Astrophysical Quantities 4th Ed.).
bsipocz: @stargaser - could you please rebase?

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 9cc0c28c2d..5d5bf3f541 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -10,6 +10,9 @@ astropy.config
 astropy.constants
 ^^^^^^^^^^^^^^^^^
 
+- The version of constants can be specified via ScienceState in a way
+  that ``constants`` and ``units`` will be consistent. [#8517]
+
 astropy.convolution
 ^^^^^^^^^^^^^^^^^^^
 
diff --git a/astropy/__init__.py b/astropy/__init__.py
index ab8d1736ad..edfea0da94 100644
--- a/astropy/__init__.py
+++ b/astropy/__init__.py
@@ -157,6 +157,80 @@ class Conf(_config.ConfigNamespace):
 
 conf = Conf()
 
+
+# Define a base ScienceState for configuring constants and units
+from .utils.state import ScienceState
+class base_constants_version(ScienceState):
+    """
+    Base class for the real version-setters below
+    """
+    _value = 'test'
+
+    _versions = dict(test='test')
+
+    @classmethod
+    def validate(cls, value):
+        if value not in cls._versions:
+            raise ValueError('Must be one of {}'
+                             .format(list(cls._versions.keys())))
+        return cls._versions[value]
+
+    @classmethod
+    def set(cls, value):
+        """
+        Set the current constants value.
+        """
+        import sys
+        if 'astropy.units' in sys.modules:
+            raise RuntimeError('astropy.units is already imported')
+        if 'astropy.constants' in sys.modules:
+            raise RuntimeError('astropy.constants is already imported')
+
+        class _Context:
+            def __init__(self, parent, value):
+                self._value = value
+                self._parent = parent
+
+            def __enter__(self):
+                pass
+
+            def __exit__(self, type, value, tb):
+                self._parent._value = self._value
+
+            def __repr__(self):
+                return ('<ScienceState {0}: {1!r}>'
+                        .format(self._parent.__name__, self._parent._value))
+
+        ctx = _Context(cls, cls._value)
+        value = cls.validate(value)
+        cls._value = value
+        return ctx
+
+
+class physical_constants(base_constants_version):
+    """
+    The version of physical constants to use
+    """
+    # Maintainers: update when new constants are added
+    _value = 'codata2014'
+
+    _versions = dict(codata2018='codata2018', codata2014='codata2014',
+                     codata2010='codata2010', astropyconst40='codata2018',
+                     astropyconst20='codata2014', astropyconst13='codata2010')
+
+
+class astronomical_constants(base_constants_version):
+    """
+    The version of astronomical constants to use
+    """
+    # Maintainers: update when new constants are added
+    _value = 'iau2015'
+
+    _versions = dict(iau2015='iau2015', iau2012='iau2012',
+                     astropyconst40='iau2015', astropyconst20='iau2015',
+                     astropyconst13='iau2012')
+
+
 # Create the test() function
 from .tests.runner import TestRunner
 test = TestRunner.make_test_runner_in(__path__[0])
@@ -320,7 +394,8 @@ def online_help(query):
 
 __dir_inc__ = ['__version__', '__githash__', '__minimum_numpy_version__',
                '__bibtex__', 'test', 'log', 'find_api_page', 'online_help',
-               'online_docs_root', 'conf']
+               'online_docs_root', 'conf', 'physical_constants',
+               'astronomical_constants']
 
 
 from types import ModuleType as __module_type__
diff --git a/astropy/constants/__init__.py b/astropy/constants/__init__.py
index 65350120f9..7607e642de 100644
--- a/astropy/constants/__init__.py
+++ b/astropy/constants/__init__.py
@@ -13,21 +13,21 @@ A typical use case might be::
     <Quantity 0.510998927603161 MeV>
 
 """
+import warnings
 from contextlib import contextmanager
 
 from astropy.utils import find_current_module
 
 # Hack to make circular imports with units work
-try:
-    from astropy import units
-    del units
-except ImportError:
-    pass
+from astropy import units
+del units
 
+# These lines import some namespaces into the top level
 from .constant import Constant, EMConstant  # noqa
 from . import si  # noqa
 from . import cgs  # noqa
-from . import codata2014, iau2015  # noqa
+from .config import codata, iaudata
+
 from . import utils as _utils
 
 # for updating the constants module docstring
@@ -38,9 +38,11 @@ _lines = [
     '========== ============== ================ =========================',
 ]
 
-# NOTE: Update this when default changes.
-_utils._set_c(codata2014, iau2015, find_current_module(),
-              not_in_module_only=True, doclines=_lines, set_class=True)
+# Catch warnings about "already has a definition in the None system"
+with warnings.catch_warnings():
+    warnings.filterwarnings('ignore', 'Constant .*already has a definition')
+    _utils._set_c(codata, iaudata, find_current_module(),
+                  not_in_module_only=False, doclines=_lines, set_class=True)
 
 _lines.append(_lines[1])
 
@@ -65,38 +67,42 @@ def set_enabled_constants(modname):
     """
 
     # Re-import here because these were deleted from namespace on init.
+    import importlib
     import warnings
     from astropy.utils import find_current_module
     from . import utils as _utils
 
-    # NOTE: Update this when default changes.
-    if modname == 'astropyconst13':
-        from .astropyconst13 import codata2010 as codata
-        from .astropyconst13 import iau2012 as iaudata
-    else:
-        raise ValueError(
-            'Context manager does not currently handle {}'.format(modname))
+    try:
+        modmodule = importlib.import_module('.constants.' + modname, 'astropy')
+        codata_context = modmodule.codata
+        iaudata_context = modmodule.iaudata
+    except ImportError as exc:
+        exc.args += ('Context manager does not currently handle {}'
+                     .format(modname),)
+        raise
 
     module = find_current_module()
 
     # Ignore warnings about "Constant xxx already has a definition..."
     with warnings.catch_warnings():
-        warnings.simplefilter('ignore')
-        _utils._set_c(codata, iaudata, module,
+        warnings.filterwarnings('ignore',
+                                'Constant .*already has a definition')
+        _utils._set_c(codata_context, iaudata_context, module,
                       not_in_module_only=False, set_class=True)
 
     try:
         yield
     finally:
         with warnings.catch_warnings():
-            warnings.simplefilter('ignore')
-            # NOTE: Update this when default changes.
-            _utils._set_c(codata2014, iau2015, module,
+            warnings.filterwarnings('ignore',
+                                    'Constant .*already has a definition')
+            _utils._set_c(codata, iaudata, module,
                           not_in_module_only=False, set_class=True)
 
 
 # Clean up namespace
 del find_current_module
+del warnings
 del contextmanager
 del _utils
 del _lines
diff --git a/astropy/constants/astropyconst13.py b/astropy/constants/astropyconst13.py
index 4e19a0f883..b22fca4ccd 100644
--- a/astropy/constants/astropyconst13.py
+++ b/astropy/constants/astropyconst13.py
@@ -8,7 +8,10 @@ from astropy.utils import find_current_module
 from . import utils as _utils
 from . import codata2010, iau2012
 
-_utils._set_c(codata2010, iau2012, find_current_module())
+codata = codata2010
+iaudata = iau2012
+
+_utils._set_c(codata, iaudata, find_current_module())
 
 # Clean up namespace
 del find_current_module
diff --git a/astropy/constants/astropyconst20.py b/astropy/constants/astropyconst20.py
index 5a04ec7d61..050db5d0af 100644
--- a/astropy/constants/astropyconst20.py
+++ b/astropy/constants/astropyconst20.py
@@ -7,7 +7,11 @@ from astropy.utils import find_current_module
 from . import utils as _utils
 from . import codata2014, iau2015
 
-_utils._set_c(codata2014, iau2015, find_current_module())
+
+codata = codata2014
+iaudata = iau2015
+
+_utils._set_c(codata, iaudata, find_current_module())
 
 # Clean up namespace
 del find_current_module
diff --git a/astropy/constants/cgs.py b/astropy/constants/cgs.py
index adca1b3411..0c7aeb5367 100644
--- a/astropy/constants/cgs.py
+++ b/astropy/constants/cgs.py
@@ -3,14 +3,13 @@
 Astronomical and physics constants in cgs units.  See :mod:`astropy.constants`
 for a complete listing of constants defined in Astropy.
 """
-
 import itertools
 
 from .constant import Constant
-from . import codata2014, iau2015
+from .config import codata, iaudata
 
-for _nm, _c in itertools.chain(sorted(vars(codata2014).items()),
-                               sorted(vars(iau2015).items())):
+for _nm, _c in itertools.chain(sorted(vars(codata).items()),
+                               sorted(vars(iaudata).items())):
     if (isinstance(_c, Constant) and _c.abbrev not in locals()
-         and _c.system in ['esu', 'gauss', 'emu']):
+            and _c.system in ['esu', 'gauss', 'emu']):
         locals()[_c.abbrev] = _c
diff --git a/astropy/constants/config.py b/astropy/constants/config.py
new file mode 100644
index 0000000000..c163b820ec
--- /dev/null
+++ b/astropy/constants/config.py
@@ -0,0 +1,14 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+Configures the codata and iaudata used, possibly using user configuration.
+"""
+# Note: doing this in __init__ causes import problems with units,
+# as si.py and cgs.py have to import the result.
+import importlib
+import astropy
+
+phys_version = astropy.physical_constants.get()
+astro_version = astropy.astronomical_constants.get()
+
+codata = importlib.import_module('.constants.' + phys_version, 'astropy')
+iaudata = importlib.import_module('.constants.' + astro_version, 'astropy')
diff --git a/astropy/constants/si.py b/astropy/constants/si.py
index 65e5a2d13e..c8f83a8411 100644
--- a/astropy/constants/si.py
+++ b/astropy/constants/si.py
@@ -3,16 +3,13 @@
 Astronomical and physics constants in SI units.  See :mod:`astropy.constants`
 for a complete listing of constants defined in Astropy.
 """
-
-
-
 import itertools
 
 from .constant import Constant
-from . import codata2014, iau2015
+from .config import codata, iaudata
 
-for _nm, _c in itertools.chain(sorted(vars(codata2014).items()),
-                               sorted(vars(iau2015).items())):
+for _nm, _c in itertools.chain(sorted(vars(codata).items()),
+                               sorted(vars(iaudata).items())):
     if (isinstance(_c, Constant) and _c.abbrev not in locals()
-         and _c.system == 'si'):
+            and _c.system == 'si'):
         locals()[_c.abbrev] = _c
diff --git a/docs/constants/index.rst b/docs/constants/index.rst
index de4e4709f2..d31aa1656a 100644
--- a/docs/constants/index.rst
+++ b/docs/constants/index.rst
@@ -118,9 +118,57 @@ Union (IAU) are collected in modules with names like ``iau2012`` or ``iau2015``:
       Reference = IAU 2015 Resolution B 3
 
 The astronomical and physical constants are combined into modules with
-names like ``astropyconst13`` and ``astropyconst20``. To temporarily set
-constants to an older version (e.g., for regression testing), a context
-manager is available as follows:
+names like ``astropyconst13`` and ``astropyconst20`` for different versions.
+However, importing these prior version modules directly will lead to
+inconsistencies with other subpackages that have already imported
+`astropy.constants`. Notably, `astropy.units` will have already used
+the default version of constants. When using prior versions of the constants
+in this manner, quantities should be constructed with constants instead of units.
+
+To ensure consistent use of a prior version of constants in other Astropy
+packages (such as `astropy.units`) that import constants, the physical and
+astronomical constants versions should be set via ScienceState classes.
+These must be set before the first import of either `astropy.constants` or
+`astropy.units`.  For example, you can use the CODATA2010 physical constants and the
+IAU 2012 astronomical constants:
+
+    >>> from astropy import physical_constants, astronomical_constants
+    >>> physical_constants.set('codata2010')  # doctest: +SKIP
+    <ScienceState physical_constants: 'codata2010'>
+    >>> physical_constants.get()  # doctest: +SKIP
+    'codata2010'
+    >>> astronomical_constants.set('iau2012')  # doctest: +SKIP
+    <ScienceState astronomical_constants: 'iau2012'>
+    >>> astronomical_constants.get()  # doctest: +SKIP
+    'iau2012'
+
+Then all other packages that import `astropy.constants` will self-consistently
+initialize with that prior version of constants.
+
+The versions may also be set using values referring to the version modules:
+
+    >>> from astropy import physical_constants, astronomical_constants
+    >>> physical_constants.set('astropyconst13')  # doctest: +SKIP
+    <ScienceState physical_constants: 'codata2010'>
+    >>> physical_constants.get()  # doctest: +SKIP
+    'codata2010'
+    >>> astronomical_constants.set('astropyconst13')  # doctest: +SKIP
+    <ScienceState astronomical_constants: 'iau2012'>
+    >>> astronomical_constants.get()  # doctest: +SKIP
+    'iau2012'
+
+If either `astropy.constants` or `astropy.units` have already been imported, a
+``RuntimeError`` will be raised.
+
+    >>> import astropy.units
+    >>> from astropy import physical_constants, astronomical_constants
+    >>> astronomical_constants.set('astropyconst13')
+    Traceback (most recent call last):
+        ...
+    RuntimeError: astropy.units is already imported
+
+To temporarily set constants to an older version (e.g.,
+for regression testing), a context manager is available, as follows:
 
     >>> from astropy import constants as const
     >>> with const.set_enabled_constants('astropyconst13'):
@@ -137,11 +185,10 @@ manager is available as follows:
       Unit  = J s
       Reference = CODATA 2014
 
-.. warning::
+The context manager may be used at any time in a Python session, but it
+uses the prior version only for `astropy.constants`, and not for any
+other subpackage such as `astropy.units`.
 
-    Units such as ``u.M_sun`` will use the current version of the
-    corresponding constant. When using prior versions of the constants,
-    quantities should be constructed with constants instead of units.
 
 .. note that if this section gets too long, it should be moved to a separate
    doc page - see the top of performance.inc.rst for the instructions on how to
diff --git a/docs/units/constants_versions.rst b/docs/units/constants_versions.rst
new file mode 100644
index 0000000000..4663b3257f
--- /dev/null
+++ b/docs/units/constants_versions.rst
@@ -0,0 +1,29 @@
+Using prior versions of constants
+*********************************
+
+By default, `astropy.units` are initialized upon first import to use
+the current versions of `astropy.constants`. For units to initialize
+properly to a prior version of constants, the constants versions must
+be set before the first import of `astropy.units` or `astropy.constants`.
+
+This is accomplished using ScienceState classes in the top-level package.
+Setting the prior versions at the start of a Python session will allow
+consistent units, as follows:
+
+>>> import astropy
+>>> astropy.physical_constants.set('codata2010')  # doctest: +SKIP
+<ScienceState physical_constants: 'codata2010'>
+>>> astropy.astronomical_constants.set('iau2012')  # doctest: +SKIP
+<ScienceState astronomical_constants: 'iau2012'>
+>>> import astropy.units as u
+>>> import astropy.constants as const
+>>> (const.M_sun / u.M_sun).to(u.dimensionless_unscaled) - 1  # doctest: +SKIP
+<Quantity 0.>
+>>> const.M_sun  # doctest: +SKIP
+  Name   = Solar mass
+  Value  = 1.9891e+30
+  Uncertainty  = 5e+25
+  Unit  = kg
+  Reference = Allen's Astrophysical Quantities 4th Ed.
+
+If `astropy.units` has already been imported, a RuntimeError is raised.
diff --git a/docs/units/index.rst b/docs/units/index.rst
index 4267de4234..82a7d412d6 100644
--- a/docs/units/index.rst
+++ b/docs/units/index.rst
@@ -165,6 +165,7 @@ Using `astropy.units`
    logarithmic_units
    format
    equivalencies
+   constants_versions
    conversion
 
 See Also


Test patch:
diff --git a/astropy/constants/tests/test_prior_version.py b/astropy/constants/tests/test_prior_version.py
index 5d13c66a7b..c7f54ef75b 100644
--- a/astropy/constants/tests/test_prior_version.py
+++ b/astropy/constants/tests/test_prior_version.py
@@ -163,6 +163,6 @@ def test_context_manager():
 
     assert const.h.value == 6.626070040e-34  # CODATA2014
 
-    with pytest.raises(ValueError):
+    with pytest.raises(ImportError):
         with const.set_enabled_constants('notreal'):
             const.h
diff --git a/astropy/constants/tests/test_sciencestate.py b/astropy/constants/tests/test_sciencestate.py
new file mode 100644
index 0000000000..1fa32b6909
--- /dev/null
+++ b/astropy/constants/tests/test_sciencestate.py
@@ -0,0 +1,21 @@
+import pytest
+
+from astropy import physical_constants, astronomical_constants
+import astropy.constants as const
+
+
+def test_version_match():
+    pversion = physical_constants.get()
+    refpversion = const.h.__class__.__name__.lower()
+    assert pversion == refpversion
+    aversion = astronomical_constants.get()
+    refaversion = const.M_sun.__class__.__name__.lower()
+    assert aversion == refaversion
+
+
+def test_previously_imported():
+    with pytest.raises(RuntimeError):
+        physical_constants.set('codata2018')
+
+    with pytest.raises(RuntimeError):
+        astronomical_constants.set('iau2015')