diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 9167ffcd5aafd01f047c046dd093b7928fcc0ad8..49e2c11b77e2704c84297600a21f5e1933594227 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -34,6 +34,7 @@
#include <boost/type_index.hpp>
using boost::typeindex::type_id_with_cvr;
+#include <fenv.h>
#include <iostream>
#include <limits>
#include <typeinfo>
@@ -48,25 +49,26 @@ using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
-using namespace corsika::units::hep;
+using namespace corsika::units::si;
class ProcessCut : public corsika::process::ContinuousProcess<ProcessCut> {
- EnergyType fECut;
+ HEPEnergyType fECut;
- EnergyType fEnergy = 0_GeV;
- EnergyType fEmEnergy = 0_GeV;
+ HEPEnergyType fEnergy = 0_GeV;
+ HEPEnergyType fEmEnergy = 0_GeV;
int fEmCount = 0;
- EnergyType fInvEnergy = 0_GeV;
+ HEPEnergyType fInvEnergy = 0_GeV;
int fInvCount = 0;
public:
- ProcessCut(const EnergyType v)
- : fECut(v) {}
+ ProcessCut(const HEPEnergyType cut)
+ : fECut(cut) {}
+
template <typename Particle>
bool isBelowEnergyCut(Particle& p) const {
// FOR NOW: center-of-mass energy hard coded
- const EnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
+ const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
if (p.GetEnergy() < fECut || Ecm < 10_GeV)
return true;
else
@@ -143,7 +145,7 @@ public:
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) {
const Code pid = p.GetPID();
- EnergyType energy = p.GetEnergy();
+ HEPEnergyType energy = p.GetEnergy();
cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy
<< ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl;
EProcessReturn ret = EProcessReturn::eOk;
@@ -188,15 +190,16 @@ public:
<< " ******************************" << endl;
}
- EnergyType GetInvEnergy() const { return fInvEnergy; }
- EnergyType GetCutEnergy() const { return fEnergy; }
- EnergyType GetEmEnergy() const { return fEmEnergy; }
+ HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
+ HEPEnergyType GetCutEnergy() const { return fEnergy; }
+ HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
//
// The example main program for a particle cascade
//
int main() {
+ feenableexcept(FE_INVALID);
// initialize random number sequence(s)
corsika::random::RNGManager::GetInstance().RegisterRandomStream("cascade");
@@ -237,14 +240,14 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const hep::EnergyType E0 = 100_TeV;
+ const HEPEnergyType E0 = 100_TeV;
double theta = 0.;
double phi = 0.;
{
auto particle = stack.NewParticle();
particle.SetPID(Code::Proton);
- hep::MomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass());
- auto momentumComponents = [](double theta, double phi, MomentumType& ptot) {
+ HEPMomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass());
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
diff --git a/Documentation/Examples/stack_example.cc b/Documentation/Examples/stack_example.cc
index fbe4cd9da69440c48d345c17a1affc52d6c0f150..33ebc230af244932375c904747b127b0b5c240ab 100644
--- a/Documentation/Examples/stack_example.cc
+++ b/Documentation/Examples/stack_example.cc
@@ -30,7 +30,7 @@ void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
void read(corsika::stack::super_stupid::SuperStupidStack& s) {
assert(s.GetSize() == 11); // stack has 11 particles
- EnergyType total_energy;
+ HEPEnergyType total_energy;
int i = 0;
for (auto& p : s) {
total_energy += p.GetEnergy();
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index e1b742624010bdf549008c82669ff637de54443b..12885054b18cafdb1312818418f0bf46e95c9a88 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -81,7 +81,7 @@ public:
template <typename Particle, typename T, typename Stack>
EProcessReturn DoContinuous(Particle& p, T&, Stack& s) const {
- EnergyType E = p.GetEnergy();
+ HEPEnergyType E = p.GetEnergy();
if (E < 85_MeV) {
p.Delete();
fCount++;
@@ -123,7 +123,7 @@ TEST_CASE("Cascade", "[Cascade]") {
stack.Clear();
auto particle = stack.NewParticle();
- EnergyType E0 = 100_GeV;
+ HEPEnergyType E0 = 100_GeV;
particle.SetPID(particles::Code::Electron);
particle.SetEnergy(E0);
particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
@@ -138,7 +138,7 @@ TEST_CASE("Cascade", "[Cascade]") {
for (int i = 0; i < 0; ++i) {
stack.Clear();
auto particle = stack.NewParticle();
- EnergyType E0 = 100_GeV * pow(10, i);
+ HEPEnergyType E0 = 100_GeV * pow(10, i);
particle.SetEnergy(E0);
EAS.Init();
EAS.Run();
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index 44b9b55ed70cfb00e6ed34e3d5a0ceb55da8907b..959c27a6be0028c607d9114d4577adcd6e791c3c 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -44,7 +44,7 @@ namespace corsika::particles {
// forward declarations to be used in GeneratedParticleProperties
int16_t constexpr GetElectricChargeNumber(Code const);
corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const);
- corsika::units::hep::MassType constexpr GetMass(Code const);
+ corsika::units::si::HEPMassType constexpr GetMass(Code const);
PDGCodeType constexpr GetPDG(Code const);
constexpr std::string const& GetName(Code const);
corsika::units::si::TimeType constexpr GetLifetime(Code const);
@@ -58,7 +58,7 @@ namespace corsika::particles {
/*!
* returns mass of particle
*/
- corsika::units::hep::MassType constexpr GetMass(Code const p) {
+ corsika::units::si::HEPMassType constexpr GetMass(Code const p) {
return detail::masses[static_cast<CodeIntType const>(p)];
}
diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py
index b1a11ae8ecedb9898348dfa27c26a24a8b389221..2e463a0e1d7c9a227d9f49160d4916ecf3390761 100755
--- a/Framework/Particles/pdxml_reader.py
+++ b/Framework/Particles/pdxml_reader.py
@@ -290,9 +290,9 @@ def gen_properties(particle_db):
string += "\n"
# particle masses table
- string += "static constexpr std::array<corsika::units::hep::MassType const, size> masses = {\n"
+ string += "static constexpr std::array<corsika::units::si::HEPMassType const, size> masses = {\n"
for p in particle_db.values():
- string += " {mass:e} * (1e9 * corsika::units::constants::eV), // {name:s}\n".format(mass = p['mass'], name = p['name'])
+ string += " {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(mass = p['mass'], name = p['name'])
string += "};\n\n"
# PDG code table
@@ -343,7 +343,7 @@ def gen_properties(particle_db):
string += "};\n"
# nucleus mass number A
- string += "static constexpr std::array<int, size> nucleusA = {\n"
+ string += "static constexpr std::array<int16_t, size> nucleusA = {\n"
for p in particle_db.values():
A = 0
if p['isNucleus']:
@@ -352,7 +352,7 @@ def gen_properties(particle_db):
string += "};\n"
# nucleus charge number Z
- string += "static constexpr std::array<int, size> nucleusZ = {\n"
+ string += "static constexpr std::array<int16_t, size> nucleusZ = {\n"
for p in particle_db.values():
Z = 0
if p['isNucleus']:
@@ -396,14 +396,14 @@ def gen_classes(particle_db):
string += "class " + cname + " {\n"
string += " public:\n"
string += " static constexpr Code GetCode() { return Type; }\n"
- string += " static constexpr corsika::units::hep::MassType GetMass() { return corsika::particles::GetMass(Type); }\n"
+ string += " static constexpr corsika::units::si::HEPMassType GetMass() { return corsika::particles::GetMass(Type); }\n"
string += " static constexpr corsika::units::si::ElectricChargeType GetCharge() { return corsika::particles::GetElectricCharge(Type); }\n"
string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetElectricChargeNumber(Type); }\n"
string += " static std::string const& GetName() { return corsika::particles::GetName(Type); }\n"
string += " static constexpr Code GetAntiParticle() { return AntiType; }\n"
string += " static constexpr bool IsNucleus() { return corsika::particles::IsNucleus(Type); }\n"
- string += " static constexpr int GetNucleusA() { return corsika::particles::GetNucleusA(Type); }\n"
- string += " static constexpr int GetNucleusZ() { return corsika::particles::GetNucleusZ(Type); }\n"
+ string += " static constexpr int16_t GetNucleusA() { return corsika::particles::GetNucleusA(Type); }\n"
+ string += " static constexpr int16_t GetNucleusZ() { return corsika::particles::GetNucleusZ(Type); }\n"
string += " static constexpr Code Type = Code::" + cname + ";\n"
string += " static constexpr Code AntiType = Code::" + antiP + ";\n"
string += " private:\n"
@@ -463,7 +463,7 @@ if __name__ == "__main__":
print("usage: {:s} <Pythia8.xml> <Nuclei.xml> <ClassNames.xml>".format(sys.argv[0]), file=sys.stderr)
sys.exit(1)
- print("\n pdxml_reader.py: Automatically produce particle-properties from input files\n")
+ print("\n pdxml_reader.py: automatically produce particle properties from input files\n")
names = class_names(sys.argv[3])
particle_db = OrderedDict()
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index f211260541665a426fd641eb506ef444869f92d1..bdd07d8370592da981e59672c8b7a50ef92df9e4 100644
--- a/Framework/Particles/testParticles.cc
+++ b/Framework/Particles/testParticles.cc
@@ -17,7 +17,7 @@
#include <catch2/catch.hpp>
using namespace corsika::units;
-using namespace corsika::units::hep;
+using namespace corsika::units::si;
using namespace corsika::particles;
TEST_CASE("ParticleProperties", "[Particles]") {
diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h
index 9aef04ee6fa63d6818d5d471f342d7a8e2b3a729..95436959107ebee231450d5fb0b237a08b0886e6 100644
--- a/Framework/Units/PhysicalConstants.h
+++ b/Framework/Units/PhysicalConstants.h
@@ -43,7 +43,7 @@ namespace corsika::units::constants {
constexpr quantity<electric_charge_d> e{Rep(1.6021766208e-19L) * coulomb};
// electronvolt
- constexpr quantity<energy_d> eV{e / coulomb * joule};
+ // constexpr quantity<hepenergy_d> eV{e / coulomb * joule};
// Planck constant
constexpr quantity<dimensions<2, 1, -1> > h{Rep(6.62606876e-34L) * joule * second};
@@ -55,9 +55,6 @@ namespace corsika::units::constants {
// unified atomic mass unit
constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
- // barn moved to PhysicalUnits
- // constexpr quantity<area_d> barn{Rep(1.e-28L) * meter * meter};
-
// etc.
} // namespace corsika::units::constants
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index fe4e284acd7f039147acebb39512390bcef1431c..262e55d1c52070e36add9404f0d2ba67bae347bc 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -23,7 +23,7 @@ namespace phys::units {
* optimal coding style.
*
*/
-
+/*
namespace corsika::units::hep {
using namespace phys::units;
using namespace phys::units::literals;
@@ -39,7 +39,7 @@ namespace corsika::units::hep {
using EnergyType = phys::units::quantity<energy_hep_d, double>;
} // namespace corsika::units::hep
-
+*/
namespace corsika::units::si {
using namespace phys::units;
using namespace phys::units::literals;
@@ -48,15 +48,11 @@ namespace corsika::units::si {
/// defining momentum you suckers
/// dimensions, i.e. composition in base SI dimensions
- using momentum_d = phys::units::dimensions<1, 1, -1>;
- // defining the unit of momentum, so far newton-meter, maybe go to HEP?
- constexpr phys::units::quantity<momentum_d> newton_second{
- phys::units::meter * phys::units::kilogram / phys::units::second};
+ using hepmomentum_d = phys::units::hepenergy_d;
+ using hepmass_d = phys::units::hepenergy_d;
/// defining cross section
- using sigma_d = phys::units::dimensions<2, 0, 0>;
- constexpr phys::units::quantity<sigma_d> barn{phys::units::Rep(1.e-28L) *
- phys::units::meter * phys::units::meter};
+ using sigma_d = phys::units::area_d;
/// add the unit-types
using LengthType = phys::units::quantity<phys::units::length_d, double>;
@@ -65,11 +61,12 @@ namespace corsika::units::si {
using FrequencyType = phys::units::quantity<phys::units::frequency_d, double>;
using ElectricChargeType =
phys::units::quantity<phys::units::electric_charge_d, double>;
- using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
+ using HEPEnergyType = phys::units::quantity<phys::units::hepenergy_d, double>;
using MassType = phys::units::quantity<phys::units::mass_d, double>;
+ using HEPMassType = phys::units::quantity<hepmass_d, double>;
using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
using GrammageType = phys::units::quantity<phys::units::dimensions<-2, 1, 0>, double>;
- using MomentumType = phys::units::quantity<momentum_d, double>;
+ using HEPMomentumType = phys::units::quantity<hepmomentum_d, double>;
using CrossSectionType = phys::units::quantity<area_d, double>;
using InverseLengthType =
phys::units::quantity<phys::units::dimensions<-1, 0, 0>, double>;
@@ -83,28 +80,19 @@ namespace corsika::units::si {
* @file PhysicalUnits
*
* Define _XeV literals, alowing 10_GeV in the code.
- * Define _meter literal
* Define _barn literal
- * Define _newton_second literal for SI momenta
*/
namespace phys {
namespace units {
namespace literals {
- QUANTITY_DEFINE_SCALING_LITERALS(eV, energy_d,
- magnitude(corsika::units::constants::eV))
-
- // QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::area_d,
- // magnitude(corsika::units::si::constants::barn))
+ QUANTITY_DEFINE_SCALING_LITERALS(eV, hepenergy_d, 1)
QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::sigma_d,
magnitude(corsika::units::constants::barn))
- QUANTITY_DEFINE_SCALING_LITERALS(meter, length_d,
- magnitude(corsika::units::constants::meter))
-
- QUANTITY_DEFINE_SCALING_LITERALS(Ns, corsika::units::si::momentum_d,
- magnitude(1_m * 1_kg / 1_s))
+ // QUANTITY_DEFINE_SCALING_LITERALS(Ns, corsika::units::si::momentum_d,
+ // magnitude(1_m * 1_kg / 1_s))
} // namespace literals
} // namespace units
diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc
index 4dc6316e1cdaf63170475d417199a2736b11d1a9..38b3b299e151274bc8c03a99ebf750ec4ee46bcd 100644
--- a/Framework/Units/testUnits.cc
+++ b/Framework/Units/testUnits.cc
@@ -1,4 +1,3 @@
-
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -16,6 +15,7 @@
#include <corsika/units/PhysicalUnits.h>
#include <array>
+#include <sstream>
using namespace corsika;
using namespace corsika::units::si;
@@ -39,12 +39,12 @@ TEST_CASE("PhysicalUnits", "[Units]") {
[[maybe_unused]] LengthType arr1[2] = {{1_mm}, {2_cm}};
- [[maybe_unused]] std::array<EnergyType, 4> arr2; // empty array
+ [[maybe_unused]] std::array<HEPEnergyType, 4> arr2; // empty array
- [[maybe_unused]] std::array<EnergyType, 4> arr3 = {1_GeV, 1_eV, 5_MeV};
+ [[maybe_unused]] std::array<HEPEnergyType, 4> arr3 = {1_GeV, 1_eV, 5_MeV};
- auto p1 = 10_Ns;
- REQUIRE(p1 == 10_Ns);
+ auto p1 = 10_s * newton;
+ REQUIRE(p1 == 10_s * newton);
}
SECTION("Powers in literal units") {
@@ -79,13 +79,13 @@ TEST_CASE("PhysicalUnits", "[Units]") {
}
SECTION("Formulas") {
- const EnergyType E2 = 20_GeV * 2;
+ const HEPEnergyType E2 = 20_GeV * 2;
REQUIRE(E2 == 40_GeV);
REQUIRE(E2 / 1_GeV == Approx(40));
const MassType m = 1_kg;
const SpeedType v = 1_m / 1_s;
- REQUIRE(m * v == 1_Ns);
+ REQUIRE(m * v == 1_s * newton);
const double lgE = log10(E2 / 1_GeV);
REQUIRE(lgE == Approx(log10(40.)));
@@ -94,21 +94,19 @@ TEST_CASE("PhysicalUnits", "[Units]") {
REQUIRE(E3 == 180_GeV);
}
- SECTION("Unit system conversion") {
-
- const units::hep::MassType m_hep = 3_GeV;
- std::cout << m_hep << std::endl;
-
- const units::si::MassType m_hep2 = 3_kg;
- std::cout << m_hep2 << std::endl;
-
- REQUIRE(m_hep == 3_GeV); // hep::mass identical to si::energy
- auto type_check = m_hep / units::constants::cSquared;
- REQUIRE(dynamic_cast<units::si::MassType*>(&type_check)); // hep::mass*c2 is mass unit
-
- const units::hep::EnergyType e_hep = 4_GeV;
-
- REQUIRE(sqrt(m_hep * m_hep + e_hep * e_hep) == 5_GeV);
+ SECTION("Output") {
+ {
+ const HEPEnergyType E = 5_eV;
+ std::stringstream stream;
+ stream << E;
+ REQUIRE(stream.str() == std::string("5 eV"));
+ }
+ {
+ const HEPEnergyType E = 5_EeV;
+ std::stringstream stream;
+ stream << E;
+ REQUIRE(stream.str() == std::string("5e+18 eV"));
+ }
}
SECTION("Special") {
diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index 42d5995d2e54c01ef0ca5468dd98b1005d70fb3a..1bf8f839f1818fb12aa5c813ee166504e3351bd1 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/Framework/Utilities/COMBoost.cc
@@ -13,8 +13,8 @@
using namespace corsika::utl;
using namespace corsika::units::si;
-COMBoost::COMBoost(EnergyType eProjectile, COMBoost::MomentumVector const& pProjectile,
- MassType mTarget)
+COMBoost::COMBoost(HEPEnergyType eProjectile, COMBoost::MomentumVector const& pProjectile,
+ HEPMassType mTarget)
: fRotation(Eigen::Matrix3d::Identity())
, fCS(pProjectile.GetCoordinateSystem()) {
// calculate matrix for rotating pProjectile to z-axis first
@@ -37,8 +37,7 @@ COMBoost::COMBoost(EnergyType eProjectile, COMBoost::MomentumVector const& pProj
}
// calculate boost
- double const x = pProjNorm * units::constants::c /
- (eProjectile + mTarget * units::constants::cSquared);
+ double const x = pProjNorm / (eProjectile + mTarget);
/* Accurracy matters here, x = 1 - epsilon for ultra-relativistic boosts */
double const coshEta = 1 / std::sqrt((1 + x) * (1 - x));
@@ -50,34 +49,34 @@ COMBoost::COMBoost(EnergyType eProjectile, COMBoost::MomentumVector const& pProj
fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
}
-std::tuple<EnergyType, corsika::geometry::QuantityVector<momentum_d>> COMBoost::toCoM(
- EnergyType E, COMBoost::MomentumVector p) const {
- corsika::geometry::QuantityVector<momentum_d> pComponents = p.GetComponents(fCS);
+std::tuple<HEPEnergyType, corsika::geometry::QuantityVector<hepmomentum_d>>
+COMBoost::toCoM(HEPEnergyType E, COMBoost::MomentumVector p) const {
+ corsika::geometry::QuantityVector<hepmomentum_d> pComponents = p.GetComponents(fCS);
Eigen::Vector3d eVecRotated = fRotation * pComponents.eVector;
Eigen::Vector2d lab;
- lab << (E * (1 / 1_GeV)), (eVecRotated(2) * (units::constants::c / 1_GeV).magnitude());
+ lab << (E * (1 / 1_GeV)), (eVecRotated(2) * (1 / 1_GeV).magnitude());
auto const boostedZ = fBoost * lab;
auto const E_CoM = boostedZ(0) * 1_GeV;
- eVecRotated(2) = boostedZ(1) * (1_GeV / units::constants::c).magnitude();
+ eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
return std::make_tuple(E_CoM,
- corsika::geometry::QuantityVector<momentum_d>{eVecRotated});
+ corsika::geometry::QuantityVector<hepmomentum_d>{eVecRotated});
}
-std::tuple<EnergyType, COMBoost::MomentumVector> COMBoost::fromCoM(
- EnergyType E, corsika::geometry::QuantityVector<units::si::momentum_d> pCoM) const {
+std::tuple<HEPEnergyType, COMBoost::MomentumVector> COMBoost::fromCoM(
+ HEPEnergyType E,
+ corsika::geometry::QuantityVector<units::si::hepmomentum_d> pCoM) const {
Eigen::Vector2d com;
- com << (E * (1 / (units::constants::c * 1_Ns))),
- (pCoM.eVector(2) * (1 / 1_Ns).magnitude());
+ com << (E * (1 / 1_GeV)), (pCoM.eVector(2) * (1 / 1_GeV).magnitude());
auto const boostedZ = fInverseBoost * com;
- auto const E_CoM = boostedZ(0) * (1_Ns * units::constants::c);
+ auto const E_CoM = boostedZ(0) * 1_GeV;
auto pLab = pCoM;
- pLab.eVector(2) = boostedZ(1) * (1_Ns).magnitude();
+ pLab.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
pLab.eVector = fRotation.transpose() * pLab.eVector;
return std::make_tuple(E_CoM, MomentumVector(fCS, pLab));
diff --git a/Framework/Utilities/COMBoost.h b/Framework/Utilities/COMBoost.h
index 6250d422dc256c4cad102bb93bd596e9a0df2a30..fa57d8f3fa56c6c34432ed2dbdbbec3db14c5a14 100644
--- a/Framework/Utilities/COMBoost.h
+++ b/Framework/Utilities/COMBoost.h
@@ -23,21 +23,22 @@ namespace corsika::utl {
Eigen::Matrix3d fRotation;
Eigen::Matrix2d fBoost, fInverseBoost;
corsika::geometry::CoordinateSystem const& fCS;
- using MomentumVector = corsika::geometry::Vector<units::si::momentum_d>;
+ using MomentumVector = corsika::geometry::Vector<units::si::hepmomentum_d>;
public:
//! construct a COMBoost given energy and momentum of projectile and mass of target
- COMBoost(units::si::EnergyType eProjectile, MomentumVector const& pProjectile,
- units::si::MassType mTarget);
+ COMBoost(units::si::HEPEnergyType eProjectile, MomentumVector const& pProjectile,
+ units::si::HEPMassType mTarget);
//! transforms a 4-momentum from lab frame to the center-of-mass frame
- std::tuple<units::si::EnergyType, geometry::QuantityVector<units::si::momentum_d>>
- toCoM(units::si::EnergyType E, MomentumVector p) const;
+ std::tuple<units::si::HEPEnergyType,
+ geometry::QuantityVector<units::si::hepmomentum_d>>
+ toCoM(units::si::HEPEnergyType E, MomentumVector p) const;
//! transforms a 4-momentum from the center-of-mass frame back to lab frame
- std::tuple<units::si::EnergyType, MomentumVector> fromCoM(
- units::si::EnergyType E,
- geometry::QuantityVector<units::si::momentum_d> pCoM) const;
+ std::tuple<units::si::HEPEnergyType, MomentumVector> fromCoM(
+ units::si::HEPEnergyType E,
+ geometry::QuantityVector<units::si::hepmomentum_d> pCoM) const;
};
} // namespace corsika::utl
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index 3af2585630be330acd547c1087b55c174732fcc0..6b3b950518bd5a488910bbbd3538f93f1cda8721 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/Framework/Utilities/testCOMBoost.cc
@@ -30,57 +30,52 @@ TEST_CASE("boosts") {
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// relativistic energy
- auto energy = [](MassType m, Vector<momentum_d> const& p) {
- return sqrt(m * m * cSquared * cSquared + p.squaredNorm() * cSquared);
+ auto energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
+ return sqrt(m * m + p.squaredNorm());
};
// mandelstam-s
- auto s = [](EnergyType E, QuantityVector<momentum_d> const& p) {
- return E * E / cSquared - p.squaredNorm();
+ auto s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
+ return E * E - p.squaredNorm();
};
// define projectile kinematics in lab frame
- MassType const projectileMass = 1._GeV / cSquared;
- std::vector<Vector<momentum_d>> labProjectiles{
- {rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}}, // standard case
- {rootCS, {0_GeV / c, 0_GeV / c, -1_GeV / c}}}; // "special" case
-
- for (auto const& pProjectileLab : labProjectiles) {
- EnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
-
- // define target kinematics in lab frame
- MassType const targetMass = 1_GeV / cSquared;
- Vector<momentum_d> pTargetLab{rootCS, {0_Ns, 0_Ns, 0_Ns}};
- EnergyType const eTargetLab = energy(targetMass, pTargetLab);
-
- // define boost to com frame
- COMBoost boost(eProjectileLab, pProjectileLab, targetMass);
-
- // boost projecticle
- auto const [eProjectileCoM, pProjectileCoM] =
- boost.toCoM(eProjectileLab, pProjectileLab);
-
- // boost target
- auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
-
- // sum of momenta in CoM, should be 0
- auto const sumPCoM = pProjectileCoM + pTargetCoM;
- CHECK(sumPCoM[2] / (1_GeV / c) == Approx(0).margin(absMargin));
-
- // mandelstam-s should be invariant under transformation
- CHECK(s(eProjectileLab + eTargetLab,
- pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
- (1_GeV / c) / (1_GeV / c) ==
- Approx(s(eProjectileCoM + eTargetCoM, pProjectileCoM + pTargetCoM) /
- (1_GeV / c) / (1_GeV / c)));
-
- // boost back...
- auto const [eProjectileBack, pProjectileBack] =
- boost.fromCoM(eProjectileCoM, pProjectileCoM);
-
- // ...should yield original values before the boosts
- CHECK(eProjectileBack / eProjectileLab == Approx(1));
- CHECK((pProjectileBack - pProjectileLab).norm() / pProjectileLab.norm() ==
- Approx(0).margin(absMargin));
- }
+ HEPMassType const projectileMass = 1._GeV;
+ Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 1_PeV, 0_GeV}};
+ HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+
+ // define target kinematics in lab frame
+ HEPMassType const targetMass = 1_GeV;
+ Vector<hepmomentum_d> pTargetLab{rootCS, {0_eV, 0_eV, 0_eV}};
+ HEPEnergyType const eTargetLab = energy(targetMass, pTargetLab);
+
+ // define boost to com frame
+ COMBoost boost(eProjectileLab, pProjectileLab, targetMass);
+
+ // boost projecticle
+ auto const [eProjectileCoM, pProjectileCoM] =
+ boost.toCoM(eProjectileLab, pProjectileLab);
+
+ // boost target
+ auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
+
+ // sum of momenta in CoM, should be 0
+ auto const sumPCoM = pProjectileCoM + pTargetCoM;
+ CHECK(sumPCoM[2] / 1_GeV == Approx(0).margin(absMargin));
+
+ // mandelstam-s should be invariant under transformation
+ CHECK(s(eProjectileLab + eTargetLab,
+ pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
+ 1_GeV / 1_GeV ==
+ Approx(s(eProjectileCoM + eTargetCoM, pProjectileCoM + pTargetCoM) / 1_GeV /
+ 1_GeV));
+
+ // boost back...
+ auto const [eProjectileBack, pProjectileBack] =
+ boost.fromCoM(eProjectileCoM, pProjectileCoM);
+
+ // ...should yield original values before the boosts
+ CHECK(eProjectileBack / eProjectileLab == Approx(1));
+ CHECK((pProjectileBack - pProjectileLab).norm() / pProjectileLab.norm() ==
+ Approx(0).margin(absMargin));
}
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 13c6a3ebc094b7618ca694fb7b2b74786571cdcc..a798f318db53e994ba7a1edff828e9010dec979d 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -136,8 +136,8 @@ namespace corsika::process {
using std::endl;
using namespace corsika::units::si;
- corsika::units::hep::EnergyType E = p.GetEnergy();
- corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID());
+ corsika::units::si::HEPEnergyType E = p.GetEnergy();
+ corsika::units::si::HEPMassType m = corsika::particles::GetMass(p.GetPID());
const double gamma = E / m;
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index fef98eb924f2a0a2d537766eebc3cfa309742aa1..fac0dd578b77cfea776b1b942b6a42b9e899f809 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -53,7 +53,6 @@ namespace corsika::process::sibyll {
corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) {
using namespace corsika::units;
- using namespace corsika::units::hep;
using namespace corsika::units::si;
using namespace corsika::geometry;
using std::cout;
@@ -79,26 +78,27 @@ namespace corsika::process::sibyll {
// FOR NOW: assume target is oxygen
const int kTarget = corsika::particles::Oxygen::GetNucleusA();
- const hep::MassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
- hep::EnergyType Etot = p.GetEnergy() + nucleon_mass;
+ const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
+ corsika::particles::Neutron::GetMass());
+ HEPEnergyType const Elab = p.GetEnergy();
+ HEPEnergyType const Etot = Elab + nucleon_mass;
MomentumVector Ptot(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
Ptot += p.GetMomentum();
Ptot += pTarget;
// calculate cm. energy
- const hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
+ const HEPEnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
const double Ecm = sqs / 1_GeV;
std::cout << "Interaction: LambdaInt: \n"
- << " input energy: " << p.GetEnergy() / 1_GeV << endl
+ << " input energy: " << Elab / 1_GeV << endl
<< " beam can interact:" << kBeam << endl
<< " beam XS code:" << kBeam << endl
<< " beam pid:" << p.GetPID() << endl
<< " target mass number:" << kTarget << std::endl;
- if (kInteraction) {
+ if (kInteraction && Elab >= 8.5_GeV && sqs >= 10_GeV) {
double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
double dumdif[3];
@@ -114,12 +114,10 @@ namespace corsika::process::sibyll {
std::cout << "Interaction: "
<< "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
- const si::MassType nucleon_mass =
- 0.93827_GeV / corsika::units::constants::cSquared;
std::cout << "Interaction: "
<< "nucleon mass " << nucleon_mass << std::endl;
// calculate interaction length in medium
- GrammageType int_length = kTarget * nucleon_mass / sig;
+ GrammageType int_length = kTarget * corsika::units::constants::u / sig;
std::cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length << std::endl;
@@ -134,7 +132,6 @@ namespace corsika::process::sibyll {
using namespace corsika::units;
using namespace corsika::utl;
- using namespace corsika::units::hep;
using namespace corsika::units::si;
using namespace corsika::geometry;
using std::cout;
@@ -188,9 +185,9 @@ namespace corsika::process::sibyll {
GetNucleusA(); // env.GetTargetParticle().GetPID();
// FOR NOW: target is always at rest
- const hep::MassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
- corsika::particles::Neutron::GetMass());
- const EnergyType Etarget = 0_GeV + nucleon_mass;
+ auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
+ corsika::particles::Neutron::GetMass());
+ const auto Etarget = 0_GeV + nucleon_mass;
const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
@@ -210,12 +207,12 @@ namespace corsika::process::sibyll {
*/
// total energy: E_beam + E_target
// in lab. frame: E_beam + m_target*c**2
- EnergyType E = p.GetEnergy();
- EnergyType Etot = E + Etarget;
+ HEPEnergyType E = p.GetEnergy();
+ HEPEnergyType Etot = E + Etarget;
// total momentum
MomentumVector Ptot = p.GetMomentum();
// invariant mass, i.e. cm. energy
- EnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
+ HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
/*
get transformation between Stack-frame and SibStack-frame
for EAS Stack-frame is lab. frame, could be different for CRMC-mode
@@ -229,25 +226,25 @@ namespace corsika::process::sibyll {
std::cout << "Interaction: "
<< " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
- Vector<si::momentum_d> pProjectileLab = p.GetMomentum() / constants::c;
+ auto const pProjectileLab = p.GetMomentum();
//{rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}};
- EnergyType const eProjectileLab = p.GetEnergy();
+ HEPEnergyType const eProjectileLab = p.GetEnergy();
//energy(projectileMass, pProjectileLab);
// define target kinematics in lab frame
- si::MassType const targetMass = nucleon_mass / constants::cSquared;
+ HEPMassType const targetMass = nucleon_mass;
// define boost to com frame
- COMBoost boost(eProjectileLab, pProjectileLab, targetMass);
+ COMBoost const boost(eProjectileLab, pProjectileLab, targetMass);
cout << "Interaction: new boost: ebeam lab: " << eProjectileLab / 1_GeV << endl
- << "Interaction: new boost: pbeam lab: " << pProjectileLab.GetComponents() / ( 1_GeV / constants::c ) << endl;
+ << "Interaction: new boost: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV << endl;
// boost projecticle
auto const [eProjectileCoM, pProjectileCoM] =
boost.toCoM(eProjectileLab, pProjectileLab);
cout << "Interaction: new boost: ebeam com: " << eProjectileCoM / 1_GeV << endl
- << "Interaction: new boost: pbeam com: " << pProjectileCoM / ( 1_GeV / constants::c ) << endl;
+ << "Interaction: new boost: pbeam com: " << pProjectileCoM / 1_GeV << endl;
int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
@@ -283,7 +280,7 @@ namespace corsika::process::sibyll {
// momentum array in Sibyll
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- EnergyType E_final = 0_GeV, Ecm_final = 0_GeV;
+ HEPEnergyType E_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
// skip particles that have decayed in Sibyll
@@ -302,10 +299,10 @@ namespace corsika::process::sibyll {
// rotatate to rootCS
const auto pSibyllComponents = SibVector.GetComponents(rootCS);
// boost to lab. frame
- hep::EnergyType en_lab = 0. * 1_GeV;
- hep::MomentumType p_lab_components[3];
+ HEPEnergyType en_lab = 0. * 1_GeV;
+ HEPMomentumType p_lab_components[3];
en_lab = psib.GetEnergy() * gamma;
- hep::EnergyType pnorm = 0. * 1_GeV;
+ HEPEnergyType pnorm = 0. * 1_GeV;
for (int j = 0; j < 3; ++j)
pnorm += (pSibyllComponents[j] * gammaBetaComponents[j]) / (gamma + 1.);
pnorm += psib.GetEnergy();
@@ -320,7 +317,7 @@ namespace corsika::process::sibyll {
auto pnew = s.NewParticle();
pnew.SetEnergy(en_lab);
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
- corsika::geometry::QuantityVector<energy_hep_d> p_lab_c{
+ corsika::geometry::QuantityVector<hepmomentum_d> p_lab_c{
p_lab_components[0], p_lab_components[1], p_lab_components[2]};
pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
pnew.SetPosition(pOrig);
diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h
index bc3dc4c0fba66fadba4e277c0e305f01b65b1131..83d902f5c2b030ed455afbbe38867fa8b50dbf3c 100644
--- a/Processes/Sibyll/SibStack.h
+++ b/Processes/Sibyll/SibStack.h
@@ -21,7 +21,7 @@
namespace corsika::process::sibyll {
- typedef corsika::geometry::Vector<corsika::units::hep::energy_hep_d> MomentumVector;
+ typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
class SibStackData {
@@ -35,30 +35,29 @@ namespace corsika::process::sibyll {
int GetCapacity() const { return 8000; }
void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
- void SetEnergy(const int i, const corsika::units::hep::EnergyType v) {
- using namespace corsika::units::hep;
+ void SetEnergy(const int i, const corsika::units::si::HEPEnergyType v) {
+ using namespace corsika::units::si;
s_plist_.p[3][i] = v / 1_GeV;
}
- void SetMass(const int i, const corsika::units::hep::MassType v) {
- using namespace corsika::units::hep;
+ void SetMass(const int i, const corsika::units::si::HEPMassType v) {
+ using namespace corsika::units::si;
s_plist_.p[4][i] = v / 1_GeV;
}
void SetMomentum(const int i, const MomentumVector& v) {
- using namespace corsika::units;
- using namespace corsika::units::hep;
+ using namespace corsika::units::si;
auto tmp = v.GetComponents();
for (int idx = 0; idx < 3; ++idx) s_plist_.p[idx][i] = tmp[idx] / 1_GeV;
}
int GetId(const int i) const { return s_plist_.llist[i]; }
- corsika::units::hep::EnergyType GetEnergy(const int i) const {
- using namespace corsika::units::hep;
+ corsika::units::si::HEPEnergyType GetEnergy(const int i) const {
+ using namespace corsika::units::si;
return s_plist_.p[3][i] * 1_GeV;
}
- corsika::units::hep::EnergyType GetMass(const int i) const {
- using namespace corsika::units::hep;
+ corsika::units::si::HEPEnergyType GetMass(const int i) const {
+ using namespace corsika::units::si;
return s_plist_.p[4][i] * 1_GeV;
}
@@ -66,10 +65,10 @@ namespace corsika::process::sibyll {
using corsika::geometry::CoordinateSystem;
using corsika::geometry::QuantityVector;
using corsika::geometry::RootCoordinateSystem;
- using namespace corsika::units::hep;
+ using namespace corsika::units::si;
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- QuantityVector<energy_hep_d> components = {
+ QuantityVector<hepmomentum_d> components = {
s_plist_.p[0][i] * 1_GeV, s_plist_.p[1][i] * 1_GeV, s_plist_.p[2][i] * 1_GeV};
MomentumVector v1(rootCS, components);
return v1;
@@ -94,29 +93,33 @@ namespace corsika::process::sibyll {
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
public:
- void SetEnergy(const corsika::units::hep::EnergyType v) {
+ void SetEnergy(const corsika::units::si::HEPEnergyType v) {
GetStackData().SetEnergy(GetIndex(), v);
}
- corsika::units::hep::EnergyType GetEnergy() const {
+
+ corsika::units::si::HEPEnergyType GetEnergy() const {
return GetStackData().GetEnergy(GetIndex());
}
- void SetMass(const corsika::units::hep::MassType v) {
+ bool HasDecayed() const { return abs(GetStackData().GetId(GetIndex())) > 100; }
+
+ void SetMass(const corsika::units::si::HEPMassType v) {
GetStackData().SetMass(GetIndex(), v);
}
- corsika::units::hep::EnergyType GetMass() const {
+
+ corsika::units::si::HEPEnergyType GetMass() const {
return GetStackData().GetMass(GetIndex());
}
- bool HasDecayed() const {
- return abs(GetStackData().GetId(GetIndex())) > 100 ? true : false;
- }
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
+
corsika::process::sibyll::SibyllCode GetPID() const {
return static_cast<corsika::process::sibyll::SibyllCode>(
GetStackData().GetId(GetIndex()));
}
+
MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); }
+
void SetMomentum(const MomentumVector& v) {
GetStackData().SetMomentum(GetIndex(), v);
}
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 8c171cdce00edf619cd0fbd968d3495233ca66f9..3d7d425d9c63aaacf5ffcef5723dfac730283b7e 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -106,9 +106,9 @@ TEST_CASE("SibyllInterface", "[processes]") {
setup::Stack stack;
auto particle = stack.NewParticle();
{
- const hep::EnergyType E0 = 10_GeV;
+ const HEPEnergyType E0 = 10_GeV;
particle.SetPID(particles::Code::Proton);
- hep::MomentumType P0 = sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
+ HEPMomentumType P0 = sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
particle.SetEnergy(E0);
particle.SetMomentum(plab);
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 39b59505eaba16853a06807d42ed91a56dcb9750..184eae5b543b9779696ccbd5b45036ed95edcfe4 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -41,10 +41,10 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
if (!fReport) return process::EProcessReturn::eOk;
[[maybe_unused]] int i = 0;
- EnergyType Etot = 0_GeV;
+ HEPEnergyType Etot = 0_GeV;
for (auto& iterP : s) {
- EnergyType E = iterP.GetEnergy();
+ HEPEnergyType E = iterP.GetEnergy();
Etot += E;
geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem(); // for printout
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index 285f275c0ae53ab2b6b128be270860819326cac6..845d3268b619b919e2197298c13ea625bf8113d8 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -33,14 +33,14 @@ using namespace corsika::geometry;
using namespace std;
using namespace corsika::units::si;
-typedef corsika::units::hep::energy_hep_d MOMENTUM;
+typedef corsika::units::si::hepmomentum_d MOMENTUM;
struct DummyParticle {
- EnergyType fEnergy;
+ HEPEnergyType fEnergy;
Vector<MOMENTUM> fMomentum;
Point fPosition;
- DummyParticle(EnergyType pEnergy, Vector<MOMENTUM> pMomentum, Point pPosition)
+ DummyParticle(HEPEnergyType pEnergy, Vector<MOMENTUM> pMomentum, Point pPosition)
: fEnergy(pEnergy)
, fMomentum(pMomentum)
, fPosition(pPosition) {}
@@ -89,7 +89,7 @@ TEST_CASE("TrackingLine") {
//~ std::cout << env.GetUniverse().get() << std::endl;
- DummyParticle p(1_J, Vector<MOMENTUM>(cs, 0_GeV, 0_GeV, 1_GeV),
+ DummyParticle p(1_GeV, Vector<MOMENTUM>(cs, 0_GeV, 0_GeV, 1_GeV),
Point(cs, 0_m, 0_m, 0_m));
auto const radius = 20_m;
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 481f8d48dc1abab991086e63d1f3b0962313844f..044ae32a1ae5adcabe5adaf9b5bc5b149d85aca1 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -29,7 +29,7 @@ namespace corsika::stack {
namespace super_stupid {
- typedef corsika::geometry::Vector<corsika::units::hep::energy_hep_d> MomentumVector;
+ typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
/**
* Example of a particle object on the stack.
@@ -46,7 +46,7 @@ namespace corsika::stack {
GetStackData().SetPID(GetIndex(), id);
}
- void SetEnergy(const corsika::units::hep::EnergyType& e) {
+ void SetEnergy(const corsika::units::si::HEPEnergyType& e) {
GetStackData().SetEnergy(GetIndex(), e);
}
@@ -66,7 +66,7 @@ namespace corsika::stack {
return GetStackData().GetPID(GetIndex());
}
- corsika::units::hep::EnergyType GetEnergy() const {
+ corsika::units::si::HEPEnergyType GetEnergy() const {
return GetStackData().GetEnergy(GetIndex());
}
@@ -82,9 +82,9 @@ namespace corsika::stack {
return GetStackData().GetTime(GetIndex());
}
- corsika::geometry::Vector<corsika::units::si::SpeedType::dimension_type>
- GetDirection() const {
- return GetMomentum() / GetEnergy() * corsika::units::constants::c;
+ corsika::geometry::Vector<corsika::units::si::dimensionless_d> GetDirection()
+ const {
+ return GetMomentum() / GetEnergy();
}
};
@@ -112,7 +112,7 @@ namespace corsika::stack {
void SetPID(const int i, const corsika::particles::Code id) { fDataPID[i] = id; }
- void SetEnergy(const int i, const corsika::units::hep::EnergyType e) {
+ void SetEnergy(const int i, const corsika::units::si::HEPEnergyType e) {
fDataE[i] = e;
}
void SetMomentum(const int i, const MomentumVector& v) { fMomentum[i] = v; }
@@ -125,7 +125,7 @@ namespace corsika::stack {
corsika::particles::Code GetPID(const int i) const { return fDataPID[i]; }
- corsika::units::hep::EnergyType GetEnergy(const int i) const { return fDataE[i]; }
+ corsika::units::si::HEPEnergyType GetEnergy(const int i) const { return fDataE[i]; }
MomentumVector GetMomentum(const int i) const { return fMomentum[i]; }
@@ -160,13 +160,14 @@ namespace corsika::stack {
using corsika::geometry::Point;
using corsika::particles::Code;
fDataPID.push_back(Code::Unknown);
- fDataE.push_back(0 * corsika::units::si::joule);
+ fDataE.push_back(0 * corsika::units::si::electronvolt);
//#TODO this here makes no sense: see issue #48
geometry::CoordinateSystem& dummyCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
fMomentum.push_back(MomentumVector(
- dummyCS, {0 * corsika::units::si::joule, 0 * corsika::units::si::joule,
- 0 * corsika::units::si::joule}));
+ dummyCS,
+ {0 * corsika::units::si::electronvolt, 0 * corsika::units::si::electronvolt,
+ 0 * corsika::units::si::electronvolt}));
fPosition.push_back(
Point(dummyCS, {0 * corsika::units::si::meter, 0 * corsika::units::si::meter,
0 * corsika::units::si::meter}));
@@ -187,7 +188,7 @@ namespace corsika::stack {
/// the actual memory to store particle data
std::vector<corsika::particles::Code> fDataPID;
- std::vector<corsika::units::hep::EnergyType> fDataE;
+ std::vector<corsika::units::si::HEPEnergyType> fDataE;
std::vector<MomentumVector> fMomentum;
std::vector<corsika::geometry::Point> fPosition;
std::vector<corsika::units::si::TimeType> fTime;
diff --git a/ThirdParty/phys/units/physical_constants.hpp b/ThirdParty/phys/units/physical_constants.hpp
index 8a5669d57908974c8e3d850892299ee6935ced87..131dfb03c1ee8497254f7d39f50bc3373136331d 100644
--- a/ThirdParty/phys/units/physical_constants.hpp
+++ b/ThirdParty/phys/units/physical_constants.hpp
@@ -24,35 +24,34 @@
#include "phys/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
diff --git a/ThirdParty/phys/units/quantity.hpp b/ThirdParty/phys/units/quantity.hpp
index 5f764409d6939ea43c393e254f2f62208c192f87..890061b950a0bb84903a271115eed9e89016b3e8 100644
--- a/ThirdParty/phys/units/quantity.hpp
+++ b/ThirdParty/phys/units/quantity.hpp
@@ -1,9 +1,8 @@
/**
* \file quantity.hpp
*
- * \brief Zero-overhead dimensional analysis and unit/quantity manipulation and conversion.
- * \author Michael S. Kenniston, Martin Moene
- * \date 7 September 2013
+ * \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.
@@ -31,36 +30,36 @@
#include <cmath>
#include <cstdlib>
-#include <utility> // std::declval
+#include <utility> // std::declval
/// namespace phys.
namespace phys {
-/// namespace units.
+ /// namespace units.
-namespace units {
+ namespace units {
#ifdef PHYS_UNITS_REP_TYPE
- using Rep = PHYS_UNITS_REP_TYPE;
+ using Rep = PHYS_UNITS_REP_TYPE;
#else
- using Rep = double;
+ using Rep = double;
#endif
-/*
- * declare now, define later.
- */
-template< typename Dims, typename T = Rep >
-class quantity;
+ /*
+ * 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 >
-struct dimensions
-{
- enum
- {
+ /**
+ * 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,
@@ -68,897 +67,835 @@ struct dimensions
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);
+ }
+ };
- is_all_zero =
- D1 == 0 && D2 == 0 && D3 == 0 && D4 == 0 && D5 == 0 && D6 == 0 && D7 == 0,
+ /// 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 Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+ };
+
+ 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 Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+ };
+
+ 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 Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+ };
+
+ 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 Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+ };
+
+ template <typename D, int N, typename T>
+ using Power = typename detail::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 Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+ };
+
+ 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
- is_base =
- 1 == (D1 != 0) + (D2 != 0) + (D3 != 0) + (D4 != 0) + (D5 != 0) + (D6 != 0) + (D7 != 0) &&
- 1 == D1 + D2 + D3 + D4 + D5 + D6 + D7,
- };
+ /**
+ * \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;
+ }
- template< int R1, int R2, int R3, int R4, int R5, int R6, int R7 >
- constexpr bool operator==( dimensions<R1, R2, R3, R4, R5, R6, R7> const & ) const
- {
- return D1==R1 && D2==R2 && D3==R3 && D4==R4 && D5==R5 && D6==R6 && D7==R7;
- }
+ /**
+ * the quantity's magnitude.
+ */
+ constexpr value_type magnitude() const { return m_value; }
- template< int R1, int R2, int R3, int R4, int R5, int R6, int R7 >
- constexpr bool operator!=( dimensions<R1, R2, R3, R4, R5, R6, R7> const & rhs ) const
- {
- return !( *this == rhs );
- }
-};
+ /**
+ * the quantity's dimensions.
+ */
+ constexpr dimension_type dimension() const { return dimension_type{}; }
-/// demensionless 'dimension'.
+ /**
+ * 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 ) };
-typedef dimensions< 0, 0, 0 > dimensionless_d;
+ private:
+ /**
+ * private initializing constructor.
+ */
+ constexpr explicit quantity(value_type x)
+ : m_value{x} {}
-/// namespace detail.
+ private:
+ value_type m_value;
-namespace detail {
+ enum { has_dimension = !Dims::is_all_zero };
-/**
- * \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;
-};
+ // static_assert( has_dimension, "quantity dimensions must not all be zero" ); //
+ // MR: removed
-template< typename T >
-struct collapse< dimensionless_d, T >
-{
- typedef T type;
-};
+ private:
+ // friends:
-template< typename D, typename T >
-using Collapse = typename collapse<D,T>::type;
+ // arithmetic
-// promote types of expression to result type.
+ template <typename D, typename X, typename Y>
+ friend constexpr quantity<D, X>& operator+=(quantity<D, X>& x,
+ quantity<D, Y> const& y);
-template < typename X, typename Y >
-using PromoteAdd = decltype( std::declval<X>() + std::declval<Y>() );
+ template <typename D, typename X>
+ friend constexpr quantity<D, X> operator+(quantity<D, X> const& x);
-template < typename X, typename Y >
-using PromoteMul = decltype( std::declval<X>() * std::declval<Y>() );
+ 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);
-/*
- * The following batch of structs are type generators to calculate
- * the correct type of the result of various operations.
- */
+ template <typename D, typename X, typename Y>
+ friend constexpr quantity<D, X>& operator-=(quantity<D, X>& x,
+ quantity<D, Y> const& y);
-/**
- * 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,
- };
+ template <typename D, typename X>
+ friend constexpr quantity<D, X> operator-(quantity<D, X> const& x);
- typedef Collapse< dimensions< d1, d2, d3, d4, d5, d6, d7 >, T > type;
-};
+ 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 DX, typename DY, typename X, typename Y>
-using Product = typename product<DX, DY, PromoteMul<X,Y>>::type;
+ template <typename D, typename X, typename Y>
+ friend constexpr quantity<D, X>& operator*=(quantity<D, X>& x, const Y& y);
-/**
- * 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,
- };
+ template <typename D, typename X, typename Y>
+ friend constexpr quantity<D, detail::PromoteMul<X, Y>> operator*(
+ quantity<D, X> const& x, const Y& y);
- typedef Collapse< dimensions< d1, d2, d3, d4, d5, d6, d7 >, T > type;
-};
+ 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>
-using Quotient = typename quotient<DX, DY, PromoteMul<X,Y>>::type;
+ 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);
-/**
- * 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,
- };
+ template <typename D, typename X, typename Y>
+ friend constexpr quantity<D, X>& operator/=(quantity<D, X>& x, const Y& y);
- typedef Collapse< dimensions< d1, d2, d3, d4, d5, d6, d7 >, T > type;
-};
+ 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>
-using Reciprocal = typename reciprocal<D, PromoteMul<X,Y>>::type;
+ template <typename D, typename X, typename Y>
+ friend constexpr detail::Reciprocal<D, X, Y> operator/(const X& x,
+ quantity<D, Y> const& y);
-/**
- * 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,
- };
+ 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);
- typedef Collapse< dimensions< d1, d2, d3, d4, d5, d6, d7 >, T > type;
-};
+ // absolute value.
-template< typename D, int N, typename T >
-using Power = typename detail::power< D, N, T >::type;
+ template <typename D, typename X>
+ friend constexpr quantity<D, X> abs(quantity<D, X> const& x);
-/**
- * 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
- };
+ // powers and roots
- 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
- };
+ template <int N, typename D, typename X>
+ friend detail::Power<D, N, X> nth_power(quantity<D, X> const& x);
- typedef Collapse< dimensions< d1, d2, d3, d4, d5, d6, d7 >, T > type;
-};
+ template <typename D, typename X>
+ friend constexpr detail::Power<D, 2, X> square(quantity<D, X> const& x);
-template< typename D, int N, typename T >
-using Root = typename detail::root< D, N, T >::type;
+ template <typename D, typename X>
+ friend constexpr detail::Power<D, 3, X> cube(quantity<D, X> const& x);
-/**
- * tag to construct a quantity from a magnitude.
- */
-constexpr struct magnitude_tag_t{} magnitude_tag{};
+ template <int N, typename D, typename X>
+ friend detail::Root<D, N, X> nth_root(quantity<D, X> const& x);
-} // namespace detail
+ template <typename D, typename X>
+ friend detail::Root<D, 2, X> sqrt(quantity<D, X> const& x);
-/**
- * \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;
+ // comparison
- typedef T value_type;
+ template <typename D, typename X, typename Y>
+ friend constexpr bool operator==(quantity<D, X> const& x, quantity<D, Y> const& y);
- typedef quantity<Dims, T> this_type;
+ template <typename D, typename X, typename Y>
+ friend constexpr bool operator!=(quantity<D, X> const& x, quantity<D, Y> const& y);
- constexpr quantity() : m_value{} { }
+ template <typename D, typename X, typename Y>
+ friend constexpr bool operator<(quantity<D, X> const& x, quantity<D, Y> const& y);
- /**
- * 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 ) { }
+ template <typename D, typename X, typename Y>
+ friend constexpr bool operator<=(quantity<D, X> const& x, quantity<D, Y> const& y);
- /**
- * converting copy-assignment constructor.
- */
- template <typename X >
- constexpr quantity( quantity<Dims, X> const & x )
- : m_value( x.magnitude() ) { }
+ template <typename D, typename X, typename Y>
+ friend constexpr bool operator>(quantity<D, X> const& x, quantity<D, Y> const& y);
-// /**
-// * convert to compatible unit, for example: (3._dm).to(meter) gives 0.3;
-// */
-// constexpr value_type to( quantity const & x ) const { return *this / x; }
+ template <typename D, typename X, typename Y>
+ friend constexpr bool operator>=(quantity<D, X> const& x, quantity<D, Y> const& y);
+ };
- /**
- * 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;
- }
+ // Give names to the seven fundamental dimensions of physical reality.
- /**
- * the quantity's magnitude.
- */
- constexpr value_type magnitude() const { return m_value; }
+ 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;
- /**
- * the quantity's dimensions.
- */
- constexpr dimension_type dimension() const { return dimension_type{}; }
+ // Addition operators
- /**
- * 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 ) };
+ /// quan += quan
-private:
- /**
- * private initializing constructor.
- */
- constexpr explicit quantity( value_type x ) : m_value{ x } { }
+ 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;
+ }
-private:
- value_type m_value;
+ /// + quan
+
+ template <typename D, typename X>
+ constexpr quantity<D, X> operator+(quantity<D, X> const& x) {
+ return quantity<D, X>(+x.m_value);
+ }
- enum { has_dimension = ! Dims::is_all_zero };
+ /// quan + quan
- // static_assert( has_dimension, "quantity dimensions must not all be zero" ); // MR: removed
+ 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);
+ }
-private:
- // friends:
+ // Subtraction operators
- // arithmetic
+ /// quan -= quan
template <typename D, typename X, typename Y>
- friend constexpr quantity<D, X> &
- operator+=( quantity<D, X> & x, quantity<D, Y> const & 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>
- friend constexpr quantity<D, X>
- operator+( quantity<D, X> const & x );
+ constexpr quantity<D, X> operator-(quantity<D, X> const& x) {
+ return quantity<D, X>(-x.m_value);
+ }
- 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 );
+ /// quan - quan
template <typename D, typename X, typename Y>
- friend constexpr quantity<D, X> &
- operator-=( quantity<D, X> & x, quantity<D, Y> const & 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);
+ }
- template <typename D, typename X>
- friend constexpr quantity<D, X>
- operator-( quantity<D, X> const & x );
+ // Multiplication operators
+
+ /// quan *= num
- 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>
+ constexpr quantity<D, X>& operator*=(quantity<D, X>& x, const Y& y) {
+ return x.m_value *= y, x;
+ }
- template< typename D, typename X, typename Y>
- friend constexpr quantity<D, X> &
- operator*=( quantity<D, X> & x, const Y & y );
+ /// quan * num
template <typename D, typename X, typename Y>
- friend constexpr quantity<D, detail::PromoteMul<X,Y>>
- operator*( quantity<D, X> const & x, const Y & 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>
- friend constexpr quantity< D, detail::PromoteMul<X,Y> >
- operator*( const X & x, quantity<D, Y> const & 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>
- friend constexpr detail::Product<DX, DY, X, Y>
- operator*( quantity<DX, X> const & lhs, quantity< DY, Y > const & rhs );
+ 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
- template< typename D, typename X, typename Y>
- friend constexpr quantity<D, X> &
- operator/=( quantity<D, X> & x, const Y & y );
+ /// quan /= num
template <typename D, typename X, typename Y>
- friend constexpr quantity<D, detail::PromoteMul<X,Y>>
- operator/( quantity<D, X> const & x, const Y & 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>
- friend constexpr detail::Reciprocal<D, X, Y>
- operator/( const X & x, quantity<D, Y> const & 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);
+ }
- 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 );
+ /// num / quan
- // absolute value.
+ 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);
+ }
- template <typename D, typename X>
- friend constexpr quantity<D,X>
- abs( quantity<D,X> const & x );
+ /// quan / quan:
- // powers and roots
+ 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);
+ }
- template <int N, typename D, typename X>
- friend detail::Power<D, N, X>
- nth_power( quantity<D, X> const & x );
+ /// absolute value.
template <typename D, typename X>
- friend constexpr detail::Power<D, 2, X>
- square( quantity<D, X> const & x );
+ constexpr quantity<D, X> abs(quantity<D, X> const& x) {
+ return quantity<D, X>(std::abs(x.m_value));
+ }
- template <typename D, typename X>
- friend constexpr detail::Power<D, 3, X>
- cube( quantity<D, X> const & x );
+ // General powers
- template <int N, typename D, typename X>
- friend detail::Root<D, N, X>
- nth_root( quantity<D, X> const & x );
+ /// N-th power.
- template <typename D, typename X>
- friend detail::Root< D, 2, X >
- sqrt( quantity<D, X> const & x );
+ template <int N, typename D, typename X>
+ detail::Power<D, N, X> nth_power(quantity<D, X> const& x) {
+ return detail::Power<D, N, X>(std::pow(x.m_value, X(N)));
+ }
- // comparison
+ // Low powers defined separately for efficiency.
- template <typename D, typename X, typename Y>
- friend constexpr bool operator==( quantity<D, X> const & x, quantity<D, Y> const & y );
+ /// square.
- 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>
+ constexpr detail::Power<D, 2, X> square(quantity<D, X> const& x) {
+ return x * x;
+ }
- template <typename D, typename X, typename Y>
- friend constexpr bool operator<( quantity<D, X> const & x, quantity<D, Y> const & y );
+ /// cube.
- 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>
+ constexpr detail::Power<D, 3, X> cube(quantity<D, X> const& x) {
+ return x * x * x;
+ }
- template <typename D, typename X, typename Y>
- friend constexpr bool operator>( quantity<D, X> const & x, quantity<D, Y> const & y );
+ // General root
- template <typename D, typename X, typename Y>
- friend constexpr bool operator>=( quantity<D, X> const & x, quantity<D, Y> const & y );
-};
+ /// n-th root.
-// Give names to the seven fundamental dimensions of physical reality.
-
-typedef dimensions< 1, 0, 0, 0, 0, 0, 0 > length_d;
-typedef dimensions< 0, 1, 0, 0, 0, 0, 0 > mass_d;
-typedef dimensions< 0, 0, 1, 0, 0, 0, 0 > time_interval_d;
-typedef dimensions< 0, 0, 0, 1, 0, 0, 0 > electric_current_d;
-typedef dimensions< 0, 0, 0, 0, 1, 0, 0 > thermodynamic_temperature_d;
-typedef dimensions< 0, 0, 0, 0, 0, 1, 0 > amount_of_substance_d;
-typedef dimensions< 0, 0, 0, 0, 0, 0, 1 > luminous_intensity_d;
+ template <int N, typename D, typename X>
+ detail::Root<D, N, X> nth_root(quantity<D, X> const& x) {
+ static_assert(detail::root<D, N, X>::all_even_multiples,
+ "root result dimensions must be integral");
-// Addition operators
+ return detail::Root<D, N, X>(std::pow(x.m_value, X(1.0) / N));
+ }
-/// 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 );
-}
+ // Low roots defined separately for convenience.
-// 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 );
-}
+ /// square root.
-/// 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 D, typename X>
+ detail::Root<D, 2, X> sqrt(quantity<D, X> const& x) {
+ static_assert(detail::root<D, 2, X>::all_even_multiples,
+ "root result dimensions must be integral");
-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
+ return detail::Root<D, 2, X>(std::pow(x.m_value, X(1.0) / 2));
+ }
-/// quan /= num
+ // Comparison operators
-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;
-}
+ /// equality.
-/// quan / num
+ 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;
+ }
-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 );
-}
+ /// inequality.
-/// num / quan
+ 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;
+ }
-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 );
-}
+ /// less-than.
-/// quan / quan:
+ 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;
+ }
-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>
-constexpr quantity<D,X> abs( quantity<D,X> const & x )
-{
- return quantity<D,X>( std::abs( x.m_value ) );
-}
+ /// less-equal.
-// General powers
-
-/// N-th power.
-
-template <int N, typename D, typename X>
-detail::Power<D, N, X>
-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, typename Y>
+ constexpr bool operator<=(quantity<D, X> const& x, quantity<D, Y> const& y) {
+ return x.m_value <= y.m_value;
+ }
-template <typename D, typename X>
-constexpr detail::Power<D, 2, X>
-square( quantity<D, X> const & x )
-{
- return x * x;
-}
+ /// greater-than.
-/// cube.
+ 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;
+ }
-template <typename D, typename X>
-constexpr detail::Power<D, 3, X>
-cube( quantity<D, X> const & x )
-{
- return x * x * x;
-}
+ /// greater-equal.
-// General root
+ 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;
+ }
-/// n-th root.
+ /// quantity's dimension.
-template <int N, typename D, typename X>
-detail::Root<D, N, X>
-nth_root( quantity<D, X> const & x )
-{
- static_assert( detail::root<D, N, X>::all_even_multiples, "root result dimensions must be integral" );
+ template <typename DX, typename X>
+ inline constexpr DX dimension(quantity<DX, X> const& q) {
+ return q.dimension();
+ }
- return detail::Root<D, N, X>( std::pow( x.m_value, X( 1.0 ) / N ) );
-}
+ /// quantity's magnitude.
-// Low roots defined separately for convenience.
-
-/// square root.
-
-template <typename D, typename X>
-detail::Root< D, 2, X >
-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::pow( x.m_value, X( 1.0 ) / 2 ) );
-}
-
-// 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 };
-
-// The standard SI prefixes.
-
-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;
-
-// 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 ); \
+ template <typename DX, typename X>
+ inline constexpr X magnitude(quantity<DX, X> const& q) {
+ return q.magnitude();
}
-#define QUANTITY_DEFINE_SCALING_LITERALS( pfx, dim, fact ) \
- 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 )
-
-
-#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 phys::units
+ // 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 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;
+
+ // 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(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)
+
+#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
diff --git a/ThirdParty/phys/units/quantity_io.hpp b/ThirdParty/phys/units/quantity_io.hpp
index 75df1ba64674bfdd407621c86227457419c9cddc..e94ffab36c1d0c84556cd326353aa1c6eb5d6881 100644
--- a/ThirdParty/phys/units/quantity_io.hpp
+++ b/ThirdParty/phys/units/quantity_io.hpp
@@ -131,6 +131,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();
}