function [F,z,coordbott,coordtop,l_top_comp,l_vert_comp] = ...
           trusscheck(theta,l_vert,l_bott,l_top,iprint)

% function [F,z,coordbott,coordtop,l_top_comp] = ...
%            trusscheck(theta,l_vert,l_bott,l_top,iprint)
% Check that the solution to the trussproblem
% satisfies all of the conditions.
%
% Input:
%  theta = angles between horizontal and the three
%          sides of the truss that have two vertices
%          on the bottom triangle and one on the top.
%  l_vert = length of the 6 vertical trusses (scalar).
%  l_bott = length of the 3 bottom trusses (scalar).
%  l_top  = desired lengths of the 3 top trusses (1x3).
%  iprint = 0 if no printing desired
%         > 0 to print a 3-line report
%
% Output:
%  F:     the residuals for the nonlinear system to be solved
%            the differences between the squares of the lengths 
%            of the top triangle and the squares of the
%            prescribed lengths.
%  z:     -log_10 of the max error in truss length
%              (=30 if error is zero)
%  coordbott: The rows contain the coordinates of the
%         three vertices of the bottom triangle.
%  coordtop : The rows contain the coordinates of the
%         three vertices of the top triangle.
%  l_top_comp  = computed lengths of the 3 top trusses (1x3).
%  l_vert_comp= computed lengths of the 6 vertical trusses (1x6).
%
% Dianne P. O'Leary 05/2008

% Compute the bottom vertices.  
% The bottom is an % equilateral triangle with vertices
%     (0,        0,        0)
%     (0,        1,        0)*l_bott
%     (cos(pi/6),sin(pi/6),0)*l_bott

c = cos(pi/6);
s = sin(pi/6);

coordbott = [0,        0,        0
             0,        1,        0
             c,        s,        0]*l_bott;

% Compute the top vertices. 
% Each has the coordinates of the point halfway inbetween
% the other two vertices on the bottom, shifted
% by h*sin(theta) vertically and d=h*cos(theta) horizontally,
% where h is determined so that the distance to one of the
% bottom vertices is equal to l_vert.
% Note that the horizontal shift must be applied with angle
%
%    pi/3 for the 1st vertex
%   -pi/3 for the 2nd vertex
%    0    for the 3rd vertex

midpt = (coordbott(2,:)+coordbott(3,:))/2;
hA    = sqrt(l_vert^2 - norm(midpt-coordbott(2,:))^2);
dA = hA*cos(theta(1));
coordtop(1,:) = midpt + [dA*cos(pi/3), dA*sin(pi/3),   hA*sin(theta(1))];

midpt = (coordbott(1,:)+coordbott(3,:))/2;
hB    = sqrt(l_vert^2 - norm(midpt-coordbott(1,:))^2);
dB = hB*cos(theta(2));
coordtop(2,:) = midpt + [dB*cos(-pi/3), dB*sin(-pi/3), hB*sin(theta(2))];

midpt = (coordbott(1,:)+coordbott(2,:))/2;
hC    = sqrt(l_vert^2 - norm(midpt-coordbott(1,:))^2);
dC = hC*cos(theta(3));
coordtop(3,:) = midpt + [dC,             0,            hC*sin(theta(3))];

% Check the lengths of the six vertical trusses.

l(1) = norm(coordbott(1,:)-coordtop(3,:)) - l_vert;
l(2) = norm(coordbott(2,:)-coordtop(3,:)) - l_vert;

l(3) = norm(coordbott(1,:)-coordtop(2,:)) - l_vert;
l(4) = norm(coordbott(3,:)-coordtop(2,:)) - l_vert;

l(5) = norm(coordbott(2,:)-coordtop(1,:)) - l_vert;
l(6) = norm(coordbott(3,:)-coordtop(1,:)) - l_vert;

% Check the lengths at the top triangle.

l_top_comp(1) = norm(coordtop(1,:)-coordtop(2,:));
l_top_comp(2) = norm(coordtop(2,:)-coordtop(3,:));
l_top_comp(3) = norm(coordtop(1,:)-coordtop(3,:));

for k=1:3,
    t(k) = l_top_comp(k)   - l_top(k);
    F(k) = l_top_comp(k)^2 - l_top(k)^2;
end

% Print a report (if desired)
% and set z=0 if the relative
% error in any length is greater than 10*eps_mach.

z = max(max(abs(t)./l_top),max(abs(l)/l_vert));

if (iprint > 0)
   sprintf('Vertical   errors in length: %e \n',l)
   sprintf('Horizontal errors in length: %e \n',t)
   sprintf('Maximum rel error in length: %e',   z)
end

if (z == 0)
   z = 30;
else 
   z = round(-log10(z));
end

l_vert_comp = l;

