updated the code

3 files changed, 317 insertions, 294 deletions

diff --git a/1040/CH1/EX1.4/Ex1_4.sce b/1040/CH1/EX1.4/Ex1_4.sce
index c4f83141c..57de822b4 100644
--- a/1040/CH1/EX1.4/Ex1_4.sce
+++ b/1040/CH1/EX1.4/Ex1_4.sce
@@ -1,118 +1,116 @@
-//Harriot P.,2003,Chemical Reactor Design (I-Edition) Marcel Dekker,Inc.,USA,pp 436. 

-//Chapter-1 Ex1.4  Pg No. 23

-//Title: Activation energy from packed bed data

-//=========================================================================================================

-clear

-clc

-clf

-// COMMON INPUT

-L= [0 1 2 3 4 5 6 9];//Bed length in feet(ft)

-T=[330 338 348 361 380 415 447 458 ] //Temperature Corresponding the bed length given (°C) 

-R=1.98587E-3;//Gas constant (kcal/mol K)

-

-//CALCLATION (Ex1.4.a)

-//Basis is 1mol of feed A(Furfural) X moles reacted to form Furfuran and CO 

-x=(T-330)./130;//Conversion based on fractional temperature rise

-n=length (T);//6 moles of steam per mole of Furfural is used to decrease temperature rise in the bed

-P_mol=x+7;//Total No. of moles in product stream

-for i=1:(n-1)

-    T_avg(i)= (T(i)+T(i+1))/2

-    P_molavg(i)= (P_mol(i)+P_mol(i+1))/2

-    delta_L(i)=L(i+1)-L(i)

-    k_1(i)=((P_molavg(i))/delta_L(i))*log((1-x(i))/(1-x(i+1)))

-    u1(i)=(1/(T_avg(i)+273.15));

-end

-v1=(log(k_1));

-i=length(u1);

-X1=[u1 ones(i,1) ];

-result1= X1\v1;

-k_1_dash=exp(result1(2,1));

-E1=(-R)*(result1(1,1));

-

-//OUTPUT (Ex1.4.a)

-//Console Output

-mprintf('\n OUTPUT Ex1.4.a');

-mprintf('\n========================================================================================\n')

-mprintf('L \t \t T \t\t x \t\t T_average \t(7+x)ave \tk_1')

-mprintf('\n(ft) \t \t (°C) \t\t   \t\t (°C) \t  ')

-mprintf('\n========================================================================================')

-for i=1:n-1

-mprintf('\n%f \t %f \t %f ',L(i+1),T(i+1),x(i+1))

-mprintf('\t %f \t %f \t %f',T_avg(i),P_molavg(i),k_1(i))

-end

-mprintf('\n\nThe activation energy from the slope =%f kcal/mol',E1 );

-//=====================================================================================================

-

-

-//Title: II Order Reaction 

-//=========================================================================================================

-//CALCULATION (Ex 1.4.b)

-for i=1:(n-1)

-    T_avg(i)= (T(i)+T(i+1))/2

-    P_molavg(i)= (P_mol(i)+P_mol(i+1))/2

-    delta_L(i)=L(i+1)-L(i)

-    k_2(i)=((P_molavg(i))/delta_L(i))*((x(i+1)-x(i))/((1-x(i+1))*(1-x(i))))

-    u2(i)=(1/(T_avg(i)+273.15));

-end

-v2=(log(k_2));

-plot(u1.*1000,v1,'o',u2.*1000,v2,'*');

-xlabel("1000/T (K^-1)");

-ylabel("ln k_1 or ln k_2");

-xtitle("ln k vs 1000/T ");

-legend('ln k_1','ln k_2');

-j=length(u2);

-X2=[u2 ones(j,1) ];

-result2= X2\v2;

-k_2_dash=exp(result2(2,1));

-E2=(-R)*(result2(1,1));

-

-//OUTPUT (Ex 1.4.b)

-mprintf('\n OUTPUT Ex1.4.b');

-mprintf('\n========================================================================================\n')

-mprintf('L \t \t T \t\t x \t\t T_average \t(7+x)ave \tk_2')

-mprintf('\n(ft) \t \t (°C) \t\t   \t\t (°C) \t  ')

-mprintf('\n========================================================================================')

-for i=1:n-1

-mprintf('\n%f \t %f \t %f ',L(i+1),T(i+1),x(i+1))

-mprintf('\t %f \t %f \t %f',T_avg(i),P_molavg(i),k_2(i))

