global myl_vert myl_bott myl_top mylambda mythetatrue myscale myl_top_easy myff myhh myfangles myhangles

% This program solves the problems in
% ``Case Study: Variable-Geometry Trusses: What's Your Angle?"
%
% It determines the configurations of various platforms, given
% the lengths of the 12 trusses that form them.
%
% We compare the use of fsolve with the use of a homotopy method.
% We compare two formulations, one geometric and one polynomial.
%
% Dianne P. O'Leary 05/2008

% The data for the four problems come from the thesis of 
% V. Arun, Virginia Polytech, 1990:
% 
%                       myl_vert = 34.5;
%                       myl_bott = 48;
%
%           Fig 6.4:    myl_top = [35, 30, 40];  8 solutions
%           Fig 6.6:    myl_top = [30, 30, 30];  8 solutions
%           Fig 6.8:    myl_top = [35, 40, 45];  4 solutions
%           Fig 6.9-10: myl_top = [10, 10, 10]; 16 solutions

% Turn on a diary.

s = date;
c = clock;
t = c(4)+c(5)/100;
diaryname = strcat('diary',s,'.',num2str(t));

diary(diaryname)

arun = [35, 30, 40
        30, 30, 30
        35, 40, 45
        10, 10, 10];
arunname = ['aruna';'arunb';'arunc';'arund'];

myl_vert = 34.5;
myl_bott = 48;

% Run the methods on the four problems.

for karun = 1:4

        
   probname=arunname(karun,:);
   myl_top = arun(karun,:);

   disp(sprintf('ARUN TRUSS PROBLEM %d',karun))

   [thetaf{karun}{1},thetah{karun}{1}] = methodgeom(probname,karun,0,0);
   [thetaf{karun}{2},thetah{karun}{2}] = methodpoly(probname,karun,0,0);

end % for karun problem number


% Print summary

disp('yyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyy')

for karun = 1:4

   myl_top = arun(karun,:);

   for meth=1:2,

     topf = [];
     toph = [];
     fcoords = myff{karun}{meth};
     hcoords = myhh{karun}{meth};
     [mf,nf] = size(fcoords);
     [mh,nh] = size(hcoords);
     disp(sprintf('%d fsolve coordinates for problem %d method %d',mf,karun,meth))
     fcoords
     disp(sprintf('%d fsolve angles for problem %d method %d',mf,karun,meth))
     myfangles{karun}{meth} 
     disp(sprintf('%d homotopy coordinates for problem %d method %d',mh,karun,meth))
     hcoords
     disp(sprintf('%d homotopy angles for problem %d method %d',mh,karun,meth))
     myhangles{karun}{meth} 
     mmh(karun,meth) = mh;
     mmf(karun,meth) = mf;
   end %for meth
end    % for karun

disp('Number of fsolve solutions: rows correspond to problem number,')
disp('       columns correspond to geometric and polynomial methods,')
mmf
disp('Number of homotopy solutions: rows correspond to problem number,')
disp('       columns correspond to geometric and polynomial methods,')
mmh

% Paranoia check

disp('We now check each solution, computing z = log of max error')
disp('in truss length.  If z < 5, the solution should be rejected.')
arun = [35, 30, 40
        30, 30, 30
        35, 40, 45
        10, 10, 10];

myl_vert = 34.5;
myl_bott = 48;

for karun=1:4,
 myl_top = arun(karun,:);
 for meth = 1:2,
    for k=1:2,
       if (k==1)
          thetas = myfangles{karun}{meth};
          [m,n] = size(thetas);
          disp(sprintf('Checking the %d fsolve solutions for problem %s, method %d',...
           n, arunname(karun,:), meth))
       else
          thetas = myhangles{karun}{meth};
          [m,n] = size(thetas);
          disp(sprintf('Checking the %d homotopy solutions for problem %s, method %d',...
           n, arunname(karun,:), meth))
       end
       zcheck = [];
       for j=1:n
           [F,z,coordbott,coordtop,l_top_comp,l_vert_comp] = ...
                    trusscheck(thetas(:,j),myl_vert,myl_bott,arun(karun,:),0);
           zcheck = [zcheck, z];
       end
       disp(zcheck)
    end
 end
end

% Turn the diary off

diary off
