from phaser import *
from cctbx import xray
i = InputEP_DAT()
HKLIN = "S-insulin.mtz"
xtalid = "insulin"
waveid = "sad"
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)
  Fmean = r.getF(xtalid,waveid)
  Smean = r.getSIGF(xtalid,waveid)
  i = InputANO()
  i.setSPAC_HALL(r.getSpaceGroupHall())
  i.setCELL(r.getUnitCell())
  i.setREFL(hkl,Fmean,Smean)
  i.addCOMP_PROT_FASTA_NUM("S-insulin.seq",1.)
  i.setROOT("insulin_ano")
  i.setMUTE(True)
  del(r)
  r = runANO(i)
  if r.Success():
    print "Anisotropy Correction"
    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 "Principal components = " , r.getEigenBs()
    print "Range of principal components = " , r.getAnisoDeltaB()
    print "Wilson Scale = " , r.getWilsonK()
    print "Wilson B-factor = " , r.getWilsonB()
    Fpos = r.getScaledCorrected(hkl,Fpos)
    Spos = r.getScaledCorrected(hkl,Spos)
    Fneg = r.getScaledCorrected(hkl,Fneg)
    Sneg = r.getScaledCorrected(hkl,Sneg)
#   print r.logfile()
    i = InputEP_SAD()
    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_CRYS_COMPLETE(xtalid,True)
    i.addLLGC_CRYS_SCAT_ELEMENT(xtalid,"S")
    i.setHKLO(True)
    i.setSCRI(False)
    i.setXYZO(False)
    i.setMUTE(True)
    i.setROOT("insulin_sad_split")
    del(r)
    r = runEP_SAD(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(xtalid)
      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()
      phib = r.getPHIB()
      fom = r.getFOM()
      nrefl = min(10,hkl.size())
      print "First ", nrefl , " reflections"
      print "%4s %4s %4s %10s %10s %10s %10s" % \
        ("H","K","L","FWT","PHWT","PHIB","FOM")
      for i in range(0,nrefl):
        print "%4d %4d %4d %10.4f %10.4f %10.4f %10.4f" % \
          (hkl[i][0],hkl[i][1],hkl[i][2],fwt[i],phwt[i],phib[i],fom[i])
      print "statistics:    %6s %6s %10s %10s" % \
        ("lores","hires","fig. merit","num refls")
      for i in range(0,r.stats_numbins()):
        print "by bin:     %d %6.2f-%6.2f %10.4f %10d" % \
          (i+1,r.stats_by_bin_lores(i),r.stats_by_bin_hires(i),r.stats_by_bin_fom(i),r.stats_by_bin_num(i))
      print "total         %6.2f-%6.2f %10.4f %10d" % \
        (r.stats_lores(),r.stats_hires(),r.stats_fom(),r.stats_num())
#     print r.logfile()
    else:
      print "Job exit status FAILURE"
      print r.ErrorName(), "ERROR :", r.ErrorMessage()
  else:
    print "Job exit status FAILURE"
    print r.ErrorName(), "ERROR :", r.ErrorMessage()
else:
  print "Job exit status FAILURE"
  print r.ErrorName(), "ERROR :", r.ErrorMessage()