Madelung Constant

import math


proc simulateIon(x: float, y: float, z: float): float =
  let charge = math.pow(-1.0, x + y + z)  # NaCl - FCC 
  let distance = math.sqrt(x * x + y * y + z * z) 

  return charge / distance


proc madelung(N): float =
  var M: float = 0

  for x in -N..N + 1:
    for y in -N..N + 1:
      for z in -N..N + 1:
        if x != 0 or y != 0 or z != 0:
          M += simulateIon(float(x), float(y), float(z))

  return M


const N = 80
var mValues: seq[float] = @[]

# for i in N - 3 .. N + l:
for i in 0..N:
  if i == 0:
    mValues.add(0.0)
  else:
    mValues.add(madelung(i))
    
    let mean = (mValues[i] + mValues[i-1]) / 2
    
    # echo("M(", i, ") = ", mValues[i], "  M(", i-1, ") = ", mValues[i-1], "  Mean = ", mean)
    echo(i, ", ", mValues[i], ", ", mean)