from __future__ import division
from libtbx.utils import Sorry, Usage
import libtbx.phil
import time
import sys

def run (args, out=None) :
  if (out is None) :
    out = sys.stdout
  if (len(args) == 0) or ("--help" in args) :
    if ("--help" in args) : args.remove("--help")
    raise Usage("phaser.find_anomalous_substructure model.pdb data.mtz wavelength=0.9792|energy=12141")
  import phaser.substructure
  import mmtbx.utils
  master_phil = libtbx.phil.parse("""
  include scope mmtbx.utils.cmdline_input_phil_str
  output_file = phaser_AX.pdb
    .type = path
  wavelength = None
    .type = float
  energy = None
    .type = float
  """, process_includes=True)
  cmdline = mmtbx.utils.cmdline_load_pdb_and_data(
    args=args,
    master_phil=master_phil,
    out=out,
    process_pdb_file=False,
    create_fmodel=True,
    prefer_anomalous=True)
  params = cmdline.params
  if ([params.energy, params.wavelength].count(None) in [0,2]) :
    raise Sorry("Please specify either the wavelength (in Angstroms) or the "+
      "energy (in eV).")
  if (not cmdline.fmodel.f_obs().anomalous_flag()) :
    raise Sorry("Anomalous data are required to run this program.")
  if (params.energy is not None) :
    params.wavelength = params.energy / 12398.424468024265
  t1 = time.time()
  result = phaser.substructure.find_anomalous_scatterers(
    fmodel=cmdline.fmodel,
    pdb_hierarchy=cmdline.pdb_hierarchy,
    wavelength=params.wavelength,
    verbose=True,
    log=out)
  t2 = time.time()
  print >> out, "elapsed time: %.1fs" % (t2 - t1)
  f = open(cmdline.params.output_file, "w")
  f.write(result.as_pdb_string())
  f.close()
  print >> out, "Wrote %s" % cmdline.params.output_file

if (__name__ == "__main__") :
  run(sys.argv[1:])