Ps=0.4;   % R$/kg-soluto
Px=0.01;  % R$/kg-solvente
F=1e4;    % kg-agua/h
xo=0.02;  % kg-soluto/kg-agua
%s=7e-4;   % kg-solvente/kg-agua
s=0;      % kg-solvente/kg-agua
m=4;      % kg-agua/kg-solvente

a=Ps*F*xo+2*Px*F/m-Px*F*s;
b=Px*F*xo/m;
c=Ps*F;
d=Px*F/m;

%a=130;
%b=0.5;
%c=4000;
%d=25;

x1=[0.0001:0.0001:0.035]';
x2=[0.0001:0.0001:0.025]';
%x1=[0.01:0.001:0.02]';
%x2=[0.006:0.001:0.015]';

[X1,X2]=meshgrid(x1',x2');
[n,m]=size(X1);
L=a-b./X1-c*X2-d*X1./X2;
r1=0.02*ones(n,1);
r2=1e-4*ones(m,1);  % 1e-4 (e nao zero) para aparecer na figura
cs=contour(x1,x2,L,[-5 0 5 10 15 18 19 19.4]);
%cs=contour(x1,x2,L,[-5 0 1 2 3 4 5]);

xlabel('X1');
ylabel('X2');
clabel(cs);
hold on;
plot(x1(1:m),r2,'b:');
plot(r1,x2(1:n),'b:');
p=min(n,m);
plot(x1(1:p),x2(1:p),'b:');
figure(gcf);
hold off;
% cinvert
%figure(2); mesh(X1,X2,L);

[X1s,X2s]=meshgrid(0.005:0.0001:0.022,0.005:0.0001:0.014);
Ls=a-b./X1s-c*X2s-d*X1s./X2s;
figure(3); mesh(X1s,X2s,Ls)
xlabel('X1');
ylabel('X2');
