c	Simplified TE model in .MEX format

c	N. L. Ricker, U. of Washington

c	STATE VARIABLES are molar holdups [kmol] 
c	in the specified location unless stated otherwise.

c	 1	A in reactor
c	 2	B in reactor
c	 3	C in reactor
c	 4	D in reactor
c	 5	E in reactor
c	 6	F in reactor
c	 7	G in reactor
c	 8	H in reactor
c	 9	A in separator
c	10	B in separator
c	11	C in separator
c	12	D in separator
c	13	E in separator
c	14	F in separator
c	15	G in separator
c	16	H in separator
c	17	A in feed zone
c	18	B in feed zone
c	19	C in feed zone
c	20	D in feed zone
c	21	E in feed zone
c	22	F in feed zone
c	23	G in feed zone
c	24	H in feed zone
c	25	D in stripper bottoms
c	26	E in stripper bottoms
c	27	F in stripper bottoms
c	28	G in stripper bottoms
c	29	H in stripper bottoms
c	30	Separator internal energy [MJ/kmol]
c	31	Stripper internal energy [MJ/kmol]

c	INPUTS are flows [kmol/h] unless noted otherwise:

c	 1	Feed stream 1 (pure A)
c	 2	Feed stream 2 (pure D)
c	 3	Feed stream 3 (pure E)
c	 4	Feed stream 4 (A & C)
c	 5	Compressor valve [Percent]
c	 6	Purge (stream 9)
c	 7	Separator underflow (stream 10)
c	 8	Stripper underflow (stream 11)
c	 9	Steam flow [kg/h]
c	10	Reactor temperature [C]
c	11	Condenser valve [%]
c	12	A in stream 4 [Mole Percent]
c	13	B in stream 4 [Mole Percent]
c	14	Reaction 1 activity factor [Percent].
c	15	Reaction 2 activity factor [Percent].
c	16	Stream 10 bias flow [kmol/h].
c	17	Reactor/Sep flow parameter [Percent].
c	18	Feed/Reactor flow parameter [Percent].
c	19	NOT USED in this version, i.e., ignored.
c	20	Adjustment to VLE of D to G in separator [Percent]
c	21	Adjustment to H VLE in separator [Percent]
c	22	C bias flow to feed zone [kmol/h]
c	23	D bias flow to feed zone [kmol/h]
c	24	E bias flow to feed zone [kmol/h]
c	25	Reaction 3 activity factor [Percent]
c	26	Adjustment to VLE in reactor [Percent]
c	27	Bias correction to compressor performance curve
c	28	Bias correction to compressor power prediction
c	29	Product density correction [Percent]
c	30	# equil. stages for D in stripper
c	31	# equil. stages for E in stripper
c	32	# equil. stages for F in stripper
c	33	Condenser cooling water inlet temperature [C]
c	34	UA in condenser [kW/K]
c	35	Stripper heat loss [kW]

c	OUTPUTS are mole % in stream unless noted otherwise:

c	 1	Reactor pressure [kPa]
c	 2	Reactor liq. holdup [%]
c	 3	Separator pressure [kPa]
c	 4	Separator liq. holdup [%]
c	 5	Product liq. holdup [%]
c	 6	Feed zone pressure [kPa]
c	 7	Total feed entering reactor (stream 6) [kscmh]
c	 8	A in reactor feed (stream 6)
c	 9	B in reactor feed (stream 6)
c	10	C in reactor feed (stream 6)
c	11	D in reactor feed (stream 6)
c	12	E in reactor feed (stream 6)
c	13	F in reactor feed (stream 6)
c	14	A in purge (stream 9)
c	15	B in purge (stream 9)
c	16	C in purge (stream 9)
c	17	D in purge (stream 9)
c	18	E in purge (stream 9)
c	19	F in purge (stream 9)
c	20	G in purge (stream 9)
c	21	H in purge (stream 9)
c	22	D in product (stream 11)
c	23	E in product (stream 11)
c	24	F in product (stream 11)
c	25	G in product (stream 11)
c	26	H in product (stream 11)
c	27	Volumetric production rate [m^3/h].
c	28	Recycle rate (stream 8) [kscmh]
c	29	Compressor power [kW]
c	30	Stripper Temperature [C]
c	31	Separator Temperature [C]
c	32	Condenser coolant outlet temp [C]
c	33	Cost [cents/kmol product]
c	34	Rate of reaction 1 [kmol G produced/h]
c	35	Rate of reaction 2 [kmol H produced/h]
c	36	Rate of reaction 3 [kmol F produced/h]
c	37	Partial pressure of A in reactor [kPa]
c	38	Partial pressure of C in reactor [kPa]
c	39	Partial pressure of D in reactor [kPa]
c	40	Partial pressure of E in reactor [kPa]
c
C---------------------------------------------------------------------

      SUBROUTINE mexFunction(NLHS, PLHS, NRHS, PRHS)
      INTEGER PLHS(*), PRHS(*)
      INTEGER NLHS, NRHS
