Weighing the Earth

Overview

This example determines the mass of the Earth from two benchtop measurements, the way the quantity was first obtained after Cavendish measured G. A pendulum is timed to get the local gravitational acceleration, which is then combined with the Earth's radius and the gravitational constant:

\[g = \frac{4\pi^2 L}{T^2} \qquad M = \frac{g R^2}{G}\]

Every input carries a measurement uncertainty, G carries the one CODATA publishes for it, and the result reports what all of them together imply. The representation type doing the work is uncertain<T>, which propagates uncertainty through arithmetic automatically.

Measuring the local gravity

A pendulum is timed over many swings so that the reaction-time error is divided by their count. Solving \(T = 2\pi\sqrt{L/g}\) for the acceleration gives the first stage:

weighing_the_earth.cpp

// A pendulum measures the magnitude of the free-fall acceleration, which is a scalar, while
// `isq::acceleration` is a vector quantity with no scalar counterpart yet (V3 will model the
// magnitude of a vector).
QUANTITY_SPEC(gravity_magnitude, isq::speed / isq::duration);

[[nodiscard]] QuantityOf<gravity_magnitude> auto local_gravity(const QuantityOf<isq::distance> auto& pendulum_length,
                                                               const QuantityOf<isq::duration> auto& total_time,
                                                               int swings)
{
  const QuantityOf<isq::duration> auto period = total_time / swings;
  return gravity_magnitude(1 * (mag<4> * pow<2>(π)) * pendulum_length / pow<2>(period));
}

The π here is a symbolic constant, so it is exact and contributes nothing to the uncertainty of the result. Only the two measurements do.

Note the quantity being defined. isq::acceleration is a vector quantity in ISO 80000-3, just as isq::velocity is, while a pendulum measures the magnitude of the local free-fall acceleration, which is a scalar. ISQ names no scalar counterpart to acceleration the way speed is the one for velocity, and the magnitude of a vector quantity is not modelled yet (it is planned for V3), so the example defines what it means. Being honest about the character is not pedantry here: it is what lets the second stage reduce to isq::mass on its own. Had g been left a vector, the framework would have refused to convert the result to a scalar quantity, and rightly so.

Weighing the planet

The second stage asks what mass a sphere must have to produce the observed surface gravity. G is a measured constant, so it is materialized with measurement_of rather than used as a bare unit, and its CODATA uncertainty joins the propagation:

weighing_the_earth.cpp

[[nodiscard]] QuantityOf<isq::mass> auto sphere_mass_from_surface_gravity(
  const QuantityOf<gravity_magnitude> auto& surface_gravity, const QuantityOf<isq::radius> auto& radius)
{
  const quantity G = measurement_of(iau::newtonian_constant_of_gravitation);
  return surface_gravity * pow<2>(radius) / G;
}

Running the experiment

First the measured inputs, each with the uncertainty its instrument justifies:

weighing_the_earth.cpp

  // A 1 m pendulum timed over 50 swings with a stopwatch, and the Earth's mean radius as
  // published. The uncertainties are what the instruments justify: half a millimetre on the
  // metre stick, six hundredths of a second on the total time, a hundred metres on the radius.
  const quantity pendulum_length = isq::distance(uncertain{1.0000, 0.0005} * m);
  const quantity total_time = uncertain{100.30, 0.06} * s;
  constexpr int swings = 50;
  const quantity earth_radius = isq::radius(uncertain{6371.0, 0.1} * km);

The two lengths are given strong quantity types, and that is typing which buys something. Both are lengths, so nothing dimensional would stop them being passed in each other's place, but isq::distance and isq::radius sit in different branches of the ISQ hierarchy (under path_length and width respectively), so each is rejected where the other is expected. The duration needs no such help, since nothing else in the call is a duration. A caller who does not want typed quantities is not forced into them either: a plain 6371.0 * km carries kind_of<isq::length>, which converts down and satisfies either constraint.

Then the two stages run and the results are printed against their published counterparts:

weighing_the_earth.cpp

  const quantity gravity = local_gravity(pendulum_length, total_time, swings);
  const quantity earth_mass = sphere_mass_from_surface_gravity(gravity, earth_radius);

  std::cout << MP_UNITS_STD_FMT::format("pendulum length     = {}\n", pendulum_length);
  std::cout << MP_UNITS_STD_FMT::format("time for {} swings  = {}\n", swings, total_time);
  std::cout << MP_UNITS_STD_FMT::format("local gravity       = {::N[.4f]}\n", gravity.in(m / s2));
  // a conventional value fixed by the CGPM, not a measurement, so it carries no uncertainty
  std::cout << MP_UNITS_STD_FMT::format("  standard gravity  = {::N[.4f]}\n", (1. * si::standard_gravity).in(m / s2));
  std::cout << MP_UNITS_STD_FMT::format("mass of the Earth   = {::N[.4e]}\n", earth_mass.in(kg));
  // the published mass is itself defined as (GM)⊕ᴺ/G, so requesting an uncertainty-capable
  // representation for the conversion reports how well `G` is known
  std::cout << MP_UNITS_STD_FMT::format("  accepted value    = {::N[.4e]}\n",
                                        (1. * iau::terrestrial_mass).in<uncertain<double>>(kg));

Finally the uncertainty budget, which is where the interesting part is. Each input contributes its own relative uncertainty scaled by the exponent it enters the formula with, and the contributions are combined by hand so that the total can be compared against the one the library propagated:

