#!/usr/bin/env /opt/exp_soft/enmr/CIRMMP/xplor-nih/2.21/64/bin/pyXplor

#
# calculate the solution X-ray scattering curve for one (or more) pdb files

import sys
usage=r"""
 usage: %s [option] <pdb files>

 where options are zero or more of:
   -psf=psf file
   -deltaQ=value
   -numQ=value
   -angleN=value
 
 deltaQ and numQ specify the values at which I(q) is calculated.
 
 angleN controls how many solid angle values are used in the
 approximation. 500 is the default, and probably sufficient for most
 SAXS experiments.
 
 If more than one pdb file is specified, the average is calculated.
""" % sys.argv[0]

#

deltaQ=0.01  #these parameters define the values of q for which I is calc'd
numQ = 100
angleN = 500 # this specifies the granularity of the solid-angle grid

(opts,pdbFiles) = xplor.parseArguments(('psf:1',
                                        'deltaQ:1','numQ:1',
                                        'angleN:1',
                                        'help-script:0'))
psf=None

for opt in opts:
    if opt[0]=='psf':
        psf=opt[1]
        pass
    if opt[0]=='deltaQ':
        deltaQ=float(opt[1])
        pass
    if opt[0]=='numQ':
        numQ=int(opt[1])
        pass
    if opt[0]=='angleN':
        angleN=int(opt[1])
        pass
    if opt[0]=='help-script':
        print usage
        sys.exit(0)
    pass

import protocol
if psf:
    protocol.initStruct(psf)
    protocol.initCoords(pdbFiles[0])
else:
    protocol.loadPDB(pdbFiles[0])
    pass

xplor.simulation.deleteAtoms("not known")

expt=[]
q=0
for i in range(numQ):
    q = i*deltaQ
    I = 1 #dummy value
    expt.append( (q,I) )
    pass

import solnXRayPotTools
#use this for large-angle DNA calculations
#solnXRayPotTools.solventVolume = solnXRayPotTools.solventVolumeSets['tiede']
from solnXRayPotTools import create_solnXRayPot

scat = create_solnXRayPot('scat',"not hydro",expt)
scat.setNumAngles(angleN)

#scat.setCalcType('n2')  #uncomment to get the exact result



calcd=0
for pdbFile in pdbFiles:
    protocol.initStruct(erase=True)
    protocol.loadPDB(pdbFile)
    xplor.command("delete sele=(not known) end")
    #scat.calcGlobCorrect('n2')

    #FIX: modified not working for Ne=1!
    #scat.calcEnergy()

    if not calcd:
        calcd = scat.calcd()
    else:
        calcd += scat.calcd()
        pass

    pass

calcd.scale( 1./len(pdbFiles) )


for i in range(len(scat.qValues())):
    print scat.qValues()[i], calcd[i]/scat.globCorrect[i], scat.expt()[i],\
          scat.globCorrect[i]
    pass
#trace.resume()