par = Reactor_numpar;		# Parameters
sym = Reactor_sympar;		# Parameter indices

F_s= [90:10:500];		# Range of flows

Z_a = []; Z_b = []; P = [];
for f_s=F_s
  par(sym.f_s) = f_s;
  [A,B,C,D] = Reactor_sm(par);	# Linearised system

  p = sort(eig(A));
  P = [P p];

  C_a = C([1 3],:);		# C vector for c_a and t
  D_a = D([1 3],:);		# D vector for c_a and t
  z_a = tzero(A,B,C_a,D_a);	# Transmission zeros for c_a and t
  Z_a = [Z_a z_a];

  C_b = C(2:3,:);		# C vector for c_b and t
  D_b = D(2:3,:);		# D vector for c_b and t
  z_b = tzero(A,B,C_b,D_b);	# Transmission zeros for c_b and t
  Z_b = [Z_b z_b];
endfor

grid; xlabel("f_s"); ylabel("p1,p2");
plot(F_s,P(1:2,:));
psfig("Reactor_pole_1_2");

grid; xlabel("f_s"); ylabel("p3");
plot(F_s,P(3,:));
psfig("Reactor_pole_3");

grid; xlabel("f_s"); ylabel("z_a");
plot(F_s,Z_a);
psfig("Reactor_zero_a");

grid; xlabel("f_s"); ylabel("z_b");
plot(F_s,Z_b);
psfig("Reactor_zero_b");