%*********************************************************
% The namespaces have been normalized. The following
% table shows the attribuation. 
% Normalized Name --> Full Name ---> User-defined Name
% =================================== 
% e0 --> e[0]122185 --> 
%*********************************************************

%*********************************************************
% The variables are named according to the notation
% provided in the Mosaic model.
% 
% The variable names can be read as follows:
% ==========================================
% 	e0_kB#
% 		kB0: Freguency B0
% 	 
% 	e0_CB#
% 		CB0: Feed composition B
% 	 
% 	e0_H_A
% 		H: Respective Heat
% 		Subscript
% 			A: Component A
% 	 
% 	e0_H_B
% 		H: Respective Heat
% 		Subscript
% 			B: Component B
% 	 
% 	e0_T_j
% 		T: Temperature
% 		Subscript
% 			j: Adjustable Temperature
% 	 
% 	e0_T#
% 		T0: Starting Temperature
% 	 
% 	e0_U.A
% 		U.A: Thermal conductivity
% 	 
% 	e0_cp
% 		cp: constant heat capacity
% 	 
% 	e0_rho
% 		rho: Density
% 	 
% 	e0_CA#
% 		CA0: Feed Composition A
% 	 
% 	e0_E_A
% 		E: Energy
% 		Subscript
% 			A: Component A
% 	 
% 	e0_E_B
% 		E: Energy
% 		Subscript
% 			B: Component B
% 	 
% 	e0_F
% 		F: Feed stream
% 	 
% 	e0_R
% 		R: Gas constant
% 	 
% 	e0_V
% 		V: Volume
% 	 
% 	e0_kA#
% 		kA0: Frequency A0
% 	 
% 	e0_C_A
% 		C: Component
% 		Subscript
% 			A: Component A
% 	 
% 	e0_C_B
% 		C: Component
% 		Subscript
% 			B: Component B
% 	 
% 	e0_T
% 		T: Temperature
% 	 
%*********************************************************

