corsika/framework/stack/StackIteratorInterface.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/stack/ParticleBase.hpp>
+
+#include <boost/type_index.hpp>
+
+namespace corsika::history {
+  template <typename T, template <typename> typename TParticleInterface>
+  class HistorySecondaryProducer; // forward decl.
+}
+
+namespace corsika {
+
+  template <typename TStackData, template <typename> typename TParticleInterface,
+            template <class T1, template <class> class T2> class MSecondaryProducer>
+  class Stack; // forward decl
+
+  template <typename TStackData, template <typename> typename TParticleInterface,
+            template <class T1, template <class> class T2> class MSecondaryProducer>
+  class SecondaryView; // forward decl
+
+  template <typename TStackData, template <typename> typename TParticleInterface,
+            template <class T1, template <class> class T2> class MSecondaryProducer,
+            typename StackType>
+  class ConstStackIteratorInterface; // forward decl
+
+  /**
+   * The StackIteratorInterface is the main interface to iterator over
+   * particles on a stack.
+   *
+   * At the same time StackIteratorInterface is a
+   * Particle object by itself, thus there is no difference between
+   * type and ref_type for convenience of the physicist.
+   *
+   * This allows to write code like:
+   * @code
+   * for (auto& p : theStack) {
+   *    p.setEnergy(newEnergy);
+   * }
+   * @endcode
+   *
+   * The template argument Stack determines the type of Stack object
+   * the data is stored in. A pointer to the Stack object is part of
+   * the StackIteratorInterface. In addition to Stack the iterator only knows
+   * the index index_ in the Stack data.
+   *
+   * The template argument `TParticleInterface` acts as a policy to provide
+   * readout function of Particle data from the stack. The TParticleInterface
+   * class must know how to retrieve information from the Stack data
+   * for a particle entry at any index index_.
+   *
+   * The TParticleInterface class must be written and provided by the
+   * user, it contains methods like <code> auto getData() const {
+   * return getStackData().getData(getIndex()); }</code>, where
+   * StackIteratorInterface::getStackData() return a reference to the
+   * object storing the particle data of type TStackData. And
+   * StackIteratorInterface::getIndex() provides the iterator index to
+   * be readout. The TStackData is another user-provided class to
+   * store data and must implement functions compatible with
+   * TParticleInterface, in this example TStackData::getData(const unsigned int
+   * vIndex).
+   *
+   * For two examples see stack_example.cc, or the
+   * corsika::sibyll::SibStack class.
+   */
+
+  template <typename TStackData, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer,
+            typename TStackType =
+                Stack<TStackData, TParticleInterface, MSecondaryProducer>>
+  class StackIteratorInterface
+      : public TParticleInterface<StackIteratorInterface<
+            TStackData, TParticleInterface, MSecondaryProducer, TStackType>> {
+
+  public:
+    typedef TParticleInterface<corsika::StackIteratorInterface<
+        TStackData, TParticleInterface, MSecondaryProducer, TStackType>>
+        particle_interface_type;
+
+    typedef TStackType stack_type;
+    typedef TStackData stack_data_type;
+
+    // it is not allowed to create a "dangling" stack iterator
+    StackIteratorInterface() = delete; //! \todo check rule of five
+
+    StackIteratorInterface(StackIteratorInterface&& rhs)
+        : index_(std::move(rhs.index_))
+        , data_(std::move(rhs.data_)) {}
+
+    StackIteratorInterface(StackIteratorInterface const& vR)
+        : index_(vR.index_)
+        , data_(vR.data_) {}
+
+    StackIteratorInterface& operator=(StackIteratorInterface const& vR) {
+      if (&vR != this) {
+        index_ = vR.index_;
+        data_ = vR.data_;
+      }
+      return *this;
+    }
+
+    /**
+     *  Iterator must always point to data, with an index.
+     *
+     *    @param data reference to the stack [rw]
+     *    @param index index on stack
+     */
+    StackIteratorInterface(stack_type& data, unsigned int const index)
+        : index_(index)
+        , data_(&data) {}
+
+    /**
+     *  Constructor that also sets new values on particle data object.
+     *
+     *  @param data reference to the stack [rw].
+     *  @param index index on stack.
+     *  @param args variadic list of data to initialize stack entry, this must be
+     *        consistent with the definition of the user-provided
+     *        particle_interface_type::setParticleData(...) function.
+     */
+    template <typename... TArgs>
+    StackIteratorInterface(stack_type& data, unsigned int const index,
+                           const TArgs... args)
+        : index_(index)
+        , data_(&data) {
+
+      (**this).setParticleData(args...);
+    }
+
+    /**
+     * Constructor that also sets new values on particle data object, including reference
+     * to parent particle.
+     *
+     *  @param data reference to the stack [rw]
+     *  @param index index on stack
+     *  @param parent to parent particle [rw]. This can be used for thinning, particle
+     *         counting, history, etc.
+     *   @param args variadic list of data to initialize stack entry, this must be
+     *        consistent with the definition of the user-provided
+     *         particle_interface_type::setParticleData(...) function.
+     */
+    template <typename... TArgs>
+    StackIteratorInterface(stack_type& data, unsigned int const index,
+                           StackIteratorInterface& parent, const TArgs... args)
+        : index_(index)
+        , data_(&data) {
+
+      (**this).setParticleData(*parent, args...);
+    }
+
+    bool isErased() const { return getStack().isErased(*this); }
+
+  public:
+    /**
+     * @name Iterator interface
+     * @{
+     */
+    StackIteratorInterface& operator++() {
+      do {
+        ++index_;
+      } while (
+          getStack().isErased(*this)); // this also check the allowed bounds of index_
+      return *this;
+    }
+    StackIteratorInterface operator++(int) {
+      StackIteratorInterface tmp(*this);
+      do {
+        ++index_;
+      } while (
+          getStack().isErased(*this)); // this also check the allowed bounds of index_
+      return tmp;
+    }
+    StackIteratorInterface operator+(int delta) const {
+      return StackIteratorInterface(*data_, index_ + delta);
+    }
+    bool operator==(StackIteratorInterface const& rhs) const {
+      return index_ == rhs.index_;
+    }
+    bool operator!=(StackIteratorInterface const& rhs) const {
+      return index_ != rhs.index_;
+    }
+    bool operator==(const ConstStackIteratorInterface<
+                    TStackData, TParticleInterface, MSecondaryProducer,
+                    stack_type>& rhs) const; // implemented below
+    bool operator!=(const ConstStackIteratorInterface<
+                    TStackData, TParticleInterface, MSecondaryProducer,
+                    stack_type>& rhs) const; // implemented below
+
+    /**
+     * Convert iterator to value type, where value type is the user-provided particle
+     * readout class.
+     */
+    particle_interface_type& operator*() {
+      return static_cast<particle_interface_type&>(*this);
+    }
+
+    /**
+     * Convert iterator to const value type, where value type is the user-provided
+     * particle readout class.
+     */
+    particle_interface_type const& operator*() const {
+      return static_cast<particle_interface_type const&>(*this);
+    }
+    //! @}
+
+  protected:
+    /**
+     * @name Stack data access
+     * @{
+     */
+    /// Get current particle index
+    unsigned int getIndex() const { return index_; }
+    /// Get current particle Stack object
+    stack_type& getStack() { return *data_; }
+    /// Get current particle const Stack object
+    stack_type const& getStack() const { return *data_; }
+    /// Get current user particle TStackData object
+    TStackData& getStackData() { return data_->getStackData(); }
+    /// Get current const user particle TStackData object
+    TStackData const& getStackData() const { return data_->getStackData(); }
+    /// Get data index as mapped in Stack class
+    unsigned int getIndexFromIterator() const {
+      return data_->getIndexFromIterator(index_);
+    }
+    //! @}
+
+    // friends are needed for access to protected methods
+    friend class Stack<TStackData, TParticleInterface,
+                       MSecondaryProducer>; // for access to getIndex for Stack
+    friend class Stack<TStackData&, TParticleInterface,
+                       MSecondaryProducer>;            // for access to getIndex
+                                                       // SecondaryView : public Stack
+    friend class ParticleBase<StackIteratorInterface>; // for access to getStackDataType
+
+    template <typename T1,                     // best fix this to: TStackData,
+              template <typename> typename M1, // best fix this to: TParticleInterface,
+              template <typename T, template <typename> typename T3> typename M2>
+    friend class SecondaryView; // access grant for SecondaryView
+
+    template <typename T, template <typename> typename TParticleInterface_>
+    friend class corsika::history::HistorySecondaryProducer;
+
+    friend class ConstStackIteratorInterface<TStackData, TParticleInterface,
+                                             MSecondaryProducer, stack_type>;
+
+  protected:
+    unsigned int index_ = 0;
+
+  private:
+    stack_type* data_ = 0; // info: Particles and StackIterators become invalid when
+                           // parent Stack is copied or deleted!
+
+  }; // end class StackIterator
+
+  /**
+   * This is the iterator class for const-access to stack data.
+   *
+   * The const counterpart of StackIteratorInterface, which is used
+   * for read-only iterator access on particle stack:
+   *
+   * @code
+   * for (auto const& p : theStack) { E += p.getEnergy(); }
+   * @endcode
+   *
+   * See documentation of StackIteratorInterface for more details:
+   * \sa StackIteratorInterface.
+   */
+
+  template <typename TStackData, template <typename> typename TParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer,
+            typename TStackType =
+                Stack<TStackData, TParticleInterface, MSecondaryProducer>>
+  class ConstStackIteratorInterface
+      : public TParticleInterface<ConstStackIteratorInterface<
+            TStackData, TParticleInterface, MSecondaryProducer, TStackType>> {
+
+  public:
+    typedef TParticleInterface<ConstStackIteratorInterface<
+        TStackData, TParticleInterface, MSecondaryProducer, TStackType>>
+        particle_interface_type;
+
+    typedef TStackType stack_type;
+    typedef TStackData stack_data_type;
+
+    // we don't want to allow dangling iterators to exist
+    ConstStackIteratorInterface() = delete; //! \todo check rule of five
+
+  public:
+    ConstStackIteratorInterface(ConstStackIteratorInterface&& rhs)
+        : index_(std::move(rhs.index_))
+        , data_(std::move(rhs.data_)) {}
+
+    ConstStackIteratorInterface(stack_type const& data, unsigned int const index)
+        : index_(index)
+        , data_(&data) {}
+
+    bool isErased() const { return getStack().isErased(*this); }
+
+    /**
+     * @name Iterator interface
+     * @{
+     */
+    ConstStackIteratorInterface& operator++() {
+      do {
+        ++index_;
+      } while (
+          getStack().isErased(*this)); // this also check the allowed bounds of index_
+      return *this;
+    }
+    ConstStackIteratorInterface operator++(int) {
+      ConstStackIteratorInterface tmp(*this);
+      do {
+        ++index_;
+      } while (
+          getStack().isErased(*this)); // this also check the allowed bounds of index_
+      return tmp;
+    }
+    ConstStackIteratorInterface operator+(int const delta) const {
+      return ConstStackIteratorInterface(*data_, index_ + delta);
+    }
+    bool operator==(ConstStackIteratorInterface const& rhs) const {
+      return index_ == rhs.index_;
+    }
+    bool operator!=(ConstStackIteratorInterface const& rhs) const {
+      return index_ != rhs.index_;
+    }
+    bool operator==(
+        StackIteratorInterface<stack_data_type, TParticleInterface, MSecondaryProducer,
+                               stack_type> const& rhs) const {
+      return index_ == rhs.index_;
+    }
+    bool operator!=(
+        StackIteratorInterface<stack_data_type, TParticleInterface, MSecondaryProducer,
+                               stack_type> const& rhs) const {
+      return index_ != rhs.index_;
+    }
+
+    particle_interface_type const& operator*() const {
+      return static_cast<particle_interface_type const&>(*this);
+    }
+    //! @}
+
+  protected:
+    /**
+     * @name Stack data access
+     * Only the const versions for read-only access
+     * @{
+     */
+    unsigned int getIndex() const { return index_; }
+    stack_type const& getStack() const { return *data_; }
+    stack_data_type const& getStackData() const { return data_->getStackData(); }
+    /// Get data index as mapped in Stack class
+    unsigned int getIndexFromIterator() const {
+      return data_->getIndexFromIterator(index_);
+    }
+    //! @}
+
+    // friends are needed for access to protected methods
+    friend class Stack<stack_data_type, TParticleInterface,
+                       MSecondaryProducer>; // for access to GetIndex for Stack
+    friend class Stack<stack_data_type&, TParticleInterface,
+                       MSecondaryProducer>; // for access to GetIndex
+
+    friend class ParticleBase<ConstStackIteratorInterface>; // for access to GetStackData
+
+    template <typename T1,                     // best fix to: stack_data_type,
+              template <typename> typename M1, // best fix to: TParticleInterface,
+              template <class T2, template <class> class T3> class M2>
+    friend class SecondaryView; // access for SecondaryView
+
+    friend class StackIteratorInterface<stack_data_type, TParticleInterface,
+                                        MSecondaryProducer, stack_type>;
+
+    template <typename T, template <typename> typename TParticleInterface_>
+    friend class corsika::history::HistorySecondaryProducer;
+
+  protected:
+    unsigned int index_ = 0;
+
+  private:
+    stack_type const* data_ = 0; // info: Particles and StackIterators become invalid when
+                                 // parent Stack is copied or deleted!
+
+  }; // end class ConstStackIterator
+
+} // namespace corsika

