degree of heterogeneity of a reservoir By Dykstra method

#### 
data = read.csv("karpur.csv")

sorted_k = sort(data$k.core, decreasing = TRUE)
head(sorted_k)

## [1] 15600.00 14225.31 13544.98 13033.53 11841.74 11117.40

###
n = length(sorted_k) 
k_distrib = c(1:n)/(n+1)
head(k_distrib)

## [1] 0.001219512 0.002439024 0.003658537 0.004878049 0.006097561 0.007317073

###
scatter.smooth(k_distrib, log(sorted_k, 10))

###
k50 = sorted_k[floor(0.50*n)]
k84.1 = sorted_k[floor(0.841*n)]
print(c(k50, k84.1))

## [1] 1599.2813  395.8943

###
permeability_variation = (k50 - k84.1)/k50
print(permeability_variation)

## [1] 0.7524549

###
d = data$depth
h = d[2:n]-d[1:n-1] 
head(h)

## [1] 0.5 0.5 0.5 0.5 0.5 0.5

###
h = append(h, 0.5)

df = data.frame(h) 
df['h']    = h
k = data$k.core
df['k']    = k
df['phi']  = data$phi.core
df         = df[order(k, decreasing = TRUE),] 

df['kh']   = df['k']*df['h'] 
df['sum_kh'] = cumsum(df['kh'])
df['sum_kh/total'] = df['sum_kh']/df[n, 'sum_kh'] 

df['phih']   = df['phi']*df['h']
df['sum_phih'] = cumsum(df['phih'])
df['sum_phih/total'] = df['sum_phih']/df[n, 'sum_phih']

head(df)

##       h        k     phi       kh   sum_kh sum_kh/total     phih sum_phih
## 809 0.5 15600.00 27.9000 7800.000  7800.00  0.008237623 13.95000 13.95000
## 808 0.5 14225.31 28.5071 7112.657 14912.66  0.015749339 14.25355 28.20355
## 810 0.5 13544.98 27.7563 6772.489 21685.15  0.022901802 13.87815 42.08170
## 807 0.5 13033.53 29.0333 6516.764 28201.91  0.029784192 14.51665 56.59835
## 806 0.5 11841.74 29.5596 5920.872 34122.78  0.036037257 14.77980 71.37815
## 811 0.5 11117.40 27.5865 5558.701 39681.48  0.041907832 13.79325 85.17140
##     sum_phih/total
## 809    0.001241781
## 808    0.002510583
## 810    0.003745967
## 807    0.005038189
## 806    0.006353836
## 811    0.007581664

###
scatter.smooth(
  df['sum_phih/total'][1:n,1], df['sum_kh/total'][1:n,1],
  xlab = "fraction of total volume (phih)", 
  ylab = "fraction of total flow capacity (kh)")

###
x = c(df['sum_phih/total'][1:n,1])
dx = x[2:n] - x[1:n-1]

y1 = c(df['sum_kh/total'][1:n,1])[1:n-1]
y2 = c(df['sum_kh/total'][1:n,1])[2:n]

lower_area = Reduce("+", y1*dx)
upper_area = Reduce("+", y2*dx)

area = (lower_area + upper_area)/2

print(area)

## [1] 0.7247318

###
Lorenz_coefficient = (area - 0.49)/0.49
print(Lorenz_coefficient)

## [1] 0.4790445