auto integral = []() -> std::function<double (double)> {
    struct Env {
      double a, x, y; 
    }; 
    auto env = std::make_shared<Env>();

    return [
      env,
      integrate = epa::default_integrator(0),
      fx = [
        env,
        integrate = epa::default_integrator(1),
        fy = [
          env,
          integrate = epa::default_integrator(2),
          fz = [env](double z) -> double {
            return z / sqrt(sqr(env->x) + sqr(env->y) + z*z);
          }  
        ](double y) -> double {
          env->y = y; 
          return y * integrate(fz, 0, sqrt(1 - sqr(env->x / env->a)) - y*y);
        }  
      ](double x) -> double {
        env->x = x; 
        return x * integrate(fy, 0, sqrt(1 - sqr(x / env->a)));
      }  
    ](double a) -> double {
      return (15.0 / a) * integrate(fx, 0, a);
    }; 
  };