Using iphreeqc module with MATLAB on Linux
4 views (last 30 days)
Show older comments
Dear MATLAB users,
Good day!
Has anyone used iphreeqc module with MATLAB on Linux?
We can simply use iphreeqc COM module with MATLAB on Window, but I dont' know how to do the same job using iphreeqc module with MATLAB on Linux.
I want two things: 1) to create an iphreeqc instance in MATLAB and 2) to link to iphreeqc library.
1) There are instrutions to create an iphreeqc instance in C, C++, and Fortran programs.
For exmple, we can create an iphreeqc instance in C as below:
#include "IPhreeqc.h"
int id;
id = CreateIPhreeqc();
How can I do the same job using MATLAB?
2) Below's the list of library that I installed:
Libraries:
/usr/local/lib/libiphreeqc.a (Static library)
/usr/local/lib/libiphreeqc.so (Shared object library)
/usr/local/lib/libiphreeqc.la (libtool library file)
How can I link any of them to MATLAB?
Thanks in advance!
Sincerely,
Minjeong
0 Comments
Answers (1)
Scott Smith
on 13 Dec 2019
I would like to do this as well !
I have a work around. Just use fprintf to write the PHREEQC input file and then use ! to run PHREEQC and inport the results.
this makes a simple pe pH diagram for hydrous ferric oxide. I modified a local version of the PHREEQC database to have constant pH and pe "phases". So this won't run as I've written it for you I don't think but it gives you an idea of how to call PHREEQC from matlab.
I find it is pretty fast. Code also work on Windows but much much slower. I use iPHREEQC when I am on my windows computer.
anyway, here is the code example ...
% make pH pe diagram for iron salts in water
function makepHpe_simpleFesystem
figure(1); clf
pH=2:1:12; pe=-15:2:22; FeT=1e-5;
% only perform calculations if within the stability field of water
for i=1:size(pH,2)
for j=1:size(pe,2)
HFO(i,j)=NaN;
if pe(j)>=-pH(i)
if pe(j)<=21.6-pH(i)
HFO(i,j)=returnFespeciation(FeT,pH(i),pe(j));
end
end
end
end
surf(pH,pe,HFO'./FeT)
set(gca,'Clim',[0 1])
set(gca,'linewidth',2,'fontsize',10)
axis([min(pH) max(pH) min(pe) max(pe) 0 1.1]); view([0 90]);
hold on
plot(pH,-pH,'k','linewidth',2)
plot(pH,21.6-pH,'k','linewidth',2)
title('HFO PHREEQC database')
xlabel('pH'); ylabel('pe')
colorbar
end
function [HFO]=returnFespeciation(FeT,pH_val,pe_val)
temp_val=25;
fileID=fopen('runphreeqc.txt','w');
fprintf(fileID,'SOLUTION 1 Pure water\n');
fclose(fileID);
fileID=fopen('runphreeqc.txt','a');
pHtxt=['pH ',num2str(pH_val),'\n']; fprintf(fileID,pHtxt); % the initial pH condition
petxt=['pe ',num2str(pe_val),'\n']; fprintf(fileID,petxt); % the initial pe condition
temptxt=['temp ',num2str(temp_val),'\n']; fprintf(fileID,temptxt); % the temperature
fprintf(fileID,'-units mol/L\n'); % the unit of the input; usually mol/L is used
FeTtxt=['Fe(+3) ' num2str(FeT),'\n']; fprintf(fileID,FeTtxt); % total Fe
ClTtxt=['Cl ' num2str(FeT*3),'\n']; fprintf(fileID,ClTtxt); % total Cl (assume adding FeCl3)
fprintf(fileID,'EQUILIBRIUM_PHASES 1\n');
fprintf(fileID,' Fe(OH)3(a) 0.0 10 precipitate_only\n'); % Fe(OH)3 will precipitate out when the saturaiton index is greater than 0
fixpHtxt=[' pH_Fix -',num2str(pH_val),' HCl 10.0\n']; fprintf(fileID,fixpHtxt);
fprintf(fileID,'-force_equality true\n');
fixpetxt=[' pe_Fix ',num2str(-1*pe_val),' O2 \n']; fprintf(fileID,fixpetxt);
fprintf(fileID,'-force_equality true\n');
fprintf(fileID,'SELECTED_OUTPUT\n');
fprintf(fileID,'-file selected.out\n');
fprintf(fileID,'-selected_out true\n');
fprintf(fileID,'-user_punch true\n');
fprintf(fileID,'-reset false\n');
fprintf(fileID,'-simulation false\n');
fprintf(fileID,'-state false\n');
fprintf(fileID,'-distance false\n');
fprintf(fileID,'-time false\n');
fprintf(fileID,'-step false\n');
fprintf(fileID,'-ph true\n');
fprintf(fileID,'-pe true\n');
fprintf(fileID,'-reaction false\n');
fprintf(fileID,'-temperature false\n');
fprintf(fileID,'-alkalinity false\n');
fprintf(fileID,'-ionic_strength false\n');
fprintf(fileID,'-water false\n');
fprintf(fileID,'-charge_balance false\n');
fprintf(fileID,'-percent_error false\n');
fprintf(fileID,'-equilibrium_phases Fe(OH)3(a)');
fclose(fileID);
%in LINUX
!phreeqc runphreeqc.txt out.txt PHREEQC.dat
%can see output of that one - so track error messages
%str=['!phreeqc runphreeqc.txt out.txt PHREEQC.dat']; evalc(str); % so no screen output
str=['!tail -n +2 "selected.out" > "selected.tmp" && mv "selected.tmp" "selected.out"']; evalc(str); load selected.out; % way faster only works in linux
%t= importdata('selected.out', '\t', 1); selected=t.data; % works in linux or windows but slower
selected=selected(2,:);
pH=selected(1)
pe=selected(2)
HFO=selected(4)
%check pH was actually fixed. it is a little off without pH fixed phase
pHtst=pH-pH_val
petst=pe-pe_val
end
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!