// input
double: x1,x2,x3,x4, y1,y2,y3,y4; // vertex coordinates
double: t;                        // thickness
double: E, G;                     // Young's modulus and shear modulus, G:=E/2/(1+nu)

// output
double: K[1..4,1..4];             // stiffness matrix

// local variables
double: c1,c2,c3,c4, s1,s2,s3,s4, r1,r2,r3,r4, k1,k2,k3,k4, l1,l2,l3,l4;
double: t1,t2,t3,t4, g1,g2,g3,g4, j1,j2,j3,j4, p1,p2,p3,p4, h1,h2,h3,h4;
double: A, k, B[1..4], s, D;
integer: i, j;

c1:=x2-x1; c2:=x3-x2; c3:=x4-x3; c4:=x1-x4;
s1:=y2-y1; s2:=y3-y2; s3:=y4-y3; s4:=y1-y4;
r1:=x1*y2-x2*y1; r2:=x2*y3-x3*y2; r3:=x3*y4-x4*y3; r4:=x4*y1-x1*y4;
k1:=c2*(s3*r4-s4*r3)-s2*(c3*r4-c4*r3)+r2*(c3*s4-c4*s3);
k2:=c1*(s3*r4-s4*r3)-s1*(c3*r4-c4*r3)+r1*(c3*s4-c4*s3);
k3:=c1*(s2*r4-s4*r2)-s1*(c2*r4-c4*r2)+r1*(c2*s4-c4*s2);
k4:=c1*(s2*r3-s3*r2)-s1*(c2*r3-c3*r2)+r1*(c2*s3-c3*s2);
if(G>0.0)and(E>0.0)and(k1>0.0)and(k2>0.0)and(k3>0.0)and(k4>0.0)then
begin
  l1:=sqrt(c1*c1+s1*s1); l2:=sqrt(c2*c2+s2*s2);
  l3:=sqrt(c3*c3+s3*s3); l4:=sqrt(c4*c4+s4*s4);
  A:=0.50*(r1+r2+r3+r4);
  k:=0.25*(k1+k2+k3+k4);
  B[1]:=-k1/k*l1; B[2]:=k2/k*l2;
  B[3]:=-k3/k*l3; B[4]:=k4/k*l4;
  t1:=c4-c2; t2:=s4-s2; t3:=c3-c1; t4:=s3-s1;
  g1:=(c1*t1+s1*t2)/(s1*t1-c1*t2); 
  g2:=(t3*c2+t4*s2)/(t4*c2-t3*s2);
  g3:=(c3*t1+s3*t2)/(s3*t1-c3*t2);
  g4:=(t3*c4+t4*s4)/(t4*c4-t3*s4);
  t1:=0.5*(c1*s3-c3*s1); t2:=0.5*(c2*s4-c4*s2);
  j1:=A+t2; j2:=A-t1; j3:=A-t2; j4:=A+t1;
  t1:=0.5/G; t2:=2.0/E;
  p1:=t1+t2*g1*g1; p2:=t1+t2*g2*g2;
  p3:=t1+t2*g3*g3; p4:=t1+t2*g4*g4;
  h1:=k1/k; h1:=h1*h1;
  h2:=k2/k; h2:=h2*h2;
  h3:=k3/k; h3:=h3*h3;
  h4:=k4/k; h4:=h4*h4;
  s:=h1*p1*j1+h2*p2*j2+h3*p3*j3+h4*p4*j4;
  D:=2.0*t/s;
  for i:=1 to 4 do for j:=1 to 4 do K[i,j]:=B[i]*D*B[j];
end
else
  for i:=1 to 4 do for j:=1 to 4 do K[i,j]:=0.0;