%% INTEGRATING F(X,Y) ON A REGION
%(Selcuk Bayin, METU-NCC, bayin@metu.edu.tr)
clc
disp(' Estimate the integral of xy over the upper-right quadrant of a circle')
f = inline( 'x.*y+2*x', 'x', 'y') % function to integrate
%%
% Use a rectangular mesh [a,b]x[c,d]
a = 0;
b = 1;            
c = 0;             
d = 1;
%Size of the mesh
M= 5;
N= 3;
dx = (b-a)/M;
dy = (d-c)/N;
meshProp=[a b c d N M dx dy]
%%
%Boundary of the region
%Take it as a circle with R=1
g = inline( 'x.^2 + y.^2', 'x', 'y')
%%
clc
fprintf('-------')
%Midpoints of the mesh
x_midpoints = a+dx/2:dx:b-dx/2
y_midpoints = c+dy/2:dy:d-dy/2
%%
[midpoints_x, midpoints_y] = meshgrid(x_midpoints, y_midpoints)
%%
clc
interior_midpoints_x = midpoints_x(g(midpoints_x, midpoints_y) <= 1);
interior_midpoints_y = midpoints_y(g(midpoints_x, midpoints_y) <= 1);
midpoints=[interior_midpoints_x,interior_midpoints_y]
%%
%Let us now evaluate the integral
midpoint_approx = sum(f(interior_midpoints_x, interior_midpoints_y))*dx*dy
