function basis
e1 = [7 .2 1.7 1.7 1 1; 7 .2 1.7 1.7 1 1];
e2 = [6 1 1 1 1 1; 3 1 1 1 1 1];

for i=1:1:8,
    e1(i) = 0.5+rand(1)*7;
    e2(i) = 0.5+rand(1)*7;
end

e1
e2

fname = 'good6.jpg';

drawseq(e1, e2);
print('-djpeg', fname);
end


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

warning off all

counter = 1;

ymin = 0.25;
ystep = 0.25;
ymax = ymin + 4*ystep;

xmin = 0.25;
xstep = 0.25;
xmax = xmin + 4*xstep;

for y=ymin:ystep:ymax,
    for x=xmin:xstep:xmax,

        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 = [x*e1(1:4) 1 1; x*e1(7:10) 1 1] + [y*e2(1:4) 1 1; y*e2(7:10) 1 1] ;

        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