function door(a,b,m,g)
% Function door
% Solve for reactions in a door with two hinges assuming the system is
% 1) Statically determinate
% 2) Overdetermined
% 3) Underdetermined
%
% The top hinge is 1, and bottom hinge is 2, the door width is a, the
% door height is b, the door mass is m, and gravity is g.
%
% Usage:  door(a,b,m,g)
%
% Inputs: a = door width
%         b = door height
%         m = door mass
%         g = gravity
%         The top hinge is 1 and the bottom hinge is 2.
%         X is to the right, Y is verticall upward, and Z is
%         perpendicular to the door face. The hinges are on the
%         left side of the door.
%
% Outputs: Hinge reaction loads Fx, Fy, and Mz for each of the three cases above,
%          where x = [Fx; Fy; Mz]

format short

% Case 1: Statically determinate system
% Set reaction loads on hinge 2 = 0

disp('Case 1: Statically determinate system')

A1 = eye(3)
detA1 = det(A1)
rankA1 = rank(A1)
condA1 = cond(A1)
rcondA1 = rcond(A1)
b1 = [0; m*g; m*g*a/2]

% Solve using the backslash operator
disp('Solution using backslash operator')
x1 = A1\b1
pause

% Solve using the matrix inverse
disp('Solution using matrix inverse')
invA1 = inv(A1)
x1 = invA1*b1
pause

% Case 2: Overdetermined system
% Re-solve case 1 but assume we make 3 experimental measurements of mg and
% m*g*a/2 that are all slightly different.

disp('Case 2: Overdetermined system')

A2 = [A1; A1; A1]
% detA2 = det(A2) % Won't work since A2 not square
rankA2 = rank(A2)
condA2 = cond(A2)
% rcondA2 = rcond(A2) % Won't work since A2 not square
b2 = [0;     m*g;      m*g*a/2;...
      0.01;  0.99*m*g; 1.01*m*g*a/2;...
      -0.01; 1.01*m*g; 0.99*m*g*a/2]

disp('Least-squares solution using backslash operator')
x2 = A2\b2
errors = A2*x2-b2
pause

% Case 3: Underdetermined system
% Use the 3 equations from case 1 but solve for all 6 reaction loads

disp('Case 3: Underdetermined system')

A3 = [1 0 0 1 0 0;...
      0 1 0 0 1 0;...
      0 0 1 b 0 1]
% detA3 = det(A3) % Won't work since A3 not square
rankA3 = rank(A3)
condA3 = cond(A3)
% rcondA3 = rcond(A3) % Won't work since A3 not square
b3 = b1

% Solve using the backslash operator
disp('Solution using backslash operator')
x3 = A3\b3
normx3 = norm(x3)
pause

% Solve using the pseudo inverse
disp('Solution using matrix pseudo inverse')
pinvA3 = pinv(A3)
x3 = pinvA3*b3
normx3 = norm(x3)
pause

% Question: How could we now solve this as an optimization problem? Let's
% formulate the problem to decide.
%
% Design variables x = [Fx1 Fy1 Mz1 Fx2 Fy2 Mz2]'
%
% Cost function: min sum(xi^2)
%
% Constraints: A3*x = b3 (linear)
%
% Question: What optimization problem type is this? Constrained quadratic
% programming. Matlab has a nice quadratic programming solver that was can
% use to solve this problem.
% 
% QUADPROG Quadratic programming. 
%    X=QUADPROG(H,f,A,b) solves the quadratic programming problem:
% 
%             min 0.5*x'*H*x + f'*x   subject to:  A*x <= b 
%              x    
%
%    X=QUADPROG(H,f,A,b,Aeq,beq) solves the problem above while additionally
%    satisfying the equality constraints Aeq*x = beq.

disp('Solution using quadratic programming optimization')
H = eye(6);
f = zeros(1,6);
A = []; % No inequality constraints, so
b = []; % make these empty matrices
Aeq = A3;
beq = b3;
xqp = quadprog(H,f,A,b,Aeq,beq)