-end

-mprintf('\n\nThe activation energy from the slope =%f kcal/mol',E2 );

-

-//FILE OUTPUT

-fid= mopen('.\Chapter1-Ex4-Output.txt','w');

-mfprintf(fid,'\n OUTPUT Ex1.4.a');

-mfprintf(fid,'\n========================================================================================\n')

-mfprintf(fid,'L \t \t T \t\t x \t\t T_average \t(7+x)ave \tk_1')

-mfprintf(fid,'\n(ft) \t \t (°C) \t\t   \t\t (°C) \t  ')

-mfprintf(fid,'\n========================================================================================')

-for i=1:n-1

-mfprintf(fid,'\n%f \t %f \t %f ',L(i+1),T(i+1),x(i+1))

-mfprintf(fid,'\t %f \t %f \t %f',T_avg(i),P_molavg(i),k_1(i))

-end

-mfprintf(fid,'\n\nThe activation energy from the slope =%f kcal/mol',E1 );

-mfprintf(fid,'\n\n========================================================================================\n')

-mfprintf(fid,'\n OUTPUT Ex1.4.b');

-mfprintf(fid,'\n========================================================================================\n')

-mfprintf(fid,'L \t \t T \t\t x \t\t T_average \t(7+x)ave \tk_2')

-mfprintf(fid,'\n(ft) \t \t (°C) \t\t   \t\t (°C) \t  ')

-mfprintf(fid,'\n========================================================================================')

-for i=1:n-1

-mfprintf(fid,'\n%f \t %f \t %f ',L(i+1),T(i+1),x(i+1))

-mfprintf(fid,'\t %f \t %f \t %f',T_avg(i),P_molavg(i),k_2(i))

-end

-mfprintf(fid,'\n\nThe activation energy from the slope =%f kcal/mol',E2 );

-mclose(all);

-

-//============================================================END OF PROGRAM===========================================

-//Disclaimer (Ex1.4.a):The last value of tavg and k_1  corresponding to L=9 in Table 1.6 (Pg No. 25)of the textbook is a misprint.

-// The value should be 452.5 and 4.955476 respectively instead of 455 and 18.2 as printed in the textbook.

-//Hence there is a change in the activation energy obtained from the code 

-// The answer obtained is 21.3935 kcal/mol instead of 27 kcal/mol as reported in the textbook.

-//Figure 1.8 is a plot between ln k_1 vs 1000/T instead of k_1 vs 1000/T as stated in the solution of Ex1.4.a

-//=========================================================================================================

-//Disclaimer (Ex1.4.b): There is a discrepancy between the computed value of activation energy and value reported in textbook 

-// Error could have been on similar lines as reported for example Ex.1.4.a 

-// Further, intermeidate values for Ex.1.4.b is not available/ reported in textbook and hence could not be compared. 

-//Figure 1.8 is a plot between ln k_2 vs 1000/T instead of k_2 vs 1000/T as stated in the solution of Ex1.4.b

-

-

