File: ef230-2023-08/www/examples/truss_wrs.m Download
%! A large truss analysis program with many different examples
%! Creates plots of geometry, results, FBDs, and equations
% Univeristy of Tennessee : Will Schleter

%=====================================================================================%
function main
% demonstrate modular program with graphics
% plots a truss specified by joint coordinates and a connection list

close all;
clc;

% note: global variables are the easiest way to use callback functions required by the
%       user interface routines.
global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL

% load the truss geometry (coordinates and connection list)
opts={'1','2','3','4','5','6','98','99','h',...
        'oa','ob','oc','od','oe','of',...
        'wa','wb','wc','wd','we','wf'};
        
ch = menu('Truss Data',opts);
if (ch==0), error('Cancelled'); end

wrsTrussGeom1(opts{ch});

[j x] = size(gXY);
[m x] = size(gLinks);


if (j*2 ~= m+3)
  j
  m
  error('Wrong number of joints and/or members')
end

% build coefficient list based upon truss geometry
buildcoef;

UpdateTruss;

return

%=====================================================================================%
function UpdateTruss
global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL

% plot the truss and associated loads, forces and reactions

%define load matrix
gLoad = zeros(gNE,1);
for i=1:size(gExtLoad,1)
   addload(gExtLoad(i,:));
end

