Skip to content

Commit

Permalink
uploading the first working version
Browse files Browse the repository at this point in the history
  • Loading branch information
simulkade committed Apr 29, 2015
1 parent 6d53e54 commit 67714f8
Show file tree
Hide file tree
Showing 40 changed files with 2,690 additions and 3 deletions.
73 changes: 73 additions & 0 deletions ActivityModels/Margules2.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
function [gErt, gama] = Margules2(temperature, x, component, BIP)
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

%{
Copyright (c) 2012, 2013, Ali Akbar Eftekhari
All rights reserved.
Redistribution and use in source and binary forms, with or
without modification, are permitted provided that the following
conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above
copyright notice, this list of conditions and the following
disclaimer in the documentation and/or other materials provided
with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%}


%NOT READY YET
% Margules 2 suffix activity coefficient method gives gE/RT in []
% Temperature in K
% A in J/mole (main diagonal of the matrix should be 1)

A = [BIP.NRTLcons]+[BIP.NRTLtdep]*temperature;
alfa = BIP.NRTLalfa;
R = 8.314; % J/mole/K
N = length(x);
mult1 = ones(N)-diag(ones(N,1));
taw = (A/R/temperature) .* mult1;
G = (exp(-alfa*taw)) .* mult1;

if (N==2)
gErt = x(1)*x(2)*(taw(2,1)*G(2,1)/(x(1)+x(2)*G(2,1))+ ...
taw(1,2)*G(1,2)/(x(2)+x(1)*G(1,2)));
gama(1) = exp(x(2)^2*(taw(2,1)*(G(2,1)/(x(1)+x(2)*G(2,1)))^2+ ...
taw(1,2)*G(1,2)/(x(2)+x(1)*G(1,2))^2));
gama(2) = exp(x(1)^2*(taw(1,2)*(G(1,2)/(x(2)+x(1)*G(1,2)))^2+ ...
taw(2,1)*G(2,1)/(x(1)+x(2)*G(2,1))^2));
else
gErt = x*( x*(G.*taw)./ (x*G) )';

part1 = x*(G.*taw)./ (x*G);
part2 = (x.*(1./(x*G)))*(G.*taw)';
part3 = (x.*(1./(x*G)).*(1./(x*G)).*(x*(taw.*G)))*G';
gama = exp(part1 +part2-part3);
end
73 changes: 73 additions & 0 deletions ActivityModels/NRTL.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
function [gErt, gama] = NRTL(temperature, x, component, BIP)
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

%{
Copyright (c) 2012, 2013, Ali Akbar Eftekhari
All rights reserved.
Redistribution and use in source and binary forms, with or
without modification, are permitted provided that the following
conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above
copyright notice, this list of conditions and the following
disclaimer in the documentation and/or other materials provided
with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%}


% NRTL activity coefficient method gives gE/RT in []
% Temperature in K
% A in J/mole (main diagonal of the matrix should be 1)
% alfa in []
A = [BIP.NRTLcons]+[BIP.NRTLtdep]*temperature+[BIP.NRTLtdep2]*temperature*temperature+ ...
[BIP.NRTLtdepm1]/temperature+[BIP.NRTLtdeplog]*log(temperature);
alfa = BIP.NRTLalfa;
R = 8.314; % J/mole/K
N = length(x);
mult1 = ones(N)-diag(ones(N,1));
taw = (A/R/temperature) .* mult1;
G = (exp(-alfa*taw)) .* mult1;

if (N==2)
gErt = x(1)*x(2)*(taw(2,1)*G(2,1)/(x(1)+x(2)*G(2,1))+ ...
taw(1,2)*G(1,2)/(x(2)+x(1)*G(1,2)));
gama(1) = exp(x(2)^2*(taw(2,1)*(G(2,1)/(x(1)+x(2)*G(2,1)))^2+ ...
taw(1,2)*G(1,2)/(x(2)+x(1)*G(1,2))^2));
gama(2) = exp(x(1)^2*(taw(1,2)*(G(1,2)/(x(2)+x(1)*G(1,2)))^2+ ...
taw(2,1)*G(2,1)/(x(1)+x(2)*G(2,1))^2));
else
gErt = x*( x*(G.*taw)./ (x*G) )';