weighing_the_earth.cpp

  // each contribution is scaled by the exponent its input enters the formula with
  const auto relative_uncertainty_of = [](const auto& quantity) {
    return quantity.numerical_value_in(quantity.unit).relative_uncertainty();
  };
  const auto from_length = relative_uncertainty_of(pendulum_length);
  const auto from_period = 2 * relative_uncertainty_of(total_time);
  const auto from_radius = 2 * relative_uncertainty_of(earth_radius);
  const auto from_G = get_value<double>(get_relative_standard_uncertainty(iau::newtonian_constant_of_gravitation));
  // independent contributions, so they combine in quadrature rather than by addition
  const auto combined = std::hypot(std::hypot(from_length, from_period), std::hypot(from_radius, from_G));

  std::cout << MP_UNITS_STD_FMT::format("\nrelative uncertainty contributions\n");
  std::cout << MP_UNITS_STD_FMT::format("  pendulum length        = {:.1e}\n", from_length);
  std::cout << MP_UNITS_STD_FMT::format("  period (squared)       = {:.1e}\n", from_period);
  std::cout << MP_UNITS_STD_FMT::format("  Earth radius (sq.)     = {:.1e}\n", from_radius);
  std::cout << MP_UNITS_STD_FMT::format("  G (CODATA 2018)        = {:.1e}\n", from_G);
  std::cout << MP_UNITS_STD_FMT::format("  combined in quadrature = {:.3e}\n", combined);
  std::cout << MP_UNITS_STD_FMT::format("  reported by the result = {:.3e}\n", relative_uncertainty_of(earth_mass));

pendulum length     = 1 ± 5e-04 m
time for 50 swings  = 100.3 ± 0.06 s
local gravity       = 9.8107 ± 0.0127 m/s²
  standard gravity  = 9.8066 m/s²
mass of the Earth   = 5.9663e+24 ± 7.7398e+21 kg
  accepted value    = 5.9722e+24 ± 1.3139e+20 kg

relative uncertainty contributions
  pendulum length        = 5.0e-04
  period (squared)       = 1.2e-03
  Earth radius (sq.)     = 3.1e-05
  G (CODATA 2018)        = 2.2e-05
  combined in quadrature = 1.297e-03
  reported by the result = 1.297e-03

What the numbers say

The two values agree. Everything needed is in the two mass lines printed above:

From the output	Value
this experiment	5.9663 × 10²⁴ ± 7.7398 × 10²¹ kg
published (iau::terrestrial_mass)	5.9722 × 10²⁴ ± 1.3139 × 10²⁰ kg

The gap between the two central values is 0.099% of the published mass:

\[\Delta = 5.9722 \times 10^{24}\,\mathrm{kg} - 5.9663 \times 10^{24}\,\mathrm{kg} = 5.9 \times 10^{21}\,\mathrm{kg}\]

Because both figures carry an uncertainty, the comparison uses their combination rather than either alone:

\[u_\text{combined} = \sqrt{\left(7.7398 \times 10^{21}\right)^2 + \left(1.3139 \times 10^{20}\right)^2} = 7.74 \times 10^{21}\,\mathrm{kg}\]

The gap is therefore well inside one combined standard uncertainty:

\[\frac{\Delta}{u_\text{combined}} = 0.76\]

so the two results are consistent. The published uncertainty is about sixty times smaller than the experiment's and barely moves that number, but it is included rather than assumed away. This is the question an uncertainty exists to answer: not "what number did I get" but "is the difference from the expected one larger than what my instruments can resolve".

Note where those two uncertainties come from. iau::terrestrial_mass is defined as (GM)⊕ᴺ/G, so its value in kilograms is known only as well as G is, and the example opts into reporting that with .in<uncertain<double>>(kg). The standard free-fall acceleration, by contrast, prints without a ±, because si::standard_gravity is an exact conventional value fixed by the CGPM rather than a measured one. It is a reference point, not the true local value: real free-fall acceleration ranges from about 9.78 m/s² at the equator to 9.83 m/s² at the poles, a spread far wider than this experiment's uncertainty.

The budget reconciles with the result. The contributions are independent, so they combine in quadrature rather than by addition: adding the four numbers would give 1.75 × 10⁻³, while combining them properly gives 1.297 × 10⁻³, exactly what the propagation through the two functions reported. The example prints both so the agreement can be checked rather than taken on trust.

The stopwatch dominates the budget. The period enters squared, so its relative uncertainty counts double, which puts it above the metre stick even though both readings look similarly careful. That is the practical lesson of an uncertainty budget: improving the radius figure or waiting for a better value of G would change nothing here, while timing more swings would.

G is negligible, and that is a finding. The least precisely known constant in physics contributes 2.2 × 10⁻⁵, some sixty times below the total. Being able to see that, rather than assume it, is why the constant carries its published uncertainty in its definition.

The agreement is partly luck, and the numbers show it

Everything above is random uncertainty. The model carries a systematic error too: treating the Earth as a uniform sphere with its mean radius ignores oblateness and rotation. That bias can be measured. Putting the exact standard gravity and the same mean radius through the same formula gives 5.9639 × 10²⁴ kg, which is 0.139% below the published mass.

Compare that with the experiment's own precision. The printed reported by the result = 1.297e-03 is 0.130% relative, so the modelling bias is larger than the random uncertainty. The result landed inside its error bar only because the measured g came out slightly above standard, which pulled the answer up to 0.099% low instead of 0.139%.

The consequence is worth sitting with: time more swings and the random uncertainty shrinks while the 0.139% bias does not, so a more precise experiment would start to disagree with the accepted value while looking more confident. No propagation formula can warn about this, which is why the assumptions and limitations of uncertain<T> are worth reading before trusting a narrow interval.