+//Harriot P.,2003,Chemical Reactor Design (I-Edition) Marcel Dekker,Inc.,USA,pp 436. 
+//Chapter-1 Ex1.4  Pg No. 23
+//Title: Activation energy from packed bed data
+//=========================================================================================================
+clear
+clc
+clf
+// COMMON INPUT
+L= [0 1 2 3 4 5 6 9];//Bed length in feet(ft)
+T=[330 338 348 361 380 415 447 458 ] //Temperature Corresponding the bed length given (°C) 
+R=1.98587E-3;//Gas constant (kcal/mol K)
+
+//CALCLATION (Ex1.4.a)
+//Basis is 1mol of feed A(Furfural) X moles reacted to form Furfuran and CO 
+x=(T-330)./130;//Conversion based on fractional temperature rise
+n=length (T);//6 moles of steam per mole of Furfural is used to decrease temperature rise in the bed
+P_mol=x+7;//Total No. of moles in product stream
+for i=1:(n-1)
+    T_avg(i)= (T(i)+T(i+1))/2
+    P_molavg(i)= (P_mol(i)+P_mol(i+1))/2
+    delta_L(i)=L(i+1)-L(i)
+    k_1(i)=((P_molavg(i))/delta_L(i))*log((1-x(i))/(1-x(i+1)))
+    u1(i)=(1/(T_avg(i)+273.15));
+end
+v1=(log(k_1));
+i=length(u1);
+X1=[u1 ones(i,1) ];
+result1= X1\v1;
+k_1_dash=exp(result1(2,1));
+E1=(-R)*(result1(1,1));
+
+//OUTPUT (Ex1.4.a)
+//Console Output
+mprintf('\n OUTPUT Ex1.4.a');
+mprintf('\n========================================================================================\n')
+mprintf('L \t \t T \t\t x \t\t T_average \t(7+x)ave \tk_1')
+mprintf('\n(ft) \t \t (°C) \t\t   \t\t (°C) \t  ')
+mprintf('\n========================================================================================')
+for i=1:n-1
+mprintf('\n%f \t %f \t %f ',L(i+1),T(i+1),x(i+1))
+mprintf('\t %f \t %f \t %f',T_avg(i),P_molavg(i),k_1(i))
+end
+mprintf('\n\nThe activation energy from the slope =%f kcal/mol',E1 );
+//=====================================================================================================
+
+
+//Title: II Order Reaction 
+//=========================================================================================================
+//CALCULATION (Ex 1.4.b)
+for i=1:(n-1)
+    T_avg(i)= (T(i)+T(i+1))/2
+    P_molavg(i)= (P_mol(i)+P_mol(i+1))/2
+    delta_L(i)=L(i+1)-L(i)
+    k_2(i)=((P_molavg(i))/delta_L(i))*((x(i+1)-x(i))/((1-x(i+1))*(1-x(i))))
+    u2(i)=(1/(T_avg(i)+273.15));
+end
+v2=(log(k_2));
+plot(u1.*1000,v1,'o',u2.*1000,v2,'*');
+xlabel("1000/T (K^-1)");
+ylabel("ln k_1 or ln k_2");
+xtitle("ln k vs 1000/T ");
+legend('ln k_1','ln k_2');
+j=length(u2);
+X2=[u2 ones(j,1) ];
+result2= X2\v2;
+k_2_dash=exp(result2(2,1));
+E2=(-R)*(result2(1,1));
+
+//OUTPUT (Ex 1.4.b)
+mprintf('\n OUTPUT Ex1.4.b');
+mprintf('\n========================================================================================\n')
+mprintf('L \t \t T \t\t x \t\t T_average \t(7+x)ave \tk_2')
+mprintf('\n(ft) \t \t (°C) \t\t   \t\t (°C) \t  ')
+mprintf('\n========================================================================================')
+for i=1:n-1
+mprintf('\n%f \t %f \t %f ',L(i+1),T(i+1),x(i+1))
+mprintf('\t %f \t %f \t %f',T_avg(i),P_molavg(i),k_2(i))
+end
+mprintf('\n\nThe activation energy from the slope =%f kcal/mol',E2 );
+
+//FILE OUTPUT
+fid= mopen('.\Chapter1-Ex4-Output.txt','w');
+mfprintf(fid,'\n OUTPUT Ex1.4.a');
+mfprintf(fid,'\n========================================================================================\n')
+mfprintf(fid,'L \t \t T \t\t x \t\t T_average \t(7+x)ave \tk_1')
+mfprintf(fid,'\n(ft) \t \t (°C) \t\t   \t\t (°C) \t  ')
+mfprintf(fid,'\n========================================================================================')
+for i=1:n-1
+mfprintf(fid,'\n%f \t %f \t %f ',L(i+1),T(i+1),x(i+1))
+mfprintf(fid,'\t %f \t %f \t %f',T_avg(i),P_molavg(i),k_1(i))
+end
+mfprintf(fid,'\n\nThe activation energy from the slope =%f kcal/mol',E1 );
+mfprintf(fid,'\n\n========================================================================================\n')
+mfprintf(fid,'\n OUTPUT Ex1.4.b');
+mfprintf(fid,'\n========================================================================================\n')
+mfprintf(fid,'L \t \t T \t\t x \t\t T_average \t(7+x)ave \tk_2')
+mfprintf(fid,'\n(ft) \t \t (°C) \t\t   \t\t (°C) \t  ')
+mfprintf(fid,'\n========================================================================================')
+for i=1:n-1
+mfprintf(fid,'\n%f \t %f \t %f ',L(i+1),T(i+1),x(i+1))
+mfprintf(fid,'\t %f \t %f \t %f',T_avg(i),P_molavg(i),k_2(i))
+end
+mfprintf(fid,'\n\nThe activation energy from the slope =%f kcal/mol',E2 );
+close;
+
+//============================================================END OF PROGRAM===========================================
+//Disclaimer (Ex1.4.a):The last value of tavg and k_1  corresponding to L=9 in Table 1.6 (Pg No. 25)of the textbook is a misprint.
+// The value should be 452.5 and 4.955476 respectively instead of 455 and 18.2 as printed in the textbook.
+//Hence there is a change in the activation energy obtained from the code 
+// The answer obtained is 21.3935 kcal/mol instead of 27 kcal/mol as reported in the textbook.
+//Figure 1.8 is a plot between ln k_1 vs 1000/T instead of k_1 vs 1000/T as stated in the solution of Ex1.4.a
+//=========================================================================================================
+//Disclaimer (Ex1.4.b): There is a discrepancy between the computed value of activation energy and value reported in textbook 
+// Error could have been on similar lines as reported for example Ex.1.4.a 
+// Further, intermeidate values for Ex.1.4.b is not available/ reported in textbook and hence could not be compared. 
+//Figure 1.8 is a plot between ln k_2 vs 1000/T instead of k_2 vs 1000/T as stated in the solution of Ex1.4.b
\ No newline at end of file
diff --git a/1040/CH2/EX2.1/Ex2_1.sce b/1040/CH2/EX2.1/Ex2_1.sce
index d0026aeba..17089276c 100644
--- a/1040/CH2/EX2.1/Ex2_1.sce
+++ b/1040/CH2/EX2.1/Ex2_1.sce
@@ -1,84 +1,97 @@
-//Harriot P.,2003,Chemical Reactor Design (I-Edition) Marcel Dekker,Inc. USA,pp 436.

