clear all; 
close all;

k = 0;
for i = -4 : 1 : 0
    for j = -4-i:1:i+4
        k = k+1;
        P(k) = struct('x',i,'y',j);
    end
end
for i = 1:1:4
    for j = -4+i:1:4-i
        k = k+1;
        P(k) = struct('x',i,'y',j);
    end
end

m = 0;
k = 0;
l = 0;
for i = -4 : 1 : -1
    m = m+1;
    k = k + 2*m-1;
    for j = 1 : 2 : (2*m-1)
        kk = k-2*(m-1)+j-1;
        l = l+1;
        V(l) = struct('v1',kk,'v2',kk+2*(m-1)+1,'v3',kk+2*(m-1)+2);
        l = l+1;
        V(l) = struct('v1',kk,'v2',kk+2*(m-1)+2,'v3',kk+2*(m-1)+3);
    end
    for j = 2 : 2 : (2*m-2)
        kk = k-2*(m-1)+j-1;
        l = l+1;
        V(l) = struct('v1',kk,'v2',kk-1,'v3',kk+2*(m-1)+2);
        l = l+1;
        V(l) = struct('v1',kk,'v2',kk+1,'v3',kk+2*(m-1)+2);
    end
end

for i = 0 : 1 : 3
    if (i == 0)
        m = m+1;
    else
        m = m-1;
    end
    k = k + 2*m-1;
    for j = 3 : 2 : (2*m-3)
        kk = k-2*(m-1)+j-1;
        l = l+1;
        V(l) = struct('v1',kk,'v2',kk+2*(m-1)-1,'v3',kk+2*(m-1));
        l = l+1;
        V(l) = struct('v1',kk,'v2',kk+2*(m-1),'v3',kk+2*(m-1)+1);
    end
    for j = 2 : 2 : (2*m-2)
        kk = k-2*(m-1)+j-1;
        l = l+1;
        V(l) = struct('v1',kk,'v2',kk-1,'v3',kk+2*(m-1));
        l = l+1;
        V(l) = struct('v1',kk,'v2',kk+1,'v3',kk+2*(m-1));
    end
end

Vlength = length(V);
figure(1); hold on;
for i = 1 : Vlength
    plot([ P(V(i).v1).x, P(V(i).v2).x ],[ P(V(i).v1).y, P(V(i).v2).y ]);
    plot([ P(V(i).v2).x, P(V(i).v3).x ],[ P(V(i).v2).y, P(V(i).v3).y ]);
    plot([ P(V(i).v3).x, P(V(i).v1).x ],[ P(V(i).v3).y, P(V(i).v1).y ]);
end
hold off;

for i = 1 : 5
    for j = 1 : 5
        x(i,j) = -4 + i - 1 + j - 1;
        y(i,j) = 1 - i + j - 1;
    end
end
u = cos(x).*sin(y);

figure(2); 
xI = -1 : 0.25 : 1;
h = contour(x,y,u,xI); clabel(h); hold on;
x1 = -4 : 0.5 : 0;
y1 = 4 + x1;
y2 = -4 - x1;
plot(x1,y1,'b',x1,y2,'b');
x2 = 0 : 0.5 : 4;
y3 = 4 - x2;
y4 = -4 + x2;
plot(x2,y3,'b',x2,y4,'b');

figure(3);
axis([-4,4,-4,4,-1,1]); hold on;
xlabel('x'); ylabel('y'); zlabel('u'); title('The trigonometric functions on a triangular grid');
view(-20,20);
surf(x,y,u); 
hold off;

xr = -4 : 1 : 4; yr = -4 : 1 : 4; [Xr,Yr] = meshgrid(xr,yr);
Ur = cos(Xr).*sin(Yr);
figure(4); h = contour(Xr,Yr,Ur,xI); clabel(h);
figure(5); axis([-4,4,-4,4,-1,1]); hold on;
xlabel('x'); ylabel('y'); zlabel('u'); title('The trigonometric functions on a triangular grid');
view(-20,20);
surf(Xr,Yr,Ur); 
hold off;
