How do I calculate the flux through the surface I have generate inside the barrel?

I have this code and the only thing that is missing is the total flux across the surface/ closed shape I have generated. Calculate the flux using the Divergence Theorem/ Integrations. The surface could be any closed shape. Also if possible, find the value of surface area and volume of the shape using vector analysis. Thank You

R = 10;
[x,y,z] = sphere;
% Define rows to keep as those with Z in the middle
topCap = z > .5;
bottomCap = z < -.5; 
% Max out z at the cap values where x and y are in the cap.
z(topCap) = 0.7;
z(bottomCap) = -0.7;
% Display it.
surf(R*x, R*y, R*z);
%transparency of object
alpha(.1)
axis equal
axis tight; hold on
%generate 100 random points 
x1= -10 + (10+10)*rand(100,1)
y1= -10 + (10+10)*rand(100,1)
z1= -10 + (10+10)*rand(100,1)
%plot the points
plot3(x1,y1,z1,'.','MarkerSize',15)
d = sqrt(x1.^2 + y1.^2 + z1.^2);
[dmax, maxidx] = max(d);
[dmin, minidx] = min(d);
xmin = x1(minidx); ymin = y1(minidx); zmin = z1(minidx);
xmax = x1(maxidx); ymax = y1(maxidx); zmax = z1(maxidx);
fprintf('The closest point is (%g,%g,%g) which is distance %g\n', xmin, ymin, zmin, dmin);
fprintf('The furthest point is (%g,%g,%g) which is distance %g\n', xmax, ymax, zmax, dmax);
%3d polygon
npts = 75;
pts = randn(npts,3);
scatter3(x1,y1, z1,'.')
shp = alphaShape(pts);
h = plot(shp);
volume(shp)
h.FaceAlpha = .25;
totalsurfarea = surfaceArea(shp)

[Timothy Turk]

Did you ever figure this out? I have <300 points (Lat/Long/Alt) in orbit around the Earth and I need to determine at what point does the orbiting objects reach maximum density.