third_party/phys/units/io.hpp → corsika/framework/units/io.hpp

@@ -16,9 +16,9 @@
 #ifndef PHYS_UNITS_IO_HPP_INCLUDED
 #define PHYS_UNITS_IO_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
-#include "phys/units/quantity_io_engineering.hpp"
-#include "phys/units/quantity_io_symbols.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io_engineering.hpp"
+#include "corsika/framework/units/quantity_io_symbols.hpp"
 
 #endif // PHYS_UNITS_IO_HPP_INCLUDED
 

third_party/phys/units/io_output.hpp → corsika/framework/units/io_output.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_IO_OUTPUT_HPP_INCLUDED
 #define PHYS_UNITS_IO_OUTPUT_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 #endif // PHYS_UNITS_IO_OUTPUT_HPP_INCLUDED
 

third_party/phys/units/io_output_eng.hpp → corsika/framework/units/io_output_eng.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_IO_ENG_HPP_INCLUDED
 #define PHYS_UNITS_IO_ENG_HPP_INCLUDED
 
-#include "phys/units/quantity_io_engineering.hpp"
+#include "corsika/framework/units/quantity_io_engineering.hpp"
 
 #endif // PHYS_UNITS_IO_ENG_HPP_INCLUDED
 

third_party/phys/units/io_symbols.hpp → corsika/framework/units/io_symbols.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_IO_SYMBOLS_HPP_INCLUDED
 #define PHYS_UNITS_IO_SYMBOLS_HPP_INCLUDED
 
-#include "phys/units/quantity_io_symbols.hpp"
+#include "corsika/framework/units/quantity_io_symbols.hpp"
 
 #endif // PHYS_UNITS_IO_SYMBOLS_HPP_INCLUDED
 

third_party/phys/units/other_units.hpp → corsika/framework/units/other_units.hpp

@@ -29,7 +29,7 @@
 #ifndef PHYS_UNITS_OTHER_UNITS_HPP_INCLUDED
 #define PHYS_UNITS_OTHER_UNITS_HPP_INCLUDED
 
