function [lm, nd]=pf_calclm(n, bc)
lm=zeros(n,3); // initialize lm to zero

[ns dummy]=size(bc); // ns = number of supports (nodes with constrained dof)
for i=1:ns
  nn=bc(i,1); // node number of ith support
  lm(nn, 1:3) = bc(i, 2:4); // constraint codes of ith support
end;

nd=0;
for i=1:n
  for j=1:3
    if lm(i,j) == 1 then // constrained dof
      lm(i,j) = 0;
    else // unconstrained dof
      nd = nd + 1;
      lm(i,j) = nd;
    end
  end
end

endfunction

function [l, cx, cy]=pf_calclen(imem, xy, conn)
n=conn(imem,:); // start and end nodes of member number i
p1 = xy(n(1),:); // x,y coordinates of start node
p2 = xy(n(2),:); // x,y coordinates of end node
dxdy=p2-p1; // x,y projections of member
dxdy2=dxdy.^2; // square of the projections
l=sqrt(sum(dxdy2)); // length of member

cx = dxdy(1) / l; // x-direction direction cosine of member
cy = dxdy(2) / l; // y-direction direction cosine of member

endfunction


function [r] = pf_calcrot(cx, cy)
r = zeros(6,6); // initialize rotation matrix to zero

r(1,1) = cx;
r(1,2) = cy;
r(2,1) = -r(1,2);
r(2,2) = r(1,1);
r(3,3) = 1;

r(4:6, 4:6) = r(1:3, 1:3); // copy rows and columns 1:3, 1:3 into 4:6, 4:6
endfunction

function [k] = pf_stiff(E, A, I, L)
k = zeros(6,6); // initialize k to zero

// Fill up the upper triangule of k
scm1 = E * A / L;
scm2 = 4.0 * E * I / L;
scm3 = 1.5 * scm2 / L;
scm4 = 2.0 * scm3 / L;

k(1,1) = scm1; k(1,4) = -k(1,1);
k(2,2) = scm4; k(2,3) = scm3;  k(2,5) = -scm4; k(2,6) = scm3;
k(3,3) = scm2; k(3,5) = -scm3; k(3,6) = scm2/2.0;
k(4,4) = scm1;
k(5,5) = scm4; k(5,6) = -scm3;
k(6,6) = scm2;

// copy off diagonal elements from upper triangle into lower triangle
for i=2:6;
  for j=1:i-1
    k(i,j) = k(j,i);
  end
end

endfunction


function [K] = pf_gstiff(imem, xy, conn, mprop)

[L cx cy] = pf_calclen(imem, xy, conn); // calculate length and direction cosines
r = pf_calcrot(cx, cy); // calculate rotation matrix

iprop = conn(imem,3); // property id of ith member
E = mprop(iprop,1); // Modulus of elasticity
A = mprop(iprop,2); // Area of cross section
I = mprop(iprop,3); // Second moment of area of cross section about NA

k = pf_stiff(E, A, I, L); // local stiffness matrix of ith member

K = r' * k * r; // global stiffness matrix of ith member
endfunction


function [ssm] = pf_assemssm(imem, xy, conn, mprop, lm, ssm)

K = pf_gstiff(imem, xy, conn, mprop);

nj = conn(imem,1); nk = conn(imem,2);
dof(1:3) = lm(nj,1:3);
dof(4:6) = lm(nk,1:3);

for i=1:6
  ii = dof(i);
  if ii == 0 then

  else
    for j=1:6
      jj = dof(j);

      if jj == 0 then

      else
        tmp = ssm(ii,jj) + K(i,j);
        ssm(ii,jj) = tmp;
      end
    end
  end
end
endfunction

function [K]=pf_ssm(xy, conn, mprop, lm, nd)
  K=zeros(nd,nd);
  [nmem dummy]=size(conn);
  for imem = 1:nmem
    K = pf_assemssm(imem, xy, conn, mprop, lm, K);
  end
endfunction


function [dof]=pf_getdof(imem, conn, lm)
dof=zeros(1,6);
n1 = conn(imem,1);
dof(1:3) = lm(n1,:);
n2 = conn(imem,2);
dof(4:6) = lm(n2,:);
endfunction


function [P]=pf_assemloadvec_jl(lm, jtloads, P)
[nloads dummy] = size(jtloads);

for i=1:nloads
  n=jtloads(i,1);
  dof=lm(n,:);

  for j=1:3
    jj = dof(j);
    if jj == 0 then
    
    else
      tmp = P(jj) + jtloads(i, j+1);
      P(jj) = tmp;
    end
  end
end
endfunction


function [P]=pf_assemloadvec_ml(iload, xy, conn, lm, memloads, P)
  imem = memloads(iload,1);
  [L cx cy]=pf_calclen(imem,xy,conn);
  r = pf_calcrot(cx, cy);

  am=-r' * memloads(iload,2:7)';
  dof=pf_getdof(imem, conn, lm);

  for i=1:6
    ii = dof(i);

    if ii ~= 0 then
      tmp = P(ii) + am(i);
      P(ii) = tmp;
    end
  end

endfunction


function [ssm,x,P]=pf(xy,conn,bc,mprop,jtloads,memloads)

[nodes dummy]=size(xy);

[lm ndof]=pf_calclm(nodes,bc);
x=zeros(ndof,1);
p=zeros(ndof,1);
ssm=pf_ssm(xy,conn,mprop,lm,ndof);

P=zeros(ndof,1);
[njtlds dummy]=size(jtloads);

P=pf_assemloadvec_jl(lm, jtloads, P);

[nmemlds dummy]=size(memloads);
for iload=1:nmemlds
  [P]=pf_assemloadvec_ml(iload, xy, conn, lm, memloads, P);
end

x=zeros(ndof,1);
x=inv(ssm)*P;
endfunction


function [f] = pf_memendforces(imem,xy,conn,mprop,lm,x,memloads)
iprop=conn(imem,3);
E=mprop(iprop,1); A=mprop(iprop,2); I=mprop(iprop,3);
[L cx cy]=pf_calclen(imem,xy,conn);
r=pf_calcrot(cx,cy);

k = pf_stiff(E,A,I,L);

u = zeros(6,1);
dof = pf_getdof(imem, conn, lm);
for i = 1:6
  idof = dof(i);
  if idof ~= 0 then
    u(i) = x(idof);
  end
end
uu=r*u;

f = zeros(6,1);
f = k * uu;

[nmemloads,dummy]=size(memloads);
for i=1:nmemloads
  if memloads(i,1) == imem
    f = f + memloads(i, 2:$)';
  end
end

endfunction