#!/usr/bin/env python

import apfel

# activate some options
apfel.SetPerturbativeOrder(1);
apfel.SetNumberOfGrids(3);
apfel.SetGridParameters(1,130,3,1e-9);
apfel.SetGridParameters(2,60,5,1e-1);
apfel.SetGridParameters(3,20,5,8e-1);

# Initializes integrals on the grids
apfel.InitializeAPFEL();

# Evolve PDFs on the grids
eps = 1e-10;

Q20 = float(raw_input("Enter initial scale in GeV^2: ")) 
Q2  = float(raw_input("Enter final scale in GeV^2: "))

Q0 = Q20**0.5 - eps
Q =  Q2**0.5

apfel.EvolveAPFEL(Q0,Q)

momsr = 0.0;
for i in range(-6,7):
    momsr += apfel.NPDF(i,2);
momsr += apfel.Ngamma(2);

uvsr = apfel.NPDF(2,1) - apfel.NPDF(-2,1);
dvsr = apfel.NPDF(1,1) - apfel.NPDF(-1,1);
svsr = apfel.NPDF(3,1) - apfel.NPDF(-3,1);

print "Sum rules at Q =",Q,"GeV:"
print "- Momentum sum rule         = ", momsr
print "- Up valence sum rule       = ", uvsr
print "- Down valence sum rule     = ", dvsr
print "- Strange valence sum rule  = ", svsr