function final

param1 =[

    6.5861    8.6363    5.6762    9.8048    1.0000    1.0000;
    7.9183    1.5259    8.3303    1.9186    1.0000    1.0000]

param2 =[

    6.5861    8.6363    5.6762    9.8048    1.0000    1.0000;
    7.9183    1.5259    8.3303    1.9186    1.0000    1.0000]


fname = 'final13.jpg'

drawseq(param1, param2);
print('-djpeg', fname);

end

function drawseq(param1, param2)
% param: 2*6 matrix, [m,n1,n2,n3,a,b (for r1); m,n1,n2,n3,a,b (for r2)]

counter = 1;
for t=0:0.2:1,

	delta = 0.05;
	d = -1:delta:1;
	[x0,y0,z0] = meshgrid(d,d,d);

	theta = atan(y0./x0) + (x0<0).*pi;
	phi = atan(z0./sqrt(x0.^2+y0.^2));

	param = t*param1 + (1-t)*param2;

	m = param(1,1); n1 = param(1,2); n2 = param(1,3); n3 = param(1,4); a = param(1,5); b = param(1,6);
	r1 = (abs(cos(m/4*theta)/a).^n2 + abs(sin(m/4*theta)/b).^n3).^(-1/n1);

	m = param(2,1); n1 = param(2,2); n2 = param(2,3); n3 = param(2,4); a = param(2,5); b = param(2,6);
	r2 = (abs(cos(m/4*phi)/a).^n2 + abs(sin(m/4*phi)/b).^n3).^(-1/n1);

	r = r1.*r2;

	v0 = r./max(max(max(r))) - sqrt(x0.^2+y0.^2+z0.^2);


	subplot(2,3,counter); counter=counter+1;
	p = patch(isosurface(x0,y0,z0,v0,0)); 
	set(p,'FaceColor','red','EdgeColor','none','AmbientStrength',.2,'SpecularStrength',.7,'DiffuseStrength',.4);
	isonormals(x0,y0,z0,v0,p); 
	colordef white
	daspect([1 1 1]);
	axis off
	view(3); 
	camlight right; 
	camlight left; 
	lighting phong;

end
end