-
Notifications
You must be signed in to change notification settings - Fork 2
/
HS_solve_phi.m
51 lines (41 loc) · 1.46 KB
/
HS_solve_phi.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
function [phi,phi2] = HS_solve_phi(rho_b,rho_h,rho_f)
% Function which solves for melt fraction using HS+
% given a known bulk resistivity, matrix resistivity, and fluid resistivity
%
% Usage: [phi,i] = HS_solve_phi(rho_b,rho_m,rho_f)
%
% Inputs: rho_b = bulk resistivity (Ohm m)
% rho_h = host rock resistivity (Ohm m)
% rho_f = fluid resistivity (Ohm m)
%
% Outputs: phi = melt fraction (or porosity)
%
% HS+ formula:
%
% sig_b = sig_f + (1-phi)/((1/(sigm-sigf))+phi/(3*sigf))
%
%
% The formula cannot be solved algebraically so the Newton-Raphson method
% is used to solve the equation numerically. Since 0<phi<1, a hard-coded
% starting guess of 0.5 is used.
% Check for input errors first:
flag = 0;
if rho_b > rho_h
disp('Bulk resistivity cannot be greater than matrix resistivity')
flag = 1;
end
if rho_b < rho_f
disp('Bulk resistivity cannot be less than fluid resistivity')
flag = 1;
end
if rho_f > rho_h
disp('Fluid resistivity cannot be greater than matrix resistivity')
flag = 1;
end
sigf = 1./rho_f;
sigm = 1./rho_h;
sigb = 1./rho_b;
phi = (3*sigf.*(sigm-sigb))./(3*sigf.*(sigm-sigf+(1/3)*sigf)+sigm.*sigb-sigf.*sigm-sigb.*sigf);
%phi2 = ((1-(sigb./sigf)).*(3*sigf)-3*sigf+3*sigm)./(3*sigm-3*sigf+(1-sigb./sigf).*(sigf-sigm));
%phi2 = (sigm-sigb)./(sigm-sigf+(1/3)*(sigf-sigm-sigb+(sigb.*sigm)./sigf));
phi2 = (sigm-sigb)./((2/3)*(sigm-sigf)+(1/3)*sigb*(sigm./sigf-1));