  Iterative Density Porosity (Python 3.x)   Algorithm for computing total porosity (Phi_t) in formations with variable matrix density. The "classic" density porosity algorithm presumes that matrix density, $$\rho _{\mathit{ma}}$$, remains constant within the lithofacies: $\phi _t=\frac{\rho _{\mathit{ma}}-\rho _b}{\rho _{\mathit{ma}}-\rho _f},\rho _{\mathit{ma}}=\mathit{Const}$ In many cases, however, such matrix density can be a function of total porosity: $\phi _t=\frac{\rho _{\mathit{ma}}(\phi _t)-\rho _b}{\rho _{\mathit{ma}}(\phi _t)-\rho _f}$ This algorithm resolves the above equation numerically. It must be stressed that the dependency between the matrix density and total porosity is empirical and relies on laboratory measurements performed on full-bore or side-wall cores.   Download source code: Iterative_Density_Porosity.py   import numpy as np def Iterative_Density_Porosity( Rho_b, Rho_fl, Rho_ma=lambda p: 2.65, accuracy=0.001, max_iter=100, rho_ma_initial=3.0, verbose=False): """ Solves the density iteratively into total porosity Rho_b - array-like, bulk density, in g/cc Rho_fl - array-like, density of fluid in flushed zone, in g/cc Rho_ma - equation for matrix density, dependency from total porosity Phi_t accuracy = 0.001 - desired solution accuracy, in v/v max_iter = 100 - maximum number of iterations rho_ma_initial = 3.0 - initial matrix density, g/cc verbose - print intermediate results returns: total porosity Phi_t as an array, in v/v estimated matrix density as an array, in g/cc """ density = np.array(Rho_b) l = len(density) fluid = np.array(Rho_fl) if len(fluid) != l: fluid = np.ones(l) * fluid Phi_t = np.clip( (rho_ma_initial-density) / (rho_ma_initial-fluid), 0, 1) Rho_ma_est = np.ones(l) * rho_ma_initial for i, d in enumerate(density): f = fluid[i] p = Phi_t[i] for j in range(max_iter): Rho_ma_est[i] = Rho_ma(p) Phi_t[i] = np.clip((Rho_ma_est[i] - d) / (Rho_ma_est[i] - f), 0, 1) if abs( Phi_t[i] - p) >= accuracy: p = Phi_t[i] continue if verbose: print("Rho_b={:.3f}, Rho_f={:.3f} : Phi_t={:.3f}, Rho_ma={:.3f}" " (after {:d} iterations)".format( d, f, Phi_t[i], Rho_ma_est[i], j)) break return Phi_t, Rho_ma_est if __name__ == "__main__": print("Polymictic feldspar sandstone with matrix density (from core)") print("\t2.56 + Phi_t*0.55:") print("Flushed zone fluid density is assumed 1.03 g/cc;") print("Bulk density is environmentally corrected;") print("Rho_b\t\tPhi_t(classic)\tPhi_t(iterative)") Rho_b = np.linspace(2.0, 2.6, 12) Phi_t_classic = (2.65-Rho_b) / (2.65-1) Phi_t, Rho_ma_apparent = Iterative_Density_Porosity( Rho_b, [1.03], Rho_ma=lambda p: np.clip(2.56+p*0.55, 2.56, 2.67), verbose=False) for rhob, phit_c, phit in zip(Rho_b, Phi_t_classic, Phi_t): print("{:.3f}\t\t{:.3f}\t\t{:.3f}".format(rhob, phit_c, phit)) (C) M.Yakimov, 2006 