function[X,Y]=solveEquationSystem()
    close all;
    clear all;
    clc;
    dbclear if error;
    
    u = 333;             % Init Tj, To
    %u(2) = 333;         % Init To
	% declare interval of independent variable
	X_START=0.0;
	X_END=3600.0;
	X_INTERVAL=linspace(X_START,X_END);	% e0_t

	% specify the initial values for the state variables 
	Y_INIT(1) = 8.0;  	% e0_C_A  
	Y_INIT(2) = 20.0;  	% e0_C_B  
	Y_INIT(3) = 333.0;  	% e0_T
    
    objective = @(u)objFcn(u,X_INTERVAL,Y_INIT');
    lb=273;
    ub=373;
    options=optimoptions('lsqnonlin','TolFun',1e-16,'FinDiffRelStep',1e-16,'FinDiffType','central');
    [X,Y] = lsqnonlin(objective, u, lb, ub, options);

end

function c_A = objFcn(u,X_INTERVAL,Y_INIT)
	% declare parameters 
	PARAMS(1) = 5.0;  	% e0_CA0 
	PARAMS(2) = 69000.0;  	% e0_E_A 
	PARAMS(3) = 72000.0;  	% e0_E_B 
	PARAMS(4) = 6.5E-4;  	% e0_F 
	PARAMS(5) = 8.314;  	% e0_R 
	PARAMS(6) = 1.0;  	% e0_V 
	PARAMS(7) = 5000000.0;  	% e0_kA0 
	PARAMS(8) = 1.0E7;  	% e0_kB0 
	PARAMS(9) = 15.0;  	% e0_CB0 
	PARAMS(10) = 45000.0;  	% e0_H_A 
	PARAMS(11) = -55000.0;  	% e0_H_B 
	PARAMS(12) = 333.0;  	% e0_T_j 
	PARAMS(13) = 333.0;  	% e0_T0 
	PARAMS(14) = 1.4;  	% e0_U.A 
	PARAMS(15) = 3.5;  	% e0_cp 
	PARAMS(16) = 800.0;  	% e0_rho 
    
    
    %ODE_Sol = ode45(@(t,z)updateStates(t,z,x), tSpan, z0);
    [X,Y] = ode23t(@(X,Y)calculateDifferentials(X,Y,PARAMS),X_INTERVAL,Y_INIT');
    c_A = -Y(:,1);
    
	%M = daeSystemMM();

	%OPTIONS = odeset('Mass',M);

	%[X,Y] = ode15s(@(X,Y)daeSystemLHS(X,Y,PARAMS),X_INTERVAL,Y_INIT',OPTIONS);

	displayResults(X,Y);
end
% evaluate the differential function.
function[DYDX] = calculateDifferentials(X,Y,PARAMS)

	% read out variables  
	e0_C_A = Y(1); 
	e0_C_B = Y(2); 
	e0_T = Y(3);
    %e0_C_C= Y(4);

	% read oudbstop if caught errort differential variable
	e0_t = X;

	% read out parameters 
	e0_CA0 = PARAMS(1); 
	e0_E_A = PARAMS(2); 
	e0_E_B = PARAMS(3); 
	e0_F = PARAMS(4); 
	e0_R = PARAMS(5); 
	e0_V = PARAMS(6); 
	e0_kA0 = PARAMS(7); 
	e0_kB0 = PARAMS(8); 
	e0_CB0 = PARAMS(9); 
	e0_H_A = PARAMS(10); 
	e0_H_B = PARAMS(11); 
	e0_T_j = PARAMS(12); 
	e0_T0 = PARAMS(13); 
	e0_U.A = PARAMS(14); 
	e0_cp = PARAMS(15); 
	e0_rho = PARAMS(16); 

	% evaluate the function values  
	DYDX(1) = ( e0_F )/( e0_V ) * ( e0_CA0 - e0_C_A ) - e0_kA0 * e0_C_A * exp(  - ( e0_E_A )/( e0_R * e0_T ) ) + e0_kB0 * e0_C_B * exp(  - ( e0_E_B )/( e0_R * e0_T ) ); 
	DYDX(2) = ( e0_F )/( e0_V ) * ( e0_CB0 - e0_C_B ) - e0_kB0 * e0_C_B * exp(  - ( e0_E_B )/( e0_R * e0_T ) ); 
	DYDX(3) = ( e0_F )/( e0_V ) * ( e0_T0 - e0_T ) + ( e0_U.A )/( e0_rho * e0_cp * e0_V ) * ( e0_T_j - e0_T ) + (  - e0_H_A )/( e0_rho * e0_cp ) * e0_kA0 * e0_C_A * exp(  - ( e0_E_A )/( e0_R * e0_T ) ) + (  - e0_H_B )/( e0_rho * e0_cp ) * e0_kB0 * e0_C_B * exp(  - ( e0_E_B )/( e0_R * e0_T ) ); 
	DYDX=DYDX';

end

function MASS = daeSystemMM()

	MASS = zeros(3, 3); %total number of equations
	MASS(1:3, 1:3)=eye(3,3); % number of odes

end
% evaluate the differential function.
function[DYDX] = daeSystemLHS(X,Y,PARAMS)

	% read out variables  
	e0_C_A = Y(1); 
	e0_C_B = Y(2); 
	e0_T = Y(3); 

	% read out differential variable
	e0_t = X;

	% read out parameters 
	e0_CA0 = PARAMS(1); 
	e0_E_A = PARAMS(2); 
	e0_E_B = PARAMS(3); 
	e0_F = PARAMS(4); 
	e0_R = PARAMS(5); 
	e0_V = PARAMS(6); 
	e0_kA0 = PARAMS(7); 
	e0_kB0 = PARAMS(8); 
	e0_CB0 = PARAMS(9); 
	e0_H_A = PARAMS(10); 
	e0_H_B = PARAMS(11); 
	e0_T_j = PARAMS(12); 
	e0_T0 = PARAMS(13); 
	e0_U.A = PARAMS(14); 
	e0_cp = PARAMS(15); 
	e0_rho = PARAMS(16); 

	% evaluate the function values  
    
	DYDX(1) = ( e0_F )/( e0_V ) * ( e0_CA0 - e0_C_A ) - e0_kA0 * e0_C_A * exp(  - ( e0_E_A )/( e0_R * e0_T ) ) + e0_kB0 * e0_C_B * exp(  - ( e0_E_B )/( e0_R * e0_T ) ); 
	DYDX(2) = ( e0_F )/( e0_V ) * ( e0_CB0 - e0_C_B ) - e0_kB0 * e0_C_B * exp(  - ( e0_E_B )/( e0_R * e0_T ) ); 
	DYDX(3) = ( e0_F )/( e0_V ) * ( e0_T0 - e0_T ) + ( e0_U.A )/( e0_rho * e0_cp * e0_V ) * ( e0_T_j - e0_T ) + (  - e0_H_A )/( e0_rho * e0_cp ) * e0_kA0 * e0_C_A * exp(  - ( e0_E_A )/( e0_R * e0_T ) ) + (  - e0_H_B )/( e0_rho * e0_cp ) * e0_kB0 * e0_C_B * exp(  - ( e0_E_B )/( e0_R * e0_T ) ); 

	DYDX=DYDX';

end
function[] = displayResults(X,Y)

	% decide for a plot type: 
	%   0 	-> Plot the variables into individual figures 
	%   1 	-> Plot into sub figures 
	%   2 	-> Plot all selected into one figure 
	%   other 	-> Do not plot 
	plotType = 1; 

	% set a line width: 
	linewidth = 1.5; 

	% define which dependent variables should be plotted
	%   1 	-> Plot.
	%   other 	-> Do not plot.
	plotControl=[ 
		1	% e0_C_A  
		1	% e0_C_B  
		1	% e0_T  
		];

	% decide wether to normalize the y axis
	%   1 	-> Normalized
	%   other 	-> Individual maximum scale
	axisControl = 2;

	%====================================================

	% labels of the dependent variables
	yAxisLabels=[
		'e0.C_{A}'	% e0_C_A
		'e0.C_{B}'	% e0_C_B
		'e0.T    '	% e0_T
		];
	xAxisLabel = 't';

	% plot the variables 
	figureIndex=1; 
	xMinVal = min(X); 
	xMaxVal = max(X); 
	yMinVal = min(min(Y)); 
	yMaxVal = max(max(Y)); 
	if (plotType==0) 
		% create a plot for each state variable individually 
		for i=1:length(Y(1,:))
			if (plotControl(i)==1)
				figure(i)
				plot(X,Y(:,i),'LineWidth',linewidth)
				title('Solution of the equation system');
				xlabel(xAxisLabel);
				ylabel(yAxisLabels(i,:));
				if (axisControl==1)
					axis([xMinVal xMaxVal yMinVal yMaxVal]);
				end
				legend(yAxisLabels(i,:));
				figureIndex=figureIndex+1;
			end
		end
	elseif (plotType==1)
		% use a subplot environment
		firstTime = 1;
		numberOfPlots = sum(plotControl);
		figure(1)
		for i=1:length(Y(1,:))
			if (plotControl(i)==1)
				subplot(numberOfPlots,1,figureIndex);
				plot(X,Y(:,i),'LineWidth',linewidth)
				ylabel(yAxisLabels(i,:));
				legend(yAxisLabels(i,:));
				if (axisControl==1)
					axis([xMinVal xMaxVal yMinVal yMaxVal]);
				end
				figureIndex=figureIndex+1;
				if (firstTime) 
					title('Solution of the equation system');
					firstTime = 0;
				end 
		end 
		end
		xlabel(xAxisLabel);
	elseif (plotType==2)
		% plot in one figure
		colors = [
			'r'	% -> Red
			'g' % -> Green
			'b' % -> Blue
			'c' % -> Cyan
			'm' % -> Magenta
			'y' % -> Yellow
			'k' % -> Black
			];
		colorCtr=1;
		maxColors=7;
		figure(1)
		hold on;
		for i=1:length(Y(1,:))
			if (plotControl(i)==1)
				linespec = colors(colorCtr);
				if (colorCtr<=maxColors)
					colorCtr = colorCtr+1;
				end
				plot(X,Y(:,i),linespec,'LineWidth',linewidth);
			end
		end
		title('Solution of the equation system');
		xlabel(xAxisLabel);
		ylabel(yAxisLabels(i,:));
		legend(yAxisLabels(:,:));
		hold off;
	end


end