kern.m

% Example: kern
% ~~~~~~~~~~~~~
% This program determines the kern for an 
% axially loaded compression member described
% by a polygonal cross section.
%
% Data is defined in the declaration statements
% below, where:
%
% x,y - vectors containing the coordinates for
%       the polygon nodes.  CCW fashion for
%       positive contributions, CW fashion
%       for negative contributions.
%
% User m functions required:
%    prop, shftprop, inout, intrsect, clip,
%    refpoly, circle, genprint
%----------------------------------------------

clear;
%...Input definitions
Problem=1;
if Problem == 1
  %...Symmetric angle
  x=[0 10 10 2 2 0]; y=[0 0 2 2 10 10];
elseif Problem == 2
  %...Z section
  x=[-1 7 7 5 5 1 1 -7 -7 -5 -5 -1];
  y=[-10 -10 -6 -6 -8 -8 10 10 6 6 8 8];
elseif Problem == 3
  %...Square with triangle cutout
  x=[0 10 10 0 0 5 2 8 5];
  y=[0 0 10 10 0 2 8 8 2];
elseif Problem == 4
  %...Two disconnected triangles
  x=[5 10 10 5 -10 -5 -10 -10];
  y=[0 -6 6 0 -6 0 6 -6];
elseif Problem == 5
  %...Channel 
  x=[0 10 10 9 9 1 1 0];
  y=[0 0 10 10 1 1 10 10];
else
  %...3/4 circular section
  [x,y]=circle(36,0,0,10);
  x=x(1:28); x(29)=0; y=y(1:28); y(29)=0;
end

%...Initialize
Epsilon=1e-8; NoPts=length(x);

%...Determine geometrical properties
[A,Qx,Qy,Ix,Iy,Ixy,xbar,ybar]=prop(x,y);

%...Shift to centroidal axis
[Ixbar,Iybar,Ixybar]= ...
  shftprop(A,xbar,ybar,Ix,Iy,Ixy);
xs=x-xbar; ys=y-ybar;

%...Define polygon to clip against
[NoClipPts,Xclip,Yclip]=refpoly(xs,ys);

%.......................................
%...Loop on each corner of user geometry
%...and clip reference polygon
%.......................................
Determinant=Ixbar*Iybar-Ixybar^2; alpha=1/A;
for i=1:NoPts
  beta=(Ixbar*xs(i)-Ixybar*ys(i))/Determinant;
  gamma=(Iybar*ys(i)-Ixybar*xs(i))/Determinant;
  sfact=sign(gamma);
  if sfact == 0, sfact=1; end
  d=sqrt(beta^2+gamma^2)*sfact;

  %...Note: The neutral axis for any candidate
  %...      point has the form:
  %...      
  %...      alpha + beta*x + gamma*y = 0
  %...      
  %...      The normal form is:
  %...      
  %...      alpha   beta*x   gamma*y
  %...      ----- + ------ + ------- = 0
  %...        d        d        d    
  %...      
  %...      where:
  %...      
  %...      d=SIGN(beta)*SQRT(beta^2+gamma^2)

  Nform(1)=alpha/d;  % perpendicular distance
  Nform(2)=beta/d;   % run of normal vector
  Nform(3)=gamma/d;  % rise of normal vector

  %...Make sure line is defined such that
  %...  a perpendicular (the normal vector for 
  %...  the clipping line) is directed away from
  %...  clipping zone
  DotProduct=xs(i)*Nform(2)+ys(i)*Nform(3);
  if DotProduct >= 0, Nform=-Nform; end

  %...Clip the reference polygon
  [NoNewPts,IsetClipPoly,NoInterpPts, ...
    Xclip,Yclip]=clip(Xclip,Yclip,NoClipPts, ...
    Nform,Epsilon);

  %...Clip complete, reorder arrays
  Loop=NoClipPts+NoInterpPts;
  xtmp(1:Loop)=Xclip(1:Loop);
  ytmp(1:Loop)=Yclip(1:Loop);
  clear Xclip Yclip;
  for j=1:NoNewPts
    itmp=IsetClipPoly(j);
    Xclip(j)=xtmp(itmp); Yclip(j)=ytmp(itmp);
  end
  NoClipPts=NoNewPts;
end
  
%...Convert back to original system
Xclip=Xclip+xbar; Yclip=Yclip+ybar;

fprintf( ...
  '\n\nKern of a Polygonal Cross Section');
fprintf('\n---------------------------------');
fprintf( ...
  '\n\nTable of Nodes Forming Cross Section');
fprintf( ...
    '\n\n  Node #         x             y');
for i=1:NoPts
  fprintf('\n  %3.0f %12.5g  %12.5g', ...
          i,x(i),y(i));
end 
fprintf('\n\nGeometrical Properties');
fprintf(  '\n  A:     %g',A);
fprintf(  '\n  Qx:    %g',Qx);
fprintf(  '\n  Qy:    %g',Qy);
fprintf(  '\n  Ix:    %g',Ix);
fprintf(  '\n  Iy:    %g',Iy);
fprintf(  '\n  Ixy:   %g',Ixy);
fprintf(  '\n  x-bar: %g',xbar);
fprintf(  '\n  y-bar: %g',ybar);
fprintf('\n\nGeometrical Properties About ');
fprintf(    'Centroidal Axis');
fprintf(  '\n  Ix-bar:  %g',Ixbar);
fprintf(  '\n  Iy-bar:  %g',Iybar);
fprintf(  '\n  Ixy-bar: %g',Ixybar);
fprintf('\n\nTable of Kern Boundary Nodes');
fprintf('\n\n  Node #         x             y');
for i=1:NoClipPts
  fprintf('\n  %3.0f %12.5g  %12.5g', ...
          i,Xclip(i),Yclip(i));
end 
fprintf('\n');

%...Close end of polygon for plotting
xp=[x(:);x(1)]; yp=[y(:);y(1)];
clf;
plot(xp,yp,'-',xp,yp,'o'); hold on; 
  fill(Xclip,Yclip,'y'); 
  plot(Xclip,Yclip,'o',Xclip,Yclip,'-');
  title('Kern of Cross Section');
  xlabel('x'); ylabel('y');
  axis('equal'); hold off; drawnow;
% genprint('kern');