clear
mm=10;
for fi=1:mm
filename1 = ['hc a30m' num2str(fi) '.csv'];
filename2=['a30m' num2str(fi) '.csv'];
m=20000;
z=10^2;%pred-prey body mass ratio
aa=dlmread(filename2)';%transposed interaction matrix where i eats J
hc=dlmread(filename1);% Prey average trophic level as calulated in R
a=aa;
n=length(hc);
Troph=hc-1;
k1=zeros(n,1);
k2=zeros(n,1);
k3=zeros(n,1);
k4=zeros(n,1);
om=zeros(n,1);
meta=zeros(n,1);
xw=zeros(1,n);
ee=zeros(n,n);
ff=zeros(n,n);
r=zeros(1,n);
c=zeros(1,n);
y=zeros(1,n);
cv=zeros(1,n);
BB=zeros(n,mm);

for i=1:n%calculating growth rate (r)
if hc(i)<=1.00
r(i)=1;
else
r(i)=0;
end
end


basalPosition=find(hc<=1.0); %Finds all basal species (i.e. potential foundation species
consumerPosition=find(hc>=1.0)
twobase = basalPosition(randsample(length(basalPosition),2));%random selection of basal species for foundation species (fSP)
fspPosition=twobase(1);
extPos=twobase(2);

links=cumsum(a(:,fspPosition));
fspLinks=links(n);

blinks=cumsum(a(:,extPos));
extLinks=blinks(n);

for i=1:n%claulating body mass
Mass(i)=z^Troph(i);
end

for i=1:n%calculating metabolic rate (x)
if Troph(i)<= 0
xw(i)=.138;
else
xw(i)=0.314*Mass(i)^-.25;
meta=xw';
Mnti=meta/10;
end
end


for i=1:n%determining if herbivore or predator for assimilation efficiancy (eij)
if hc(i) == 2
ee(i,1:n)=.45;
else
ee(i,1:n)=.85;
end
end

K=1;% Carrying Capacity
b1=.5;% Half saturation constant in the absence of nti
b=1;% Hill # for functional response

for i=1:n
y(i)=8;% Max consumption rate
c(i)=0;% Predator interference
end
for i=1:n%relative consumption rate Wij in Brose et al 2006.
sum=0;
for j=1:n
sum=sum+a(i,j);
end
if sum==0 ;
om(i)=0;
else
om(i)=1/sum;
end
end
for i=1:n %conversion efficiency
for j=1:n
ff(i,j)=1.0;
end
end

Initial_Biomass=((1.0-.05)/(1-0))*rand(n,1)+.05; % uniformly distributed biomass from .05 to 1.
%B(1)=.9456; B(2)=0; B(3)=.00; B(4)=.765; B(5)=.3414; B(6)=.5925; B(7)=.5925; B(8)=.9905; B(9)=.5153; B(10)=.0970;

B=Initial_Biomass;

dt=.001;%Time step 
ii=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%_No_Foundation_Species_Run_%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:m
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
t(i)=dt*i;
   k1=dt*f1(B,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k1/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k2=dt*f1(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k2/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k3=dt*f1(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k3;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k4=dt*f1(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
   BN = B + (k1 + 2*k2 + 2*k3 + k4)/6;
   t(i) = i*dt;
   B=BN;
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
if mod(i,1000) == 0
ii=ii+1
BB(:,ii)=B;%what does : mean? All
tt(ii)=i*dt;
B
i
end
end

%%%%%%%%%%%%%_Output for saving%%%%%%%%%%%%%%
for i=1:n %interaction matrix without extinc species%
if B(i)==0;
a(i,:)=0;
a(:,i)=0;
end
end


for i=1:n
Col=cumsum(a,1);
Row=cumsum(a,2);
pp=Row(:,n)';
zz=(Col(n,:));
qq=pp+zz;
if qq(i)==0;
B(i)=0;
else
B(i)=B(i);
end
end

Final=B;%final species list

for i=1:n%Calculating the -cv for population stability
if B(i)>0;
cv(i)=-std(wkeep(BB(i,:),(m*dt)/2,'r'))/mean(wkeep(BB(i,:),(m*dt)/2,'r'));
else
cv(i)=NaN;
end
end
cvAvg=nanmean(cv);
for i=1:n %turning densites in presence absence
if B(i)==0;
yb(i)=0;
else
yb(i)=1;
end
end

ich=cumsum(yb);%calculating species richness
rich=ich(n);
sprich=[rich,cvAvg,fspPosition,fspLinks];
newMatrix=a;

filename1=['Pops_nofsp_' num2str(n) '_' num2str(fi) '.csv'];
filename2=['cv_nofsp_' num2str(n) '_' num2str(fi) '.csv'];
filename3=['Matrix_nofsp_' num2str(n) '_' num2str(fi) '.csv'];
filename4=['fPop_nofsp_' num2str(n) '_' num2str(fi) '.csv'];
dlmwrite(filename1,BB, ',')
dlmwrite(filename2,sprich,',')
dlmwrite(filename3,a,',')
dlmwrite(filename4,Final,',')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%_No_Foundation_Species_Run_Extinction of FSP%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ii=0;
B=Final;
B(fspPosition)=0;
a=newMatrix;

for i = 1:m
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
t(i)=dt*i;
   k1=dt*f1(B,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k1/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k2=dt*f1(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k2/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k3=dt*f1(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k3;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k4=dt*f1(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
   BN = B + (k1 + 2*k2 + 2*k3 + k4)/6;
   t(i) = i*dt;
   B=BN;
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
if mod(i,1000) == 0
ii=ii+1
BB(:,ii)=B;%what does : mean? All
tt(ii)=i*dt;
B
i
end
end

%%%%%%%%%%%%%_Output for saving%%%%%%%%%%%%%%
for i=1:n %interaction matrix without extinc species%
if B(i)==0;
a(i,:)=0;
a(:,i)=0;
end
end

links=cumsum(a(:,fspPosition));
xfspLinks=links(n);

blinks=cumsum(a(:,extPos));
xextLinks=blinks(n);


for i=1:n
Col=cumsum(a,1);
Row=cumsum(a,2);
pp=Row(:,n)';
zz=(Col(n,:));
qq=pp+zz;
if qq(i)==0;
B(i)=0;
else
B(i)=B(i);
end
end

Finalx=B;%final species list

for i=1:n%Calculating the -cv for population stability
if B(i)>0;
cv(i)=-std(wkeep(BB(i,:),(m*dt)/2,'r'))/mean(wkeep(BB(i,:),(m*dt)/2,'r'));
else
cv(i)=NaN;
end
end
cvAvg=nanmean(cv);
for i=1:n %turning densites in presence absence
if B(i)==0;
yb(i)=0;
else
yb(i)=1;
end
end

ich=cumsum(yb);%calculating species richness
rich=ich(n);
sprich=[rich,cvAvg,fspPosition,xfspLinks];

filename1=['Pops_nofspx_' num2str(n) '_' num2str(fi) '.csv'];
filename2=['cv_nofspx_' num2str(n) '_' num2str(fi) '.csv'];
filename3=['Matrix_nofspx_' num2str(n) '_' num2str(fi) '.csv'];
filename4=['fPop_nofspx_' num2str(n) '_' num2str(fi) '.csv'];
dlmwrite(filename1,BB, ',')
dlmwrite(filename2,sprich,',')
dlmwrite(filename3,a,',')
dlmwrite(filename4,Finalx,',')




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%_Basal_Run_%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=Initial_Biomass;
ii=0;
a=aa;

for i=1:n%calculating metabolic rate (x)
if Troph(i)<= 0
xw(i)=.138;
else
xw(i)=0.314*Mass(i)^-.25;
meta=xw';
metaC=meta/10;
Mnti=meta;
end
end
for i=1:length(basalPosition)
Mnti(basalPosition(i))=metaC(i);
end


for i = 1:m
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
t(i)=dt*i;
   k1=dt*flow(B,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k1/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k2=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k2/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k3=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k3;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k4=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
   BN = B + (k1 + 2*k2 + 2*k3 + k4)/6;
   t(i) = i*dt;
   B=BN;
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
if mod(i,1000) == 0
ii=ii+1
BB(:,ii)=B;
tt(ii)=i*dt;
B
i
end
end

%%%%%%%%%%%%%_Output for saving%%%%%%%%%%%%%%
for i=1:n %interaction matrix without extinc species%
if B(i)==0;
a(i,:)=0;
a(:,i)=0;
end
end

for i=1:n
Col=cumsum(a,1);
Row=cumsum(a,2);
pp=Row(:,n)';
zz=(Col(n,:));
qq=pp+zz;
if qq(i)==0;
B(i)=0;
else
B(i)=B(i);
end
end

Final=B;%final species list

for i=1:n%Calculating the -cv for population stability
if B(i)>0;
cv(i)=-std(wkeep(BB(i,:),(m*dt)/2,'r'))/mean(wkeep(BB(i,:),(m*dt)/2,'r'));
else
cv(i)=NaN;
end
end
cvAvg=nanmean(cv);
for i=1:n %turning densites in presence absence
if B(i)==0;
yb(i)=0;
else
yb(i)=1;
end
end

ich=cumsum(yb);%calculating species richness
rich=ich(n);
sprich=[rich,cvAvg,fspPosition,fspLinks];
newMatrix=a;


filename1=['Pops_Lower_' num2str(n) '_' num2str(fi) '.csv'];
filename2=['cv_Lower_' num2str(n) '_' num2str(fi) '.csv'];
filename3=['Matrix_Lower_' num2str(n) '_' num2str(fi) '.csv'];
filename4=['fPop_Lower_' num2str(n) '_' num2str(fi) '.csv'];
dlmwrite(filename1,BB, ',')
dlmwrite(filename2,sprich, ',')
dlmwrite(filename3,a, ',')
dlmwrite(filename4,Final, ',')


for j5=1:n
semilogy(tt,BB(j5,:))
hold on
end
hold off



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%_Lower_Run_FSP_Extinction%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ii=0;
B=Final;
B(fspPosition)=0;
a=newMatrix;

for i = 1:m
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
t(i)=dt*i;
   k1=dt*flow(B,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k1/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k2=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k2/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k3=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k3;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k4=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
   BN = B + (k1 + 2*k2 + 2*k3 + k4)/6;
   t(i) = i*dt;
   B=BN;
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
if mod(i,1000) == 0
ii=ii+1
BB(:,ii)=B;
tt(ii)=i*dt;
B
i
end
end

%%%%%%%%%%%%%_Output for saving%%%%%%%%%%%%%%
for i=1:n %interaction matrix without extinc species%
if B(i)==0;
a(i,:)=0;
a(:,i)=0;
end
end

links=cumsum(a(:,fspPosition));
xfspLinks=links(n);

blinks=cumsum(a(:,extPos));
xextLinks=blinks(n);

for i=1:n
Col=cumsum(a,1);
Row=cumsum(a,2);
pp=Row(:,n)';
zz=(Col(n,:));
qq=pp+zz;
if qq(i)==0;
B(i)=0;
else
B(i)=B(i);
end
end

Finalx=B;%final species list

for i=1:n%Calculating the -cv for population stability
if B(i)>0;
cv(i)=-std(wkeep(BB(i,:),(m*dt)/2,'r'))/mean(wkeep(BB(i,:),(m*dt)/2,'r'));
else
cv(i)=NaN;
end
end
cvAvg=nanmean(cv);
for i=1:n %turning densites in presence absence
if B(i)==0;
yb(i)=0;
else
yb(i)=1;
end
end

ich=cumsum(yb);%calculating species richness
rich=ich(n);
sprich=[rich,cvAvg,fspPosition,xfspLinks];



filename1=['Pops_Lowerx_' num2str(n) '_' num2str(fi) '.csv'];
filename2=['cv_Lowerx_' num2str(n) '_' num2str(fi) '.csv'];
filename3=['Matrix_Lowerx_' num2str(n) '_' num2str(fi) '.csv'];
filename4=['fPop_Lowerx_' num2str(n) '_' num2str(fi) '.csv'];
dlmwrite(filename1,BB, ',')
dlmwrite(filename2,sprich, ',')
dlmwrite(filename3,a, ',')
dlmwrite(filename4,Finalx, ',')
%semilogy(tt,BB(1,:),tt,BB(2,:),tt,BB(3,:),tt,BB(4,:),tt,BB(5,:),tt,BB(6,:),tt,BB(7,:),tt,BB(8,:),tt,BB(9,:),tt,BB(10,:),tt,BB(11,:),tt,BB(12,:),tt,BB(13,:),tt,BB(14,:),tt,BB(15,:),tt,BB(16,:),tt,BB(17,:),tt,BB(18,:),tt,BB(19,:),tt,BB(20,:))

for j5=1:n
semilogy(tt,BB(j5,:))
hold on
end
hold off





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%_Upper_Run_%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=Initial_Biomass;
ii=0;
a=aa;

for i=1:n%calculating metabolic rate (x)
if Troph(i)<= 0
xw(i)=.138;
else
xw(i)=0.314*Mass(i)^-.25;
meta=xw';
metaC=meta/10;
Mnti=meta;
end
end
for i=1:length(consumerPosition)
Mnti(consumerPosition(i))=metaC(i);
end

for i = 1:m
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
t(i)=dt*i;
   k1=dt*flow(B,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k1/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k2=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k2/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k3=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k3;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k4=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
   BN = B + (k1 + 2*k2 + 2*k3 + k4)/6;
   t(i) = i*dt;
   B=BN;
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
if mod(i,1000) == 0
ii=ii+1
BB(:,ii)=B;
tt(ii)=i*dt;
B
i
end
end

%%%%%%%%%%%%%_Output for saving%%%%%%%%%%%%%%
for i=1:n %interaction matrix without extinc species%
if B(i)==0;
a(i,:)=0;
a(:,i)=0;
end
end

for i=1:n
Col=cumsum(a,1);
Row=cumsum(a,2);
pp=Row(:,n)';
zz=(Col(n,:));
qq=pp+zz;
if qq(i)==0;
B(i)=0;
else
B(i)=B(i);
end
end

Final=B;%final species list

for i=1:n%Calculating the -cv for population stability
if B(i)>0;
cv(i)=-std(wkeep(BB(i,:),(m*dt)/2,'r'))/mean(wkeep(BB(i,:),(m*dt)/2,'r'));
else
cv(i)=NaN;
end
end
cvAvg=nanmean(cv);
for i=1:n %turning densites in presence absence
if B(i)==0;
yb(i)=0;
else
yb(i)=1;
end
end

ich=cumsum(yb);%calculating species richness
rich=ich(n);
sprich=[rich,cvAvg,fspPosition,fspLinks];
newMatrix=a;


filename1=['Pops_Upper_' num2str(n) '_' num2str(fi) '.csv'];
filename2=['cv_Upper_' num2str(n) '_' num2str(fi) '.csv'];
filename3=['Matrix_Upper_' num2str(n) '_' num2str(fi) '.csv'];
filename4=['fPop_Upper_' num2str(n) '_' num2str(fi) '.csv'];
dlmwrite(filename1,BB, ',')
dlmwrite(filename2,sprich, ',')
dlmwrite(filename3,a, ',')
dlmwrite(filename4,Final, ',')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%_Upper_Run_fsp_Extinction%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ii=0;
B=Final;
B(fspPosition)=0;
a=newMatrix;

for i = 1:m
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
t(i)=dt*i;
   k1=dt*flow(B,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k1/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k2=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k2/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k3=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k3;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k4=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
   BN = B + (k1 + 2*k2 + 2*k3 + k4)/6;
   t(i) = i*dt;
   B=BN;
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
if mod(i,1000) == 0
ii=ii+1
BB(:,ii)=B;
tt(ii)=i*dt;
B
i
end
end

%%%%%%%%%%%%%_Output for saving%%%%%%%%%%%%%%
for i=1:n %interaction matrix without extinc species%
if B(i)==0;
a(i,:)=0;
a(:,i)=0;
end
end

links=cumsum(a(:,fspPosition));
xfspLinks=links(n);

blinks=cumsum(a(:,extPos));
xextLinks=blinks(n);

for i=1:n
Col=cumsum(a,1);
Row=cumsum(a,2);
pp=Row(:,n)';
zz=(Col(n,:));
qq=pp+zz;
if qq(i)==0;
B(i)=0;
else
B(i)=B(i);
end
end

Finalx=B;%final species list

for i=1:n%Calculating the -cv for population stability
if B(i)>0;
cv(i)=-std(wkeep(BB(i,:),(m*dt)/2,'r'))/mean(wkeep(BB(i,:),(m*dt)/2,'r'));
else
cv(i)=NaN;
end
end
cvAvg=nanmean(cv);
for i=1:n %turning densites in presence absence
if B(i)==0;
yb(i)=0;
else
yb(i)=1;
end
end

ich=cumsum(yb);%calculating species richness
rich=ich(n);
sprich=[rich,cvAvg,fspPosition,xfspLinks];



filename1=['Pops_Upperx_' num2str(n) '_' num2str(fi) '.csv'];
filename2=['cv_Upperx_' num2str(n) '_' num2str(fi) '.csv'];
filename3=['Matrix_Upperx_' num2str(n) '_' num2str(fi) '.csv'];
filename4=['fPop_Upperx_' num2str(n) '_' num2str(fi) '.csv'];
dlmwrite(filename1,BB, ',')
dlmwrite(filename2,sprich, ',')
dlmwrite(filename3,a, ',')
dlmwrite(filename4,Finalx, ',')




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%_METABOLIC_RUN_%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=Initial_Biomass;
ii=0;
a=aa;

for i=1:n%calculating metabolic rate (x)
if Troph(i)<= 0
xw(i)=.138;
else
xw(i)=0.314*Mass(i)^-.25;
meta=xw';
Mnti=meta/10;
end
end

for i = 1:m
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
t(i)=dt*i;
   k1=dt*flow(B,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k1/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k2=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k2/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k3=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k3;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k4=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
   BN = B + (k1 + 2*k2 + 2*k3 + k4)/6;
   t(i) = i*dt;
   B=BN;
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
if mod(i,1000) == 0
ii=ii+1
BB(:,ii)=B;
tt(ii)=i*dt;
B
i
end
end

%%%%%%%%%%%%%_Output for saving%%%%%%%%%%%%%%
for i=1:n %interaction matrix without extinc species%
if B(i)==0;
a(i,:)=0;
a(:,i)=0;
end
end

for i=1:n
Col=cumsum(a,1);
Row=cumsum(a,2);
pp=Row(:,n)';
zz=(Col(n,:));
qq=pp+zz;
if qq(i)==0;
B(i)=0;
else
B(i)=B(i);
end
end

Final=B;%final species list

for i=1:n%Calculating the -cv for population stability
if B(i)>0;
cv(i)=-std(wkeep(BB(i,:),(m*dt)/2,'r'))/mean(wkeep(BB(i,:),(m*dt)/2,'r'));
else
cv(i)=NaN;
end
end
cvAvg=nanmean(cv);
for i=1:n %turning densites in presence absence
if B(i)==0;
yb(i)=0;
else
yb(i)=1;
end
end

ich=cumsum(yb);%calculating species richness
rich=ich(n);
sprich=[rich,cvAvg,fspPosition,fspLinks];
newMatrix=a;

filename1=['Pops_Met_' num2str(n) '_' num2str(fi) '.csv'];
filename2=['cv_Met_' num2str(n) '_' num2str(fi) '.csv'];
filename3=['Matrix_Met_' num2str(n) '_' num2str(fi) '.csv'];
filename4=['fpops_Met_' num2str(n) '_' num2str(fi) '.csv'];
dlmwrite(filename1,BB, ',')
dlmwrite(filename2,sprich, ',')
dlmwrite(filename3,a, ',')
dlmwrite(filename4,Final, ',')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%_METABOLIC_RUN_fsp_extinction%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ii=0;
B=Final;
B(fspPosition)=0;
a=newMatrix;

for i = 1:m
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
t(i)=dt*i;
   k1=dt*flow(B,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition, Mnti);
p=B+k1/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k2=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k2/2;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k3=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
p=B+k3;
for ji=1:n
if p(ji) < 10^-30;
p(ji)=0;
end
end
   k4=dt*flow(p,n,b,b1,om,c,ff,K,meta,y,ee,a,r,fspPosition,Mnti);
   BN = B + (k1 + 2*k2 + 2*k3 + k4)/6;
   t(i) = i*dt;
   B=BN;
for ji=1:n
if B(ji) < 10^-30;
B(ji)=0;
end
end
if mod(i,1000) == 0
ii=ii+1
BB(:,ii)=B;
tt(ii)=i*dt;
B
i
end
end

%%%%%%%%%%%%%_Output for saving%%%%%%%%%%%%%%
for i=1:n %interaction matrix without extinc species%
if B(i)==0;
a(i,:)=0;
a(:,i)=0;
end
end

links=cumsum(a(:,fspPosition));
xfspLinks=links(n);

blinks=cumsum(a(:,extPos));
xextLinks=blinks(n);


for i=1:n
Col=cumsum(a,1);
Row=cumsum(a,2);
pp=Row(:,n)';
zz=(Col(n,:));
qq=pp+zz;
if qq(i)==0;
B(i)=0;
else
B(i)=B(i);
end
end

Finalx=B;%final species list

for i=1:n%Calculating the -cv for population stability
if B(i)>0;
cv(i)=-std(wkeep(BB(i,:),(m*dt)/2,'r'))/mean(wkeep(BB(i,:),(m*dt)/2,'r'));
else
cv(i)=NaN;
end
end
cvAvg=nanmean(cv);
for i=1:n %turning densites in presence absence
if B(i)==0;
yb(i)=0;
else
yb(i)=1;
end
end

ich=cumsum(yb);%calculating species richness
rich=ich(n);
sprich=[rich,cvAvg,fspPosition,xfspLinks];


filename1=['Pops_Metx_' num2str(n) '_' num2str(fi) '.csv'];
filename2=['cv_Metx_' num2str(n) '_' num2str(fi) '.csv'];
filename3=['Matrix_Metx_' num2str(n) '_' num2str(fi) '.csv'];
filename4=['fpops_Metx_' num2str(n) '_' num2str(fi) '.csv'];
dlmwrite(filename1,BB, ',')
dlmwrite(filename2,sprich, ',')
dlmwrite(filename3,a, ',')
dlmwrite(filename4,Finalx, ',')

end