% Script for calculating population growth rate for a set of simulations when a given element is reduced or increased
% Can be used to simulate herbivory, variation in seedling germination and so on
function herbivory(mat,stg,disp,dbg,runs,m,basmat) 
% mat - transition matrix, or multiple matrices one below each other
% stg - number of stages
% disp - transitions to be reduced
% dbg - reduction coefficients, there can be more differetn coefficients arranged in a column 
% runs - number of values in dbg file
% m - number of matrices
% basmat - add values to transitions after multiplication, contains as many matrices as file mat

% reads transition matrix
fid=fopen(mat,'r');
aaa=fscanf(fid,'%f',[stg,stg*m]);
aaa=transpose(aaa);
fclose(fid);
aaa

for i=1:m
    for j=1:stg   
       ba(j,:,i)=aaa((j+(i-1)*stg),:);
    end
 end
ba

%reads basmat
fid=fopen(basmat,'r');
hhh=fscanf(fid,'%f',[stg,stg*m]);
hhh=transpose(hhh);
fclose(fid);
for i=1:m
    for j=1:stg   
       ha(j,:,i)=hhh((j+(i-1)*stg),:);
    end
 end

%reads dispersal matrix:
fid1=fopen(disp,'r');
b=fscanf(fid1,'%f',[stg,stg]);
b=transpose(b);
fclose(fid1);

%reads dbg number
fid2=fopen(dbg,'r');
c=fscanf(fid2,'%f',[runs]);
c=transpose(c);
fclose(fid2);
for rr=1:m
    a=ba(:,:,rr);
    h=ha(:,:,rr);
    for k = 1:runs % for each dbg value
	     for i=1:stg
   	      for j=1:stg
                if b(i,j)==1
                   aa(i,j)= a(i,j)*c(k);
                else 
                   aa(i,j)=a(i,j);
                end
            end
        end
        aa=aa+h; % add basmat
        [W,d] = eig(aa);
	     lambda = diag(d);
	     imax = find(lambda==max(lambda));
	     V=conj(inv(W));
	     w = W(:,imax);
	     v=real(V(imax,:));
	     lambda1 = lambda(imax);
	     sen = v'*w';
        ela = sen.*aa/max(eig(aa));
        fid4=fopen('lambda.out','a');
	     fprintf(fid4,'%10.9f\t%10.9f\t%10.9f\n',c(k),lambda1,rr);
	     fclose(fid4);
        fid4=fopen('sen.out','a');
        fid5=fopen('ela.out','a');
        fid6=fopen('mat.out','a');
        for i=1:stg
            for j=1:stg
	             fprintf(fid4,'%10.7f \t',sen(i,j));
                fprintf(fid5,'%10.7f \t',ela(i,j));
                fprintf(fid6,'%10.7f \t',aa(i,j));
          		 if j==stg
      	          fprintf(fid4,'\n');
                   fprintf(fid5,'\n');
                   fprintf(fid6,'\n');
                end
			   end
        end
        fclose(fid4);
        fclose(fid5);
	     fclose(fid6);
     end % end of runs
 end %end of matrices   