1、function demoDisplayclearclcclose all;x,y,z = heart;bottomZ = min(z(:);radius = abs(bottomZ);axis vis3daxis offdaspect(1.6, 1, 1.875);hold on;markerColor = 0.5, 0.5, 0.5;lineColor = 0.7, 0.7, 0.7;figWidth = 800figHeight = 600;screenSize = get(0, ScreenSize);X0 = (screenSize(3)-figWidth)/2;Y0 = (scre

2、enSize(4)-figHeight)/2;set(gcf, Position, X0,Y0, figWidth, figHeight);camtarget(0, 0, 0);set(gcf,Renderer,zbuffer);set(gcf,DoubleBuffer,on);% RGB(255 105 180) is for pink colorset(gcf, color, 255 255 255/255);startAZ = 30;satrtEL = 20;view(startAZ, satrtEL);% draw pointsfor iLoop = 1:size(x, 1)plot3

3、(x(iLoop, :), y(iLoop,:), z(iLoop, :), Marker, ., MarkerEdgeColor, markerColor, LineStyle, none);pause(0.2);end% draw wire framesfor iLoop = 1:size(x, 1)plot3(x(iLoop, :), y(iLoop,:), z(iLoop, :), color, lineColor, LineStyle, -);pause(0.1);endfor iLoop = 1:size(x, 2)plot3(x(:, iLoop), y(:, iLoop), z

4、(:, iLoop), color, lineColor, LineStyle, -);pause(0.1);end% rotate wire framefor i=1:36camorbit(10,0,data);pause(0.1);end% clean the curent axesdelete(get(gca, Children);camlight left;newLineColor = 0.7, 0.4, 0.4;% re-draw the whole heart with surf functionh = surf(x,y,z, EdgeColor, newLineColor, Fa

5、ceColor, r);for i=1:36camorbit(10,0,data);pause(0.1);endmediumColor = 0.8, 0.3, 0.3;set(h, EdgeColor, mediumColor,FaceLighting,gouraud);for i=1:36camorbit(10,0,data);drawnow;endset(h, EdgeColor, none);drawnow;% define the Z-path for heart re-shapedeltaZ = linspace(-.3, 0, 5);deltaZ = fliplr(deltaZ),

6、 deltaZ;deltaZ = repmat(deltaZ, 1, 3);zlim(manual);for iLoop = 1:length(deltaZ)% draw the bouncing heartcurDeltaZ = deltaZ(mod(iLoop,length(deltaZ)+1);ratio = (radius+curDeltaZ)*radius;newZ = z*ratio + curDeltaZ;if exist(h, var)& ishandle(h), delete(h);end;h = surf(x,y, newZ, EdgeColor, none, FaceCo

7、lor, r,FaceLighting,gouraud);camorbit(-1, 0, data);drawnow;end% draw the bouncing heartdeltaZ = generateSimulatedBouncingPath;zlim(-1, 3);for iLoop = 1:length(deltaZ)if exist(h, var)& ishandle(h), delete(h);end;curDeltaZ = deltaZ(iLoop);if curDeltaZ 0% when deltaZ is under the zero bar, its hitting

8、the groundratio = (radius+curDeltaZ)*radius;newZ = z*ratio + curDeltaZ;else% when deltaZ is bigger than zero, its a free objectnewZ = z + curDeltaZ;endh = surf(x, y, newZ, EdgeColor, none, FaceColor, r,FaceLighting,gouraud);% zoom in and out for more vivid effectif iLoop = length(deltaZ)/2camzoom(0.

9、99);elsecamzoom(1.005);endcamorbit(390/length(deltaZ), 0, data);drawnow;endtext(0.3,2.4, X,FontSize,25,FontName,Times New Roman)%text(0.3,1.8, X,FontSize,30,FontName,Times New Roman)%endfunction z = generateSimulatedBouncingPath(iterationNum, roundTimes)% define the z path of the center of ball for

10、rountTimes of bouncing along z-axisif nargin 2roundTimes = 3;endif nargin 1iterationNum = 40;endradius = 1;deltaH = 1.5; % the highest point comapred with radiusdeltaL = -0.3;% the lowest point comapred with radiusrangeY = -2, 1;% gravity accelerator constant% Note: this is smaller to the real 9.8 v

11、alue, but this makes the bouncing slowerg = 5;% define the accelerator for the z 0plusZtime = sqrt(2*deltaH*radius/g);% time for z 0zPlus = startZ - 0.5 * g * qPeriodTime(qPeriodTimeplusZtime).2;% z path when z =plusZtime)-plusZtime) + 0.5 * k* (qPeriodTime(qPeriodTime =plusZtime)-plusZtime).2;z = z

12、Plus, zMinus, fliplr(zMinus), fliplr(zPlus);% the last zPlus is for the stop actionz = repmat(z, 1, roundTimes), zPlus;y = linspace(rangeY(1), rangeY(2), numel(z);endfunction cordX, cordY, cordZ = heart(sizeTheta, sizeFai)if nargin =0sizeTheta = 30;sizeFai = 40;elseif nargin 2error(heart:NotEnoughAr

13、ugments, Please give at least two input arugments);end;theta = linspace(0, pi, sizeTheta);nudge = 0.0001; % used to avoid the operlappingfai = linspace(0 + nudge, pi/2 - nudge, round(sizeFai/4);a = 9/4;b = 9/80;A = 1+(a-1)*(sin(theta).2) * (sin(fai).2);B = (sin(theta).2.*cos(theta).3) * (1 + (b-1)*(

14、sin(fai).2);rou = zeros(size(A);for iLoop = 1:numel(A);curA = A(iLoop);curB = B(iLoop);% this is the polar coordinates version of the sextic equation found on% /HeartSurface.htmlpolyFactors = curA3, -curB, -3*curA2, 0, 3*curA, 0, -1;solutions = roots(polyFactors);realRou = real(solutions(abs(imag(solutions)0);end% x,y,z for the quater of whole heartx = rou .* (sin(theta)