 	}
 }
 
+__global__ void MeanKernel(const float* __restrict__ d_Tpot, const float* __restrict__ d_MR,
+                           const float* __restrict__ d_prevTpot, const float* __restrict__ d_prevMR,
+                           float* __restrict__ d_meanTpot, float* __restrict__ d_meanMR, float* __restrict__ d_range,
+                           size_t N)
+{
+	const int idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+	if (idx < N)
+	{
+		if (IsValid(d_Tpot[idx]))
+		{
+			// trapezoidal integration
+			d_range[idx] += 1;
+			d_meanTpot[idx] += (d_prevTpot[idx] + d_Tpot[idx]) * 0.5;
+			d_meanMR[idx] += (d_prevMR[idx] + d_MR[idx]) * 0.5;
+		}
+	}
+}
+
+__global__ void MeanFinalizeKernel(float* __restrict__ d_meanTpot, float* __restrict__ d_meanMR,
+                                   const float* __restrict__ d_range, size_t N)
+{
+	const int idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+	if (idx < N)
+	{
+		d_meanTpot[idx] = d_meanTpot[idx] / d_range[idx];
+		d_meanMR[idx] = d_meanMR[idx] / d_range[idx];
+	}
+}
+
 cape_multi_source cape_cuda::GetNHighestThetaEValuesGPU(const std::shared_ptr<const plugin_configuration>& conf,
                                                         std::shared_ptr<info<float>> myTargetInfo, int n)
 {

	CUDA_CHECK(cudaFree(d_Psurf));

	CUDA_CHECK(cudaStreamDestroy(stream));

	return std::make_tuple(T, TD, VEC(Psurf));
}


			P[i] = Psample[i][k];
		}

		if (static_cast<unsigned int>(count_if(P.begin(), P.end(), [](const float& v) { return IsMissing(v); })) ==
		    P.size())
		{
