%POMIKI vrne vektor pomikov v vozliscih
%dX=POMIKI(datoteka) vrne vektor pomikov v tockah
%za podatke opisane v datoteki 'datoteka'
%'datoteka mora imeti naslednjo stukturo (po vrsticah)
%n - stevilo vozlisc
%m - stevilo palic
%d - dimenzija problema (2,3)
%stpod - stevilo podpor
%stsil - stevilo sil, ki obremenjujejo palicje
%xi yi - n vrstic, ki opisujejo koordinate vozlisc
%zi ki Ei Si - zacetek i-te palice, konec, Youngov modul, presek
%T0 - indeksi podpornih tock
%f_Ti fxi fyi fzi - indeks, kjer deluje i-ta sila, komponente sile

function dX=pomiki(datoteka);
  %preberemo podatke
  % branje podatkov z datoteke "datoteka"

[fid,str]=fopen(datoteka,'rt'); %odpre datoteko za branje teksta
if fid==-1 
	error(str); 	%ce je napaka pri odpiranju datoteke, koncaj
			%in izpisi sporocilo o napaki
end

%bere ustrezne podatke
n=fscanf(fid,'%d',1); %st. vozlisc
m=fscanf(fid,'%d',1); %st. palic
d=fscanf(fid,'%d',1); % dimenzija
stpod=fscanf(fid,'%d',1); %stevilo podpor
stsil=fscanf(fid,'%d',1); %stevilo sil
N=n*d;	%dolzina vektorja pomikov
T=zeros(d,n); %matrika tock
T=fscanf(fid,'%f',[d,n]); % prebere koordinate vozlisc v T

P=zeros(2,m);
E=zeros(1,m);
S=zeros(1,m);
for j=1:m
	P(:,j)=fscanf(fid,'%d',[2,1]);
	%I(j,1)=fscanf(fid,'%d',1);
	%I(j,2)=fscanf(fid,'%d',1);
	E(j)=fscanf(fid,'%f',1);
	S(j)=fscanf(fid,'%f',1);
end

% prebere podpore v vektor podpor
TP=fscanf(fid,'%d',[1,stpod]);

%prebere obremenitvene sile
f=zeros(d,n);
for j=1:stsil
    indeks=fscanf(fid,'%d',1); %indeks tocke, kje je sila
    f(:,indeks)=fscanf(fid,'%f',[d,1]); %koordinate sile
end

fclose(fid); %zapre datoteko

%konec branja

%izracuna dolzine palic
  m=size(P,2);
  L=zeros(1,m);
  for j=1:m
    L(j)=norm(T(:,P(2,j))-T(:,P(1,j)),2);
  end
  %izracuna matriko smernih kosinusov
  C=zeros(d,m);
  for j=1:m
      C(:,j)=1/L(j)*[T(:,P(2,j))-T(:,P(1,j))];
  end
  x=zeros(N,1); %matriko pomikov raztegne v vektor
  %sestavi incidencno matriko A po formuli na str 115 Mat. mod.
  A=zeros(m,N);
  if d==2
  for i=1:m
    A(i,:)=([enka(P(2,i)*d-1,N)-enka(P(1,i)*d-1,N),...
             enka(P(2,i)*d,N)-enka(P(1,i)*d,N)]*C(:,i))';
  end
  else
    for i=1:m
      A(i,:)=([enka(P(2,i)*d-2,N)-enka(P(1,i)*d-2,N),...
               enka(P(2,i)*d-1,N)-enka(P(1,i)*d-1,N),...
               enka(P(2,i)*d,N)-enka(P(1,i)*d,N)]*C(:,i))';
    end
  end
  %izracuna vektor konstant vzmeti palic
  k=zeros(1,m);
  k=E.*S./L;
  D=diag(sqrt(k));
  %reducira matriko A in vektor f
  %torej v A "crta stolpce, katerih indeksi ustrezajo
  %pritrjenim komponentan, v f pa ustrezne vrstice
  TN=setdiff([1:n],TP); %indeksi nepodprtih tock
  if d==2
     stolpci=union(2*TN-1,2*TN); %stolpci, ki jih ohranimo
  else
     stolpci1=union(3*TN-2,3*TN-1);
     stolpci=union(stolpci1,3*TN); %stolpci, ki jih ohranimo
  end 
  A0=A(:,stolpci);
  f0=f(:);
  f0=f0(stolpci);
  [Q,U]=qr(D*A0);
  x0=(U'*U)\f0;
  x0=reshape(x0,d,length(TN));
  dX=zeros(d,n);
  dX(:,TN)=x0;
  T_nove=T+dX;
  clf
  hold on
    for i=1:m
      indeksi=[P(1,i),P(2,i)];
      if d==2
        plot(T(1,indeksi),T(2,indeksi),'b');
        plot(T_nove(1,indeksi),T_nove(2,indeksi),'r');
      else
        plot3(T(1,indeksi),T(2,indeksi),T(3,indeksi),'b');
        plot3(T_nove(1,indeksi),T_nove(2,indeksi),T_nove(3,indeksi),'r');
      end
    end
    axis equal;
    hold off;
    

%ENKA vrne enotski vektor
%E=enka(I,N);
%E=[0,0,...,0,1,0,...]' kjer je 1 na I-tem mestu, dolzina
%vektorja pa je N
function e=enka(i,n);
  e=zeros(n,1);
  e(i)=1;