part1 = x*(G.*taw)./ (x*G);
part2 = (x.*(1./(x*G)))*(G.*taw)';
part3 = (x.*(1./(x*G)).*(1./(x*G)).*(x*(taw.*G)))*G';
gama = exp(part1 +part2-part3);
end
81 changes: 81 additions & 0 deletions ActivityModels/UNIQUAC.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,81 @@
function [gErt, gama] = UNIQUAC(temperature, x, component, BIP)
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

%{
Copyright (c) 2012, 2013, Ali Akbar Eftekhari
All rights reserved.
Redistribution and use in source and binary forms, with or
without modification, are permitted provided that the following
conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above
copyright notice, this list of conditions and the following
disclaimer in the documentation and/or other materials provided
with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%}


N = length(x);
R = 8.314;
z = 10; %coordination number
r= [component.uniquacR];
q= [component.uniquacQ];
A = [BIP.UNIQUACcons]+[BIP.UNIQUACtdep]*temperature+[BIP.UNIQUACtdep2]*temperature*temperature+ ...
[BIP.UNIQUACtdepm1]/temperature+[BIP.UNIQUACtdeplog]*log(temperature);
mult1 = ones(N)-diag(ones(N,1));
fay = (x.*r)/(x*r');
teta = (x.*q)/(x*q');
taw = exp(-A/R/temperature).*mult1;
l = z/2*(r-q)-(r-1);
if (N==2)
lngama = zeros(1,N);
gEcrt = x(1)*log(fay(1)/x(1))+x(2)*log(fay(2)/x(2))+ ...
z/2*(q(1)*x(1)*log(teta(1)/fay(1))+q(2)*x(2)*log(teta(2)/fay(2)));
gErrt = -q(1)*x(1)*log(teta(1)+teta(2)*taw(2,1))- ...
q(2)*x(2)*log(teta(2)+teta(1)*taw(1,2));
gErt = gEcrt + gErrt;
for i=1:N
for j=1:N
if (i~=j)
lngama(i) = log(fay(i)/x(i))+z/2*q(i)*log(teta(i)/fay(i)) ...
+ fay(j)*(l(i)-r(i)/r(j)*l(j))-q(i)*log(teta(i)+teta(j)*taw(j,i)) ...
+ teta(j)*q(i)*(taw(j,i)/(teta(i)+teta(j)*taw(j,i))- ...
taw(i,j)/(teta(j)+teta(i)*taw(i,j)));
end
end
end
gama = exp(lngama);
else
gErt = x*log(fay./x)' + z/2*x*(q.*log(teta./fay))'-x*(q.*log(teta*taw'))';
lngama = log(fay./x) + z/2*q.*log(teta./fay) + l - fay./x*(x*l') ...
- q.*(teta*taw') + q - q.*((1./(teta*taw))*(teta*taw')');
gama = exp(lngama);
end
66 changes: 66 additions & 0 deletions ActivityModels/Wilson.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
function [gErt, gama] = Wilson(temperature, x, component, BIP)
%NOT READY YET
% Margules 2 suffix activity coefficient method gives gE/RT in []
% Temperature in K
% A in J/mole (main diagonal of the matrix should be 1)
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

%{
Copyright (c) 2012, 2013, Ali Akbar Eftekhari
All rights reserved.
Redistribution and use in source and binary forms, with or
without modification, are permitted provided that the following
conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above
copyright notice, this list of conditions and the following
disclaimer in the documentation and/or other materials provided
with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%}


A = [BIP.Wilsoncons]+[BIP.Wilsontdep]*temperature;
R = 8.314; % J/mole/K
N = length(x);

if (N==2)
gErt = -x(1)*log(x(1)+A(1,2)*x(2))-x(2)*log(x(2)+A(2,1)*x(1));
gama(1) = exp(-log(x(1)+A(1,2)*x(2))+x(2)*(A(1,2)/(x(1)+A(1,2)*x(2)) ...
-A(2,1)/(A(2,1)*x(1)+x(2))));
gama(2) = exp(-log(x(2)+A(2,1)*x(1))-x(1)*(A(1,2)/(x(1)+A(1,2)*x(2)) ...
-A(2,1)/(A(2,1)*x(1)+x(2))));
else
gErt = x*log(x*A')';

part1 = -log(x*A')+1;
part2 = -(x.*(1./(x*A')))*A;
gama = exp(part1+part2);
end
Loading

0 comments on commit 67714f8

Please sign in to comment.