
function [Tvh,Tivh]=van_Herden(Tinit)
% Newton-Raphson iteration to find steady-state temperature profiles
%   fcorr - correction for NR-iteration; reduce if convergense problems
%   Tss = steady-state temperature profile
global mc mh Cp Cpc Q1 Q2 Q3 Npoint Tf cf dm1 dm2 dm3 dHrx u

fcorr=0.2; epsnr=1.e-4;

% Start from some steady-state (Tinit is estimate)
% Tinit=linspace(340,460,30);  % gives middle steady-state ;  Ti=382.6 K
Tss=ssnr(Tinit,fcorr,epsnr); 
u=hx(Tss(Npoint),Tf);
Check_Inlet_temperature=u % should be about 367
% Linearize at this steady-state to find Jacobi matrix
delta=0.001;
Tit=Tss;
for i=1:Npoint
 perturb=zeros(Npoint,1);
 perturb(i)=delta;
 pert1 = righthopenl(Tit+perturb')';
 pert2 = righthopenl(Tit-perturb')';
 A0(:,i)=( pert1 - pert2 )/(2*delta);
end
 J=inv(A0);

% Now start generating S-curve
Tit=Tss;
fcorr=0.5;
for j=1:51
  u=200+(j-1)*5;  % from 200 to 450C
  Ti=u
  error=1;
  k=0;
  while error>1.e-4 
    func=righthopenl(Tit);
    Tit=(Tit'-fcorr*J*func')';
    error=norm(func);
	k=k+1; if k>50, erro2=error; error=-1; end
  end % while
     if error<0, erro2, end
To=Tit(30)
Tivh(j)=Ti;
Tvh(j,:)=Tit;
end %for

% Plotting
%Tovh=Tvh(:,30);
%figure(7)
%epsi=0.6285;
%T0=linspace(200,550,100);
%Ti=epsi*T0+(1-epsi)*Tf;
%plot(Ti,T0,Tivh,Tovh);
%axis([200 450 200 550]);
%xlabel('Inlet temp. reactor, Ti, [deg C]');
%ylabel('Outlet temp. third bed, To, [deg C]');

%% For checking only: Sample Temp profile through reactor
%bl=linspace(1,3,30);
%figure(8)
%idat=43
%Ti1=Tivh(idat), Ti2=epsi*Tvh(idat,30)+(1-epsi)*Tf
%plot(bl,Tvh(idat,:))

