#!/usr/bin/python

import apfel

xlha = [1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 
		1e-1, 3e-1, 5e-1, 7e-1, 9e-1]

# activate some options
apfel.SetPerturbativeOrder(2);
apfel.SetPDFSet("MRST2004qed");

# initializes integrals on the grids
apfel.InitializeAPFEL()

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)

# tabulate PDFs for the LHA x values
print "alpha_QCD(muF) =", apfel.AlphaQCD(Q)
print "alpha_QED(muF) =", apfel.AlphaQED(Q)


print "x\t","u-ubar\t","d-dbar\t","2(ubr+dbr)\t","c+cbar\t","gluon\t","photon\t\n"

for i in range(2,11):
    print '%e' % xlha[i],
    print '%e' % float(apfel.xPDF(2,xlha[i])-apfel.xPDF(-2,xlha[i])),    
    print '%e' % float(apfel.xPDF(1,xlha[i])-apfel.xPDF(-1,xlha[i])),   
    print '%e' % float(2*(apfel.xPDF(-1,xlha[i])+apfel.xPDF(-2,xlha[i]))),
    print '%e' % float(apfel.xPDF(4,xlha[i])+apfel.xPDF(-4,xlha[i])),
    print '%e' % float(apfel.xPDF(0,xlha[i])),
    print '%e' % float(apfel.xgamma(xlha[i]))