-//Chapter-2 Ex2.1 Pg No.52

-//Title: Effectiveness factor for solid catalyzed reaction

-//======================================================================================================================

-clear

-clc

-clf

-//INPUT

-// Case: I  constant hydrogen pressure: P_H2= 2110 torr  

-P_B=[70 185 286];// Benzene Pressure (torr)

-r_1=1E-3 *[4.27 5.4 6.12];//(mol/hr g) observed rates

-P_H2_const=2110;//Constant Hydrogen Pressure (torr)

-

-

-// Case: II  Constant benzene pressure P_B_const=70 torr

-P_H2=[1050 2105 2988];// Hydrogen Pressure (torr)

-r_2=1E-3 * [3.81 4.27 4.5];//(mol/hr g) observed rates

-P_B_const=70;//Constant Benzene Pressure (torr)

-

-//CALCULATION

-// Case: I  constant hydrogen pressure: P_H2= 2110 torr 

-

-n=length(P_B)

-for i=1:n

-    Y_1(i)=(P_B(i)*P_H2_const/r_1(i))^(1/3);

-    X_1(i)=P_B(i);

-end

-coefs_I=regress(X_1',Y_1');

-intercept_1=coefs_I(1)

-slope_1=coefs_I(2)

-

-// Case: II  Constant benzene pressure P_B_const=70 torr

-m=length(P_H2)

-for i=1:n

-    Y_2(i)=(P_B_const*P_H2(i)/r_2(i))^(1/3);

-    X_2(i)=(P_H2(i))^0.5;

-end

-coefs_II=regress(X_2',Y_2');

-intercept_2=coefs_II(1);

-slope_2=coefs_II(2);

-coef_1=(intercept_1)^0.5;

-coef_2=(slope_1*slope_2)^(1/2)*(slope_1/slope_2)*intercept_1;

-

-function y=funct1(K_H2)

-    y=coef_2*K_H2^0.5-coef_1*K_H2^(4/3)-1

-endfunction

-

-[K_H2_res]=fsolve(0,funct1);

-

-K_B=K_H2_res^(4/3)*(slope_1/slope_2);

-

-k=(0.635)^(-1/3)*K_B^2/K_H2_res;

-scf(0)

-plot(X_1,Y_1,'-*-')

-xtitle('Benzene Hydrogenation(a)Variable benzene pressure')

-xlabel("P_B (torr)");

-ylabel("(P_H2 P_B/10^3 r)^(1/3)");

-legend('T=67.6 °C');

-

-scf(1)

-plot(X_2,Y_2,'-*-')

-xtitle('Benzene Hydrogenation(b)Variable hydrogen pressure')

-xlabel("P_H2 (torr)");

-ylabel("(P_H2 P_B/10^3 r)^(1/3)");

