function retro(mat,k,m)
% mat = set of transition matrices, one bellow each other, without any headings 
% k=number of stages
% m - number of matrices
% calculates standard deviation of each matrix element
% and the contribution of each matrix element to observed variation in population growth rate
% the contribution is calculated as SD*sen^2

% Reads the transition matrices
fid=fopen(mat,'r');
bb=fscanf(fid,'%f',[k,k*m]);
l=0;
for i=1:m
    for j=1:k
      for p=1:k
   	    a (j,p,i)= bb(j,p+k*(i-1));
       end
   end
a(:,:,i)=transpose(a(:,:,i));   
end
fclose(fid);

var=zeros(k,k);
for i=1:k
   for j=1:k
       meana(i,j)=mean(a(i,j,:));
       if meana(i,j)>0
       for mm=1:m
          var(i,j)=var(i,j)+ ((meana(i,j)-a(i,j,mm))^2)); % calculates standart deviation of each matrix element
       end
       var(i,j)=var(i,j)/(m-1);
       else var(i,j)=0;
       end
   end
end   

for i=1:k
   for j=1:k
      b(i,j)=mean(a(i,j,:));
   end
end   
[W,d] = eig(b);
lambda = diag(d);
imax = find(lambda==max(lambda));
V=conj(inv(W));
% w = stable stage distribution
w = W(:,imax);
% v = reproductive value
v=real(V(imax,:));
% gets population growth rate - the dominant eigenvalue of the matrix
lambda1 = lambda(imax);
% calculates sensitivity of the matrix
sen = v'*w';
% calculates elasticity of the matrix
ela = sen.*b/max(eig(b));
% prints the results
fid=fopen('lambda','a');
fprintf(fid,'%10.7f\t%10.7f\n',lambda1,i);
fclose(fid);
var
x=0;
for i=1:k
   for j=1:k
      contr(i,j)=var(i,j)*(sen(i,j)*sen(i,j)); % contribution of each stage to total variation, in population growth rate
      x=x+1;
       
      slsen(x)=sen(i,j);
      slela(x)=ela(i,j);
      slvar(x)=var(i,j);
      for mm=1:m
         sla(x,mm)=a(i,j,mm);
      end   
   end
end
contr
c=0;
covariance=zeros(k*k*k*k,1);
for i=1:k*k
   for j=i:k*k
      c=c+1;
      for mm=1:m
          covariance(c)=covariance(c)+((mean(sla(i,:)))-sla(i,mm))*((mean(sla(j,:)))-sla(j,mm));
      end
   end 
end
c=0;
for i=1:k*k
   for j=i:k*k
      c=c+1;
      covmat(i,j)=covariance(c);
   end 
end
covmat;
figure;
plot(slela,slvar,'o','MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',10);
xlabel('elasticity');
ylabel('standart deviation');
figure;
plot(slsen,slvar,'o','MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',10);
xlabel('sensitivity');
ylabel('standart deviation');
ela
sen
contr