diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index c6f2dbe6111e7141f2880e320c83a4bdead4a030..7068bb26a0c88f42558a8e16f0f2b049c0a6adbd 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -15,9 +15,12 @@
 #include <corsika/cascade/Cascade.h>
 #include <corsika/environment/Environment.h>
 #include <corsika/environment/FlatExponential.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/IMagneticFieldModel.h>
 #include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
 #include <corsika/environment/NuclearComposition.h>
 #include <corsika/environment/ShowerAxis.h>
+#include <corsika/environment/UniformMagneticField.h>
 #include <corsika/geometry/Plane.h>
 #include <corsika/geometry/Sphere.h>
 #include <corsika/logging/Logging.h>
@@ -143,8 +146,8 @@ int main(int argc, char** argv) {
   cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
        << endl;
 
-  auto const observationHeight = 0_km + builder.getEarthRadius();
-  auto const injectionHeight = 112.75_km + builder.getEarthRadius();
+  auto const observationHeight = 1.4_km + seaLevel;
+  auto const injectionHeight = 112.75_km + seaLevel;
   auto const t = -observationHeight * cos(thetaRad) +
                  sqrt(-units::static_pow<2>(sin(thetaRad) * observationHeight) +
                       units::static_pow<2>(injectionHeight));
diff --git a/Setup/SetupEnvironment.h b/Setup/SetupEnvironment.h
index f4265e3ddf7342b2972e6317a14bbf946c664ff1..9c893beb2c79800601c739096f5fe56ffb0c12b5 100644
--- a/Setup/SetupEnvironment.h
+++ b/Setup/SetupEnvironment.h
@@ -13,6 +13,7 @@
 #include <corsika/environment/IMediumModel.h>
 #include <corsika/environment/IMediumPropertyModel.h>
 #include <corsika/environment/IRefractiveIndexModel.h>
+#include <corsika/environment/IMagneticFieldModel.h>
 
 namespace corsika::setup {