Ported Code

BeginPackage["KKP`"];
(* Written by William Horowitz (horowitz atsym phys . columbia . edu), 12/23/07. *) (* Ported from kkp.f found at http://www.desy.de/~poetter/kkp.html. *)
(* Corrected for the error in the LO kaon contribution in kkp.f on lines 658-661 in the original code in which the up and down quark contributions were switched. *)
KKP::"usage"="KKP(ih,iset,x,qs).
ih, iset, x, qs are input; KKP[[dh+1]] is output.
ih   = 1 : (pi^+ + pi^-)  /2
ih   = 2 : (K^+ + K^-)    /2   
ih   = 3 : (K^0 + K^0_bar)/2   
ih   = 4 : (p + p_bar)    /2
ih   = 5 : (pi^0)
ih   = 6 : (n + n_bar)    /2
ih   = 7 : (h^+ + h^-)         [as sum of pions, kaons and protons]

iset = 0 : LO         
iset = 1 : NLO

x    = longitudinal-momentum fraction
qs   = fragmentation scale (in GeV)
c
Parton label:
0    1    2    3    4    5    6    7    8     9    10
g    u   ubar  d   dbar  s   sbar  c   cbar   b   bbar

Lambda_QCD (in GeV):
0.088 in LO
0.213 in NLO

Note that qs is in GeV.";

KKP::xarg="The input `1` is outside the recommended range of .1≤x≤.8.  Result may be off by more than 5%.";
KKP::Qsarg="The input `1` is outside the recommended range of (Subscript[Q, 0])=`2`≤(Subscript[Q, s])≤100 GeV for dh=`3`.  Result may be off by more than 5%.";
KKP::input1="Incorrect input.  Either ih, iset, or dh is not an integer or x or Qs is not numeric.";
KKP::input2="Incorrect input.";

Begin["privateKKP`"];
(* Mass-thresholds: *)
rmcc=2.9788; (* GeV *)
rmbb=9.46037; (* GeV *)

Dform[blist_,alist_,x_,s_]:=(blist[[1]]+alist[[1]] s+alist[[2]] (s^2)+alist[[3]] (s^3)) (x^(blist[[2]]+alist[[4]] s+alist[[5]] (s^2)+alist[[6]] (s^3))) (((1.-x))^(blist[[3]]+alist[[7]] s+alist[[8]] (s^2)+alist[[9]] (s^3))) (1.+((alist[[10]] s+alist[[11]] (s^2))/x));
(* needed only for the LO proton gluon contribution *)
DformprgLO[blist_,alist_,x_,s_]:=(blist[[1]]+alist[[1]] s+alist[[2]] (s^2)+alist[[3]] (s^3)) (x^(blist[[2]]+alist[[4]] s+alist[[5]] (s^2)+alist[[6]] (s^3))) (((1.-x))^(blist[[3]]+alist[[7]] s+alist[[8]] (s^2)+alist[[9]] (s^3))) (1.+((alist[[10]] s+alist[[11]] (s^2)+alist[[12]] (s^3))/x));

ΛqcdLO = 0.088; (* GeV *)
ΛqcdNLO = 0.213; (* GeV *)
Q0[dh_]:=Which[
MemberQ[Range[0,6],dh],(Sqrt[2.]),
MemberQ[{7,8},dh],rmcc,
MemberQ[{9,10},dh],rmbb
]; (* GeV *)
Q[Qs_,dh_]:=If[Qs>Q0[dh],Qs,Q0[dh]]; (* Q in GeV *)
(* LO *)
sLO[Qs_,dh_]:=Log[((Log[(((Q[Qs,dh])^2)/((ΛqcdLO)^2))])/(Log[(((Q0[dh])^2)/((ΛqcdLO)^2))]))];
(* NLO *)
sNLO[Qs_,dh_]:=Log[((Log[(((Q[Qs,dh])^2)/((ΛqcdNLO)^2))])/(Log[(((Q0[dh])^2)/((ΛqcdNLO)^2))]))];
(********************)
(*********LO*********)
(********************)
(* LO PION *)
(* gluon *)
bpgLO={6.0451, -0.71378, 2.92133};
apgLO={-6.61523, -1.64978, 2.68223, 0.14705, -1.08423, -0.43182, 1.48429, 1.32887, -1.78696, 0.23086, -0.29182};
dpgLO[x_,Qs_]:=Dform[bpgLO,apgLO,x,sLO[Qs,0]];
(* up *)
bpuLO={0.5461, -1.46616, 1.01864};
apuLO={-0.22946, -0.22594, 0.21119, -0.45404, -0.12684, 0.27646, 0.95367, -1.09835, 0.74657, -0.01877, 0.02949};
dpuLO[x_,Qs_]:=Dform[bpuLO,apuLO,x,sLO[Qs,1]];
(* strange *)
bpsLO={22.2815, 0.12732, 6.13697};
apsLO={-20.8125, -11.5725, 15.5372, 0.23075, -2.71424, 1.72456, 2.18849, -5.04475, 3.29117, 0.09044, -0.07589};
dpsLO[x_,Qs_]:=Dform[bpsLO,apsLO,x,sLO[Qs,5]];
(* charm *)
bpcLO={8.755, -0.38611, 5.61846};
apcLO={-9.32277, 1.806, 2.02179, -0.4119, -0.48496, 0.42525, 0.74035, -0.64929, 0.66788, 0.06652, -0.05531};
dpcLO[x_,Qs_]:=Dform[bpcLO,apcLO,x,sLO[Qs,7]];
(* bottom *)
bpbLO={0.31147, -1.92993, 3.47086};
apbLO={-0.19319, -0.10487, 0.18824, -0.44692, -0.08271, 0.30441, 0.79775, -0.28091, 0.39504, -0.04887, 0.03212};
dpbLO[x_,Qs_]:=Dform[bpbLO,apbLO,x,sLO[Qs,9]];
(* LO KAON *)
(* gluon *)
bkgLO={0.02862, -2.94091, 2.73474};
akgLO={-0.02113, 0.00389, 0.00901, 0.66881, -0.2967, 0.20574, -0.58222, 0.04329, 0.78033, 0.03586, -0.0122};
dkgLO[x_,Qs_]:=Dform[bkgLO,akgLO,x,sLO[Qs,0]];
(* up *)
bkuLO={0.25937, -0.61925, 0.85946};
akuLO={-0.10502, 0.00572, -0.00269, 0.09956, 0.07389, -0.0007, 0.57965, 0.26397, -0.12764, 0.15303, 0.14807};
dkuLO[x_,Qs_]:=Dform[bkuLO,akuLO,x,sLO[Qs,1]];
(* down *)
bkdLO={5.38115, -0.00321, 3.07632};
akdLO={-3.05084, -1.10056, 1.31207, -0.25889, -0.18494, 0.13994, 1.13745, -0.90413, 0.56581, 0.05141, -0.00697};
dkdLO[x_,Qs_]:=Dform[bkdLO,akdLO,x,sLO[Qs,3]];
(* charm *)
bkcLO={5.18266, -0.17751, 4.30306};
akcLO={-3.48519, -1.00982, 1.17996, 0.02309, -0.61327, -0.03532, 1.00547, -0.51779, 0.20683, 0.13514, -0.17778};
dkcLO[x_,Qs_]:=Dform[bkcLO,akcLO,x,sLO[Qs,7]];
(* bottom *)
bkbLO={1.57044, -0.84143, 6.01488};
akbLO={-1.7834, 0.571, 0.15469, -0.43448, -0.05314, -0.36621, 0.72953, -0.64433, 0.92351, 0.01024, -0.0616};
dkbLO[x_,Qs_]:=Dform[bkbLO,akbLO,x,sLO[Qs,9]];
(* LO PROTON *)
(* gluon *)
bprgLO={0.73953, -0.76986, 7.69079};
aprgLO={-1.64519, 1.01189, -0.10175, -3.58787, 13.8025, -13.8902, -2.8447, -0.36719, -2.21825, 1.26515, -1.96117,0.54769};
dprgLO[x_,Qs_]:=DformprgLO[bprgLO,aprgLO,x,sLO[Qs,0]];
(* up *)
bpruLO={0.40211, -0.85973, 2.8016};
apruLO={-0.21633, -0.07045, 0.07831, 0.13987, -0.82412, 0.43114, 0.78923, -0.05344, 0.0146, 0.05198, -0.04623};
dpruLO[x_,Qs_]:=Dform[bpruLO,apruLO,x,sLO[Qs,1]];
(* strange *)
bprsLO={4.07885, -0.09735, 4.99191};
aprsLO={-2.97392, -0.92973, 1.23517, 0.25834, -1.52246, 0.7706, 1.14379, -0.8532, 0.45607, 0.07174, -0.08321};
dprsLO[x_,Qs_]:=Dform[bprsLO,aprsLO,x,sLO[Qs,5]];
(* charm *)
bprcLO={0.11061, -1.5434, 2.20681};
aprcLO={-0.07726, 0.05422, -0.03364, -0.20804, 0.29038, -0.23662, 0.62274, 0.29713, -0.21861, 0.00831, 0.00065};
dprcLO[x_,Qs_]:=Dform[bprcLO,aprcLO,x,sLO[Qs,7]];
(* bottom *)
bprbLO={40.0971, 0.74249, 12.3729};
aprbLO={-123.531, 128.666, -29.1808, -1.29639, -3.65003, 3.0534, -1.04932, 0.34662, -1.34412, -0.0429, -0.30359};
dprbLO[x_,Qs_]:=Dform[bprbLO,aprbLO,x,sLO[Qs,9]];
(*********************)
(*********NLO*********)
(*********************)
(* NLO PION *)
(* gluon *)
bpgNLO={3.73331, -0.74159, 2.33092};
apgNLO={-3.16946, -0.47683, 0.7027, -0.51377, -0.19705, -0.17917, 2.03394, -0.50764, -0.08565, 0.09466, -0.10222};
dpgNLO[x_,Qs_]:=Dform[bpgNLO,apgNLO,x,sNLO[Qs,0]];
(* up *)
bpuNLO={0.44809, -1.47598, 0.91338};
apuNLO={-0.13828, -0.06951, 0.01354, -0.30498, -0.01863, -0.12529, 0.64145, 0.0727, -0.16989, 0.07396, -0.07757};
dpuNLO[x_,Qs_]:=Dform[bpuNLO,apuNLO,x,sNLO[Qs,1]];
(* strange *)
bpsNLO={16.5987, 0.13345, 5.89903};
apsNLO={-18.3856, 2.44225, 2.13225, 0.22712, -0.83625, 0.38526, -0.16911, 0.59886, -0.2563, -0.18619, 0.87362};
dpsNLO[x_,Qs_]:=Dform[bpsNLO,apsNLO,x,sNLO[Qs,5]];
(* charm *)
bpcNLO={6.17173, -0.53618, 5.60108};
apcNLO={-4.8245, -1.30844, 1.95527, -0.27879, -0.51337, 0.109, 0.83571, -1.15141, 0.77027, 0.09268, -0.11267};
dpcNLO[x_,Qs_]:=Dform[bpcNLO,apcNLO,x,sNLO[Qs,7]];
(* bottom *)
bpbNLO={0.25944, -1.98713, 3.52857};
apbNLO={-0.11449, 0.03733, -0.18028, -0.35858, 0.22277, -0.66413, 0.72303, 0.4626, -0.99235, -0.02701, -0.02089};
dpbNLO[x_,Qs_]:=Dform[bpbNLO,apbNLO,x,sNLO[Qs,9]];
(* NLO KAON *)
(* gluon *)
bkgNLO={0.2314, -1.364, 1.79761};
akgNLO={-0.33644, 0.16204, -0.02598, 0.97182, -0.02908, -0.43195, 1.57116, 0.71847, -0.68331, 0.36906, 2.3906};
dkgNLO[x_,Qs_]:=Dform[bkgNLO,akgNLO,x,sNLO[Qs,0]];
(* up *)
bkuNLO={0.17806, -0.53733, 0.7594};
akuNLO={-0.10988, -0.02524, 0.03142, -0.60058, 0.07863, 0.13276, 0.61356, -0.43886, 0.23942, 0.10742, 0.128};
dkuNLO[x_,Qs_]:=Dform[bkuNLO,akuNLO,x,sNLO[Qs,1]];
(* down *)
bkdNLO={4.96269, 0.05562, 2.79926};
akdNLO={1.54098, -9.06376, 4.94791, 1.8866, -2.9435, 1.04227, 3.02991, -4.14807, 1.91494, 0.8545, -0.61016};
dkdNLO[x_,Qs_]:=Dform[bkdNLO,akdNLO,x,sNLO[Qs,3]];
(* charm *)
bkcNLO={4.25954, -0.24144, 4.21265};
akcNLO={-5.44309, 6.11031, -3.13973, -1.07757, 1.52364, -0.74308, 0.2559, 0.98423, -0.52839, -0.04, 0.08695};
dkcNLO[x_,Qs_]:=Dform[bkcNLO,akcNLO,x,sNLO[Qs,7]];
(* bottom *)
bkbNLO={1.32443, -0.88351, 6.15221};
akbNLO={-1.41156, -0.04809, 0.79066, -0.44818, -0.60073, 0.45526, 0.46679, -0.50792, 0.67006, -0.00477, -0.05503};
dkbNLO[x_,Qs_]:=Dform[bkbNLO,akbNLO,x,sNLO[Qs,9]];
(* NLO PROTON *)
(* gluon *)
bprgNLO={1.56255, 0.01567, 3.57583};
aprgNLO={-1.48158, -0.39439, 0.51249, -2.16232, 2.47127, -0.93259, 3.33958, -3.05265, 1.21042, -0.84816, 1.23583};
dprgNLO[x_,Qs_]:=Dform[bprgNLO,aprgNLO,x,sNLO[Qs,0]];
(* up *)
bpruNLO={1.25946, 0.07124, 4.12795};
apruNLO={-1.17505, 0.3755, -0.01416, -0.29533, -0.2454, 0.16543, 0.98867, -0.46846, 0.2075, 0.18957, -0.01116};
dpruNLO[x_,Qs_]:=Dform[bpruNLO,apruNLO,x,sNLO[Qs,1]];
(* strange *)
bprsNLO={4.01135, 0.17258, 5.20766};
aprsNLO={8.67124, -22.7888, 11.472, 4.57608, -9.64835, 4.61792, 7.25144, -12.6313, 6.07314, 0.16931, -0.09541};
dprsNLO[x_,Qs_]:=Dform[bprsNLO,aprsNLO,x,sNLO[Qs,5]];
(* charm *)
bprcNLO={0.0825, -1.6129, 2.01255};
aprcNLO={-0.04512, -0.00565, 0.009, -0.38012, -0.0684, 0.08888, 0.63782, -0.14146, 0.06083, -0.02958, 0.0113};
dprcNLO[x_,Qs_]:=Dform[bprcNLO,aprcNLO,x,sNLO[Qs,7]];
(* bottom *)
bprbNLO={24.2916, 0.57939, 12.1207};
aprbNLO={-88.3524, 93.1056, -17.4089, -0.80783, -5.072, -2.45377, -3.2737, 1.21188, -5.50374, 0.14628, -0.78634};
dprbNLO[x_,Qs_]:=Dform[bprbNLO,aprbNLO,x,sNLO[Qs,9]];

ihtestlist=Range[7];
startdh=0;enddh=10;
dhtestlist=Range[0,enddh];
isettestlist={0,1};

KKP[ih_Integer,iset_Integer,x_?NumericQ,Qs_?NumericQ]:=
If[
(* check for incorrect inputs *)
Not[MemberQ[ihtestlist,ih]]||Not[MemberQ[isettestlist,iset]]||Re[N[x]]≠N[x]||x≤0||Re[N[Qs]]≠N[Qs]||Qs<0,
(* then *)
Message[KKP::input2];
False,
(* else *)
If[x<.1||x>.8,Message[KKP::xarg,x]];
For[i=startdh,i≤enddh,i++,If[Qs<Q0[i]||Qs>100,Message[KKP::Qsarg,Qs,Q0[i],i]]];
If[
iset==0,
(* then LO *)
dpg[tx_,tQs_]:=dpgLO[tx,tQs];
dpu[tx_,tQs_]:=dpuLO[tx,tQs];
dps[tx_,tQs_]:=dpsLO[tx,tQs];
dpc[tx_,tQs_]:=dpcLO[tx,tQs];
dpb[tx_,tQs_]:=dpbLO[tx,tQs];
dkg[tx_,tQs_]:=dkgLO[tx,tQs];
dku[tx_,tQs_]:=dkuLO[tx,tQs];
dkd[tx_,tQs_]:=dkdLO[tx,tQs];
dkc[tx_,tQs_]:=dkcLO[tx,tQs];
dkb[tx_,tQs_]:=dkbLO[tx,tQs];
dprg[tx_,tQs_]:=dprgLO[tx,tQs];
dpru[tx_,tQs_]:=dpruLO[tx,tQs];
dprs[tx_,tQs_]:=dprsLO[tx,tQs];
dprc[tx_,tQs_]:=dprcLO[tx,tQs];
dprb[tx_,tQs_]:=dprbLO[tx,tQs],
(* then NLO *)
dpg[tx_,tQs_]:=dpgNLO[tx,tQs];
dpu[tx_,tQs_]:=dpuNLO[tx,tQs];
dps[tx_,tQs_]:=dpsNLO[tx,tQs];
dpc[tx_,tQs_]:=dpcNLO[tx,tQs];
dpb[tx_,tQs_]:=dpbNLO[tx,tQs];
dkg[tx_,tQs_]:=dkgNLO[tx,tQs];
dku[tx_,tQs_]:=dkuNLO[tx,tQs];
dkd[tx_,tQs_]:=dkdNLO[tx,tQs];
dkc[tx_,tQs_]:=dkcNLO[tx,tQs];
dkb[tx_,tQs_]:=dkbNLO[tx,tQs];
dprg[tx_,tQs_]:=dprgNLO[tx,tQs];
dpru[tx_,tQs_]:=dpruNLO[tx,tQs];
dprs[tx_,tQs_]:=dprsNLO[tx,tQs];
dprc[tx_,tQs_]:=dprcNLO[tx,tQs];
dprb[tx_,tQs_]:=dprbNLO[tx,tQs]
];
dpd[tx_,tQs_]:=dpu[tx,tQs];
dprd[tx_,tQs_]:=.5 dpru[tx,tQs];
dks[tx_,tQs_]:=dku[tx,tQs];
If[Qs<rmbb,dpb[tx_,tQs_]:=0.;dkb[tx_,tQs_]:=0.;dprb[tx_,tQs_]:=0.];
If[Qs<rmcc,dpc[tx_,tQs_]:=0.;dkc[tx_,tQs_]:=0.;dprc[tx_,tQs_]:=0.];
Which[
ih==1||ih==5,
Return[{
dpg[x,Qs]/2.,
dpu[x,Qs]/2.,
dpu[x,Qs]/2.,
dpd[x,Qs]/2.,
dpd[x,Qs]/2.,
dps[x,Qs]/2.,
dps[x,Qs]/2.,
dpc[x,Qs]/2.,
dpc[x,Qs]/2.,
dpb[x,Qs]/2.,
dpb[x,Qs]/2.
}],
ih==2||ih==3,
Return[{
dkg[x,Qs]/2.,
dku[x,Qs]/2.,
dku[x,Qs]/2.,
dkd[x,Qs]/2.,
dkd[x,Qs]/2.,
dks[x,Qs]/2.,
dks[x,Qs]/2.,
dkc[x,Qs]/2.,
dkc[x,Qs]/2.,
dkb[x,Qs]/2.,
dkb[x,Qs]/2.
}],
ih==4,
Return[{
dprg[x,Qs]/2.,
dpru[x,Qs]/2.,
dpru[x,Qs]/2.,
dprd[x,Qs]/2.,
dprd[x,Qs]/2.,
dprs[x,Qs]/2.,
dprs[x,Qs]/2.,
dprc[x,Qs]/2.,
dprc[x,Qs]/2.,
dprb[x,Qs]/2.,
dprb[x,Qs]/2.
}],
ih==6,
Return[{
dprg[x,Qs]/2.,
dpru[x,Qs]/4.,
dpru[x,Qs]/4.,
dprd[x,Qs],
dprd[x,Qs],
dprs[x,Qs]/2.,
dprs[x,Qs]/2.,
dprc[x,Qs]/2.,
dprc[x,Qs]/2.,
dprb[x,Qs]/2.,
dprb[x,Qs]/2.
}],
ih==7,
Return[{
dpg[x,Qs]+dkg[x,Qs]+dprg[x,Qs],
dpu[x,Qs]+dku[x,Qs]+dpru[x,Qs],
dpu[x,Qs]+dku[x,Qs]+dpru[x,Qs],
dpd[x,Qs]+dkd[x,Qs]+dprd[x,Qs],
dpd[x,Qs]+dkd[x,Qs]+dprd[x,Qs],
dps[x,Qs]+dks[x,Qs]+dprs[x,Qs],
dps[x,Qs]+dks[x,Qs]+dprs[x,Qs],
dpc[x,Qs]+dkc[x,Qs]+dprc[x,Qs],
dpc[x,Qs]+dkc[x,Qs]+dprc[x,Qs],
dpb[x,Qs]+dkb[x,Qs]+dprb[x,Qs],
dpb[x,Qs]+dkb[x,Qs]+dprb[x,Qs]
}]
];
];

KKP[___]:=(Message[KKP::input1];False);
End[];
EndPackage[];

testing with kkpREVISED.f

correct what looks to be an error in the KKP.f code

results from testkkp4.f

ih=Range[7];
ilist={0,1};
xlist={.1,.2,.3,.4,.5,.6,.7,.8};
qslist={2.,5.,10.,50.,100.};
testlist=Tuples[{ih,ilist,xlist,qslist}];
Length[testlist]

xlistF={".1d0",".2d0",".3d0",".4d0",".5d0",".6d0",".7d0",".8d0"};
qslistF={"2.d0","5.d0","1.d1","5.d1","1.d2"};
testlistF=Tuples[{ih,ilist,xlistF,qslistF}];

SetDirectory["G:\cygwin\home\George P. Burdell\gyulassy\Fragmentation Functions\KKP\KKPrevised"]

Export["testlistF.dat",testlistF]

kkpoutput=KKP@@@testlist;

kkpoutputFORTRAN=Import["testresultsREVISEDF.dat"];
kkpoutputF3=Join[kkpoutputFORTRAN[[3 (#-1)+1]],kkpoutputFORTRAN[[3 (#-1)+2]],kkpoutputFORTRAN[[3 (#-1)+3]]]&/@Range[560];

kkpoutputdiff=kkpoutput-kkpoutputF3;

Max[Abs[kkpoutputdiff]]

percentdiff=Table[
If[kkpoutputF3[[iter,#]]≠0.,
Abs[kkpoutputdiff[[iter,#]]]/kkpoutputF3[[iter,#]],
If[kkpoutput[[iter,#]]≠0.,1.,0.]
]&/@Range[11],{iter,1,560}];

Max[percentdiff]*100

Maximum difference between the two outputs is .002%.