C
      INTEGER mxCreateFull, mxGetPr, mxGetPi, mexGetMatrixPtr
C
      DOUBLE PRECISION mxGetScalar
C
C THE ABOVE SUBROUTINE, ARGUMENT, AND FUNCTION DECLARATIONS ARE STANDARD
C IN ALL FORTRAN MEX FILES.
C---------------------------------------------------------------------
C

c	This is a MEX interface routine for Matlab version 4.x

c	It is designed to emulate a Simulink S-function.

c	The following will hold pointers:

	integer t,x,u,flag

c		The number of states, inputs, and outputs are specified
c		by the parameters nx, nu, and ny, respectively.  
c		Also define corresponding character constants
c		lenx, lenu, to use in error messages, and double precision
c		constants for the SIZES vector.

	integer iflag,m,n,nx,nu,ny
	character*(*) lenu,lenx
	double precision dnx, dnu, dny
	parameter (nx=31, nu=35, ny=40, lenx='31', lenu='35', dnx=31.,
     +		dnu=35.,dny=40.)
	
	
	double precision 	rx(nx),rx0(nx),dxdt(nx),sizes(6),rt,ru(nu),
     +			ry(ny)

	character*80 message	! For debug messages to screen
	
c		The following is the appropriate definition of the SIZES
c		vector for SIMULINK system functions when there are no
c		discrete-time states and no direct feed-through from input
c		to output.

	data sizes/dnx, 0., dny, dnu, 0., 0./
	
	data rx0/nx*0.d0/

c		Check the number of input arguments

	if (nrhs .gt. 0 .and. nrhs .ne. 4) then
		call mexErrMsgTxt('TEest5 needs 4 input variables.')
	else if (nrhs .le. 0) then
		iflag=0
	else