-#include "phys/units/quantity.hpp"
+#include "corsika/framework/units/quantity.hpp"
 
 namespace phys { namespace units {
 

third_party/phys/units/physical_constants.hpp → corsika/framework/units/physical_constants.hpp

@@ -22,37 +22,36 @@
 #ifndef PHYS_UNITS_PHYSICAL_CONSTANTS_HPP_INCLUDED
 #define PHYS_UNITS_PHYSICAL_CONSTANTS_HPP_INCLUDED
 
-#include "phys/units/quantity.hpp"
+#include "corsika/framework/units/quantity.hpp"
 
-namespace phys { namespace units {
+namespace phys {
+  namespace units {
 
-// acceleration of free-fall, standard
-constexpr quantity< acceleration_d >
-                                g_sub_n { Rep( 9.80665L ) * meter / square( second ) };
+    // acceleration of free-fall, standard
+    constexpr quantity<acceleration_d> g_sub_n{Rep(9.80665L) * meter / square(second)};
 
-// Avogadro constant
-constexpr quantity< dimensions< 0, 0, 0, 0, 0, -1 > >
-                                N_sub_A { Rep( 6.02214199e+23L ) / mole };
-// electronvolt
-constexpr quantity< energy_d >  eV { Rep( 1.60217733e-19L ) * joule };
+    // Avogadro constant
+    constexpr quantity<dimensions<0, 0, 0, 0, 0, -1> > N_sub_A{Rep(6.02214199e+23L) /
+                                                               mole};
+    // electronvolt
+    // constexpr quantity< energy_d >  eV { Rep( 1.60217733e-19L ) * joule };
 
-// elementary charge
-constexpr quantity< electric_charge_d >
-                                e { Rep( 1.602176462e-19L ) * coulomb };
+    // elementary charge
+    constexpr quantity<electric_charge_d> e{Rep(1.602176462e-19L) * coulomb};
 
-// Planck constant
-constexpr quantity< dimensions< 2, 1, -1 > >
-                                h { Rep( 6.62606876e-34L ) * joule * second };
+    // Planck constant
+    constexpr quantity<dimensions<2, 1, -1> > h{Rep(6.62606876e-34L) * joule * second};
 
-// speed of light in a vacuum
-constexpr quantity< speed_d >   c { Rep( 299792458L ) * meter / second };
+    // speed of light in a vacuum
+    constexpr quantity<speed_d> c{Rep(299792458L) * meter / second};
 
-// unified atomic mass unit
-constexpr quantity< mass_d >    u { Rep( 1.6605402e-27L ) * kilogram };
+    // unified atomic mass unit
+    constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
 
-// etc.
+    // etc.
 
-}} // namespace phys { namespace units {
+  } // namespace units
+} // namespace phys
 
 #endif // PHYS_UNITS_PHYSICAL_CONSTANTS_HPP_INCLUDED
 

corsika/framework/units/quantity.hpp

+/**
+ * \file quantity.hpp
+ *
+ * \brief   Zero-overhead dimensional analysis and unit/quantity manipulation and
+ * conversion. \author  Michael S. Kenniston, Martin Moene \date    7 September 2013
+ * \since   0.4
+ *
+ * Copyright 2013 Universiteit Leiden. All rights reserved.
+ *
+ * Copyright (c) 2001 by Michael S. Kenniston.  For the most
+ * recent version check www.xnet.com/~msk/quantity.  Permission is granted
+ * to use this code without restriction so long as this copyright
+ * notice appears in all source files.
+ *
+ * This code is provided as-is, with no warrantee of correctness.
+ *
+ * Distributed under the Boost Software License, Version 1.0. (See accompanying
+ * file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+/*
+ * Unless otherwise specified, the definitions of all units in this
+ * file are from NIST Special Publication 811, found online at
+ * http://physics.nist.gov/Document/sp811.pdf
+ * Other sources: OED = Oxford English Dictionary
+ */
+
+#ifndef PHYS_UNITS_QUANTITY_HPP_INCLUDED
+#define PHYS_UNITS_QUANTITY_HPP_INCLUDED
+
+#include <cmath>
+#include <cstdlib>
+#include <utility> // std::declval
+#include <type_traits>
+
+/// namespace phys.
+
+namespace phys {
+
+  /// namespace units.
+
+  namespace units {
+
+#ifdef PHYS_UNITS_REP_TYPE
+    using Rep = PHYS_UNITS_REP_TYPE;
+#else
+    using Rep = double;
+#endif
+
+    /*
+     * declare now, define later.
+     */
+    template <typename Dims, typename T = Rep>
+    class quantity;
+
+    /**
+     * We could drag dimensions around individually, but it's much more convenient to
+     * package them.
+     */
+    template <int D1, int D2, int D3, int D4 = 0, int D5 = 0, int D6 = 0, int D7 = 0,
+              int D8 = 0>
+    struct dimensions {
+      enum {
+        dim1 = D1,
+        dim2 = D2,
+        dim3 = D3,
+        dim4 = D4,
+        dim5 = D5,
+        dim6 = D6,
+        dim7 = D7,
+        dim8 = D8,
+
+        is_all_zero = D1 == 0 && D2 == 0 && D3 == 0 && D4 == 0 && D5 == 0 && D6 == 0 &&
+                      D7 == 0 && D8 == 0,
+
+        is_base = 1 == (D1 != 0) + (D2 != 0) + (D3 != 0) + (D4 != 0) + (D5 != 0) +
+                           (D6 != 0) + (D7 != 0) + (D8 != 0) &&
+                  1 == D1 + D2 + D3 + D4 + D5 + D6 + D7 + D8,
+      };
+
+      template <int R1, int R2, int R3, int R4, int R5, int R6, int R7, int R8>
+      constexpr bool operator==(dimensions<R1, R2, R3, R4, R5, R6, R7, R8> const&) const {
+        return D1 == R1 && D2 == R2 && D3 == R3 && D4 == R4 && D5 == R5 && D6 == R6 &&
+               D7 == R7 && D8 == R8;
+      }
+
+      template <int R1, int R2, int R3, int R4, int R5, int R6, int R7, int R8>
+      constexpr bool operator!=(
+          dimensions<R1, R2, R3, R4, R5, R6, R7, R8> const& rhs) const {
+        return !(*this == rhs);
+      }
+    };
+
+    /// demensionless 'dimension'.
+
+    typedef dimensions<0, 0, 0> dimensionless_d;
+
+    /// namespace detail.
+
+    namespace detail {
+
+      /**
+       * \brief The "collapse" template is used to avoid quantity< dimensions< 0, 0, 0 >
+       * >, i.e. to make dimensionless results come out as type "Rep".
+       */
+      template <typename D, typename T>
+      struct collapse {
+        typedef quantity<D, T> type;
+      };
+
+      template <typename T>
+      struct collapse<dimensionless_d, T> {
+        typedef T type;
+      };
+
+      template <typename D, typename T>
+      using Collapse = typename collapse<D, T>::type;
+
+      // promote types of expression to result type.
+
+      template <typename X, typename Y>
+      using PromoteAdd = decltype(std::declval<X>() + std::declval<Y>());
+
+      template <typename X, typename Y>
+      using PromoteMul = decltype(std::declval<X>() * std::declval<Y>());
+
+      /*
+       * The following batch of structs are type generators to calculate
+       * the correct type of the result of various operations.
+       */
+
+      /**
+       * product type generator.
+       */
+      template <typename DX, typename DY, typename T>
+      struct product {
+        enum {
+          d1 = DX::dim1 + DY::dim1,
+          d2 = DX::dim2 + DY::dim2,
+          d3 = DX::dim3 + DY::dim3,
+          d4 = DX::dim4 + DY::dim4,
+          d5 = DX::dim5 + DY::dim5,
+          d6 = DX::dim6 + DY::dim6,
+          d7 = DX::dim7 + DY::dim7,
+          d8 = DX::dim8 + DY::dim8
+        };
+
+        typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+        typedef Collapse<dim, T> type;
+      };
+
+      /**
+       * convenience type for product dimension
+       */
+
+      template <typename DX, typename DY>
+      using product_d = typename product<DX, DY, Rep>::dim;
+
+      template <typename DX, typename DY, typename X, typename Y>
+      using Product = typename product<DX, DY, PromoteMul<X, Y>>::type;
+
+      /**
+       * quotient type generator.
+       */
+      template <typename DX, typename DY, typename T>
+      struct quotient {
+        enum {
+          d1 = DX::dim1 - DY::dim1,
+          d2 = DX::dim2 - DY::dim2,
+          d3 = DX::dim3 - DY::dim3,
+          d4 = DX::dim4 - DY::dim4,
+          d5 = DX::dim5 - DY::dim5,
+          d6 = DX::dim6 - DY::dim6,
+          d7 = DX::dim7 - DY::dim7,
+          d8 = DX::dim8 - DY::dim8
+        };
+
+        typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+        typedef Collapse<dim, T> type;
+      };
+
+      /**
+       * convenience type for quotient dimension
+       */
+
+      template <typename DX, typename DY>
+      using quotient_d = typename quotient<DX, DY, Rep>::dim;
+
+      template <typename DX, typename DY, typename X, typename Y>
+      using Quotient = typename quotient<DX, DY, PromoteMul<X, Y>>::type;
+
+      /**
+       * reciprocal type generator.
+       */
+      template <typename D, typename T>
+      struct reciprocal {
+        enum {
+          d1 = -D::dim1,
+          d2 = -D::dim2,
+          d3 = -D::dim3,
+          d4 = -D::dim4,
+          d5 = -D::dim5,
+          d6 = -D::dim6,
+          d7 = -D::dim7,
+          d8 = -D::dim8
+        };
+
+        typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+        typedef Collapse<dim, T> type;
+      };
+
+      /**
+       * convenience type for reciprocal dimension
+       */
+
+      template <typename DX, typename DY>
+      using reciprocal_d = typename reciprocal<DX, Rep>::dim;
+
+      template <typename D, typename X, typename Y>
+      using Reciprocal = typename reciprocal<D, PromoteMul<X, Y>>::type;
+
+      /**
+       * power type generator.
+       */
+      template <typename D, int N, typename T>
+      struct power {
+        enum {
+          d1 = N * D::dim1,
+          d2 = N * D::dim2,
+          d3 = N * D::dim3,
+          d4 = N * D::dim4,
+          d5 = N * D::dim5,
+          d6 = N * D::dim6,
+          d7 = N * D::dim7,
+          d8 = N * D::dim8
+        };
+
+        typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+        typedef Collapse<dim, T> type;
+      };
+
+      /**
+       * convenience type for power dimension
+       */
+
+      template <typename DX, int N>
+      using power_d = typename power<DX, N, Rep>::dim;
+
+      template <typename D, int N, typename T>
+      using Power = typename power<D, N, T>::type;
+
+      /**
+       * root type generator.
+       */
+      template <typename D, int N, typename T>
+      struct root {
+        enum {
+          all_even_multiples = D::dim1 % N == 0 && D::dim2 % N == 0 && D::dim3 % N == 0 &&
+                               D::dim4 % N == 0 && D::dim5 % N == 0 && D::dim6 % N == 0 &&
+                               D::dim7 % N == 0 && D::dim8 % N == 0
+        };
+
+        enum {
+          d1 = D::dim1 / N,
+          d2 = D::dim2 / N,
+          d3 = D::dim3 / N,
+          d4 = D::dim4 / N,
+          d5 = D::dim5 / N,
+          d6 = D::dim6 / N,
+          d7 = D::dim7 / N,
+          d8 = D::dim8 / N
+        };
+
+        typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+        typedef Collapse<dim, T> type;
+      };
+
+      /**
+       * convenience type for root dimension
+       */
+      template <typename D, int N>
+      using root_d = typename root<D, N, Rep>::dim;
+
+      template <typename D, int N, typename T>
+      using Root = typename detail::root<D, N, T>::type;
+
+      /**
+       * tag to construct a quantity from a magnitude.
+       */
+      constexpr struct magnitude_tag_t { } magnitude_tag{}; } // namespace detail
+
+    /**
+     * \brief class "quantity" is the heart of the library. It associates
+     * dimensions  with a single "Rep" data member and protects it from
+     * dimensionally inconsistent use.
+     */
+    template <typename Dims, typename T /*= Rep */>
+    class quantity {
+    public:
+      typedef Dims dimension_type;
+
+      typedef T value_type;
+
+      typedef quantity<Dims, T> this_type;
+
+      constexpr quantity()
+          : m_value{} {}
+
+      /**
+       * public converting initializing constructor;
+       * requires magnitude_tag to prevent constructing a quantity from a raw magnitude.
+       */
+      template <typename X>
+      constexpr explicit quantity(detail::magnitude_tag_t, X x)
+          : m_value(x) {}
+
+      /**
+       * converting copy-assignment constructor.
+       */
+      template <typename X>
+      constexpr quantity(quantity<Dims, X> const& x)
+          : m_value(x.magnitude()) {}
+
+      //    /**
+      //     * convert to compatible unit, for example: (3._dm).to(meter) gives 0.3;
+      //     */
+      //    constexpr value_type to( quantity const & x ) const {  return *this / x; }
+
+      /**
+       * convert to given unit, for example: (3._dm).to(meter) gives 0.3;
+       */
+      template <typename DX, typename X>
+      constexpr auto to(quantity<DX, X> const& x) const
+          -> detail::Quotient<Dims, DX, T, X> {
+        return *this / x;
+      }
+
+      /**
+       * the quantity's magnitude.
+       */
+      constexpr value_type magnitude() const { return m_value; }
+
+      /**
+       * the quantity's dimensions.
+       */
+      constexpr dimension_type dimension() const { return dimension_type{}; }
+
+      /**
+       * We need a "zero" of each type -- for comparisons, to initialize running
+       * totals, etc.  Note:  0 m != 0 kg, since they are of different dimensionality.
+       * zero is really just defined for convenience, since
+       * quantity< length_d >::zero == 0 * meter, etc.
+       */
+      static constexpr quantity zero() { return quantity{value_type(0.0)}; }
+      //    static constexpr quantity zero = quantity{ value_type( 0.0 ) };
+
+      /**
+       * We also define "infinity" for each type.
+       *
+       * RU, Do 7. Jan 02:10:18 CET 2021
+       */
+      static constexpr quantity infinity() {
+        return quantity{value_type(std::numeric_limits<value_type>::infinity())};
+      }
+
+      template <typename _dim = dimension_type,
+                std::enable_if_t<std::is_same_v<_dim, dimensionless_d>, bool> = true>
+      operator double() const {
+        return m_value;
+      }
+
+    private:
+      /**
+       * private initializing constructor.
+       */
+      constexpr explicit quantity(value_type x)
+          : m_value{x} {}
+
+    private:
+      value_type m_value;
+
+      enum { has_dimension = !Dims::is_all_zero };
+
+      // static_assert( has_dimension, "quantity dimensions must not all be zero" ); //
+      // MR: removed
+
+    private:
+      // friends:
+
+      // arithmetic
+
+      template <typename D, typename X, typename Y>
+      friend constexpr quantity<D, X>& operator+=(quantity<D, X>& x,
+                                                  quantity<D, Y> const& y);
+
+      template <typename D, typename X>
+      friend constexpr quantity<D, X> operator+(quantity<D, X> const& x);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr quantity<D, detail::PromoteAdd<X, Y>> operator+(
+          quantity<D, X> const& x, quantity<D, Y> const& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr quantity<D, X>& operator-=(quantity<D, X>& x,
+                                                  quantity<D, Y> const& y);
+
+      template <typename D, typename X>
+      friend constexpr quantity<D, X> operator-(quantity<D, X> const& x);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr quantity<D, detail::PromoteAdd<X, Y>> operator-(
+          quantity<D, X> const& x, quantity<D, Y> const& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr quantity<D, X>& operator*=(quantity<D, X>& x, const Y& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr quantity<D, detail::PromoteMul<X, Y>> operator*(
+          quantity<D, X> const& x, const Y& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr quantity<D, detail::PromoteMul<X, Y>> operator*(
+          const X& x, quantity<D, Y> const& y);
+
+      template <typename DX, typename DY, typename X, typename Y>
+      friend constexpr detail::Product<DX, DY, X, Y> operator*(
+          quantity<DX, X> const& lhs, quantity<DY, Y> const& rhs);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr quantity<D, X>& operator/=(quantity<D, X>& x, const Y& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr quantity<D, detail::PromoteMul<X, Y>> operator/(
+          quantity<D, X> const& x, const Y& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr detail::Reciprocal<D, X, Y> operator/(const X& x,
+                                                             quantity<D, Y> const& y);
+
+      template <typename DX, typename DY, typename X, typename Y>
+      friend constexpr detail::Quotient<DX, DY, X, Y> operator/(quantity<DX, X> const& x,
+                                                                quantity<DY, Y> const& y);
+
+      // absolute value.
+
+      template <typename D, typename X>
+      friend constexpr quantity<D, X> abs(quantity<D, X> const& x);
+
+      // powers and roots
+
+      template <int N, typename D, typename X>
+      friend constexpr detail::Power<D, N, X> nth_power(quantity<D, X> const& x);
+
+      template <typename D, typename X>
+      friend constexpr detail::Power<D, 2, X> square(quantity<D, X> const& x);
+
+      template <typename D, typename X>
+      friend constexpr detail::Power<D, 3, X> cube(quantity<D, X> const& x);
+
+      template <int N, typename D, typename X>
+      friend detail::Root<D, N, X> constexpr nth_root(quantity<D, X> const& x);
+
+      template <typename D, typename X>
+      friend detail::Root<D, 2, X> constexpr sqrt(quantity<D, X> const& x);
+
+      template <typename D, typename X>
+      friend detail::Root<D, 3, X> constexpr cbrt(quantity<D, X> const& x);
+
+      // comparison
+
+      template <typename D, typename X, typename Y>
+      friend constexpr bool operator==(quantity<D, X> const& x, quantity<D, Y> const& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr bool operator!=(quantity<D, X> const& x, quantity<D, Y> const& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr bool operator<(quantity<D, X> const& x, quantity<D, Y> const& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr bool operator<=(quantity<D, X> const& x, quantity<D, Y> const& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr bool operator>(quantity<D, X> const& x, quantity<D, Y> const& y);
+
+      template <typename D, typename X, typename Y>
+      friend constexpr bool operator>=(quantity<D, X> const& x, quantity<D, Y> const& y);
+    }; // namespace units
+
+    // helper to check whether some type is a quantity
+    template <typename TNonQuantityType>
+    struct is_quantity : std::false_type {};
+
+    template <typename Dims, typename T>
+    struct is_quantity<quantity<Dims, T>> : std::true_type {};
+
+    template <typename T>
+    constexpr bool is_quantity_v = is_quantity<T>::value;
+
+    // Give names to the seven fundamental dimensions of physical reality.
+
+    typedef dimensions<1, 0, 0, 0, 0, 0, 0, 0> length_d;
+    typedef dimensions<0, 1, 0, 0, 0, 0, 0, 0> mass_d;
+    typedef dimensions<0, 0, 1, 0, 0, 0, 0, 0> time_interval_d;
+    typedef dimensions<0, 0, 0, 1, 0, 0, 0, 0> electric_current_d;
+    typedef dimensions<0, 0, 0, 0, 1, 0, 0, 0> thermodynamic_temperature_d;
+    typedef dimensions<0, 0, 0, 0, 0, 1, 0, 0> amount_of_substance_d;
+    typedef dimensions<0, 0, 0, 0, 0, 0, 1, 0> luminous_intensity_d;
+    typedef dimensions<0, 0, 0, 0, 0, 0, 0, 1> hepenergy_d; // this is not an SI unit !
+
+    // Addition operators
+
+    /// quan += quan
+
+    template <typename D, typename X, typename Y>
+    constexpr quantity<D, X>& operator+=(quantity<D, X>& x, quantity<D, Y> const& y) {
+      return x.m_value += y.m_value, x;
+    }
+
+    /// + quan
+
+    template <typename D, typename X>
+    constexpr quantity<D, X> operator+(quantity<D, X> const& x) {
+      return quantity<D, X>(+x.m_value);
+    }
+
+    /// quan + quan
+
+    template <typename D, typename X, typename Y>
+    constexpr quantity<D, detail::PromoteAdd<X, Y>> operator+(quantity<D, X> const& x,
+                                                              quantity<D, Y> const& y) {
+      return quantity<D, detail::PromoteAdd<X, Y>>(x.m_value + y.m_value);
+    }
+
+    // Subtraction operators
+
+    /// quan -= quan
+
+    template <typename D, typename X, typename Y>
+    constexpr quantity<D, X>& operator-=(quantity<D, X>& x, quantity<D, Y> const& y) {
+      return x.m_value -= y.m_value, x;
+    }
+
+    /// - quan
+
+    template <typename D, typename X>
+    constexpr quantity<D, X> operator-(quantity<D, X> const& x) {
+      return quantity<D, X>(-x.m_value);
+    }
+
+    /// quan - quan
+
+    template <typename D, typename X, typename Y>
+    constexpr quantity<D, detail::PromoteAdd<X, Y>> operator-(quantity<D, X> const& x,
+                                                              quantity<D, Y> const& y) {
+      return quantity<D, detail::PromoteAdd<X, Y>>(x.m_value - y.m_value);
+    }
+
+    // Multiplication operators
+
+    /// quan *= num
+
+    template <typename D, typename X, typename Y>
+    constexpr quantity<D, X>& operator*=(quantity<D, X>& x, const Y& y) {
+      return x.m_value *= y, x;
+    }
+
+    /// quan * num
+
+    template <typename D, typename X, typename Y>
+    constexpr quantity<D, detail::PromoteMul<X, Y>> operator*(quantity<D, X> const& x,
+                                                              const Y& y) {
+      return quantity<D, detail::PromoteMul<X, Y>>(x.m_value * y);
+    }
+
+    /// num * quan
+
+    template <typename D, typename X, typename Y>
+    constexpr quantity<D, detail::PromoteMul<X, Y>> operator*(const X& x,
+                                                              quantity<D, Y> const& y) {
+      return quantity<D, detail::PromoteMul<X, Y>>(x * y.m_value);
+    }
+
+    /// quan * quan:
+
+    template <typename DX, typename DY, typename X, typename Y>
+    constexpr detail::Product<DX, DY, X, Y> operator*(quantity<DX, X> const& lhs,
+                                                      quantity<DY, Y> const& rhs) {
+      return detail::Product<DX, DY, X, Y>(lhs.m_value * rhs.m_value);
+    }
+
+    // Division operators
+
+    /// quan /= num
+
+    template <typename D, typename X, typename Y>
+    constexpr quantity<D, X>& operator/=(quantity<D, X>& x, const Y& y) {
+      return x.m_value /= y, x;
+    }
+
+    /// quan / num
+
+    template <typename D, typename X, typename Y>
+    constexpr quantity<D, detail::PromoteMul<X, Y>> operator/(quantity<D, X> const& x,
+                                                              const Y& y) {
+      return quantity<D, detail::PromoteMul<X, Y>>(x.m_value / y);
+    }
+
+    /// num / quan
+
+    template <typename D, typename X, typename Y>
+    constexpr detail::Reciprocal<D, X, Y> operator/(const X& x, quantity<D, Y> const& y) {
+      return detail::Reciprocal<D, X, Y>(x / y.m_value);
+    }
+
+    /// quan / quan:
+
+    template <typename DX, typename DY, typename X, typename Y>
+    constexpr detail::Quotient<DX, DY, X, Y> operator/(quantity<DX, X> const& x,
+                                                       quantity<DY, Y> const& y) {
+      return detail::Quotient<DX, DY, X, Y>(x.m_value / y.m_value);
+    }
+
+    /// absolute value.
+
+    template <typename D, typename X>
+    quantity<D, X> constexpr abs(quantity<D, X> const& x) {
+      return quantity<D, X>(std::abs(x.m_value));
+    }
+
+    // General powers
+
+    /// N-th power.
+
+    template <int N, typename D, typename X>
+    detail::Power<D, N, X> constexpr nth_power(quantity<D, X> const& x) {
+      return detail::Power<D, N, X>(std::pow(x.m_value, X(N)));
+    }
+
+    // Low powers defined separately for efficiency.
+
+    /// square.
+
+    template <typename D, typename X>
+    constexpr detail::Power<D, 2, X> square(quantity<D, X> const& x) {
+      return x * x;
+    }
+
+    /// cube.
+
+    template <typename D, typename X>
+    constexpr detail::Power<D, 3, X> cube(quantity<D, X> const& x) {
+      return x * x * x;
+    }
+
+    // General root
+
+    /// n-th root.
+
+    template <int N, typename D, typename X>
+    detail::Root<D, N, X> constexpr nth_root(quantity<D, X> const& x) {
+      static_assert(detail::root<D, N, X>::all_even_multiples,
+                    "root result dimensions must be integral");
+
+      return detail::Root<D, N, X>(std::pow(x.m_value, X(1.0) / N));
+    }
+
+    // Low roots defined separately for convenience.
+
+    /// square root.
+
+    template <typename D, typename X>
+    detail::Root<D, 2, X> constexpr sqrt(quantity<D, X> const& x) {
+      static_assert(detail::root<D, 2, X>::all_even_multiples,
+                    "root result dimensions must be integral");
+
+      return detail::Root<D, 2, X>(std::sqrt(x.m_value));
+    }
+
+    /// cubic root.
+
+    template <typename D, typename X>
+    detail::Root<D, 3, X> constexpr cbrt(quantity<D, X> const& x) {
+      static_assert(detail::root<D, 3, X>::all_even_multiples,
+                    "root result dimensions must be integral");
+
+      return detail::Root<D, 3, X>(std::cbrt(x.m_value));
+    }
+
+    // Comparison operators
+
+    /// equality.
+
+    template <typename D, typename X, typename Y>
+    constexpr bool operator==(quantity<D, X> const& x, quantity<D, Y> const& y) {
+      return x.m_value == y.m_value;
+    }
+
+    /// inequality.
+
+    template <typename D, typename X, typename Y>
+    constexpr bool operator!=(quantity<D, X> const& x, quantity<D, Y> const& y) {
+      return x.m_value != y.m_value;
+    }
+
+    /// less-than.
+
+    template <typename D, typename X, typename Y>
+    constexpr bool operator<(quantity<D, X> const& x, quantity<D, Y> const& y) {
+      return x.m_value < y.m_value;
+    }
+
+    /// less-equal.
+
+    template <typename D, typename X, typename Y>
+    constexpr bool operator<=(quantity<D, X> const& x, quantity<D, Y> const& y) {
+      return x.m_value <= y.m_value;
+    }
+
+    /// greater-than.
+
+    template <typename D, typename X, typename Y>
+    constexpr bool operator>(quantity<D, X> const& x, quantity<D, Y> const& y) {
+      return x.m_value > y.m_value;
+    }
+
+    /// greater-equal.
+
+    template <typename D, typename X, typename Y>
+    constexpr bool operator>=(quantity<D, X> const& x, quantity<D, Y> const& y) {
+      return x.m_value >= y.m_value;
+    }
+
+    /// quantity's dimension.
+
+    template <typename DX, typename X>
+    inline constexpr DX dimension(quantity<DX, X> const& q) {
+      return q.dimension();
+    }
+
+    /// quantity's magnitude.
+
+    template <typename DX, typename X>
+    inline constexpr X magnitude(quantity<DX, X> const& q) {
+      return q.magnitude();
+    }
+
+    // The seven SI base units.  These tie our numbers to the real world.
+
+    constexpr quantity<length_d> meter{detail::magnitude_tag, 1.0};
+    constexpr quantity<mass_d> kilogram{detail::magnitude_tag, 1.0};
+    constexpr quantity<time_interval_d> second{detail::magnitude_tag, 1.0};
+    constexpr quantity<electric_current_d> ampere{detail::magnitude_tag, 1.0};
+    constexpr quantity<thermodynamic_temperature_d> kelvin{detail::magnitude_tag, 1.0};
+    constexpr quantity<amount_of_substance_d> mole{detail::magnitude_tag, 1.0};
+    constexpr quantity<luminous_intensity_d> candela{detail::magnitude_tag, 1.0};
+    constexpr quantity<hepenergy_d> electronvolt{detail::magnitude_tag, 1.0};
+
+    // The standard SI prefixes.
+
+    constexpr long double quetta = 1e+30L;
+    constexpr long double ronna = 1e+27L;
+    constexpr long double yotta = 1e+24L;
+    constexpr long double zetta = 1e+21L;
+    constexpr long double exa = 1e+18L;
+    constexpr long double peta = 1e+15L;
+    constexpr long double tera = 1e+12L;
+    constexpr long double giga = 1e+9L;
+    constexpr long double mega = 1e+6L;
+    constexpr long double kilo = 1e+3L;
+    constexpr long double hecto = 1e+2L;
+    constexpr long double deka = 1e+1L;
+    constexpr long double deci = 1e-1L;
+    constexpr long double centi = 1e-2L;
+    constexpr long double milli = 1e-3L;
+    constexpr long double micro = 1e-6L;
+    constexpr long double nano = 1e-9L;
+    constexpr long double pico = 1e-12L;
+    constexpr long double femto = 1e-15L;
+    constexpr long double atto = 1e-18L;
+    constexpr long double zepto = 1e-21L;
+    constexpr long double yocto = 1e-24L;
+    constexpr long double ronto = 1e-27L;
+    constexpr long double quecto = 1e-30L;
+
+    // Binary prefixes, pending adoption.
+
+    constexpr long double kibi = 1024;
+    constexpr long double mebi = 1024 * kibi;
+    constexpr long double gibi = 1024 * mebi;
+    constexpr long double tebi = 1024 * gibi;
+    constexpr long double pebi = 1024 * tebi;
+    constexpr long double exbi = 1024 * pebi;
+    constexpr long double zebi = 1024 * exbi;
+    constexpr long double yobi = 1024 * zebi;
+
+    // The rest of the standard dimensional types, as specified in SP811.
+
+    using absorbed_dose_d = dimensions<2, 0, -2>;
+    using absorbed_dose_rate_d = dimensions<2, 0, -3>;
+    using acceleration_d = dimensions<1, 0, -2>;
+    using activity_of_a_nuclide_d = dimensions<0, 0, -1>;
+    using angular_velocity_d = dimensions<0, 0, -1>;
+    using angular_acceleration_d = dimensions<0, 0, -2>;
+    using area_d = dimensions<2, 0, 0>;
+    using capacitance_d = dimensions<-2, -1, 4, 2>;
+    using concentration_d = dimensions<-3, 0, 0, 0, 0, 1>;
+    using current_density_d = dimensions<-2, 0, 0, 1>;
+    using dose_equivalent_d = dimensions<2, 0, -2>;
+    using dynamic_viscosity_d = dimensions<-1, 1, -1>;
+    using electric_charge_d = dimensions<0, 0, 1, 1>;
+    using electric_charge_density_d = dimensions<-3, 0, 1, 1>;
+    using electric_conductance_d = dimensions<-2, -1, 3, 2>;
+    using electric_field_strenth_d = dimensions<1, 1, -3, -1>;
+    using electric_flux_density_d = dimensions<-2, 0, 1, 1>;
+    using electric_potential_d = dimensions<2, 1, -3, -1>;
+    using electric_resistance_d = dimensions<2, 1, -3, -2>;
+    using energy_d = dimensions<2, 1, -2>;
+    using energy_density_d = dimensions<-1, 1, -2>;
+    using exposure_d = dimensions<0, -1, 1, 1>;
+    using force_d = dimensions<1, 1, -2>;
+    using frequency_d = dimensions<0, 0, -1>;
+    using heat_capacity_d = dimensions<2, 1, -2, 0, -1>;
+    using heat_density_d = dimensions<0, 1, -2>;
+    using heat_density_flow_rate_d = dimensions<0, 1, -3>;
+    using heat_flow_rate_d = dimensions<2, 1, -3>;
+    using heat_flux_density_d = dimensions<0, 1, -3>;
+    using heat_transfer_coefficient_d = dimensions<0, 1, -3, 0, -1>;
+    using illuminance_d = dimensions<-2, 0, 0, 0, 0, 0, 1>;
+    using inductance_d = dimensions<2, 1, -2, -2>;
+    using irradiance_d = dimensions<0, 1, -3>;
+    using kinematic_viscosity_d = dimensions<2, 0, -1>;
+    using luminance_d = dimensions<-2, 0, 0, 0, 0, 0, 1>;
+    using luminous_flux_d = dimensions<0, 0, 0, 0, 0, 0, 1>;
+    using magnetic_field_strength_d = dimensions<-1, 0, 0, 1>;
+    using magnetic_flux_d = dimensions<2, 1, -2, -1>;
+    using magnetic_flux_density_d = dimensions<0, 1, -2, -1>;
+    using magnetic_permeability_d = dimensions<1, 1, -2, -2>;
+    using mass_density_d = dimensions<-3, 1, 0>;
+    using mass_flow_rate_d = dimensions<0, 1, -1>;
+    using molar_energy_d = dimensions<2, 1, -2, 0, 0, -1>;
+    using molar_entropy_d = dimensions<2, 1, -2, -1, 0, -1>;
+    using moment_of_force_d = dimensions<2, 1, -2>;
+    using permittivity_d = dimensions<-3, -1, 4, 2>;
+    using power_d = dimensions<2, 1, -3>;
+    using pressure_d = dimensions<-1, 1, -2>;
+    using radiance_d = dimensions<0, 1, -3>;
+    using radiant_intensity_d = dimensions<2, 1, -3>;
+    using speed_d = dimensions<1, 0, -1>;
+    using specific_energy_d = dimensions<2, 0, -2>;
+    using specific_heat_capacity_d = dimensions<2, 0, -2, 0, -1>;
+    using specific_volume_d = dimensions<3, -1, 0>;
+    using substance_permeability_d = dimensions<-1, 0, 1>;
+    using surface_tension_d = dimensions<0, 1, -2>;
+    using thermal_conductivity_d = dimensions<1, 1, -3, 0, -1>;
+    using thermal_diffusivity_d = dimensions<2, 0, -1>;
+    using thermal_insulance_d = dimensions<0, -1, 3, 0, 1>;
+    using thermal_resistance_d = dimensions<-2, -1, 3, 0, 1>;
+    using thermal_resistivity_d = dimensions<-1, -1, 3, 0, 1>;
+    using torque_d = dimensions<2, 1, -2>;
+    using volume_d = dimensions<3, 0, 0>;
+    using volume_flow_rate_d = dimensions<3, 0, -1>;
+    using wave_number_d = dimensions<-1, 0, 0>;
+
+    // Handy values.
+
+    constexpr Rep pi{Rep(3.141592653589793238462L)};
+    constexpr Rep percent{Rep(1) / 100};
+
+    //// Not approved for use alone, but needed for use with prefixes.
+    constexpr quantity<mass_d> gram{kilogram / 1000};
+
+    // The derived SI units, as specified in SP811.
+
+    constexpr Rep radian{Rep(1)};
+    constexpr Rep steradian{Rep(1)};
+    constexpr quantity<force_d> newton{meter * kilogram / square(second)};
+    constexpr quantity<pressure_d> pascal{newton / square(meter)};
+    constexpr quantity<energy_d> joule{newton * meter};
+    constexpr quantity<power_d> watt{joule / second};
+    constexpr quantity<electric_charge_d> coulomb{second * ampere};
+    constexpr quantity<electric_potential_d> volt{watt / ampere};
+    constexpr quantity<capacitance_d> farad{coulomb / volt};
+    constexpr quantity<electric_resistance_d> ohm{volt / ampere};
+    constexpr quantity<electric_conductance_d> siemens{ampere / volt};
+    constexpr quantity<magnetic_flux_d> weber{volt * second};
+    constexpr quantity<magnetic_flux_density_d> tesla{weber / square(meter)};
+    constexpr quantity<inductance_d> henry{weber / ampere};
+    constexpr quantity<thermodynamic_temperature_d> degree_celsius{kelvin};
+    constexpr quantity<luminous_flux_d> lumen{candela * steradian};
+    constexpr quantity<illuminance_d> lux{lumen / meter / meter};
+    constexpr quantity<activity_of_a_nuclide_d> becquerel{1 / second};
+    constexpr quantity<absorbed_dose_d> gray{joule / kilogram};
+    constexpr quantity<dose_equivalent_d> sievert{joule / kilogram};
+    constexpr quantity<frequency_d> hertz{1 / second};
+
+    // The rest of the units approved for use with SI, as specified in SP811.
+    // (However, use of these units is generally discouraged.)
+
+    constexpr quantity<length_d> angstrom{Rep(1e-10L) * meter};
+    constexpr quantity<area_d> are{Rep(1e+2L) * square(meter)};
+    constexpr quantity<pressure_d> bar{Rep(1e+5L) * pascal};
+    constexpr quantity<area_d> barn{Rep(1e-28L) * square(meter)};
+    constexpr quantity<activity_of_a_nuclide_d> curie{Rep(3.7e+10L) * becquerel};
+    constexpr quantity<time_interval_d> day{Rep(86400L) * second};
+    constexpr Rep degree_angle{pi / 180};
+    constexpr quantity<acceleration_d> gal{Rep(1e-2L) * meter / square(second)};
+    constexpr quantity<area_d> hectare{Rep(1e+4L) * square(meter)};
+    constexpr quantity<time_interval_d> hour{Rep(3600) * second};
+    constexpr quantity<speed_d> knot{Rep(1852) / 3600 * meter / second};
+    constexpr quantity<volume_d> liter{Rep(1e-3L) * cube(meter)};
+    constexpr quantity<time_interval_d> minute{Rep(60) * second};
+    constexpr Rep minute_angle{pi / 10800};
+    constexpr quantity<length_d> mile_nautical{Rep(1852) * meter};
+    constexpr quantity<absorbed_dose_d> rad{Rep(1e-2L) * gray};
+    constexpr quantity<dose_equivalent_d> rem{Rep(1e-2L) * sievert};
+    constexpr quantity<exposure_d> roentgen{Rep(2.58e-4L) * coulomb / kilogram};
+    constexpr Rep second_angle{pi / 648000L};
+    constexpr quantity<mass_d> ton_metric{Rep(1e+3L) * kilogram};
+
+    // Alternate (non-US) spellings:
+
+    constexpr quantity<length_d> metre{meter};
+    constexpr quantity<volume_d> litre{liter};
+    constexpr Rep deca{deka};
+    constexpr quantity<mass_d> tonne{ton_metric};
+
+    // cooked literals for base units;
+    // these could also have been created with a script.
+
+#define QUANTITY_DEFINE_SCALING_LITERAL(sfx, dim, factor)                   \
+  constexpr quantity<dim, double> operator"" _##sfx(unsigned long long x) { \
+    return quantity<dim, double>(detail::magnitude_tag, factor * x);        \
+  }                                                                         \
+  constexpr quantity<dim, double> operator"" _##sfx(long double x) {        \
+    return quantity<dim, double>(detail::magnitude_tag, factor * x);        \
+  }
+
+#define QUANTITY_DEFINE_SCALING_LITERALS(pfx, dim, fact)    \
+  QUANTITY_DEFINE_SCALING_LITERAL(Q##pfx, dim, fact* quetta) \
+  QUANTITY_DEFINE_SCALING_LITERAL(R##pfx, dim, fact* ronna) \
+  QUANTITY_DEFINE_SCALING_LITERAL(Y##pfx, dim, fact* yotta) \
+  QUANTITY_DEFINE_SCALING_LITERAL(Z##pfx, dim, fact* zetta) \
+  QUANTITY_DEFINE_SCALING_LITERAL(E##pfx, dim, fact* exa)   \
+  QUANTITY_DEFINE_SCALING_LITERAL(P##pfx, dim, fact* peta)  \
+  QUANTITY_DEFINE_SCALING_LITERAL(T##pfx, dim, fact* tera)  \
+  QUANTITY_DEFINE_SCALING_LITERAL(G##pfx, dim, fact* giga)  \
+  QUANTITY_DEFINE_SCALING_LITERAL(M##pfx, dim, fact* mega)  \
+  QUANTITY_DEFINE_SCALING_LITERAL(k##pfx, dim, fact* kilo)  \
+  QUANTITY_DEFINE_SCALING_LITERAL(h##pfx, dim, fact* hecto) \
+  QUANTITY_DEFINE_SCALING_LITERAL(da##pfx, dim, fact* deka) \
+  QUANTITY_DEFINE_SCALING_LITERAL(pfx, dim, fact * 1)       \
+  QUANTITY_DEFINE_SCALING_LITERAL(d##pfx, dim, fact* deci)  \
+  QUANTITY_DEFINE_SCALING_LITERAL(c##pfx, dim, fact* centi) \
+  QUANTITY_DEFINE_SCALING_LITERAL(m##pfx, dim, fact* milli) \
+  QUANTITY_DEFINE_SCALING_LITERAL(u##pfx, dim, fact* micro) \
+  QUANTITY_DEFINE_SCALING_LITERAL(n##pfx, dim, fact* nano)  \
+  QUANTITY_DEFINE_SCALING_LITERAL(p##pfx, dim, fact* pico)  \
+  QUANTITY_DEFINE_SCALING_LITERAL(f##pfx, dim, fact* femto) \
+  QUANTITY_DEFINE_SCALING_LITERAL(a##pfx, dim, fact* atto)  \
+  QUANTITY_DEFINE_SCALING_LITERAL(z##pfx, dim, fact* zepto) \
+  QUANTITY_DEFINE_SCALING_LITERAL(y##pfx, dim, fact* yocto) \
+  QUANTITY_DEFINE_SCALING_LITERAL(r##pfx, dim, fact* ronto) \
+  QUANTITY_DEFINE_SCALING_LITERAL(q##pfx, dim, fact* quecto)
+
+#define QUANTITY_DEFINE_LITERALS(pfx, dim) QUANTITY_DEFINE_SCALING_LITERALS(pfx, dim, 1)
+
+    /// literals
+
+    namespace literals {
+
+      QUANTITY_DEFINE_SCALING_LITERALS(g, mass_d, 1e-3)
+
+      QUANTITY_DEFINE_LITERALS(m, length_d)
+      QUANTITY_DEFINE_LITERALS(s, time_interval_d)
+      QUANTITY_DEFINE_LITERALS(A, electric_current_d)
+      QUANTITY_DEFINE_LITERALS(K, thermodynamic_temperature_d)
+      QUANTITY_DEFINE_LITERALS(mol, amount_of_substance_d)
+      QUANTITY_DEFINE_LITERALS(cd, luminous_intensity_d)
+
+    } // namespace literals
+
+  } // namespace units
+} // namespace phys
+
+#endif // PHYS_UNITS_QUANTITY_HPP_INCLUDED
+
+/*
+ * end of file
+ */

third_party/phys/units/quantity_io.hpp → corsika/framework/units/quantity_io.hpp

@@ -22,7 +22,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_HPP_INCLUDED
 
-#include "phys/units/quantity.hpp"
+#include "corsika/framework/units/quantity.hpp"
 
 #include <algorithm>
 #include <iosfwd>
@@ -81,6 +81,10 @@ inline Rep prefix( std::string const prefix_ )
         { "da", deka },
         { "d", deci  },
         { "c", centi },
+        { "r", ronto },
+        { "q", quecto },
+        { "R", ronna },
+        { "Q", quetta },
     };
 
     auto pos = table.find( prefix_ );
@@ -131,6 +135,7 @@ struct unit_info
         emit_dim( os, "K",   Dims::dim5, first );
         emit_dim( os, "mol", Dims::dim6, first );
         emit_dim( os, "cd",  Dims::dim7, first );
+        emit_dim( os, "eV",  Dims::dim8, first );
 
         return os.str();
     }

third_party/phys/units/quantity_io_ampere.hpp → corsika/framework/units/quantity_io_ampere.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_AMPERE_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_AMPERE_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {
 

third_party/phys/units/quantity_io_becquerel.hpp → corsika/framework/units/quantity_io_becquerel.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_BECQUEREL_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_BECQUEREL_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {
 

third_party/phys/units/quantity_io_candela.hpp → corsika/framework/units/quantity_io_candela.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_CANDELA_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_CANDELA_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {
 

third_party/phys/units/quantity_io_celsius.hpp → corsika/framework/units/quantity_io_celsius.hpp

@@ -17,19 +17,19 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_CELSIUS_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_CELSIUS_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {
 
 /**
- * celsius, [C].
+ * celsius, [C].
  */
 template<>
 struct unit_info< thermodynamic_temperature_d >
 {
     static bool        single() { return true; }
     static std::string name()   { return "celsius"; }
-    static std::string symbol() { return "C"; }
+    static std::string symbol() { return "C"; }
 };
 
 namespace literals {

third_party/phys/units/quantity_io_coulomb.hpp → corsika/framework/units/quantity_io_coulomb.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_COULOMB_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_COULOMB_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {
 

third_party/phys/units/quantity_io_dimensionless.hpp → corsika/framework/units/quantity_io_dimensionless.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_DIMENSIONLESS_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_DIMENSIONLESS_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {
 

third_party/phys/units/quantity_io_engineering.hpp → corsika/framework/units/quantity_io_engineering.hpp

@@ -23,7 +23,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_ENGINEERING_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_ENGINEERING_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 #include <cmath>
 #include <iomanip>
@@ -31,7 +31,7 @@
 #include <sstream>
 
 /*
- * Note: micro, , may not work everywhere, so you can define a glyph yourself:
+ * Note: micro, may not work everywhere, so you can define a glyph yourself:
  */
 #ifndef ENG_FORMAT_MICRO_GLYPH
 # define ENG_FORMAT_MICRO_GLYPH "u"

third_party/phys/units/quantity_io_farad.hpp → corsika/framework/units/quantity_io_farad.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_FARAD_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_FARAD_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {
 

third_party/phys/units/quantity_io_gray.hpp → corsika/framework/units/quantity_io_gray.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_GRAY_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_GRAY_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {
 

third_party/phys/units/quantity_io_henry.hpp → corsika/framework/units/quantity_io_henry.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_HENRY_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_HENRY_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {
 

third_party/phys/units/quantity_io_hertz.hpp → corsika/framework/units/quantity_io_hertz.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_HERTZ_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_HERTZ_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {