import sys, os, time, math
from scitbx.array_family import flex
import iotbx.pdb

pdb_str = """\
CRYST1   10.000   10.000   10.000  90.00  90.00  90.00 P 1
HETATM  120  CA  CA  A   1       5.000   5.000   5.000  1.00  0.01          CA
TER
END
"""

def run(b_iso = 25.0, element = "Ca"):
  xrs = iotbx.pdb.input(source_info=None, lines=pdb_str).xray_structure_simple()
  xrs = xrs.set_b_iso(value=b_iso)
  xrs.scattering_type_registry(
    table = "n_gaussian",
    d_min = 0.,
    types_without_a_scattering_contribution=["?"])
  rhos = []
  #
  rho_1d_calc = flex.double()
  ed = xrs._scattering_type_registry.gaussian("Ca")
  x = 0
  dist = flex.double()
  while x<5:
    dist.append(x)
    rho_1d_calc.append( ed.electron_density(x, b_iso) )
    x+=0.01
  rhos.append(rho_1d_calc)
  #
  of = open("z","w")
  for i, r in enumerate(dist):
    rhos_i = " ".join(["%12.6f"%rho_i[i] for rho_i in rhos])
    print >> of, "%12.6f %s"%(r, rhos_i)
  of.close()

if (__name__ == "__main__"):
  run()