c			There are 4 input arguments.  Set pointers and copy the
c			inputs to local storage.

		m=mxGetM(prhs(1))				! Get size of "t" array.
		n=mxGetN(prhs(1))
		if (m.gt.0 .and. n.gt.0) then		! If it's non-empty...
			rt=mxGetScalar(prhs(1))		! gets time value (a scalar).
		else
			rt=0.0				! Otherwise set to zero.
		end if
		
		m=mxGetM(prhs(4))				! FLAG
		n=mxGetN(prhs(4))
		if (m.gt.0 .and. n.gt.0) then
			if (m.eq.1 .and. n.eq.1) then
				iflag=abs(nint(mxGetScalar(prhs(4))))
			else
				call mexErrMsgTxt('FLAG must be a scalar.')
			end if
		else
			iflag=0
		end if
			
		m=mxGetM(prhs(2))				! X
		n=mxGetN(prhs(2))
		if (m.gt.0 .and. n.gt.0) then		
		   if (min(m,n).ne.1 .or. max(m,n) .ne. nx) then	
			call mexErrMsgTxt('X must be a vector, length '//lenx)
		   end if
		   x=mxGetPr(prhs(2))
		   call mxCopyPtrToReal8(x,rx,nx)		! copy x into local storage
		else if (iflag.eq.1 .or. iflag.eq.3) then
			call mexErrMsgTxt('X must be supplied.')
		end if
				
		m=mxGetM(prhs(3))				! U
		n=mxGetN(prhs(3))
		if (m.gt.0 .and. n.gt.0) then		
			if (min(m,n).ne.1 .or. max(m,n) .ne. nu) then	
				call mexErrMsgTxt
     +			('U must be a vector, length '//lenu)
			end if
			u=mxGetPr(prhs(3))
			call mxCopyPtrToReal8(u,ru,nu)		! copy u into local storage
		else if (iflag.eq.1 .or. iflag.eq.3) then
			call mexErrMsgTxt('U must be supplied.')
		end if
		
	end if	
	
c		Now return outputs according to the value of "flag"

	if (iflag.eq.0) then
	
c			Return default initial states and size info
	
		plhs(1)=mxCreateFull(6,1,0)				! "sizes"
		call mxCopyReal8ToPtr(sizes,mxGetPr(plhs(1)),6)
		if (nlhs.gt.1) then
			plhs(2)=mxCreateFull(nx,1,0)			! "x0"
			call mxCopyReal8ToPtr(rx0,mxGetPr(plhs(2)),nx)
		end if
		if (nlhs.gt.2) then
			plhs(3)=mxCreateFull(0,0,0)			! "xstr" (empty)
		end if
		
	else if (iflag.eq.1 .or. iflag.eq.3) then
	
		call TEsimf(nx,nu,ny,rx,ru,ry,dxdt)    	! Main simulation call

		if (iflag.eq.1) then

c				Return dxdt.  If there are 2 LHS arguments
c				return y as the second one.  (non-standard!!)

			plhs(1)=mxCreateFull(nx,1,0)
			call mxCopyReal8ToPtr(dxdt,mxGetPr(plhs(1)),nx)  	! returns dxdt to Matlab
			if (nlhs.ge.2) then
				plhs(2)=mxCreateFull(ny,1,0)
				call mxCopyReal8ToPtr(ry,mxGetPr(plhs(2)),ny)  	! returns y to Matlab
			end if
		else
			plhs(1)=mxCreateFull(ny,1,0)
			call mxCopyReal8ToPtr(ry,mxGetPr(plhs(1)),ny)  	! returns y to Matlab
		end if		
	else

c			For all other flags, return an empty matrix.	

		plhs(1)=mxCreateFull(0,0,0)
	
	end if
	
	return
	end
	
C---------------------------------------------------------------------

	SUBROUTINE TEsimf(nx,nu,ny,x,u,y,dxdt)
	
c	Calculates outputs and rates of change of states.

c	This subroutine can be used independently, i.e., without
c	the MEX interface.

c	NOTE:  all non-integer variables are double precision.

c	Inputs to routine:

c	nx	number of states.  Set to integer 31.
c	nu	number of inputs.  Set to integer 35.
c	ny	number of outputs.  Set to integer 40.
c	x	current state vector, length nx.
c	u	current input vector, length nu.

c	Outputs from routine

c	y	output vector, length ny.
c	dxdt	time derivatives of states, length nx.

	IMPLICIT DOUBLE PRECISION (A-H,O-Z)

	integer nx,nu,ny,nc
	double precision cf2cm,dzero,pctGH
	
	double precision x(nx),u(nu),y(ny),dxdt(nx)
	
	parameter (nc=8,cf2cm=0.028317,TauTr=0.08333,TauTs=0.08333,
     +		TauF=0.005,TauC=0.08333,dzero=0.0d0,Rg=8.314)
	
	double precision 	rlmol(nc),slmol(nc),rvmol(nc),svmol(nc),
     +			plmol(nc),fvmol(nc),xlr(nc),xls(nc),RR(4),
     +			ppr(nc),pps(nc),xvr(nc),xvs(nc),xvf(nc),
     +			Fcmol(nc),CFbias,CWbias,dP,Cflow,Crecyc,
     +			F8mass,Kval(nc),Pvc(nc),Sfac(nc),F4i(nc),
     +			F5i(nc),pvapr(nc),pvaps(nc),pvapc(nc),
     +			x4(8),x11(8)
     
	double precision 	mwts(nc),molvol(nc),
     +			sfr(nc),gami(nc),gamr(nc)
     
	integer code
	
	data mwts/2.0,25.4,28.0,32.0,46.0,48.0,62.0,76.0/	! molec. wts.
	data molvol/3*0.0,0.1070,0.1260,0.1463,0.1013,0.1231/	! [m3/kmol]
	data sfr/6*1.0,0.07,0.04/					! stripping factors
	data gami/8*1.0/				! sep. activ. coeffs.
	data gamr/8*1.0/				! react. activ. coeffs.
	data Cpw/4.18/				! specific heat of water [kJ/kg-K]
	data Cvw/2725./				! condenser coolant coefficient  [kg/h-%]	

c		Initialization

	do i=1,8
		x4(i)=0.0
		F4i(i)=0.0
	end do
	
	Fsteam=u(9)				! steam flow [kg/h]

c		Get disturbances from input vector. Flow disturbances are scaled so
c		they're comparable to manipulated variables (0-100% basis)

	x4(1)=u(12)/100.			! mol frac A in stream 4 (A+C feed)
	x4(2)=u(13)/100.			! mol frac B in stream 4 (A+C feed)
	x4(3)=1.0-x4(1)-x4(2)		! get C in stream 4 by difference
	R1F=u(14)/100.			! Reactivity factor, reaction 1
	R2F=u(15)/100.			! Reactivity factor, reaction 2
	F10b=u(16)				! Bias adjustment for stream 10 flow [kmol/h]
	RSflow=u(17)/100.			! Adjusts flow between reactor and separator
	FRflow=u(18)/100.			! Adjusts flow between feed zone and reactor
	fracGH=u(19)/100.			! Adjusts for impurity in product (components other
						! than G and H)
	do i=4,7
		gami(i)=u(20)/100.	! Adjusts VLE of D to G in separator
	end do
	gami(8)=u(21)/100.		! Adjusts VLE of H in separator
	Cbias=u(22)				! C flow bias at feed point
	Dbias=u(23)				! D flow bias at feed point
	Ebias=u(24)				! E flow bias at feed point
	R3F=u(25)/100.			! Reactivity factor, reaction 3
	do i=4,8
		gamr(i)=u(26)/100.	! Adjusts VLE in reactor
	end do
	CFbias=u(27)			! Bias correction to compressor flow curve
	CWbias=u(28)			! Bias correction to compressor power curve
	Pvolcf=u(29)/100.			! Product density correction
	DNeq=u(30)				! Equil. stages for D
	ENeq=u(31)				! Equil. stages E
	FNeq=u(32)				! Equil. stages F
	Tw0=u(33)				! Condenser coolant inlet temp [C]
	UAcon=u(34)*3600.			! UA in condenser [kJ/h-K]
	QLoss=u(35)*3600.			! Heat loss in stripper [kJ/h]
	

c		Get molar holdups in reactor, separator, and feed vapor.  
c		Assumes A, B, C insoluble in liquid.

	rlsum=0.0
	slsum=0.0
	plsum=0.0
	rvsum=0.0
	svsum=0.0
	VLR=0.0
	VLS=0.0
	VLP=0.0
	
	do i=1,3
		rlmol(i)=0.0
		slmol(i)=0.0
		plmol(i)=0.0
	end do
	
	do i=4,nc	
		rlmol(i)=x(i)	! gets liq. mol from state vector
		slmol(i)=x(i+8)
		plmol(i)=x(i+21)
		rlsum=rlsum+rlmol(i)		! accumulate total moles
		slsum=slsum+slmol(i)
		plsum=plsum+plmol(i)
		VLR=VLR+rlmol(i)*molvol(i)	! accumulate liquid volume [m3]
		VLS=VLS+slmol(i)*molvol(i)
		VLP=VLP+plmol(i)*molvol(i)
	end do
	
	fvsum=0.0
	do i=1,nc
		fvmol(i)=x(i+16)			! gets feed vapor from state vector
		fvsum=fvsum+fvmol(i)		! total feed vapor moles
	end do

	
c	Temperatures in reactor, separator, & stripper

	TCR=u(10)				! reactor temp [C]
	TKR=TCR+273.2			!              [K]
	enint=1000.0*x(30)*slsum	! get extensive internal energy [kJ]
	call UTfunc(enint,slmol,TCS)	! get separator temp [C]
	TKS=TCS+273.2			!                [K]
	enint=1000.0*x(31)*plsum	! get extensive internal energy [kJ]
	call UTfunc(enint,plmol,TCC)	! get stripper temp [C]

c		Volumes in reactor, separator, & feed zone

	VTR=36.8				! total reactor volume [m3]
	VTS=99.1				! total separator volume [m3]
	VTV=150.0				! total feed vapor holdup volume [m3]
	VVR=VTR-VLR				! vapor in reactor [m3]
	VVS=VTS-VLS				! vapor in separator [m3]
	
c		Mole fractions, partial pressures, etc.

	PTR=0.0
	PTS=0.0
	PTV=fvsum*Rg*(273.2+86.1)/VTV		! Total pressure in feed zone [kPa]
							! (assumes constant temperature of 86.1 C)

	do i=1,3							! for A, B, C
		ppr(i)=x(i)*Rg*TKR/VVR		! pp in reactor [kPa]
		pps(i)=x(i+8)*Rg*TKS/VVS	! pp in separator [kPa]
		PTR=PTR+ppr(i)			! total P in reactor [kPa]
		PTS=PTS+pps(i)			! total P in separator [kPa]
		xlr(i)=0.0
		xls(i)=0.0
		x11(i)=0.0
	end do

	call VPfunc(tcr,pvapr)						! vap. pres. in reactor [kPa]
	call VPfunc(tcs,pvaps)						! vap. pres. in separator [kPa]
	do i=4,nc								! For D,E,F,G,H
		xlr(i)=rlmol(i)/rlsum					! reactor liquid mol frac
		ppr(i)=pvapr(i)*xlr(i)*gamr(i)			! VLE
		xls(i)=slmol(i)/slsum					! Repeat for separator
		pps(i)=pvaps(i)*xls(i)*gami(i)				! Note use of VLE adjustment
		PTR=PTR+ppr(i)						! Reactor total press [kPa]
		PTS=PTS+pps(i)						! Separator total press [kPa]
		x11(i)=plmol(i)/plsum					! Product mol fractions
	end do
	
	wtmolr=0.0
	wtmolf=0.0
	wtmols=0.0
	do i=1,8
		xvr(i)=ppr(i)/PTR			! reactor vapor mol frac
		xvs(i)=pps(i)/PTS			! sep vapor mol frac
		xvf(i)=fvmol(i)/fvsum		! feed vapor mol frac
		wtmolr=wtmolr+xvr(i)*mwts(i)	! average mol wt of reactor vapor
		wtmols=wtmols+xvs(i)*mwts(i)	! average mol wt of separator vapor (& purge)
		wtmolf=wtmolf+xvf(i)*mwts(i)	! average mol wt of feed vapor
	end do

c		Rate Laws.  All rates in [kmol/h].

      RR(1)=VVR*DEXP(44.06-42600.0/1.987/TKR)*R1F*(PPR(1)**1.08)
     +		*(PPR(3)**0.311)*(PPR(4)**0.874)
      RR(2)=VVR*DEXP(10.27-19500.0/1.987/TKR)*R2F*(PPR(1)**1.15)
     +		*(PPR(3)**0.370)*(PPR(5)**1.00)
      RR(3)=VVR*DEXP(59.50-59500.0/1.987/TKR)*R3F*PPR(1)
     +		*(0.77*PPR(4)+PPR(5))
	

c		Define total flows [kmol/h].  First get those defined directly
c		as inputs.

	F1=u(1)
	F2=u(2)
	F3=u(3)
	F4=u(4)
	F9=u(6)
	F10=u(7) - F10b
	F11=u(8)
	
c		Now use the compressor curve and valve setting to get the recycle rate.
c		dP is pressure drop across compresor (feed zone P - separator P)
c		Cflow is compressor throughput in kg/h.
c		Crecyc is back-flow through compressor valve [kg/h]
c		F8mass is net mass flow of recycle (stream 8) [kg/h]
c		F8 is corresponding molar flow.
c		Cpower is compressor power [kW]

	dP=PTV-PTS				! force non-negative pressure drop					
	Cflow=-0.05152*dP*dP - 92.52*dP + 1.17e5 + CFbias
	Crecyc=u(5) * 60.83 * sqrt(dP)
	F8mass=Cflow-Crecyc
	F8=F8mass/wtmols									! [kmol/h]
	Cpower=TKS*Rg*(-9.365e-8*dP*dP + 4.300e-4*dP - 8.667e-3 + CWbias)	! [kJ/kmol]
	Cpower=Cpower*(Cflow/wtmols/3600.0)				! [kW]

c		Stripper modeling.  Removal of D, E, & F is based on the KSB equations.
c		Product composition is calculated on basis of zero A, B, and C.
c		Note that stream 10 bias is included in liquid rate for consistency 
c		with steady-state calculations in TEest6_ss.m.

	call VPfunc(Tcc,pvapc)
	do i=4,8
		Sfac(i)=(pvapc(i)/PTV)*F4/(F10+F10b)		! Ki * V / L
	end do

	sfr(4)=1.0d0 - (1.0d0 - Sfac(4))/(1.0d0 - Sfac(4)**(DNeq+1.0d0))
	sfr(5)=1.0d0 - (1.0d0 - Sfac(5))/(1.0d0 - Sfac(5)**(ENeq+1.0d0))
	sfr(6)=1.0d0 - (1.0d0 - Sfac(6))/(1.0d0 - Sfac(6)**(FNeq+1.0d0))
	sfr(7)=Sfac(7)*0.809
	sfr(8)=Sfac(8)*0.720
	
c		Stripper balance to get molar rates in stream 5.

	F4i(1)=x4(1)*F4
	F4i(2)=x4(2)*F4
	F4i(3)=x4(3)*F4
	F5i(1)=F4i(1)
	F5i(2)=F4i(2)
	F5i(3)=F4i(3)
	
	do i=4,8
		F5i(i)=sfr(i)*F10*xls(i)
	end do

c		Stripper balance to get combined feeds (streams 1+2+3+5)
c		Also account for flow biases.

	Fcmol(1)=F5i(1)+F1
	Fcmol(2)=F5i(2)
	Fcmol(3)=F5i(3)+Cbias
	Fcmol(4)=F5i(4)+F2+Dbias
	Fcmol(5)=F5i(5)+F3+Ebias
	Fcmol(6)=F5i(6)
	Fcmol(7)=F5i(7)
	Fcmol(8)=F5i(8)
	Fcomb=0.0
	do i=1,8
		Fcomb=Fcomb+Fcmol(i)
	end do
		
c		Pressure drop vs. flow to get flows of streams 6 & 7.

	F6=(2413.7/wtmolf)*((PTV-PTR)**0.5)*FRflow
	if (ptr .gt. ptv) then
		F6=-F6
	end if
	F7=(5722.0/wtmolr)*((PTR-PTS)**0.5)*RSflow
	if (pts .gt. ptr) then
		F7=-F7
	end if

c		Balances to get rates of change of states.  First do  
c		convection terms (in - out).

	do i=1,8
		dxdt(i)=F6*xvf(i) - F7*xvr(i)					! Reactor
		dxdt(i+8)=F7*xvr(i) - (F8+F9)*xvs(i)
     +		 		- F10*xls(i)				! Separator
     		dxdt(i+16)=Fcmol(i)+F8*xvs(i)-F6*xvf(i)			! Feed zone
	end do
	
	dmdt1=0.0	! Accumulation of liquid mass in separator
	dmdt2=0.0	! Accumulation of liquid mass in stripper
	do i=4,8
		dxdt(i+21)=(1.0-sfr(i))*F10*xls(i) - F11*x11(i)		! Stripper
		dmdt1=dmdt1+dxdt(i+8)
		dmdt2=dmdt2+dxdt(i+21)
	end do
	
c		Now account for reaction.

	dxdt(1)=dxdt(1) - rr(1) - rr(2) - 0.333*rr(3)
	dxdt(2)=dxdt(2)
	dxdt(3)=dxdt(3) - rr(1) - rr(2)
	dxdt(4)=dxdt(4) - rr(1) - rr(3)
	dxdt(5)=dxdt(5) - rr(2) - 0.333*rr(3)
	dxdt(6)=dxdt(6) + rr(3)
	dxdt(7)=dxdt(7) + rr(1)
	dxdt(8)=dxdt(8) + rr(2)

c		Energy balances

	T4=45.0						! temperature of stream 4 (from base case)
	call HVfunc(T4,F4i,H4)				! enthalpy of stream 4 [kJ/h]
	call HVfunc(tcc,F5i,H5)				! enthalpy of stream 5 [kJ/h]
	call HVfunc(tcr,xvr,H7)
	H7=F7*H7						! enthalpy of reactor effluent [kJ/h]
	call HVfunc(tcs,xvs,H89)
	H89=(F8+F9)*H89					! enthalpy of vapor leaving separator
	call HLfunc(tcs,xls,H10)
	H10=F10*H10						! enthalpy of liquid leaving separator
	call HLfunc(tcc,x11,H11)
	H11=F11*H11						! enthalpy of product
	Qsteam=Fsteam*2260.				! stripper heat duty (steam) [kJ/h]
	Fw=Cvw*u(11)					! condenser coolant rate  [kg/h]
	
c		Combine heat balance on condenser coolant with heat transfer rate
c		equation to get coolant outlet temperature and duty in condenser.

	alph=Fw*Cpw/UAcon					! Dimensionless group
	Tw=((2*alph-1.)*Tw0 + 2*Tcs)/(2*alph+1.)	! Temp of coolant leaving condenser [C]
	Qcond=UAcon*(0.5*(Tw+Tw0) - Tcs)			! Condenser duty [kJ/h]

c		Accumulation of internal energy

	dxdt(30)=0.001*(Qcond+H7-H89-H10)			! separator heat balance [MJ/h]
	dxdt(31)=0.001*(Qsteam-QLoss+H4+H10-H5-H11)	! stripper heat balance [MJ/h]

c		Convert to intensive form [MJ/kmol]

	dxdt(30)=(dxdt(30) - x(30)*dmdt1)/slsum
	dxdt(31)=(dxdt(31) - x(31)*dmdt2)/plsum

c		Calculate the outputs.

	y(1)=PTR-101.					! Reactor pressure 	[kPa guage]
	y(2)=5.263*VLR - 12.105				! Reactor liq. 		[%]
	y(3)=PTS-101.					! Separator pressure	[kPa guage]
	y(4)=12.28*VLS - 10.53				! Separator liquid	[%]
	y(5)=22.58*VLP - 49.03				! Product liq. holdup 	[%]
	y(6)=PTV-101.					! Pressure in feed zone [kPa guage]
	y(7)=F6/44.79					! Stream 6 flowrate 	[kscmh]
	do i=1,6
		y(i+7)=xvf(i)*100.0			! Reactor feed 		[mol %]
		y(i+13)=xvs(i)*100.0			! Purge 			[mol %]
	end do

	y(20)=xvs(7)*100.0				! G in purge		[mol %]
	y(21)=xvs(8)*100.0				! H in purge		[mol %]
	y(22)=x11(4)*100.0				! D in product		[mol %]
	y(23)=x11(5)*100.0				! E in product		[mol %]
	y(24)=x11(6)*100.0				! F in product		[mol %]
	y(25)=x11(7)*100.0				! G in product		[mol %]
	y(26)=x11(8)*100.0				! H in product		[mol %]
	y(27)=F11*VLP/plsum*Pvolcf			! Production 		[m^3/h]
	y(28)=F8/44.79					! Recycle rate		[kscmh]
	y(29)=Cpower					! Compressor power	[kW]
	y(30)=Tcc						! Stripper Temperature	[C]
	y(31)=Tcs						! Separator temperature  [C]
	y(32)=Tw						! Condenser outlet temperature  [C]

	purgls=F9*(221.*xvs(1) 
     +		+ 618.*xvs(3) + 2210.*xvs(4)
     +		+ 1460.*xvs(5) + 1790.*xvs(6) 
     +		+ 3040.*xvs(7) + 2290.*xvs(8))	! Purge loss cost [$/h]
	prodls=F11*(2210.*x11(4) + 1460.*x11(5) + 1790.*x11(6))	! Product loss cost [$/h]
     	
	y(33)=(5.36*Cpower + 3.18*u(9) + purgls + prodls)/F11     	! [cents/kmol]

	y(34)=RR(1)						! G production		[kmol/h]
	y(35)=RR(2)						! H production		[kmol/h]
	y(36)=RR(3)						! F production		[kmol/h]
	y(37)=PPR(1)					! A PP in reactor		[kPa]
	y(38)=PPR(3)					! C PP in reactor		[kPa]
	y(39)=PPR(4)					! D PP in reactor		[kPa]
	y(40)=PPR(5)					! F PP in reactor		[kPa]

	return
	end
	
C---------------------------------------------------------------------

	subroutine HLfunc(T,x,hl)
	
c		Calculates liquid enthalpy for TE problem [kJ].
c
c		T = temperature [C]
c		x = vector of molar compositions (length 8).
c		hl= liquid enthalpy [kJ]

c		Reference is liquid at 100 C.  Components A, B,
c		and C are ignored for liquid-phase calculations.

	double precision T,x(8),hl
	
	double precision mwts(8),CpL(8)

        data mwts/2.0,25.4,28.0,32.0,46.0,48.0,62.0,76.0/       ! molec. wts.
        data CpL/3*0.0d0,7.66,4.17,4.45,2.55,2.45/              ! liquid ht cap
        
        hl=0.0
        do i=4,8  
        	hl=hl+(T-100.0)*CpL(i)*x(i)*mwts(i)
        end do
        return
        end      
	
C---------------------------------------------------------------------

	subroutine HVfunc(T,x,hv)
	
c		Calculates vapor enthalpy for TE problem [kJ].
c
c		T = temperature [C]
c		x = vector of molar compositions (length 8).
c		hv= enthalpy of vapor [kJ]

c		Reference is liquid at 100 C.

	double precision T,x(8),hv
	
	double precision mwts(8),CpV(8),dhv(8)

        data mwts/2.0,25.4,28.0,32.0,46.0,48.0,62.0,76.0/       ! molec. wts.
        data CpV/14.6,2.04,1.05,1.85,1.87,2.02,.712,.628/       ! vapor ht cap
        data dhv/3*0.0d0,202.,372.,372.,523.,486./              ! latent heats
                        
        hv=0.0
        do i=1,8  
        	hv=hv+(dhv(i)+(T-100.0)*CpV(i))*x(i)*mwts(i)
        end do
        return
        end      
        
C---------------------------------------------------------------------

	subroutine UTfunc(U,x,T)
	
c		Calculates temperature [C] for specified
c		internal energy and liquid composition.
c		For TE problem.

c		U	Internal energy of liquid [kJ]
c		x	molar composition [kmol]
c		T	calculated temperature [C]

c		The first 3 elements of x are ignored


	double precision U,x(8),T
        
	double precision mwts(8),CpL(8)
                
	data mwts/2.0,25.4,28.0,32.0,46.0,48.0,62.0,76.0/       ! molec. wts.
	data CpL/3*0.0d0,7.66,4.17,4.45,2.55,2.45/              ! liquid ht cap

c		First get heat capacity, then back-calculate
c		temperature.  Assumes heat capacity is
c		independent of temperature (which is the
c		case for the data supplied by Downs and Vogel).

	T=0.0
	do i=4,8
		T=T+x(i)*CpL(i)*mwts(i)
	end do
	T=100.0+U/T      
	return
	end
        
C---------------------------------------------------------------------

	subroutine VPfunc(T,pvap)
	
c		Calculates pure-component vapor pressure [kPa]

c		T	temperature [C]

	double precision T,pvap(8)
	integer i
	
	double precision AVP(8),BVP(8),CVP(8)

	data AVP/3*20.72,20.81,2*21.24,21.32,22.10/ 		! A in Antoine
	data BVP/3*0.0,-1444.0,2*-2114.0,-2748.0,-3318.0/ 	! B in Antoine
	data CVP/3*273.0,259.0,2*266.0,233.0,250.0/		! C in Antoine
        
	do i=1,8
        	pvap(i)=1.d-3*dexp(AVP(i)+
     +			BVP(i)/(CVP(i)+T))
     	end do
     	return
     	end
     	              	                                       
