simulate_mms.m

r0 = [2,0,0];
r1 = [-1,-1.732,0];
r2 = [-1,1.732,0];
r3 = [0,0,2.8284];
r4 = [0,0,-2.8284];
r5 = [0,2,2];

r = [r2;r5;r3;r0];

map0 = [3,4; 2,4; 3,2; 1,4; 3,1; 2,1];
map1 = [3,4; 3,2; 1,4; 2,4; 3,1; 1,2];
map2 = [2,4; 3,4; 1,4; 2,3; 3,1; 2,1];

Vmat = [ [1;r(1,:)'], [1;r(2,:)'], [1;r(3,:)'], [1;r(4,:)'] ];
V = det(Vmat)/6;

a1 = det([r(2,:)',r(3,:)',r(4,:)']);
b1 = -det([[1;r(2,2:3)'],[1;r(3,2:3)'],[1;r(4,2:3)']]);
c1 = det([[1;r(2,1:2:3)'],[1;r(3,1:2:3)'],[1;r(4,1:2:3)']]);
d1 = -det([[1;r(2,1:2)'],[1;r(3,1:2)'],[1;r(4,1:2)']]);

a2 = -det([r(1,:)',r(3,:)',r(4,:)']);
b2 = det([[1;r(1,2:3)'],[1;r(3,2:3)'],[1;r(4,2:3)']]);
c2 = -det([[1;r(1,1:2:3)'],[1;r(3,1:2:3)'],[1;r(4,1:2:3)']]);
d2 = det([[1;r(1,1:2)'],[1;r(3,1:2)'],[1;r(4,1:2)']]);

a3 = -det([r(2,:)',r(1,:)',r(4,:)']);
b3 = det([[1;r(2,2:3)'],[1;r(1,2:3)'],[1;r(4,2:3)']]);
c3 = -det([[1;r(2,1:2:3)'],[1;r(1,1:2:3)'],[1;r(4,1:2:3)']]);
d3 = det([[1;r(2,1:2)'],[1;r(1,1:2)'],[1;r(4,1:2)']]);

a4 = -det([r(2,:)',r(3,:)',r(1,:)']);
b4 = det([[1;r(2,2:3)'],[1;r(3,2:3)'],[1;r(1,2:3)']]);
c4 = -det([[1;r(2,1:2:3)'],[1;r(3,1:2:3)'],[1;r(1,1:2:3)']]);
d4 = det([[1;r(2,1:2)'],[1;r(3,1:2)'],[1;r(1,1:2)']]);

a = [a1,a2,a3,a4]; b = [b1,b2,b3,b4]; c = [c1,c2,c3,c4]; d = [d1,d2,d3,d4];

L = [a1,a2,a3,a4;
     b1,b2,b3,b4;
     c1,c2,c3,c4;
     d1,d2,d3,d4];

T = zeros(4,3,6);

for j=1:6

    j1 = map2(j,1); j2 = map2(j,2);

    T(:,1,j) = (b(j2)*L(:,j1) - b(j1)*L(:,j2))/(36*V^2);
    T(:,2,j) = (c(j2)*L(:,j1) - c(j1)*L(:,j2))/(36*V^2);
    T(:,3,j) = (d(j2)*L(:,j1) - d(j1)*L(:,j2))/(36*V^2);

end

U = [[0;0;0;1],[0;1;0;0],[0;0;1;0]];

BigU = [U(:,1);U(:,2);U(:,3)];

BigT = [[T(:,1,1);T(:,2,1);T(:,3,1)],...
        [T(:,1,2);T(:,2,2);T(:,3,2)],...
        [T(:,1,3);T(:,2,3);T(:,3,3)],...
        [T(:,1,4);T(:,2,4);T(:,3,4)],...
        [T(:,1,5);T(:,2,5);T(:,3,5)],...
        [T(:,1,6);T(:,2,6);T(:,3,6)]];
BigT\BigU