cd /home/mitran/courses/MATH661/data/EiffelTower
load EiffelPoints; X=Expression1; [NX,ndims]=size(X);
load EiffelLines; L=Expression1; [NL,nseg]=size(L);
m=3*NX; disp([NX NX^2 NL m])
K=spalloc(m,m,3*NL);
for k=1:NL
  i = L(k,1); j = L(k,2);
  p = 3*(i-1); q = 3*(j-1);
  for id=1:3
    for jd=1:3
      K(p+id,q+jd) = 0;
    end;
  end;
end

for k=1:NL
  i = L(k,1); j = L(k,2);
  lij = (X(j,:)-X(i,:))'; lij=lij/norm(lij); Pij = lij*lij';
  p = 3*(i-1); q = 3*(j-1);
  for id=1:3
    for jd=1:3
      K(p+id,p+jd) = K(p+id,p+jd) - Pij(id,jd);
      K(p+id,q+jd) = K(p+id,q+jd) + Pij(id,jd);
      K(q+id,p+jd) = K(q+id,p+jd) + Pij(id,jd);
      K(q+id,q+jd) = K(q+id,q+jd) - Pij(id,jd);
    end;
  end;
end

[rows,cols,vals]=find(K);
save -v6 KCOO.mat rows cols vals