%% Chapter 8. Multiple Integration
% This is the script for the integration chapter.

syms x y z r theta a b real
syms f1(x) f2(x) g(x, y)
int(int(g, y, f1, f2), x, a, b)

%

f = @(x) sin(x + 1); g = @(x) x.^3 - 3*x + 1;
fplot(f, [-3, 3]), hold on
fplot(g, [-3, 3], 'Linewidth', 2)

%

x1 = fzero(@(x) f(x) - g(x), -2)
x2 = fzero(@(x) f(x) - g(x), 0)
x3 = fzero(@(x) f(x) - g(x), 2)

%

xcoord = linspace(x1, x3, 10);
ycoord = [f(xcoord); g(xcoord)];
plot([xcoord;xcoord], ycoord), hold off

%

format long
leftarea = integral(@(x) g(x) - f(x), x1, x2)
rightarea = integral(@(x) f(x) - g(x), x2, x3)
area = leftarea + rightarea
integral(@(x) abs(f(x) - g(x)), x1, x3)

% Note the discrepancy starting in the seventh decimal place.
%
% Computation with symbolic double integration

leftarea = double(int(int(1, sin(x+1), x^3-3*x+1), x1, x2))
format short

%
% Plotting with verticalRegion

verticalRegion(f, g, x1, x3)

%

int(int(1, x, sin(y), cos(y)), y, -3*pi/4, pi/4)

%
% Plotting with horizontalRegion

horizontalRegion(@sin, @cos, -3*pi/4, pi/4)

%
%% Polar Regions

fimplicit3(x^2+y^2+z^2-4, ...
    [-2.1, 2.1, -2.1, 2.1, -2.1, 2.1], 'MeshDensity', 15), hold on
fimplicit3((x-1)^2+y^2-1, ...
    [-2.1, 2.1, -2.1, 2.1, -2.1, 2.1], 'MeshDensity', 15), hold off
view([3,3,3]), axis equal, axis off

%

circ1 = @(t) max(2*cos(t), zeros(size(t)));
circ2 = @(t) 2*ones(size(t));
polarRegion(circ1, circ2, 0, 2*pi)

%

int(int(r, r, 0, 2), theta, pi/2, 3*pi/2) + ...
    int(int(r, r, 2*cos(theta), 2), theta, -pi/2, pi/2)

%

vol = int(int(4*sqrt(4-r^2)*r, r, 0, 2), theta, pi/2, pi) + ...
    int(int(4*sqrt(4-r^2)*r, r, 2*cos(theta), 2), theta, 0, pi/2)

%
%% Numerical Integration

format long, integral2(@(x,t) 2*(1-x.^2).^(1/2), -1, 1, 0, 1)
ans - pi

%

integral2(@(x,t) 2*(1-(x.^2).^(1/3)).^(3/2), -1, 1, 0, 1)
ans - 3*pi/8

%
%% Solid Regions

viewSolid('x', 1, 2, 'y', @(x) x/2, @(x) x, 'z', @(x,y) zeros(size(x)), @(x,y) (x+y)/4)

%
%% Region between a parabolic cylinder and a paraboloid

cyl = 4 - y^2; par = x^2 + 3*y^2;
xycurve = simplify(cyl == par)
fimplicit(xycurve, [-2, 2, -1, 1])
axis equal

%

ysols = solve(xycurve, y)
ylim1 = matlabFunction(ysols(1)); 
ylim2 = matlabFunction(ysols(2)); 
viewSolid('x', -2, 2, 'y', ylim1, ylim2, 'z', @(x, y) ...
  x.^2 + 3*y.^2, @(x, y) 4 - y.^2)

%

firststep = int(int(1, z, par, cyl), y, ysols(1), ysols(2));
int(simplify(firststep), x, -2, 2)

%
%% A torus intersected with a cylinder

syms u v t r x y z real
torus = [(4 + cos(u))*cos(v), (4 + cos(u))*sin(v), sin(u)];
toreqn = simplify(expand(subs((r-4)^2 + z^2 == 1, r, sqrt(x^2 + y^2))))
cylinder = [u, cos(v), sin(v)];
cyleqn = (y^2 + z^2 == 1);

%

fmesh(torus(1), torus(2), torus(3), [0, 2*pi]), hold on
fmesh(cylinder(1), cylinder(2), cylinder(3), [-5.5, 5.5, 0, 2*pi])
hold off, axis equal

%
%% Cross-sections with constant x

for j=0:5
   subplot(3,2,j+1)
   fimplicit(cyleqn, 'r', [-1, 1, -1, 1]), hold on
   fimplicit(subs(toreqn, x, 3+2*j/5), 'b', [-1, 1, -1, 1]), hold off
   axis equal
end

%
%% Cross-sections with constant y

for j=0:7
   subplot(4,2,j+1)
   zval = sqrt(1-(-1+2*j/7)^2);
   plot([-5 -5; 5 5], [zval -zval; zval -zval], 'r'), hold on
   fimplicit(subs(toreqn, y, -1+2*j/7), 'b', [-5, 5, -1, 1])
   hold off, axis equal
end

%
%% Sketching the solid

figure
xsol = solve(toreqn, x)
subs(xsol, [y, z], [0, 0])
xlo = matlabFunction(xsol(1)); xhi = matlabFunction(xsol(2));
zlo = matlabFunction(-sqrt(1-y^2)); zhi = matlabFunction(sqrt(1-y^2));
viewSolid('y', -1, 1, 'z', zlo, zhi, 'x', xlo, xhi)

%
%% Computing the volume

vol = 2*integral3(@(x,y,z) ones(size(x)), -1, 1, zlo, zhi, xlo, xhi)

%
%% Cylindrical Coordinates

newtoreqn = ((r-4)^2 + z^2 == 1);
newcyleqn = subs(cyleqn, y, r*sin(theta))
8*integral3(@(r, theta, z) r, 3, 5, @(r) zeros(size(r)), ...
   @(r) asin(1./r), @(r, theta) zeros(size(r)), @(r, theta) ...
   sqrt(min(1-(r-4).^2, 1-r.^2.*sin(theta).^2)))

%
%% More General Coordinate Changes
% This shows the two hyperbolas and two parabolas.

figure
fcontour(x^2-y^2, [0, 5, 0, 5], 'LevelList', [1 4], 'LineWidth', 2), hold on
fcontour(x^2+y, [0, 5, 0, 5], 'LevelList', [9 16], 'LineWidth', 2), hold off

%

f = x^2 - y^2; g = x^2 + y;
jake = det(jacobian([f, g], [x, y]))
backsoly = solve(g - f == v - u, y)
backsoly = backsoly(2); backsolx = sqrt(v - backsoly)
jacobianuv = subs(jake^(-1), [x, y], [backsolx, backsoly]);
area = integral2(matlabFunction(jacobianuv), 1, 4, 9, 16)