Skip to content

Commit 4b37d7a

Browse files
committed
Release MultipleIntegrals and Applications
1 parent 50c5f53 commit 4b37d7a

File tree

72 files changed

+1952
-12
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

72 files changed

+1952
-12
lines changed

FunctionLibrary/AddArrow.m

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@ function AddArrow(ax,pt1,pt2,opts)
2323
dim = length(pt1);
2424
% Ensure both points have the same length
2525
assert(length(pt1)==length(pt2))
26+
assert(dim < 4)
2627

2728
% Calculate the scale based on the axes limits
2829
xScale = range(ax.XLim);
@@ -70,7 +71,7 @@ function AddArrow(ax,pt1,pt2,opts)
7071
elseif dim == 2
7172
perpdir = [-slope(2) slope(1)];
7273
else
73-
error("Unexpected dimension: " + dim)
74+
error("AddArrow:IncorrectDimension","Unexpected dimension: " + dim)
7475
end
7576

7677
% Normalize the perpendicular direction

FunctionLibrary/CreatePumpkin.m

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
%[text] Original pumpkin code from Eric Ludlam, [Gourds to Graphics: The MATLAB Pumpkin](http://blogs.mathworks.com/graphics-and-apps/2023/10/31/gourds-to-graphics-the-matlab-pumpkin/blogs.mathworks.com/graphics-and-apps/2023/10/31/gourds-to-graphics-the-matlab-pumpkin/)
2+
function [Cavity,Punkin,StemData,Cpunkin] = CreatePumpkin()
3+
4+
% Draw seed cavity
5+
xrCavity = 0.5;
6+
yrCavity = 0.5;
7+
zrCavity = 0.35;
8+
9+
[Xcavity,Ycavity,Zcavity] = ellipsoid(0,0,0,xrCavity,yrCavity,zrCavity,100);
10+
11+
% Add irregularity
12+
Xcavity = Xcavity + rand(size(Xcavity))/60;
13+
Ycavity = Ycavity + rand(size(Ycavity))/60;
14+
Zcavity = Zcavity + rand(size(Zcavity))/100;
15+
16+
numPrimaryBumps = 10;
17+
totalBumps = numPrimaryBumps*2; % Add in-between creases.
18+
vertsPerBump = 12; % Originally 10
19+
numVerts = totalBumps*vertsPerBump+1;
20+
rPrimary = linspace(0,numPrimaryBumps*2,numVerts);
21+
rSecondary = linspace(0,totalBumps*2,numVerts);
22+
crease_depth = .04;
23+
crease_depth2 = .01;
24+
25+
Rxy_primary = 0 - (1-mod(rPrimary,2)).^2*crease_depth;
26+
Rxy_secondary = 0 - (1-mod(rSecondary,2)).^2*crease_depth2;
27+
Rxy = Rxy_primary + Rxy_secondary;
28+
29+
[Xsphere,Ysphere,Zsphere] = sphere(numVerts-1); % Sphere creates +1 verts
30+
Xpunkin = (1+Rxy).*Xsphere;
31+
Ypunkin = (1+Rxy).*Ysphere;
32+
33+
dimple = .2; % Fraction to dimple into top/bottom
34+
35+
rho = linspace(-1,1,numVerts)';
36+
Rz_dimple = (0-rho.^4)*dimple;
37+
38+
39+
HeightRatio = .8;
40+
Zpunkin = (1+Rxy).*Zsphere.*(HeightRatio+Rz_dimple);
41+
42+
Rstem = (1-(1-mod(rPrimary+1,2)).^2)*.05;
43+
44+
thetac = linspace(0,2,numVerts);
45+
Xcyl = cospi(thetac);
46+
Ycyl = sinpi(thetac);
47+
48+
Zcyl = linspace(0,1,11)'; % column vector
49+
Rstemz = .7+(1-Zcyl).^2*.6;
50+
51+
Xstem = (.1+Rstem).*Xcyl.*Rstemz;
52+
Ystem = (.1+Rstem).*Ycyl.*Rstemz;
53+
Zstem = repmat(Zcyl*.15,1,numVerts);
54+
55+
StemData.Xstem = Xstem;
56+
StemData.Ystem = Ystem;
57+
StemData.Zstem = Zstem;
58+
StemData.heightratio = HeightRatio;
59+
StemData.numVerts = numVerts;
60+
StemData.Zsphere = Zsphere;
61+
62+
Cpunkin = hypot(hypot(Xpunkin,Ypunkin),(1+Rxy).*Zsphere); % As if pumpkin were round with no dimples
63+
Punkin = struct("X",Xpunkin,"Y",Ypunkin,"Z",Zpunkin);
64+
Cavity = struct("X",Xcavity,"Y",Ycavity,"Z",Zcavity);
65+
DrawPumpkin(Cavity,Punkin,StemData,Cpunkin);
66+
end
67+
68+
%[appendix]{"version":"1.0"}
69+
%---
Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
function DrawCartesianVolumeElement(Copt,dx,dy,dz,Lx,Ly,Lz)
2+
arguments
3+
Copt = "summer"
4+
dx = 0.6
5+
dy = 0.8
6+
dz = 0.5
7+
Lx = 1
8+
Ly = 1
9+
Lz = 1
10+
end
11+
CData = colormap(Copt);
12+
% Define vertices of a rectangular prism centered at origin
13+
vx = Lx + [0 dx dx 0 0 dx dx 0];
14+
vy = Ly + [0 0 dy dy 0 0 dy dy];
15+
vz = Lz + [0 0 0 0 dz dz dz dz];
16+
% Faces (each row is indices into vertices)
17+
faces = [1 2 3 4; 5 6 7 8; 1 2 6 5; 2 3 7 6; 3 4 8 7; 4 1 5 8];
18+
patch(Vertices=[vx' vy' vz'],Faces=faces,FaceColor=CData(1,:),EdgeColor="none",FaceAlpha=0.8);
19+
20+
xlim([0 3])
21+
ylim([0 3])
22+
zlim([0 3])
23+
xlabel("x"), ylabel("y"), zlabel("z")
24+
title("Cartesian Volume Element")
25+
CurVol = dx*dy*dz;
26+
subtitle("V = " + CurVol)
27+
view(35,20)
28+
grid on
29+
end
30+
31+
%[appendix]{"version":"1.0"}
32+
%---
Lines changed: 82 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,82 @@
1+
function DrawCylindricalVolumeElement(Copt,dr,dtheta,dz,r0,theta0,z0)
2+
arguments
3+
Copt = "summer"
4+
dr = 0.6
5+
dtheta = pi/6
6+
dz = 0.5
7+
r0 = sqrt(2)
8+
theta0 = pi/4
9+
z0 = 1
10+
end
11+
12+
CData = colormap(Copt);
13+
% Create surface patch representing area element as curved rectangle on cylinder
14+
theta = linspace(theta0, theta0+dtheta, 30);
15+
z = [z0 z0+dz];
16+
[Th, Z] = meshgrid(theta, z);
17+
X = r0 * cos(Th);
18+
Y = r0 * sin(Th);
19+
surf(X, Y, Z,FaceColor="interp",EdgeColor="none",FaceAlpha=0.9)
20+
hold on
21+
% Show radial thickness by plotting inner and outer curved faces
22+
X2 = (r0+dr) * cos(Th);
23+
Y2 = (r0+dr) * sin(Th);
24+
surf(X2, Y2, Z,FaceColor="interp",EdgeColor="none",FaceAlpha=0.9)
25+
% Connect edges to form the small curved shell patch
26+
for k = [1 size(Th,2)]% 1:size(Th,2)
27+
patch([X(1,k) X2(1,k) X2(2,k) X(2,k)], [Y(1,k) Y2(1,k) Y2(2,k) Y(2,k)], [Z(1,k) Z(1,k) Z(2,k) Z(2,k)],...
28+
sqrt(([X(1,k) X2(1,k) X2(2,k) X(2,k)]-X(1,1)).^2+([Y(1,k) Y2(1,k) Y2(2,k) Y(2,k)]-Y(1,1)).^2),...
29+
EdgeColor="none",FaceAlpha=0.9)
30+
end
31+
32+
fill3([X(1,:) fliplr(X2(1,:))],[Y(1,:) fliplr(Y2(1,:))],z0*ones(1,2*size(Y,2)),z0*ones(1,2*size(Y,2)))
33+
fill3([X(end,:) fliplr(X2(end,:))],[Y(end,:) fliplr(Y2(end,:))],(z0+dz)*ones(1,2*size(Y,2)),(z0+dz)*ones(1,2*size(Y,2)))
34+
35+
% Plot the twelve edges of the volume element using SeriesIndex="none" to set color
36+
% Define the 8 corner points of the small curved shell (r0/r0+dr) x (theta0/theta0+dtheta) x (0/dz)
37+
thetas = [theta0, theta0+dtheta];
38+
rs = [r0, r0+dr];
39+
zs = [z0, z0+dz];
40+
Vertices = zeros(8,3); idx = 1;
41+
for ir = 1:2
42+
for it = 1:2
43+
for iz = 1:2
44+
th = thetas(it); r = rs(ir); zc = zs(iz);
45+
Vertices(idx,:) = [r*cos(th), r*sin(th), zc];
46+
idx = idx+1;
47+
end
48+
end
49+
end
50+
51+
for EdgeIdx = 1:4
52+
% Draw edge from (r,theta,z) -> (r,theta,z+dz)
53+
plot3(Vertices(2*EdgeIdx-1:2*EdgeIdx,1),Vertices(2*EdgeIdx-1:2*EdgeIdx,2),Vertices(2*EdgeIdx-1:2*EdgeIdx,3),SeriesIndex="none")
54+
% Draw edge from (r,theta,z) -> (r+dr,theta,z)
55+
plot3([Vertices(EdgeIdx,1) Vertices(EdgeIdx+4,1)],[Vertices(EdgeIdx,2) Vertices(EdgeIdx+4,2)],[Vertices(EdgeIdx,3) Vertices(EdgeIdx+4,3)],SeriesIndex="none")
56+
end
57+
58+
plot3(X(1,:),Y(1,:),Z(1,:),SeriesIndex="none")
59+
plot3(X(end,:),Y(end,:),Z(end,:),SeriesIndex="none")
60+
plot3(X(:,1),Y(:,1),Z(:,1),SeriesIndex="none")
61+
plot3(X(:,end),Y(:,end),Z(:,end),SeriesIndex="none")
62+
plot3(X2(1,:),Y2(1,:),Z(1,:),SeriesIndex="none")
63+
plot3(X2(end,:),Y2(end,:),Z(end,:),SeriesIndex="none")
64+
plot3(X2(:,1),Y2(:,1),Z(:,1),SeriesIndex="none")
65+
plot3(X2(:,end),Y2(:,end),Z(:,end),SeriesIndex="none")
66+
67+
xlim([0 3])
68+
ylim([0 3])
69+
zlim([0 3])
70+
xlabel("x")
71+
ylabel("y")
72+
zlabel("z")
73+
title("Cylindrical Volume Element")
74+
CurVol = r0*dr*dtheta*dz;
75+
subtitle("V = " + CurVol)
76+
view(35,20)
77+
grid on
78+
hold off
79+
end
80+
81+
%[appendix]{"version":"1.0"}
82+
%---

FunctionLibrary/DrawKnife.m

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
function DrawKnife(xShift,yShift,zShift,Orientation)
2+
arguments
3+
xShift (1,1) double = 0
4+
yShift (1,1) double = 0
5+
zShift (1,1) double = 0
6+
Orientation (1,1) {mustBeMember(Orientation,["horizontal","vertical"])} ="horizontal"
7+
end
8+
% Paring knife blade (3-D surface)
9+
L = 1.2; % blade length (m)
10+
w0 = 0.2; % max width at heel (m)
11+
t0 = 0.003; % max thickness at heel (m)
12+
nx = 200; ny = 60;
13+
14+
x = linspace(0,L,nx); % 0 = heel, L = tip
15+
% width tapers from w0 to nearly zero at tip (use smooth polynomial)
16+
w = w0*(1 - (x./L).^2);
17+
[X,Yu] = meshgrid(x,linspace(-1,1,ny)); % Yu in [-1,1] lateral coord
18+
W = interp1(x,w,X(1,:)); % width along x
19+
20+
% Scale Yu to physical lateral coordinate per column
21+
Y = bsxfun(@times, Yu, W/2);
22+
23+
% Thickness along length (tapers faster than width)
24+
t = t0*(1 - (x./L).^1.9);
25+
T = repmat(t,ny,1);
26+
27+
% Blade geometry:
28+
% % spine (upper surface) is nearly flat with slight taper/curve
29+
% Z_top = 0.0003*(X/L) + 0.0001*(X/L).^2; % tiny curvature for realism
30+
31+
% belly (lower surface) — convex cross-section; sharper near edge
32+
% use normalized lateral coordinate s = 2*y/w in [-1,1]
33+
S = bsxfun(@rdivide, Y, W/2);
34+
% convex profile: z_drop = thickness * (1 - (1 - |s|^p) ), p controls edge sharpness
35+
p = 2.5;
36+
Z_bottom = -T .* (1 + (1 - abs(S).^p));
37+
Z_top = T.* (1 + (1 - abs(S).^p));
38+
39+
SharpenedWidth = 12;
40+
for idx = SharpenedWidth-1:-1:0
41+
Z_bottom(idx+1,:) = Z_bottom(SharpenedWidth,:)*idx/SharpenedWidth;
42+
Z_top(idx+1,:) = Z_top(SharpenedWidth,:)*idx/SharpenedWidth;
43+
end
44+
Z_bottom(60,:) = 0;
45+
Z_top(60,:) = 0;
46+
47+
% Combine to make full blade volume (plot top and bottom as surfaces)
48+
hold on
49+
% top surface
50+
hTop = surf(X+xShift, Y+yShift, Z_top+zShift, FaceColor=[0.8 0.8 0.82], EdgeColor="none");
51+
% bottom surface
52+
hBottom = surf(X+xShift, Y+yShift, Z_bottom+zShift, FaceColor=[0.75 0.75 0.78], EdgeColor="none");
53+
54+
55+
% tip edge (closing at tip)
56+
tipIdx = nx;
57+
hTip = fill3( X(1,tipIdx)*ones(1,ny)+xShift, Y(:,tipIdx)'+yShift, Z_bottom(:,tipIdx)'+zShift, [0.72 0.72 0.76], EdgeColor="none");
58+
59+
%make it look nice
60+
axis equal off; view(30,20);
61+
camlight headlight; material dull;
62+
lighting gouraud;
63+
64+
65+
xHandMax = 0+xShift;
66+
xHandMin = -0.8+xShift;
67+
yHandMin = -0.05+yShift;
68+
yHandMax = 0.1+yShift;
69+
zHandMin = -0.03+zShift;
70+
zHandMax = 0.03+zShift;
71+
72+
HandleColor = [129 65 65]./255;
73+
74+
% Vertices (8 x 3)
75+
Vhand = [ xHandMin, yHandMin, zHandMin;
76+
xHandMax, yHandMin, zHandMin;
77+
xHandMax, yHandMax, zHandMin;
78+
xHandMin, yHandMax, zHandMin;
79+
xHandMin, yHandMin, zHandMax;
80+
xHandMax, yHandMin, zHandMax;
81+
xHandMax, yHandMax, zHandMax;
82+
xHandMin, yHandMax, zHandMax ];
83+
84+
% Faces as index into V (each row is one face)
85+
Fhand = [1 2 3 4; % bottom (zMin)
86+
5 6 7 8; % top (zMax)
87+
1 2 6 5; % y = yMin side
88+
2 3 7 6; % x = xMax side
89+
3 4 8 7; % y = yMax side
90+
4 1 5 8]; % x = xMin side
91+
92+
KnifeHandle = patch(Vertices=Vhand,Faces=Fhand, ...
93+
FaceColor=HandleColor, FaceAlpha=1, EdgeColor=[0.8 0.8 0.82]);
94+
hold off
95+
96+
if Orientation == "vertical"
97+
% Rotate the entire graphic by 90 degrees about the x-axis.
98+
% Create rotation matrix for 90 degrees (pi/2) about x and apply to all plotted objects.
99+
theta = pi/2;
100+
R = [1 0 0; 0 cos(theta) -sin(theta); 0 sin(theta) cos(theta)];
101+
102+
% Collect all relevant graphics handles
103+
hObjs = [hTop; hBottom; hTip; KnifeHandle];
104+
105+
for h = hObjs.'
106+
try
107+
if isa(h,'matlab.graphics.chart.primitive.Surface')
108+
Xdata = get(h,'XData'); Ydata = get(h,'YData'); Zdata = get(h,'ZData');
109+
pts = [Xdata(:)'; Ydata(:)'; Zdata(:)'];
110+
pts = R * pts;
111+
set(h,'XData',reshape(pts(1,:),size(Xdata)), ...
112+
'YData',reshape(pts(2,:),size(Ydata)), ...
113+
'ZData',reshape(pts(3,:),size(Zdata)));
114+
elseif isa(h,'matlab.graphics.primitive.Patch')
115+
V = get(h,'Vertices');
116+
V = (R * V')';
117+
set(h,'Vertices',V);
118+
else
119+
% for other object types (e.g. patch returned by fill3)
120+
try
121+
Xdata = get(h,'XData'); Ydata = get(h,'YData'); Zdata = get(h,'ZData');
122+
pts = [Xdata(:)'; Ydata(:)'; Zdata(:)'];
123+
pts = R * pts;
124+
set(h,'XData',reshape(pts(1,:),size(Xdata)), ...
125+
'YData',reshape(pts(2,:),size(Ydata)), ...
126+
'ZData',reshape(pts(3,:),size(Zdata)));
127+
catch
128+
% ignore objects that don't support transformation
129+
end
130+
end
131+
catch
132+
% ignore any that error
133+
end
134+
end
135+
136+
else
137+
xlim([xHandMin L+xShift])
138+
end
139+
end
140+
141+
%[appendix]{"version":"1.0"}
142+
%---

FunctionLibrary/DrawPumpkin.m

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
function ax = DrawPumpkin(Cavity,Punkin,StemData,Cpunkin)
2+
% Plot the pumpkin
3+
Cavity = surf(Cavity.X,Cavity.Y,Cavity.Z,FaceColor="#fbbf4a",EdgeColor="none",FaceAlpha=0.9);
4+
5+
hold on
6+
Spunkin = surf(Punkin.X,Punkin.Y,Punkin.Z,FaceColor="interp",EdgeColor="none");
7+
8+
colormap(validatecolor({'#da8e26' '#dfc727'},'multiple'));
9+
10+
Sstem = surf(StemData.Xstem,StemData.Ystem,StemData.Zstem+StemData.heightratio^2,FaceColor="#3d6766",EdgeColor="none");
11+
12+
Pstem = patch(Vertices=[StemData.Xstem(end,:)' StemData.Ystem(end,:)' StemData.Zstem(end,:)'+StemData.heightratio^2],...
13+
Faces=1:StemData.numVerts, ...
14+
FaceColor="#b1cab5",EdgeColor="none");
15+
16+
%daspect([1 1 1])
17+
camlight
18+
19+
Cavity.CData = [];
20+
% Spunkin.CData = Cpunkin; % Pumpkin CData
21+
Sstem.CData = [];
22+
Spunkin.CData = Cpunkin+randn(StemData.numVerts)*0.022; % Pumpkin CData
23+
% set(Cavity,"CData",[])
24+
% set(Spunkin,"CData",Cpunkin); % Pumpkin CData
25+
% set(Sstem,"CData",[]); % Make sure the stem doesn't contribute to auto Color Limits
26+
% set(Spunkin,"CData",Cpunkin+randn(StemData.numVerts)*0.022); % orig 0.015
27+
daspect([1 1 1])
28+
axis off
29+
camzoom(1.8)
30+
lighting gouraud
31+
material([Spunkin Sstem Pstem],[ .6 .9 .3 2 .6 ])
32+
hold off
33+
ax = gca;
34+
end
35+
36+
%[appendix]{"version":"1.0"}
37+
%---

0 commit comments

Comments
 (0)