-legend('T=67.6 °C');

-

-//OUTPUT

-mprintf('\n Solving for the three parameters gives');

-mprintf('\n K_H2 = %f torr^-1',K_H2_res);

-mprintf('\n K_B = %f torr^-1',K_B);

-mprintf('\n k = %E ',k);

-

-//FILE OUTPUT

-fid= mopen('.\Chapter2-Ex1-Output.txt','w');

-mfprintf(fid,'\n Solving for the three parameters gives');

-mfprintf(fid,'\n K_H2 = %f torr^-1',K_H2_res);

-mfprintf(fid,'\n K_B = %f torr^-1',K_B);

-mfprintf(fid,'\n k = %E ',k);

-mclose(fid);

-

-//===============================================================================================================================================================

-//Disclaimer: Page 53 There is a typo in the equation for Y obtained for Model case I: Constant hydrogen pressure and variable benzene pressure formulation

-// From Fig 2.7(a), It is evident that for P_H2 = 2110 torr, three experimental points are considered for linear regression. However, from table 2.1, only two points corresponds to P_H2 = 2110 torr. In comparison with Fig. 2.7(a), the table value corresponding to P_H2 = 2105 is also read as P_H2 = 2110. 

-//Therefore the values of the constants are different from that obtained in the textbook. Also regression is used to obtain the values of slopes and intercept whereas the textbook considers graphical method for the computation of the codes

+//Harriot P.,2003,Chemical Reactor Design (I-Edition) Marcel Dekker,Inc. USA,pp 436.
+//Chapter-2 Ex2.1 Pg No.52
+//Title: Effectiveness factor for solid catalyzed reaction
+//======================================================================================================================
+clear
+clc
+function [coefs]=regress(x,y)
+ coefs=[]
+  if (type(x) <> 1)|(type(y)<>1) then error(msprintf(gettext("%s: Wrong type for input arguments: Numerical expected.\n"),"regress")), end
+  lx=length(x)
+  if lx<>length(y) then error(msprintf(gettext("%s: Wrong size for both input arguments: same size expected.\n"),"regress")), end
+  if lx==0 then error(msprintf(gettext("%s: Wrong size for input argument #%d: Must be > %d.\n"),"regress", 1, 0)), end
+  x=matrix(x,lx,1)
+  y=matrix(y,lx,1)
+  xbar=sum(x)/lx
+  ybar=sum(y)/lx
+  coefs(2)=sum((x-xbar).*(y-ybar))/sum((x-xbar).^2)
+  coefs(1)=ybar-coefs(2)*xbar
+endfunction
+clf
+//INPUT
+// Case: I  constant hydrogen pressure: P_H2= 2110 torr  
+P_B=[70 185 286];// Benzene Pressure (torr)
+r_1=1E-3 *[4.27 5.4 6.12];//(mol/hr g) observed rates
+P_H2_const=2110;//Constant Hydrogen Pressure (torr)
+
+
+// Case: II  Constant benzene pressure P_B_const=70 torr
+P_H2=[1050 2105 2988];// Hydrogen Pressure (torr)
+r_2=1E-3 * [3.81 4.27 4.5];//(mol/hr g) observed rates
+P_B_const=70;//Constant Benzene Pressure (torr)
+
+//CALCULATION
+// Case: I  constant hydrogen pressure: P_H2= 2110 torr 
+
+n=length(P_B)
+for i=1:n
+    Y_1(i)=(P_B(i)*P_H2_const/r_1(i))^(1/3);
+    X_1(i)=P_B(i);
+end
+coefs_I=regress(X_1',Y_1');
+intercept_1=coefs_I(1)
+slope_1=coefs_I(2)
+
+// Case: II  Constant benzene pressure P_B_const=70 torr
+m=length(P_H2)
+for i=1:n
+    Y_2(i)=(P_B_const*P_H2(i)/r_2(i))^(1/3);
+    X_2(i)=(P_H2(i))^0.5;
+end
+coefs_II=regress(X_2',Y_2');
+intercept_2=coefs_II(1);
+slope_2=coefs_II(2);
+coef_1=(intercept_1)^0.5;
+coef_2=(slope_1*slope_2)^(1/2)*(slope_1/slope_2)*intercept_1;
+
+function y=funct1(K_H2)
+    y=coef_2*K_H2^0.5-coef_1*K_H2^(4/3)-1
+endfunction
+
+[K_H2_res]=fsolve(0,funct1);
+
+K_B=K_H2_res^(4/3)*(slope_1/slope_2);
+
+k=(0.635)^(-1/3)*K_B^2/K_H2_res;
+scf(0)
+plot(X_1,Y_1,'-*-')
+xtitle('Benzene Hydrogenation(a)Variable benzene pressure')
+xlabel("P_B (torr)");
+ylabel("(P_H2 P_B/10^3 r)^(1/3)");
+legend('T=67.6 °C');
+
+scf(1)
+plot(X_2,Y_2,'-*-')
+xtitle('Benzene Hydrogenation(b)Variable hydrogen pressure')
+xlabel("P_H2 (torr)");
+ylabel("(P_H2 P_B/10^3 r)^(1/3)");
+legend('T=67.6 °C');
+
+//OUTPUT
+mprintf('\n Solving for the three parameters gives');
+mprintf('\n K_H2 = %f torr^-1',K_H2_res);
+mprintf('\n K_B = %f torr^-1',K_B);
+mprintf('\n k = %E ',k);
+
+//FILE OUTPUT
+fid= mopen('.\Chapter2-Ex1-Output.txt','w');
+mfprintf(fid,'\n Solving for the three parameters gives');
+mfprintf(fid,'\n K_H2 = %f torr^-1',K_H2_res);
+mfprintf(fid,'\n K_B = %f torr^-1',K_B);
+mfprintf(fid,'\n k = %E ',k);
+mclose(fid);
+
+//===============================================================================================================================================================
+//Disclaimer: Page 53 There is a typo in the equation for Y obtained for Model case I: Constant hydrogen pressure and variable benzene pressure formulation
+// From Fig 2.7(a), It is evident that for P_H2 = 2110 torr, three experimental points are considered for linear regression. However, from table 2.1, only two points corresponds to P_H2 = 2110 torr. In comparison with Fig. 2.7(a), the table value corresponding to P_H2 = 2105 is also read as P_H2 = 2110. 
+//Therefore the values of the constants are different from that obtained in the textbook. Also regression is used to obtain the values of slopes and intercept whereas the textbook considers graphical method for the computation of the codes
\ No newline at end of file
diff --git a/1040/CH3/EX3.5/Ex3_5.sce b/1040/CH3/EX3.5/Ex3_5.sce
index 6e81ca290..f3a104e27 100644
--- a/1040/CH3/EX3.5/Ex3_5.sce
+++ b/1040/CH3/EX3.5/Ex3_5.sce
@@ -1,92 +1,104 @@
-//Harriot P.,2003,Chemical Reactor Design (I-Edition) Marcel Dekker,Inc.,USA,pp 436.

