% What this is modelling:
% psiN(x) = aN + bNixi
% aN = 1/6C * alternating(NRST) * alternating(ijk)*xRi*xSj*xTk
% bNi = -1(6C) * alternating(NRST) *
%         alternating(mjk)*betaRm(i)*betaSj(i)*betaTk(i)



% Triangle points
p1 = [4,0]
p2 = [8,0]
p3 = [4,7]

tArea = .5*((7-0)*(8-4))

points = buildPoints(p1,p2,p3)


% Not in use yet, these will be values at nodes
% c1 = 255
% c2 = 180
% c3 = 240
% c4 = 170

coords = [7,.2]
nodalVals = [20,40,30]

coeffs = buildShapeCoeff(tArea, coords, points)


% Alternating tensor testing examples
% a = [ 4 3 2 1] % 1
% b = [ 1 3 2 ] % -1
% c = [ 2 3 2 ] % 0
% res = alternating(a)
% res = alternating(b)
% res = alternating(c)

final = nodalVals(1)*coeffs(1) + nodalVals(2)*coeffs(2)+nodalVals(3)*coeffs(3)

% Nonsense line so you can set a breakpoint here when debugging
setDebuggingBreakpointHere = 1;



% Alternating tensor, takes array as input
function res = alternating(input)
    I = speye(length(input));
    res = det(I(:,input));
end

% build i x j C shape coefficients
function shapeCoeffRes = buildShapeCoeff(tArea, coords, points)
    shapeCoeffRes = zeros(3) % row count, col count
    shapeCoeffRes(1) = (1/(2*tArea))*((points(2,1)*points(3,2))-(points(3,1)*points(2,2))+(points(2,2)-points(3,2))*coords(1)+(points(3,1)-points(2,1))*coords(2))
    shapeCoeffRes(2) = (1/(2*tArea))*((points(3,1)*points(1,2))-(points(1,1)*points(3,2))+(points(3,2)-points(1,2))*coords(1)+(points(1,1)-points(3,1))*coords(2))
    shapeCoeffRes(3) = (1/(2*tArea))*((points(1,1)*points(2,2))-(points(2,1)*points(1,2))+(points(1,2)-points(2,2))*coords(1)+(points(2,1)-points(1,1))*coords(2))
end

% i x j points
function pointsRes = buildPoints(point1, point2, point3)
    pointsRes = zeros(3,2)
    pointsRes(1,:) = point1
    pointsRes(2,:) = point2
    pointsRes(3,:) = point3

end