%param1 = [7 .2 1.7 1.7 1 1; 7 .2 1.7 1.7 1 1]
%param2 = [6 1 1 1 1 1; 3 1 1 1 1 1]
%param3 = [5.6 0.28 1.56 1.96 1 1; 1.2 0.64 0.74 0.74 1 1]

function shapetwo

param1t = [7 .2 1.7 1.7 1 1 7 .2 1.7 1.7 1 1]
param2t = [6 1 1 1 1 1 3 1 1 1 1 1]
param3t = [5.6 0.28 1.56 1.96 1 1 1.2 0.64 0.74 0.74 1 1]

%checking linear independence
rref([param1t; param2t; param3t])
%should get 3 leading 1s

param1 = [7 .2 1.7 1.7 1 1; 7 .2 1.7 1.7 1 1]
param2 = [6 1 1 1 1 1; 3 1 1 1 1 1]
param3 = [5.6 0.28 1.56 1.96 1 1; 1.2 0.64 0.74 0.74 1 1]

%set arbitrary values for a, b, and c
a = 0.5;
b = 0.25;
c = 0.25;

param4 = a * param1 + b * param2 + c * param3

%convert from 8 parameters to 2 parameters
%map from R8 to R2

%define transformation matrix A
%A = [a11 a12; a21 a 22; a31 a32; a41 a42; a51 a52; a61 a62; a71 a72; a81 a82];

%define A transpose as At
%At = A';

fname = 'test.jpg'

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

end

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

counter = 1;
for i=0:0.25:1,
    for t=0:0.25: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 = (1-i) * ((1-t)*param1 + t*param2) + i * ((1-t)*param3 + t*param4);

        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(5,5,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
end