-//Chapter-3 Ex3.5 Pg No. 104

-//Title: Rate Equation to fit Initial Rate data

-//==========================================================================================================

-clear

-clc

-clf()

-//INPUT (Ex3.5.1)

-//Initial Rate Data

-B_by_A= [5 7 10 20 37];//B/A Mol Ratio

-r_0=[75 65 50 33 18];//Rate (mol/hr g)

-

-//CALCULATION (Ex3.5.1)

-//Assuming Eley Rideal Mechanism for the benzene alkylation with propylene 

-for i=1:5

-    C_B(i)= (B_by_A(i)/(1+B_by_A(i)));//In terms of Mol Fraaction

-    C_A(i)= (1/(1+B_by_A(i)));

-    CA_CB(i)=C_B(i)*C_A(i);

-    C_by_r(i)=CA_CB(i)/r_0(i);

-end

-coefs=regress(C_A,C_by_r);//The equation ((C_B*C_A)/r_0)= 1/(k*K_A) + (C_A/k)

-scf(0)

-plot(C_A,C_by_r,'*');

-xtitle('Test of Eley-Rideal model for benzene alkylation');

-xlabel(' CA ,Mol Fraction');

-ylabel('CA CB/r_0');

-intercept=coefs(1);

-slope=coefs(2);

-K_A=slope/intercept;

-k=1/(slope);

-K_A_k=k*K_A;

-

-//OUTPUT (Ex3.5.1)

-mprintf('\n OUTPUT Ex3.5.1');

-mprintf('\n=================================================')

-mprintf('\nThe rate equation for Eley-Ridely Mechanism is:\n    r= %0.0fC_A C_B/(1+%0.2fC_A)',K_A_k,K_A);

-//=========================================================================================================

-

-//Title:Conversion as a function of Space velocity

