Newton-Raphson routine

MatrixOp /O mDelta = Inv(mJcb)*Diagonal(Cf)

  delta_C	= dif_C[p]*mDelta[p][p]*FxC 

Function NewtonRaphson(m, k , Ct , Ci, fxC, It)

// m		= model in matrix form
// k		= binding constants for species (i.e. K=[HX]/([H][X]) , by definitition K=1 for the free components, (i.e. K= [H]/[H]=1)
// Ct		= total concentration for Components
// Ci		= Initial guess values for Components
// FxC		=wave used to fix concentration of components in Co (i.e. fx={1,1,1,0,1}, fourth component's concentration is fixed)
  		
wave m, k , Ct , Ci, fxC
variable	It

Variable n_comp 	= numpnts (Ct) //number of componets 
Variable n_spec  	= numpnts (k)	 //numer of species
Variable ii = 0
Variable jj,kk=0
	//Check for zero values in Co and replace them with 1e-15
	Ct=Ct[p]==0?1e-15:Ct[p]
	
Duplicate /Free k	 	C_spec, C_Jcb
Duplicate /Free Ct	 	calc_Ct, dif_C, sq_dif, delta_C,absD_C
Duplicate /Free Ci		Cf
Duplicate /Free m		wM						 //Dummy matrix of similar dimentions to model used to perform matrix calculations
Make       /Free /N=(n_comp,n_comp), mJcb, mDelta  //square Jacobian and delta matrix

do
	ii+=1
	//Concentration for species
	wM		= Cf[p]
	wM		= wM^m
	C_spec	= Prod(wM,1)*K
	//Calculated total concentration for components
	wM		= C_spec[q]	
	wM		= m*wM
	MatrixOp /O calc_Ct  = sumRows(wM)
	//Calculate difference between total concentration and calculated concentration for components
	dif_C 	= Ct-calc_Ct
	absD_C	= abs(dif_C)
	
		if (all(absD_C, "<", 1e-15))
			print "small dif"
			return Cf
		else
			//Calculate Jacobian from model matrix and the concentration of species
			for(jj=0;jj<n_comp;jj+=1)
				for (kk=jj;kk<n_comp;kk+=1)
					C_Jcb=m[jj][p]*m[kk][p]*C_spec
					mJcb[jj][kk]=sum(C_Jcb)
					mJcb[kk][jj]=mJcb[jj][kk]
				endfor
			endfor
		endif
	//Calculate shift in concentrations 
	//at the moment the next two lines are not calculating the shift wave properly	
	MatrixOp /O mDelta = Inv(mJcb)*Diagonal(Cf)
	delta_C	= dif_C[p]*mDelta[p][p]*FxC
	Cf		+= delta_C

		if (any(Cf,"<=",0))
			do
				delta_C *=	0.5
				Cf		=	Cf-delta_C
				absD_C	=	abs(delta_C)
				if (all(absD_C, "<", 1e-15))
					break
				endif				
			while (any(Cf,"<=",0))				
		endif
		
while (ii < It)
return Cf
End

<br />
function c_spec=NewtonRaphson(Model, beta, c_tot, c,i)<br />
<br />
ncomp=length(c_tot); % number of components<br />
nspec=length(beta); % number of species<br />
c_tot(c_tot==0)=1e-15; % numerical difficulties if c_tot=0<br />
<br />
it=0;<br />
while it<=99<br />
   it=it+1;<br />
   c_spec =beta.*prod(repmat(c',1,nspec).^Model,1); %species conc<br />
   c_tot_calc=sum(Model.*repmat(c_spec,ncomp,1),2)'; %comp ctot calc<br />
   d =c_tot-c_tot_calc; % diff actual and calc total conc<br />
<br />
   if all(abs(d) <1e-15) % return if all diff small<br />
      return<br />
   end<br />
   for j=1:ncomp % Jacobian (J*)<br />
      for k=j:ncomp<br />
         J_s(j,k)=sum(Model(j,:).*Model(k,:).*c_spec);<br />
         J_s(k,j)=J_s(j,k); % J_s is symmetric<br />
      end<br />
   end<br />
   delta_c=(d/J_s)*diag(c); % equation (2.43)<br />
   c=c+delta_c;<br />
<br />
   while any(c <= 0) % take shift back if conc neg.<br />
      delta_c=0.5*delta_c;<br />
      c=c-delta_c;<br />
      if all(abs(delta_c)<1e-15)<br />
         break<br />
      end<br />
   end<br />
end<br />

Function all (w, cond , n)
	wave 	w
	string	cond
	variable 	n

	variable 	np 	= numpnts (w)
	variable	ii	= 0
	variable 	boolean=0
		strswitch(cond)	// string switch
			case "<":		
					do	
						if (w[ii] < n)
							boolean = 1
						else
							boolean = 0 
						endif
					  	  ii +=1	 
					while (ii < np && boolean==1)
					return boolean
				break		
			case ">":		
					do	
						if (w[ii] > n)
							boolean = 1
						else
							boolean = 0 
						endif
					  	  ii +=1	 
					while (ii < np && boolean==1)
					return boolean
			break
			default:
				return boolean
		endswitch
return boolean

End

Function any ( w,cond,n)
	wave 	w
	string	cond
	variable 	n

	variable 	np 	= numpnts (w)
	variable	ii	= 0
	variable	boolean = 0
		strswitch(cond)	// string switch
			case "<=":		
					do
						if (w[ii] <= n)
							boolean=1 
						endif
						ii +=1
					while (ii < np && boolean==0)
					return boolean
				break				
			case ">=":
					do
						if (w[ii] >= n)
							boolean = 1 
						endif
						ii +=1								
					while (ii < np && boolean==0)						
					return boolean				
				break
			default:
				return boolean
		endswitch
return boolean
End

Function Prod (w,dim)
wave w
variable dim
variable n  = DimSize (w,dim)
variable ii, jj

Make/O/N=(n)/Free  w1=1
	if (dim==0)
		variable col=	DimSize (w,1)
		for(ii=0; ii<n;ii+=1)
				for(jj=0; jj<col;jj+=1)	
					w1[ii] *= w[ii][jj]
				endfor					
		endfor												

	elseif (dim==1)
	variable row=	DimSize (w,0)
			for(ii=0; ii<n;ii+=1)	
					for(jj=0; jj<row;jj+=1)	
						w1[ii] *= w[jj][ii]
					endfor
			endfor	
	endif
	
return w1

End