We cannot test without your data, but the approach seems fine. You could use the PDIST2 function and extract everything in one shot, but as a first approach this is fine.
I will assume that the question is about saving the points, so here is small update of your code that should (if I understand your example well) generate an Excel file with 4 columns that store Circle ID, Point ID, Point X, Point Y, Point Z:
% Prealloc array of counts, and cell array for storing relevant point parameters per circle.
counts = zeros( row, 1 ) ;
output = cell( row, 1 ) ;
% Initialize figure.
figure() ;
hold on ;
% Process circles.
for circleId = 1 : row
% Plot circle.
pos = [x0(circleId)-R, y0(circleId)-R, 2*R, 2*R] ;
rectangle( 'Position', pos, 'Curvature', [1 1] ) ;
% Compute distance from all points to current circle center.
distances = sqrt( (x - x0(circleId)).^2 + (y - y0(circleId)).^2 ) ;
% Flag and count relevant points.
isInsideCircle = distances <= R ;
counts(circleId) = sum( isInsideCircle ) ;
% Extract relevant parameters.
pointsId = find( isInsideCircle ) ;
pointsX = x(isInsideCircle) ; % Selected Lon
pointsY = y(isInsideCircle) ; % Selected Lat
pointsZ = z(isInsideCircle) ; % Selected h
% Store in relevant output cell. We build the first column as a vector of elements equal
% to circleId, whose size matches the size of the other columns (e.g. pointsId).
output{circleId} = [circleId * ones( size( pointsId )), pointsId, pointsX, pointsY, pointsZ] ;
end
% Build data to export to Excel by concatenating all cells of the output cell array vertically
% into a large numeric array, converting it into a cell array, and concatenating to a first
% "header" row.
output = num2cell( vertcat( output{:} )) ;
output = [{'CircleID', 'PointID', 'PointX', 'PointY', 'PointZ'}; output] ;
% Export to Excel.
xlswrite( 'PointsPerCircle.xlsx', output ) ;
Note that I didn't test it, but it illustrates a possible approach. Let me know if you have any question.