-//==========================================================================================================

-//INPUT (Ex3.5.2)

-x= [0.16 0.31 0.40 0.75];

-Exp_Inverse_WHSV=(10^-3)*[4 8.2 17 39];//Weight Hourly Space Velocity

-Feed_ratio=10;

-

-//CALCULATION (Ex3.5.2)

-//The integrated rate equation in terms of conversion ln(1/(1-X))+0.236X= 60.4/WHSV (Page no. 106)

-function [y]=integrated_rate_eqn(x0)

-    y=log(1 ./(1-x0))+ 0.236.*x0 - 60.4.*Exp_Inverse_WHSV

-endfunction

-

-n=length(x)

-x0=0.9*ones(1,n);  // Provide guess value for conversion

-[x_predicted]=fsolve(x0,integrated_rate_eqn,1d-15); // Using fsolve to determine conversion from integrated rate expression for each operating WHSV

-

-scf(1)

-plot(Exp_Inverse_WHSV,x,'*',Exp_Inverse_WHSV,x_predicted,'--')

-xtitle('Integral analysis','Inverse of WHSV','Conversion')

-legend('Experimental','Predicted')

-

-//OUTPUT (Ex3.5.2)

-//Console Output

-mprintf('\n=================================================\n');

-mprintf('\n OUTPUT Ex3.5.2');

-mprintf('\n Predicted and Experimental Conversion Values')

-mprintf('\n=================================================')

-mprintf('\n10^3/WHSV\tX_experimental\tX_predicted')

-mprintf('\n=================================================')

-for i=1:n

-    mprintf('\n %0.2f\t\t%0.2f\t\t%0.2f ',Exp_Inverse_WHSV(i)*10^3,x(i),x_predicted(i))

-end

-

-//FILE OUTPUT

-fid= mopen('.\Chapter3-Ex5-Output.txt','w');

-mfprintf(fid,'\n OUTPUT Ex3.5.1');

-mprintf('\n=================================================')

-mfprintf(fid,'\nThe rate equation for Eley-Ridely Mechanism is:\n    r= %0.0fC_A C_B/(1+%0.2fC_A)',K_A_k,K_A);

-mfprintf(fid,'\n=================================================\n')

-mfprintf(fid,'\n OUTPUT Ex3.5.2');

-mfprintf(fid,'\n Predicted and Experimental Conversion Values')

-mfprintf(fid,'\n=================================================')

-mfprintf(fid,'\n10^3/WHSV\tX_experimental\tX_predicted')

-mfprintf(fid,'\n=================================================')

-for i=1:n

-  mfprintf(fid,'\n %0.2f\t\t%0.2f\t\t%0.2f ',Exp_Inverse_WHSV(i)*10^3,x(i),x_predicted(i))

-end

-mclose(fid)

-

-//===========================================END OF PROGRAM=================================

-//Disclaimer:Regression method is used to find the slope and intercept in Ex3.5.2 .

-// Hence the rate equation differ from the graphically obtained values of slope and intercept in the textbook.

-