%checks for "correct" coefficient matrix
%condition number should be < 100
fprintf('Coefficient condition number: %f\n',cond(gCoef)');

%each column should add up to 0, last three columns should add up to 1 or -1
nr = size(gCoef,2);
for ii=1:nr
    tmpsum = sum(gCoef(:,ii));
    cmp = 0;
    if ii > (nr-3); cmp = 1; end
	if (abs(tmpsum) - cmp) > .01
        fprintf('Problem: summation of column %u = %5.2f\n',ii,tmpsum);
    end
end


gForce = gCoef \ (-gLoad);

initfig(gNJ);

plottruss(4); %the system FBD
plottruss(1); %the results
drawfbd;

return

%=====================================================================================%
function r1=jn(n)
r1 = n-'@';
return


%=====================================================================================%
function r1=d2r(x)
%convert degrees to radians
r1=x*pi/180;
return;

%=====================================================================================%
function r1=r2d(x)
%convert radians to degrees
r1=x*180/pi;
return;
  
%=====================================================================================%
function r1=cosd(x)
%cos function with input in degrees
r1=cos(d2r(x));
return

%=====================================================================================%
function r1=sind(x)
%sin function with input in degrees
r1=sin(d2r(x));
return;

%=====================================================================================%
function r1=ptoff(xy,d,a)
%generate an offset of point xy at an angle (a) and distance (d)
r1(1)=xy(1)+d*cosd(a);
r1(2)=xy(2)+d*sind(a);
return;

%=====================================================================================%
function vecplot(xy,d,a,col)
%plot representation of a vector, complete with arrowhead and text value
xyf(1,:) = xy;
% vv = d*vsf; %vector scale factor
vv = 2;
if (d<0)
   % reverse direction for negative magnitude
   d = abs(d);
   if (a>0) a=a-180;
   else a=a+180;
   end
end
xyf(2,:) = ptoff(xy,vv,a);
plot(xyf(:,1),xyf(:,2),col);

% draw the arrowhead
xyf(1,:) = ptoff(xyf(2,:),vv/5,a+180-20);
xyf(3,:) = ptoff(xyf(2,:),vv/5,a+180+20);
fill(xyf(:,1),xyf(:,2),col);

% draw text label
tal = 'left';
xoff = .1;
if (abs(a)>90)
   tal = 'right';
   xoff = -.1;
end
txt = sprintf('%.0f',d);
text(xyf(2,1) + xoff,xyf(2,2),txt,'HorizontalAlignment',tal,'FontSize',9);
%text((xyf(1,1)+xyf(2,1))/2,(xyf(1,2)+xyf(2,2))/2,sprintf('%.0f',d),'HorizontalAlignment','center');
return

%=====================================================================================%
function r1=forcename(p1,p2)
if (p1>p2)
   tmp=p1;
   p1=p2;
   p2=tmp;
end

r1 = sprintf('F_{%c%c}',p1+64,p2+64);
return

%=====================================================================================%
function r1=forcename1(n)
global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL
switch n-gNL
   case 1
      r1 = 'Px';
   case 2
      r1 = 'Py';
   case 3
      r1 = 'R';
   otherwise
      p1=gLinks(n,1);
      p2=gLinks(n,2);
      if (p1>p2)
         tmp=p1;
         p1=p2;
         p2=tmp;
      end
      r1 = sprintf('%c%c',p1+64,p2+64);
end
return

%=====================================================================================%
function r1 = pointname(p1)
r1 = sprintf('%c',p1+64);
return

%=====================================================================================%
function tbox(p1,dx,dy,col)
xy(1,:) =p1 + [dx dy];
xy(2,:) =p1 + [-dx dy];
xy(3,:) =p1 + [-dx -dy];
xy(4,:) =p1 + [dx -dy];
fill(xy(:,1),xy(:,2),col);
return

%=====================================================================================%
function vecplot1(p1,p2,org,d,txt,col)
%plot representation of a vector, complete with arrowhead and text value

a = r2d(atan2(p2(2)-p1(2),p2(1)-p1(1)));
xyv(1,:) = org;
xyv(2,:) = ptoff(org,d,a); 
plot(xyv(:,1),xyv(:,2),col);

% draw the arrowhead
xyv(1,:) = ptoff(xyv(2,:),.25,a+180-20);
xyv(3,:) = ptoff(xyv(2,:),.25,a+180+20);
fill(xyv(:,1),xyv(:,2),col);

% draw text label
tal = 'left';
xoff = .1;
if (abs(a)>90)
   tal = 'right';
   xoff = -.1;
end
txt = sprintf('%s _{%.0f\\circ}',txt,a);
text(xyv(2,1) + xoff,xyv(2,2),txt,'HorizontalAlignment',tal,'FontSize',9);
return

%=====================================================================================%
function plottruss(fig)
%plot the given truss
%loop thru the connectivity matrix (one loop per member)
global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL

figure(fig);
for i=1:gNL
   %extract the members joint coordinates
   j1 = gLinks(i,1);
   j2 = gLinks(i,2);
   xym(1,:)=gXY(j1,:);
   xym(2,:)=gXY(j2,:);
   if (fig==1)
      %determine color of member based upon the internal load
      if (gForce(i)>0.5)
         col = 'bo-';
      elseif (gForce(i)<-0.5)
         col = 'ro-';
      else
         col = 'ko-';
      end
      plot(xym(:,1),xym(:,2),col);
      %display internal load values
      aa = membang(gXY,j1,j2);
      if (aa>90) aa=aa-180; end;
      if (aa<-90) aa=aa+180; end;
      text((xym(1,1)+xym(2,1))/2,(xym(1,2)+xym(2,2))/2,sprintf('%.0f',abs(gForce(i))),'HorizontalAlignment','center','rotation',aa);
   else
      plot(xym(:,1),xym(:,2),'k-','LineWidth',3);
   end
end

% display all loads
for i=1:2:length(gLoad)
    ix=i;
    iy=ix+1;
    ij=iy/2;
    if (gLoad(ix,1)~=0 | gLoad(iy,1)~=0)
        vecplot(gXY(ij,:),sqrt(gLoad(ix).^2 + gLoad(iy).^2),r2d(atan2(gLoad(iy),gLoad(ix))),'m');
    end
end

% pin reaction force
ix=gNL+1;
iy=gNL+2;
if (fig==1)
   %vecplot(gXY(gPJ,:),sqrt(gForce(ix).^2 + gForce(iy).^2),r2d(atan2(gForce(iy),gForce(ix))),'g');
   vecplot(gXY(gPJ,:),gForce(ix),0,'g');
   vecplot(gXY(gPJ,:),gForce(iy),90,'g');
else
   pp = gXY(gPJ,:);
   vecplot1(pp,ptoff(pp,1,0),pp,2,'P_X','g');
   vecplot1(pp,ptoff(pp,1,90),pp,2,'P_Y','g');
end

% roller reaction force
ix=gNL+3;
if (fig==1)
   vecplot(gXY(gRJ,:),gForce(ix),gRA,'g');
else
   pp = gXY(gRJ,:);
   vecplot1(pp,ptoff(pp,1,gRA),pp,2,'R','g');
end

% display the pointnames
for i=1:gNJ
   tbox(gXY(i,:),.2,.2,'w');
   text(gXY(i,1),gXY(i,2),pointname(i),'horiz','center','vert','middle');
end

return

function drawfbd
%drawfbd for each truss joint
%loop thru the xy matrix (i.e. each joint)
global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL

org = [0 0];
d=1;
for j=1:gNJ
   %loop thru the connectivity matrix (one loop per member)
   %fprintf('Joint %u\n',j);
   org = gXY(j,:);
   sumx = '\SigmaF_X = 0 = ';
   sumy = '\SigmaF_Y = 0 = ';
   

   for i=1:gNL
      k=0;
      if (gLinks(i,1) == j) k=gLinks(i,2); end;
      if (gLinks(i,2) == j) k=gLinks(i,1); end;
      if (k>0)
         fname = forcename(j,k);
         %extract the members joint coordinates
	      xym(1,:)=gXY(j,:);
         xym(2,:)=gXY(k,:);
         figure(2)
         vecplot1(xym(1,:),xym(2,:),org,1,fname,'k');
         figure(jfig(j));
         vecplot1(xym(1,:),xym(2,:),org,1,fname,'k');
         figure(3);
         aa = membang(gXY,j,k);
         text(i,j*2-1,sprintf('C(%.0f)',aa),'fontsize',8,'horiz','center');
         text(i,j*2,sprintf('S(%.0f)',aa),'fontsize',8,'horiz','center');
         sumx = sprintf('%s + %scos(%.0f)',sumx,fname,aa);
         sumy = sprintf('%s + %ssin(%.0f)',sumy,fname,aa);
      end
   end
   
   % display external loads
   org = gXY(j,:) - .1;
   for i=1:size(gExtLoad,1)
      if (gExtLoad(i,1)==j)
         xym(1,:)=gXY(j);
         aa = gExtLoad(i,3);
         xym(2,:) = ptoff(xym(1,:),gExtLoad(i,2),aa);
         fname = sprintf('%.0f',gExtLoad(i,2));
         figure(2);
        	vecplot1(xym(1,:),xym(2,:),org,2,fname,'m');
         figure(jfig(j));
        	vecplot1(xym(1,:),xym(2,:),org,2,fname,'m');
         sumx = sprintf('%s + %scos(%.0f)',sumx,fname,aa);
         sumy = sprintf('%s + %ssin(%.0f)',sumy,fname,aa);
         figure(3);
         text(gNL+7,j*2-1,sprintf('%sC(%.0f)',fname,aa),'fontsize',8,'horiz','center');
         text(gNL+7,j*2,sprintf('%sS(%.0f)',fname,aa),'fontsize',8,'horiz','center');
      end
   end
   
   % display pin joint reactions
   org = gXY(j,:) - .1;
   if (j == gPJ)
      % X direction
      xym(1,:) = gXY(j);
      xym(2,:) = ptoff(xym(1,:),1.3,0);
      fname = 'P_X';
      aa = 0;
      figure(2);
   	vecplot1(xym(1,:),xym(2,:),org,2,fname,'g');
      figure(jfig(j));
   	vecplot1(xym(1,:),xym(2,:),org,2,fname,'g');
      sumx = sprintf('%s + %scos(%.0f)',sumx,fname,aa);
      sumy = sprintf('%s + %ssin(%.0f)',sumy,fname,aa);
      figure(3);
      text(gNL+1,j*2-1,sprintf('C(%.0f)',aa),'fontsize',8,'horiz','center');
      text(gNL+1,j*2,sprintf('S(%.0f)',aa),'fontsize',8,'horiz','center');

      % Y direction
      xym(1,:) = gXY(j);
      xym(2,:) = ptoff(xym(1,:),2,90);
      fname = 'P_Y';
      aa = 90;
    	figure(2);
		vecplot1(xym(1,:),xym(2,:),org,2,fname,'g');
      figure(jfig(j));
   	vecplot1(xym(1,:),xym(2,:),org,2,fname,'g');
      sumx = sprintf('%s + %scos(%.0f)',sumx,fname,aa);
      sumy = sprintf('%s + %ssin(%.0f)',sumy,fname,aa);
      figure(3);
      text(gNL+2,j*2-1,sprintf('C(%.0f)',aa),'fontsize',8,'horiz','center');
      text(gNL+2,j*2,sprintf('S(%.0f)',aa),'fontsize',8,'horiz','center');
   end

   % display roller joint reactions
   if (j == gRJ)
      xym(1,:) = gXY(j);
      xym(2,:) = ptoff(xym(1,:),2,gRA);
      fname = 'R';
      aa = gRA;
    	figure(2);
      vecplot1(xym(1,:),xym(2,:),org,2,fname,'g');
      figure(jfig(j));
      vecplot1(xym(1,:),xym(2,:),org,2,fname,'g');
      sumx = sprintf('%s + %scos(%.0f)',sumx,fname,aa);
      sumy = sprintf('%s + %ssin(%.0f)',sumy,fname,aa);

      figure(3);
      text(gNL+3,j*2-1,sprintf('C(%.0f)',aa),'fontsize',8,'horiz','center');
      text(gNL+3,j*2,sprintf('S(%.0f)',aa),'fontsize',8,'horiz','center');
       
   end
    
   % draw the point name
   org = gXY(j,:);
   figure(2);
   tbox(org,.2,.2,'w');
   text(org(1),org(2),pointname(j),'horiz','center','vert','middle');

   figure(jfig(j));
   tbox(org,.2,.2,'w');
   text(org(1),org(2),pointname(j),'horiz','center','vert','middle');

   % draw the equilibrium equations
   figure(jfig(j));
   text(org(1),org(2)-3,strvcat(sumx,sumy),'HorizontalAlignment','center');
   axis([org(1)-5,org(1)+5,org(2)-4,org(2)+2]);
   
end

figure(3);
for i=1:gNL
    text(gNL+5,i,forcename(gLinks(i,1),gLinks(i,2)),'fontsize',8,'horiz','center');
end
text(gNL+5,gNL+1,'P_X','fontsize',8,'horiz','center');
text(gNL+5,gNL+2,'P_Y','fontsize',8,'horiz','center');
text(gNL+5,gNL+3,'R','fontsize',8,'horiz','center');
 
return

%=====================================================================================%
function initfig(n)
global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL
%initialize the base figure
axlim = [ min(gXY(:,1))-4 max(gXY(:,1))+4 min(gXY(:,2))-4 max(gXY(:,2))+4 ] ;
figure(4);
clf reset;
set(gcf,'position',[100 100 600 600]);
axis equal;
axis ( axlim );
set(gcf,'Name','System FBD','NumberTitle','off');
hold on;


figure(2);
clf reset;
axis equal;
axis ( axlim );
set(gcf,'Name','All FBDs','NumberTitle','off');
hold on;

for j=jfig(1):jfig(n)
   figure(j);
   clf reset;
   axis equal;
   set(gcf,'Name',['FBD ' j] ,'NumberTitle','off');
   hold on;
end

figure(3);
clf reset;
axis equal;
set(gca,'xticklabelmode','manual');
set(gca,'yticklabelmode','manual');
set(gca,'xtickmode','manual');
set(gca,'ytickmode','manual');
set(gca,'XAxisLocation','top');
axis([0 2*n+5 0 2*n+1])
set(gcf,'Name','Matrices','NumberTitle','off');
xlabel('Coefficients');
ylabel('Equations');
set(gca,'YDir','reverse');
hold on;
x1=[.5 .25 .25 .5];
x2=[.5 .75 .75 .5];
y=[.5 .5 2*n+.5 2*n+.5];
plot(x1,y,'-');
plot(x2+2*n,y,'-');
plot(x1+1+2*n,y,'-');
plot(x2+2+2*n,y,'-');
plot(x1+3+2*n,y,'-');
plot(x2+4+2*n,y,'-');
text(2*n+1,n+.5,'*','fontsize',12,'fontweight','bold','horiz','center','vert','middle');
text(2*n+3,n+.5,'+','fontsize',12,'fontweight','bold','horiz','center','vert','middle');
for (i=1:n)
   ylabs{2*i-1}=sprintf('%cx',i+64);
   ylabs{2*i}=  sprintf('%cy',i+64);
end
set(gca,'yticklabel',ylabs);
nn = 2*n;
set(gca,'ytick',1:nn);

for (i=1:2*n)
   xlabs{i}=forcename1(i);
end
xlabs{nn+1} = 'Unknowns';
xlabs{nn+2} = 'Ext Loads';
set(gca,'xticklabel',xlabs);
set(gca,'xtick',[1:nn nn+2 nn+4]);

figure(1);
clf reset;
set(gcf,'position',[100 100 600 600]);
axis equal;
axis ( axlim );
set(gcf,'Name','Results','NumberTitle','off');
hold on;

return

%=====================================================================================%
function buildcoef
%build truss coefficient matrix
global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL

gNL = size(gLinks,1);
gNJ = size(gXY,1);
gNE = gNL + 3;

clen = gNL + 3;
gCoef = zeros(clen,clen);

%add coefficients for all members
for i=1:1:clen-3
    j1=gLinks(i,1);
    j2=gLinks(i,2);
    addmember(i,j1,j2,membang(gXY,j1,j2));
end

%add coefficients for all reactions
addreaction(clen-2,gPJ,0);
addreaction(clen-1,gPJ,90);
addreaction(clen, gRJ, gRA);

return

%=====================================================================================%
function addreaction(fn,nn,ang)
%adds a reaction to the matrix used to solve a truss problem
%mn - member number
%nn - node number
%ang - angle  from n1 to n2, measure relative to +x axis

global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL

if (nn<=0 | nn>size(gCoef,1)), return; end;
row = nn*2 - 1;
gCoef(row,  fn) = cos(d2r(ang));
gCoef(row+1,fn) = sin(d2r(ang));
return

%=====================================================================================%
function addmember(mn,n1,n2,ang)
%adds a member to the matrix used to solve a truss problem
%mn - member number
%n1 - first node number
%n2 - second node number
%ang - angle  from n1 to n2, measure relative to +x axis

global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL

addreaction(mn,n1,ang);
addreaction(mn,n2,ang+180);
return

%=====================================================================================%
function addload(load)
%adds a load to the matrix used to solve a truss problem
%load - [ joint #, load magnitude, load angle ]

global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL

nn = load(1);
ff = load(2);
ang = load(3);

if (nn<=0 | nn>length(gLoad)), return; end;

row = nn*2 - 1;
gLoad(row) = ff*cos(d2r(ang));
gLoad(row+1) = ff*sin(d2r(ang));
return

%=====================================================================================%
function r1=jfig(j)
r1 = j+'A'-1;
return

%=====================================================================================%
function r1 = membang(xy,j1,j2)
% calculate a member's angle, based upon joint coordinates
r1 = r2d(atan2(xy(j2,2)-xy(j1,2),xy(j2,1)-xy(j1,1)));
return

%=====================================================================================%
function wrsTrussGeom1(cn)
global gXY gLinks gCoef gLoad gForce gPJ gRJ gRA gExtLoad gNE gNJ gNL

switch cn
case '98'
    t = d2r(60);
    h = 3;
    d2 = h / tan(t);
    gXY = [0 0
            d2 h
            2*d2 0
            3*d2 h
            4*d2 0
            5*d2 h
            6*d2 0
            7*d2 h
            8*d2 0
            -d2 h
            -2*d2 0 ];
     gLinks = [1 2
                2 3
                3 4
                4 5
                5 6
                6 7
                7 5
                3 1
                2 4
                4 6
                7 8
                8 9
                7 9
                1 10
                1 11
                10 11
                2 10
                6 8
                3 5 ];
     gPJ = 11;
     gRJ = 9;
     gRA = 90;
     gExtLoad = [ 3 50 -90
                 5 50 -90 ];
     return
case '99'
    t = d2r(60);
    h = 3;
    d2 = h / tan(t);
    gXY = [0 0
            d2 h
            2*d2 0
            3*d2 h
            4*d2 0
            5*d2 h
            6*d2 0
            7*d2 h
            8*d2 0
            -d2 h
            -2*d2 0
            2*d2 -h/2
            4*d2 -h/2 ];
     gLinks = [1 2
                2 3
                3 4
                4 5
                5 6
                6 7
                7 5
                3 1
                2 4
                4 6
                7 8
                8 9
                7 9
                1 10
                1 11
                10 11
                2 10
                6 8
                1 12
                3 12
                12 13
                5 13
                7 13 ];
     gPJ = 11;
     gRJ = 9;
     gRA = 90;
     gExtLoad = [ 3 50 -90
                 5 50 -90 ];
     return
case '3'
    gXY=[0  0
          4  0
          8  0
         12  0
         16  0
         12  4
          8  6
          4  4 ];
     gLinks = [1 2;2 3;3 4;4 5;5 6;6 7;7 8;1 8;2 8;3 8;3 7;3 6;4 6];
     gPJ = 1;
     gRJ = 5;
     gRA = 90;
     gExtLoad = [ 3 5000 -90
                 4 3000 -90 ];
     return
case '4'
     gXY=[-5 15
        0 15
        5 15
       10 15
       10 10
        5 10
        0 10
        0  5
        0  0
     -5/3  5
     -10/3 10];
	
	gLinks = [1 2;2 3;3 4;4 5;5 6;6 7;7 8;8 9;9 10;10 11;11 1;1 7;7 2;2 6;6 3;3 5;7 11;7 10;8 10];
	
	gPJ = 9;
	gRJ = 5;
	gRA = 90;
	gExtLoad = [1 5000 -90
       			6 3000 -135 ];
	return

case '2'
    
gXY=[0 8
     6 8
     16 8
     6 0
     16 0 ];

gLinks = [1 2
      2 3
      3 5
      4 5
      1 4
      2 4
      3 4];
   
gPJ = 5;
gRJ = 3;
gRA = 0;
gExtLoad = [1 2000 -90
   			2 1000 -90 ];
return
case '5'
gXY=[0 0
     2 0
     4 0
     6 0
     8 0
     10 0
     12 0
     14 0
     16 0
     14 2
     12 4
     8 4
     4 4
     2 2
     6 2
     10 2];

gLinks = [1 2
      1 14
      2 3
      2 14
      3 4
      3 13
      3 14
      3 15
      4 5
      4 15
      5 6
      5 12
      5 15
      5 16
      6 7
      6 16
      7 8
      7 10
      7 11
      7 16
      8 9
      8 10
      9 10
      10 11
      11 12
      11 16
      12 13
      13 14
      13 15 ];
   
gPJ = 1;
gRJ = 9;
gRA = 90;
gExtLoad = [3 2000 -90
   			5 5000 -90 
            6 3000 -90
            7 2000 -90];         
         
return

case '1'
    
gXY=[0  sqrt(48)
      4  0
      0  -4/sqrt(3) ];
 gLinks = ['AB';'AC';'BC'];
 gLinks = gLinks - 'A' + 1;
 gPJ = 1;
 gRJ = 3;
 gRA = 0;
 gExtLoad = [ 2 400*9.81 -90 ];
 return;
 
case '6'
gXY=[0 4; 4 4; 8 4; 12 4; 16 4; 4 2; 12 2; 8 0];
gLinks = ['AB';'AF';'BC';'BF';'CD';'CG';'CF';'CH';'DE';'DG';'EG';'FH';'GH'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('E');
gRA = 90;
gExtLoad = [ jn('B') 1000 -90 ];
return;


case 'oa'
[pa,pb,pc] = getAS4input;
gXY=[0 0; pa 0; pa pa; 3*pa pa];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('D');
gRA = 0;
gExtLoad = [ jn('C') 400 pb
             jn('B') 300 -pc];
return

case 'ob'
[pa,pb,pc] = getAS4input;
gXY=[0 0; pa 0; pa pa; 3*pa 0];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('D');
gRA = 90;
gExtLoad = [ jn('C') 400 pb
             jn('B') 300 -pc];
return

case 'oc'
[pa,pb,pc] = getAS4input;
gXY=[0 0; pa 0; pa pa; 4*pa 0];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('D');
gRA = 90;
gExtLoad = [ jn('C') 400 pb
             jn('B') 500 -pc];
return

case 'od'
[pa,pb,pc] = getAS4input;
gXY=[0 0; pa 0; pa pa; 4*pa pa];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('D');
gRA = 180;
gExtLoad = [ jn('C') 400 pb
             jn('B') 600 -pc];
return

case 'oe'
[pa,pb,pc] = getAS4input;
gXY=[0 0; 2*pa 0; 2*pa pa; 3*pa pa];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('D');
gRA = 180;
gExtLoad = [ jn('C') 400 pb
             jn('B') 600 -pc];
return
case 'of'
[pa,pb,pc] = getAS4input;
gXY=[0 0; 2*pa 0; 2*pa pa; 4*pa pa];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('D');
gRA = 180;
gExtLoad = [ jn('C') 400 pb
             jn('B') 450 -pc];
return

case 'wa'
[pa,pb,pc] = getAS4input;
gXY=[0 0; pa 0; pa pa; 3*pa pa];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('B');
gRA = 90;
gExtLoad = [ jn('C') 400 pb
             jn('D') 300 -pc];
return

case 'wb'
[pa,pb,pc] = getAS4input;
gXY=[0 0; pa 0; pa pa; 3*pa 0];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('B');
gRA = 90;
gExtLoad = [ jn('C') 400 pb
             jn('D') 300 -pc];
return

case 'wc'
[pa,pb,pc] = getAS4input;
gXY=[0 0; pa 0; pa pa; 4*pa 0];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('B');
gRA = 90;
gExtLoad = [ jn('C') 400 pb
             jn('D') 500 -pc];
return

case 'wd'
[pa,pb,pc] = getAS4input;
gXY=[0 0; pa 0; pa pa; 4*pa pa];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('B');
gRA = 90;
gExtLoad = [ jn('C') 400 pb
             jn('D') 600 -pc];
return

case 'we'
[pa,pb,pc] = getAS4input;
gXY=[0 0; 2*pa 0; 2*pa pa; 4*pa pa];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('B');
gRA = 90;
gExtLoad = [ jn('C') 400 pb
             jn('D') 450 -pc];
return

case 'wf'
[pa,pb,pc] = getAS4input;
gXY=[0 0; 2*pa 0; 2*pa pa; 3*pa pa];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('B');
gRA = 90;
gExtLoad = [ jn('C') 400 pb
             jn('D') 600 -pc];
return

case 'h'
[pa,pb,pc] = getAS4input;   
gXY=[0 0; 2*pa 0; 2*pa pa; 4*pa pa/2];
gXY=gXY./4;
gLinks = ['AB';'AC';'BC';'BD';'CD'];
gLinks = jn(gLinks);
gPJ = jn('A');
gRJ = jn('B');
gRA = 90;
gExtLoad = [ jn('C') 400 pb
             jn('D') 450 -pc];
return

end

return

%%%%%%%%%%%%%%%%

function r1=jn(n)
r1 = n-'@';
return

function [a,b,c] = getAS4input()
mob = 5;
dob = 30;
a = 3*mob;
b = 4*dob;
c = 5*dob;
return

function r1=d2r(x)
%convert degrees to radians
r1=x*pi/180;
return;