fn1<-function(C0,k2,r){
m<-matrix(nrow=26,ncol=length(C0))
Ct<-vector("numeric")
C1<-vector("numeric")
dC<-vector("numeric")
j=1
  repeat{
  C1[j]<-C0[j]/2;dC[j]<-C0[j]/4
    repeat{
    m[1,j]=C1[j];m[2,j]=k2*C1[j]^2
    m[26,j]=26*(k2/2)^25*r*C1[j]^26   
	i=3
          repeat{
          m[i,j]<-i*(k2/2)^(i-1)*C1[j]^i
          i<-i+1
          if(i>=26) break}

     Ct[j]<-sum(m[,j])
     if(abs(Ct[j]-C0[j])<10^-6*C0[j]) break
    
     if(Ct[j]>C0[j]) {
     C1[j]=C1[j]-dC[j];dC[j]<-dC[j]/2}
     else{
     C1[j]=C1[j]+dC[j];dC[j]<-dC[j]/2}
    }

         
  j=j+1
  if (j>length(C0)) break }  

C1
}


fn2<-function(C1,C0,k2,r){
   m<-matrix(nrow=26,ncol=length(C0))
   m[1,]<-5.8*(1-0.018*C0)*C1;
   m[2,]<-2^(2/3)*5.8*(1-0.018*C0)*k2*C1^2
   m[26,]<-42*(1-0.019*C0)*26*(k2/2)^25*r*C1^26
     
      j=3
         repeat{
         m[j,]<-j^(5/3)*5.8*(1-0.018*C0)*(k2/2)^(j-1)*C1^j
         j=j+1
         if(j>=26)break   
         }
  
  
     p<-vector("numeric")
      
     k=1
       repeat{
       p[k]<-sum(m[,k])/C0[k]
       k=k+1
       if (k>length(C0))break
       }
p
   }

# load data

ass<-read.table("data.txt",header=T)
plot(ass$C0,ass$Sav,cex=1.5, xlab="C0(mg/ml)",ylab="Sav(S)")

# give initial guess of k2 and r
k2=0.3;r=10^6
fit<-nls(Sav~fn2(fn1(ass$C0,k2fit,rfit),ass$C0,k2fit,rfit),data=ass,start=list(k2fit=k2,rfit=r))