+//Harriot P.,2003,Chemical Reactor Design (I-Edition) Marcel Dekker,Inc.,USA,pp 436.
+//Chapter-3 Ex3.5 Pg No. 104
+//Title: Rate Equation to fit Initial Rate data
+//==========================================================================================================
+clear
+clc
+clf()
+function [coefs]=regress(x,y)
+coefs=[]
+  if (type(x) <> 1)|(type(y)<>1) then error(msprintf(gettext("%s: Wrong type for input arguments: Numerical expected.\n"),"regress")), end
+  lx=length(x)
+  if lx<>length(y) then error(msprintf(gettext("%s: Wrong size for both input arguments: same size expected.\n"),"regress")), end
+  if lx==0 then error(msprintf(gettext("%s: Wrong size for input argument #%d: Must be > %d.\n"),"regress", 1, 0)), end
+  x=matrix(x,lx,1)
+  y=matrix(y,lx,1)
+  xbar=sum(x)/lx
+  ybar=sum(y)/lx
+  coefs(2)=sum((x-xbar).*(y-ybar))/sum((x-xbar).^2)
+  coefs(1)=ybar-coefs(2)*xbar
+endfunction
+//INPUT (Ex3.5.1)
+//Initial Rate Data
+B_by_A= [5 7 10 20 37];//B/A Mol Ratio
+r_0=[75 65 50 33 18];//Rate (mol/hr g)
+
+//CALCULATION (Ex3.5.1)
+//Assuming Eley Rideal Mechanism for the benzene alkylation with propylene 
+for i=1:5
+    C_B(i)= (B_by_A(i)/(1+B_by_A(i)));//In terms of Mol Fraaction
+    C_A(i)= (1/(1+B_by_A(i)));
+    CA_CB(i)=C_B(i)*C_A(i);
+    C_by_r(i)=CA_CB(i)/r_0(i);
+end
+coefs=regress(C_A,C_by_r);//The equation ((C_B*C_A)/r_0)= 1/(k*K_A) + (C_A/k)
+scf(0)
+plot(C_A,C_by_r,'*');
+xtitle('Test of Eley-Rideal model for benzene alkylation');
+xlabel(' CA ,Mol Fraction');
+ylabel('CA CB/r_0');
+intercept=coefs(1);
+slope=coefs(2);
+K_A=slope/intercept;
+k=1/(slope);
+K_A_k=k*K_A;
+
+//OUTPUT (Ex3.5.1)
+mprintf('\n OUTPUT Ex3.5.1');
+mprintf('\n=================================================')
+mprintf('\nThe rate equation for Eley-Ridely Mechanism is:\n    r= %0.0fC_A C_B/(1+%0.2fC_A)',K_A_k,K_A);
+//=========================================================================================================
+
+//Title:Conversion as a function of Space velocity
+//==========================================================================================================
+//INPUT (Ex3.5.2)
+x= [0.16 0.31 0.40 0.75];
+Exp_Inverse_WHSV=(10^-3)*[4 8.2 17 39];//Weight Hourly Space Velocity
+Feed_ratio=10;
+
+//CALCULATION (Ex3.5.2)
+//The integrated rate equation in terms of conversion ln(1/(1-X))+0.236X= 60.4/WHSV (Page no. 106)
+function [y]=integrated_rate_eqn(x0)
+    y=log(1 ./(1-x0))+ 0.236.*x0 - 60.4.*Exp_Inverse_WHSV
+endfunction
+
+n=length(x)
+x0=0.9*ones(1,n);  // Provide guess value for conversion
+[x_predicted]=fsolve(x0,integrated_rate_eqn,1d-15); // Using fsolve to determine conversion from integrated rate expression for each operating WHSV
+
+scf(1)
+plot(Exp_Inverse_WHSV,x,'*',Exp_Inverse_WHSV,x_predicted,'--')
+xtitle('Integral analysis','Inverse of WHSV','Conversion')
+legend('Experimental','Predicted')
+
+//OUTPUT (Ex3.5.2)
+//Console Output
+mprintf('\n=================================================\n');
+mprintf('\n OUTPUT Ex3.5.2');
+mprintf('\n Predicted and Experimental Conversion Values')
+mprintf('\n=================================================')
+mprintf('\n10^3/WHSV\tX_experimental\tX_predicted')
+mprintf('\n=================================================')
+for i=1:n
+    mprintf('\n %0.2f\t\t%0.2f\t\t%0.2f ',Exp_Inverse_WHSV(i)*10^3,x(i),x_predicted(i))
+end
+
+//FILE OUTPUT
+fid= mopen('.\Chapter3-Ex5-Output.txt','w');
+mfprintf(fid,'\n OUTPUT Ex3.5.1');
+mprintf('\n=================================================')
+mfprintf(fid,'\nThe rate equation for Eley-Ridely Mechanism is:\n    r= %0.0fC_A C_B/(1+%0.2fC_A)',K_A_k,K_A);
+mfprintf(fid,'\n=================================================\n')
+mfprintf(fid,'\n OUTPUT Ex3.5.2');
+mfprintf(fid,'\n Predicted and Experimental Conversion Values')
+mfprintf(fid,'\n=================================================')
+mfprintf(fid,'\n10^3/WHSV\tX_experimental\tX_predicted')
+mfprintf(fid,'\n=================================================')
+for i=1:n
+  mfprintf(fid,'\n %0.2f\t\t%0.2f\t\t%0.2f ',Exp_Inverse_WHSV(i)*10^3,x(i),x_predicted(i))
+end
+mclose(fid)
+
+//===========================================END OF PROGRAM=================================
+//Disclaimer:Regression method is used to find the slope and intercept in Ex3.5.2 .
+// Hence the rate equation differ from the graphically obtained values of slope and intercept in the textbook.
\ No newline at end of file