from phaser import *
from cctbx import xray
i = InputEP_DAT()
HKLIN = "S-insulin.mtz"
xtalid = "insulin"
waveid = "cuka"
i.setHKLI(HKLIN)
i.addCRYS_ANOM_LABI(xtalid,waveid,"F(+)","SIGF(+)","F(-)","SIGF(-)")
i.setMUTE(True)
r = runEP_DAT(i)
if r.Success():
  hkl = r.getMiller()
  Fpos = r.getFpos(xtalid,waveid)
  Spos = r.getSIGFpos(xtalid,waveid)
  Ppos = r.getPpos(xtalid,waveid)
  Fneg = r.getFneg(xtalid,waveid)
  Sneg = r.getSIGFneg(xtalid,waveid)
  Pneg = r.getPneg(xtalid,waveid)
  i = InputEP_AUTO()
  i.setSPAC_HALL(r.getSpaceGroupHall())
  i.setCELL6(r.getUnitCell())
  i.setCRYS_MILLER(hkl)
  i.addCRYS_ANOM_DATA(xtalid,waveid,Fpos,Spos,Ppos,Fneg,Sneg,Pneg)
  i.setATOM_PDB(xtalid,"S-insulin_hyss.pdb")
  i.setLLGC_COMP(True)
  i.addLLGC_SCAT("S")
  i.addCOMP_PROT_SEQ_NUM("S-insulin.seq",1.)
  i.setHKLO(False)
  i.setXYZO(False)
  i.setMUTE(True)
  i.setHAND("both")
  i.setWAVE(1.54)
  r = runEP_AUTO(i)
  if r.Success():
    print "SAD phasing"
    print "Data read from: " , HKLIN
    print "Data output to : " , r.getMtzFile()
    print "Spacegroup Name (Hall symbol) = %s (%s)" % \
      (r.getSpaceGroupName(), r.getSpaceGroupHall())
    print "Unitcell = " , r.getUnitCell()
    print "LogLikelihood = " , r.getLogLikelihood()
    atom = r.getAtoms()
    print atom.size(), " refined atoms"
    print "%5s %10s %10s %10s %10s %10s" % \
      ("atom","x","y","z","occupancy","u-iso")
    for i in range(0,atom.size()):
      print "%5s %10.4f %10.4f %10.4f %10.4f %10.4f" % \
        (atom[i].scattering_type,atom[i].site[0],atom[i].site[1],atom[i].site[2],atom[i].occupancy,atom[i].u_iso)
    hkl = r.getMiller();
    fwt = r.getFWT()
    phwt = r.getPHWT()
    fom = r.getFOM()
    nrefl = min(10,hkl.size())
    print "First ", nrefl , " reflections"
    print "%4s %4s %4s %10s %10s %10s" % \
      ("H","K","L","FWT","PHWT","FOM")
    for i in range(0,nrefl):
      print "%4d %4d %4d %10.4f %10.4f %10.4f" % \
        (hkl[i][0],hkl[i][1],hkl[i][2],fwt[i],phwt[i],fom[i])
    print "statistics:    %6s %6s %10s %10s" % \
      ("lores","hires","fig. merit","num refls")
    for hand in [False, True]:
      print "hand:", hand
      print "               %6.2f-%6.2f %10.4f %10d" % \
        (r.stats_lores(hand),
         r.stats_hires(hand),
         r.stats_fom(hand),
         r.stats_num(hand))
      print
  else:
    print "Job exit status FAILURE"
    print r.ErrorName(), "ERROR :", r.ErrorMessage()
else:
  print "Job exit status FAILURE"
  print r.ErrorName(), "ERROR :", r.ErrorMessage()