(* ::Package:: *)

BeginPackage["tmd1p5`"]

PDFCollinear::usage=
"The function PDFCollinear[parton,xA,bT] gives the collinear part of the TMD PDF for a quark of type parton with bmax=1.5 GeV^-1. parton can be up, upbar, down, downbar. The valid range for xA and bT are 0.001 <= xA <= 0.69 and 0.001 <= bT <= 14.99."

FFCollinear::usage=
"The function FFCollinear[parton,zA,bT] gives the collinear part of the TMD FF for a quark of type parton with bmax=1.5 GeV^-1. parton can be up,upbar,down,downbar. The valid range for zA and bT are 0.2 <= zA <= 0.99 and 0.001 <= bT <= 14.99."
    
    
    

dummytable={};
PDFxGrid={};
PDFbGrid={};
PDFup=Table[.0,{113},{1189}];
PDFdown=Table[.0,{113},{1189}];
PDFupbar=Table[.0,{113},{1189}];
PDFdownbar=Table[.0,{113},{1189}];

FFxGrid={};
FFbGrid={};
FFup=Table[.0,{140},{1149}];
FFdown=Table[.0,{140},{1149}];
FFupbar=Table[.0,{140},{1149}];
FFdownbar=Table[.0,{140},{1149}];


PDFxincrement=0.0001;
PDFbincrement=0.0001;
PDFnx=113;
PDFnb=1189;
FFxincrement=0.005;
FFbincrement=0.0001;
FFnx=140;
FFnb=1149;





For[xx=0.001,xx<0.7,xx=xx+PDFxincrement,
PDFxGrid=AppendTo[PDFxGrid,xx];
If[xx>0.003,PDFxincrement=0.0003];
If[xx>0.0098,PDFxincrement=0.01];
];


For[btt=0.001,btt<=15,btt=btt+PDFbincrement,
PDFbGrid=AppendTo[PDFbGrid,btt];
If[btt>0.09,PDFbincrement=0.05];
];

For[xx=0.2,xx<1.0,xx=xx+FFxincrement,
FFxGrid=AppendTo[FFxGrid,xx];
If[xx>=0.8,FFxincrement=0.01];
];


For[btt=0.005,btt<=15,btt=btt+FFbincrement,
FFbGrid=AppendTo[FFbGrid,btt];
If[btt>0.09,FFbincrement=0.05];
];



PDFGridcount=1;
dummytable=AppendTo[dummytable,Import["./PDFGrids/PDFup1p5.m"]];

For[i=1,i<=PDFnx,i=i+1,
For[j=1,j<=PDFnb,j=j+1,
PDFup[[i,j]]=dummytable[[1,PDFGridcount]];
PDFGridcount++;
];
];

dummytable={};

PDFGridcount=1;
dummytable=AppendTo[dummytable,Import["./PDFGrids/PDFdown1p5.m"]];

For[i=1,i<=PDFnx,i=i+1,
For[j=1,j<=PDFnb,j=j+1,
PDFdown[[i,j]]=dummytable[[1,PDFGridcount]];
PDFGridcount++;
];
];

dummytable={};



PDFGridcount=1;
dummytable=AppendTo[dummytable,Import["./PDFGrids/PDFupbar1p5.m"]];

For[i=1,i<=PDFnx,i=i+1,
For[j=1,j<=PDFnb,j=j+1,
PDFupbar[[i,j]]=dummytable[[1,PDFGridcount]];
PDFGridcount++;
];
];

dummytable={};

PDFGridcount=1;
dummytable=AppendTo[dummytable,Import["./PDFGrids/PDFdownbar1p5.m"]];

For[i=1,i<=PDFnx,i=i+1,
For[j=1,j<=PDFnb,j=j+1,
PDFdownbar[[i,j]]=dummytable[[1,PDFGridcount]];
PDFGridcount++;
];
];

dummytable={};

FFGridcount=1;
dummytable=AppendTo[dummytable,Import["./FFGrids/FFup1p5.m"]];

For[i=1,i<=FFnx,i=i+1,
For[j=1,j<=FFnb,j=j+1,
FFup[[i,j]]=dummytable[[1,FFGridcount]];
FFGridcount++;
];
];

dummytable={};

FFGridcount=1;
dummytable=AppendTo[dummytable,Import["./FFGrids/FFdown1p5.m"]];

For[i=1,i<=FFnx,i=i+1,
For[j=1,j<=FFnb,j=j+1,
FFdown[[i,j]]=dummytable[[1,FFGridcount]];
FFGridcount++;
];
];

dummytable={};

FFGridcount=1;
dummytable=AppendTo[dummytable,Import["./FFGrids/FFup1p5.m"]];

For[i=1,i<=FFnx,i=i+1,
For[j=1,j<=FFnb,j=j+1,
FFupbar[[i,j]]=dummytable[[1,FFGridcount]];
FFGridcount++;
];
];

dummytable={};

FFGridcount=1;
dummytable=AppendTo[dummytable,Import["./FFGrids/FFdown1p5.m"]];

For[i=1,i<=FFnx,i=i+1,
For[j=1,j<=FFnb,j=j+1,
FFdownbar[[i,j]]=dummytable[[1,FFGridcount]];
FFGridcount++;
];
];





Clear[dummytable,xx,btt,PDFGridcount,i,j,PDFxincrement,PDFbincrement,PDFnx,PDFnb,FFGridcount,FFxincrement,FFbincrement,FFnx,FFnb];


Begin["`Private`"]


makefilenamePDF[string_]:=ToExpression["PDF"<>ToString[string]];
makefilenameFF[string_]:=ToExpression["FF"<>ToString[string]];



locate[xx_,n_,x_]:=Module[{ju,jm,jl=0,j},
ju=n+1;
While[ju-jl>1,
jm=IntegerPart[(ju+jl)/2];
If[x>=xx[[jm]],jl=jm,ju=jm]
];
If[x==xx[[1]],j=1];
If[x==xx[[n]],j=n-1];
If[x!= xx[[1]] && x!=xx[[n]],j=jl];
j
]
polderivative1[x1_,x2_,x3_,y1_,y2_,y3_]:=(x3*x3*(y1-y2)+2.*x1*(x3*(-y1+y2)+x2*(y1-y3))+x2*x2*(-y1+y3)+x1*x1*(-y2+y3))/((x1-x2)*(x1-x3)*(x2-x3));
polderivative2[x1_,x2_,x3_,y1_,y2_,y3_]:=(x3*x3*(y1-y2)-2.*x2*(x3*(y1-y2)+x1*(y2-y3))+x2*x2*(y1-y3)+x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3));
polderivative3[x1_,x2_,x3_,y1_,y2_,y3_]:=(x3*x3*(-y1+y2)+2.*x2*x3*(y1-y3)+x1*x1*(y2-y3)+x2*x2*(-y1+y3)+2.*x1*x3*(-y2+y3))/((x1-x2)*(x1-x3)*(x2-x3));
Ccoeffs[y_,y1_,y2_,y12_,d1_,d2_]:=Module[
{l,k,j,i,xx,d1d2,cl=Table[0.,{16}],x=Table[0.,{16}],wt=Table[0,{16},{16}],c=Table[0.,{4},{4}]},

wt={{1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0},{-3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0},{2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0},{0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0},{0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1},{0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1},{-3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0},{9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2},{-6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2},{2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0},{-6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1},{4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1}};

d1d2=d1*d2;

For[i=1,i<=4,i++,
x[[i]]=y[[i]];
x[[i+4]]=y1[[i]]*d1;
x[[i+8]]=y2[[i]]*d2;
x[[i+12]]=y12[[i]]*d1d2;
];

For[i=1,i<=16,i++,
xx=0.0;
For[k=1,k<=16,k++,
xx=xx+wt[[i,k]]*x[[k]];
];
cl[[i]]=xx;
];



l=1;

For[i=1,i<=4,i++,
For[j=1,j<=4,j++,
c[[i,j]]=cl[[l]];
l++;
];
];
c
]
bicubint[y_,y1_,y2_,y12_,x1l_,x1u_,x2l_,x2u_,x1_,x2_]:=Module[
{i,t,u,d1,d2,c=Table[0.,{4},{4}],ansy=0.},

d1=x1u-x1l;
d2=x2u-x2l;


c=Ccoeffs[y,y1,y2,y12,d1,d2];

t=(x1-x1l)/d1;
u=(x2-x2l)/d2;

For[i=4,i>=1,i--,
ansy=t*ansy+((c[[i,4]]*u+c[[i,3]])*u+c[[i,2]])*u+c[[i,1]];
];

ansy
]
PDFCollinear[parton_,xA_,bT_]:=Module[
{i,j,k,l,xx,btt,PDFnx=113,PDFnb=1189,posx,posb,F=Table[0.0,{113},{1189}],f=Table[0.,{4}],f1=Table[0.,{4}],f2=Table[0.,{4}],f12=Table[0.,{4}],f21=Table[0.,{4}]},


If[ToString[parton]!="up" && ToString[parton]!="down"&&ToString[parton]!="upbar" && ToString[parton]!="downbar",
Print["You haven't entered a valid parton"];
Print["Choices are: up, upbar, down, downbar"];
Abort[];
];

If[xA<0.001 || xA>0.69,
Print["You haven't entered a valid range for xA. You must have 0.001<=xA<=0.69."];
Abort[];
];

If[bT<0.001 || bT>14.99,
Print["You haven't entered a valid range for bT. You must have 0.001<=bT<=14.99."];
Abort[];
];

F=makefilenamePDF[parton];




posx=locate[PDFxGrid,PDFnx,xA];
posb=locate[PDFbGrid,PDFnb,bT];

f[[1]]=F[[posx,posb]];
f[[2]]=F[[posx+1,posb]];
f[[3]]=F[[posx+1,posb+1]];
f[[4]]=F[[posx,posb+1]];



If[posx==1 && posb==1,

f1[[1]]=polderivative1 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb]],F[[posx+1,posb]],F[[posx+2,posb]]];f1[[2]]=polderivative2 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb]],F[[posx+1,posb]],F[[posx+2,posb]]];
f1[[3]]=polderivative2 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb+1]],F[[posx+1,posb+1]],F[[posx+2,posb+1]]];
f1[[4]]=polderivative1 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb+1]],F[[posx+1,posb+1]],F[[posx+2,posb+1]]];

f2[[1]]=polderivative1 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx,posb]],F[[posx,posb+1]],F[[posx,posb+2]]];
f2[[2]]=polderivative1 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx+1,posb]],F[[posx+1,posb+1]],F[[posx+1,posb+2]]];f2[[3]]=polderivative2 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx+1,posb]],F[[posx+1,posb+1]],F[[posx+1,posb+2]]];f2[[4]]=polderivative2 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx,posb]],F[[posx,posb+1]],F[[posx,posb+2]]];

f12[[1]]=polderivative1 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],f2[[1]],f2[[2]],polderivative1[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx+2,posb]],F[[posx+2,posb+1]],F[[posx+2,posb+2]]]];
f12[[2]]=polderivative2 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],f2[[1]],f2[[2]],polderivative1[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx+2,posb]],F[[posx+2,posb+1]],F[[posx+2,posb+2]]]];f12[[3]]=polderivative2 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],f2[[4]],f2[[3]],polderivative2[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx+2,posb]],F[[posx+2,posb+1]],F[[posx+2,posb+2]]]];f12[[4]]=polderivative1 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],f2[[4]],f2[[3]],polderivative2[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx+2,posb]],F[[posx+2,posb+1]],F[[posx+2,posb+2]]]];

f21[[1]]=polderivative1 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],f1[[1]],f1[[4]],polderivative1[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb+2]],F[[posx+1,posb+2]],F[[posx+2,posb+2]]]];f21[[2]]=polderivative1 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],f1[[2]],f1[[3]],polderivative2[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb+2]],F[[posx+1,posb+2]],F[[posx+2,posb+2]]]];
f21[[3]]=polderivative2 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],f1[[2]],f1[[3]],polderivative2[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb+2]],F[[posx+1,posb+2]],F[[posx+2,posb+2]]]];f21[[4]]=polderivative2 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],f1[[1]],f1[[4]],polderivative1[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb+2]],F[[posx+1,posb+2]],F[[posx+2,posb+2]]]];



f12=0.5*(f12+f21);


];


If[posx==PDFnx-1 && posb==PDFnb-1,

f1[[1]]=polderivative2 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb]],F[[posx,posb]],F[[posx+1,posb]]];f1[[2]]=polderivative3 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb]],F[[posx,posb]],F[[posx+1,posb]]];f1[[3]]=polderivative3 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb+1]],F[[posx,posb+1]],F[[posx+1,posb+1]]];f1[[4]]=polderivative2 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb+1]],F[[posx,posb+1]],F[[posx+1,posb+1]]];

f2[[1]]=polderivative2 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx,posb-1]],F[[posx,posb]],F[[posx,posb+1]]];f2[[2]]=polderivative2 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx+1,posb-1]],F[[posx+1,posb]],F[[posx+1,posb+1]]];f2[[3]]=polderivative3 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx+1,posb-1]],F[[posx+1,posb]],F[[posx+1,posb+1]]];f2[[4]]=polderivative3 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx,posb-1]],F[[posx,posb]],F[[posx,posb+1]]];

f12[[1]]=polderivative2 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],polderivative2[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx-1,posb-1]],F[[posx-1,posb]],F[[posx-1,posb+1]]],f2[[1]],f2[[2]]];f12[[2]]=polderivative3 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],polderivative2[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx-1,posb-1]],F[[posx-1,posb]],F[[posx-1,posb+1]]],f2[[1]],f2[[2]]];f12[[3]]=polderivative3 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],polderivative3[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx-1,posb-1]],F[[posx-1,posb]],F[[posx-1,posb+1]]],f2[[4]],f2[[3]]];f12[[4]]=polderivative2 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],polderivative3[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx-1,posb-1]],F[[posx-1,posb]],F[[posx-1,posb+1]]],f2[[4]],f2[[3]]];

f21[[1]]=polderivative2 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],polderivative2[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb-1]],F[[posx,posb-1]],F[[posx+1,posb-1]]],f1[[1]],f1[[4]]];f21[[2]]=polderivative2 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],polderivative3[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb-1]],F[[posx,posb-1]],F[[posx+1,posb-1]]],f1[[2]],f1[[3]]];f21[[3]]=polderivative3 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],polderivative3[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb-1]],F[[posx,posb-1]],F[[posx+1,posb-1]]],f1[[2]],f1[[3]]];f21[[4]]=polderivative3 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],polderivative2[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb-1]],F[[posx,posb-1]],F[[posx+1,posb-1]]],f1[[1]],f1[[4]]];

f12 = 0.5*(f12+f21);
];


If[posx==1 && posb==PDFnb-1,

f1[[1]]=polderivative1 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb]],F[[posx+1,posb]],F[[posx+2,posb]]];f1[[2]]=polderivative2 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb]],F[[posx+1,posb]],F[[posx+2,posb]]];f1[[3]]=polderivative2 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb+1]],F[[posx+1,posb+1]],F[[posx+2,posb+1]]];f1[[4]]=polderivative1 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb+1]],F[[posx+1,posb+1]],F[[posx+2,posb+1]]];

f2[[1]]=polderivative2 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx,posb-1]],F[[posx,posb]],F[[posx,posb+1]]];f2[[2]]=polderivative2 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx+1,posb-1]],F[[posx+1,posb]],F[[posx+1,posb+1]]];f2[[3]]=polderivative3 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx+1,posb-1]],F[[posx+1,posb]],F[[posx+1,posb+1]]];f2[[4]]=polderivative3 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx,posb-1]],F[[posx,posb]],F[[posx,posb+1]]];

f12[[1]]=polderivative1 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],f2[[1]],f2[[2]],polderivative2[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx+2,posb-1]],F[[posx+2,posb]],F[[posx+2,posb+1]]]];f12[[2]]=polderivative2 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],f2[[1]],f2[[2]],polderivative2[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx+2,posb-1]],F[[posx+2,posb]],F[[posx+2,posb+1]]]];f12[[3]]=polderivative2 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],f2[[4]],f2[[3]],polderivative3[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx+2,posb-1]],F[[posx+2,posb]],F[[posx+2,posb+1]]]];f12[[4]]=polderivative1 [PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],f2[[4]],f2[[3]],polderivative3[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx+2,posb-1]],F[[posx+2,posb]],F[[posx+2,posb+1]]]];

f21[[1]]=polderivative2 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],polderivative1[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb-1]],F[[posx+1,posb-1]],F[[posx+2,posb-1]]],f1[[1]],f1[[4]]];f21[[2]]=polderivative2 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],polderivative2[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb-1]],F[[posx+1,posb-1]],F[[posx+2,posb-1]]],f1[[2]],f1[[3]]];f21[[3]]=polderivative3 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],polderivative2[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb-1]],F[[posx+1,posb-1]],F[[posx+2,posb-1]]],f1[[2]],f1[[3]]];f21[[4]]=polderivative3 [PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],polderivative1[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb-1]],F[[posx+1,posb-1]],F[[posx+2,posb-1]]],f1[[1]],f1[[4]]];

f12=0.5*(f12+f21);


];


If[posx==PDFnx-1 && posb==1,

f1[[1]]=polderivative2 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb]],F[[posx,posb]],F[[posx+1,posb]]];f1[[2]]=polderivative3 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb]],F[[posx,posb]],F[[posx+1,posb]]];f1[[3]]=polderivative3 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb+1]],F[[posx,posb+1]],F[[posx+1,posb+1]]];f1[[4]]=polderivative2 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb+1]],F[[posx,posb+1]],F[[posx+1,posb+1]]];

f2[[1]]=polderivative1 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx,posb]],F[[posx,posb+1]],F[[posx,posb+2]]];f2[[2]]=polderivative1 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx+1,posb]],F[[posx+1,posb+1]],F[[posx+1,posb+2]]];f2[[3]]=polderivative2 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx+1,posb]],F[[posx+1,posb+1]],F[[posx+1,posb+2]]];f2[[4]]=polderivative2 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx,posb]],F[[posx,posb+1]],F[[posx,posb+2]]];

f12[[1]]=polderivative2 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],polderivative1[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx-1,posb]],F[[posx-1,posb+1]],F[[posx-1,posb+2]]],f2[[1]],f2[[2]]];f12[[2]]=polderivative3 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],polderivative1[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx-1,posb]],F[[posx-1,posb+1]],F[[posx-1,posb+2]]],f2[[1]],f2[[2]]];f12[[3]]=polderivative3 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],polderivative2[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx-1,posb]],F[[posx-1,posb+1]],F[[posx-1,posb+2]]],f2[[4]],f2[[3]]];f12[[4]]=polderivative2 [PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],polderivative2[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx-1,posb]],F[[posx-1,posb+1]],F[[posx-1,posb+2]]],f2[[4]],f2[[3]]];

f21[[1]]=polderivative1 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],f1[[1]],f1[[4]],polderivative2[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb+2]],F[[posx,posb+2]],F[[posx+1,posb+2]]]];f21[[2]]=polderivative1 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],f1[[2]],f1[[3]],polderivative3[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb+2]],F[[posx,posb+2]],F[[posx+1,posb+2]]]];f21[[3]]=polderivative2 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],f1[[2]],f1[[3]],polderivative3[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb+2]],F[[posx,posb+2]],F[[posx+1,posb+2]]]];f21[[4]]=polderivative2 [PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],f1[[1]],f1[[4]],polderivative2[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb+2]],F[[posx,posb+2]],F[[posx+1,posb+2]]]];

f12 = 0.5*(f12+f21);


];

If[posx!= 1 && posx!= PDFnx-1 && posb != 1 && posb != PDFnb-1,

f1[[1]]=polderivative2[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb]],F[[posx,posb]],F[[posx+1,posb]]];f1[[2]]=polderivative2[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb]],F[[posx+1,posb]],F[[posx+2,posb]]];f1[[3]]=polderivative2[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb+1]],F[[posx+1,posb+1]],F[[posx+2,posb+1]]];f1[[4]]=polderivative2[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb+1]],F[[posx,posb+1]],F[[posx+1,posb]]];

f2[[1]]=polderivative2[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx,posb-1]],F[[posx,posb]],F[[posx,posb+1]]];f2[[2]]=polderivative2[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx+1,posb-1]],F[[posx+1,posb]],F[[posx+1,posb+1]]];f2[[3]]=polderivative2[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx+1,posb]],F[[posx+1,posb+1]],F[[posx+1,posb+2]]];f2[[4]]=polderivative2[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx,posb]],F[[posx,posb+1]],F[[posx,posb+2]]];

f12[[1]]=polderivative2[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],polderivative2[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx-1,posb-1]],F[[posx-1,posb]],F[[posx-1,posb+1]]],f2[[1]],f2[[2]]];f12[[2]]=polderivative2[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],f2[[1]],f2[[2]],polderivative2[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],F[[posx+2,posb-1]],F[[posx+2,posb]],F[[posx+2,posb+1]]]];f12[[3]]=polderivative2[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],f2[[4]],f2[[3]],polderivative2[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx+2,posb]],F[[posx+2,posb+1]],F[[posx+2,posb+2]]]];f12[[4]]=polderivative2[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],polderivative2[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],F[[posx-1,posb]],F[[posx-1,posb+1]],F[[posx-1,posb+2]]],f2[[4]],f2[[3]]];

f21[[1]]=polderivative2[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],polderivative2[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb-1]],F[[posx,posb-1]],F[[posx+1,posb-1]]],f1[[1]],f1[[4]]];f21[[2]]=polderivative2[PDFbGrid[[posb-1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],polderivative2[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb-1]],F[[posx+1,posb-1]],F[[posx+2,posb-1]]],f1[[2]],f1[[3]]];f21[[3]]=polderivative2[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],f1[[2]],f1[[3]],polderivative2[PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFxGrid[[posx+2]],F[[posx,posb+2]],F[[posx+1,posb+2]],F[[posx+2,posb+2]]]];f21[[4]]=polderivative2[PDFbGrid[[posb]],PDFbGrid[[posb+1]],PDFbGrid[[posb+2]],f1[[1]],f1[[4]],polderivative2[PDFxGrid[[posx-1]],PDFxGrid[[posx]],PDFxGrid[[posx+1]],F[[posx-1,posb+2]],F[[posx,posb+2]],F[[posx+1,posb+2]]]];

f12= 0.5*(f12+f21);

];


bicubint[f,f1,f2,f12,PDFxGrid[[posx]],PDFxGrid[[posx+1]],PDFbGrid[[posb]],PDFbGrid[[posb+1]],xA,bT]
]


FFCollinear[parton_,xA_,bT_]:=Module[
{i,j,k,l,xx,btt,FFnx=140,FFnb=1149,posx,posb,F=Table[0.0,{140},{1149}],f=Table[0.,{4}],f1=Table[0.,{4}],f2=Table[0.,{4}],f12=Table[0.,{4}],f21=Table[0.,{4}]},


If[ToString[parton]!="up" && ToString[parton]!="down" && ToString[parton]!="upbar" && ToString[parton]!="downbar",
Print["You haven't entered a valid parton"];
Print["Choices are: up,down,upbar,downbar"];
Abort[];
];

If[xA<0.2 || xA>0.99,
Print["You haven't entered a valid range for zA. You must have 0.2<=zA<=0.99."];
Abort[];
];

If[bT<0.005 || bT>14.99,
Print["You haven't entered a valid range for bT. You must have 0.005<=bT<=14.99."];
Abort[];
];

F=makefilenameFF[parton];




posx=locate[FFxGrid,FFnx,xA];
posb=locate[FFbGrid,FFnb,bT];

f[[1]]=F[[posx,posb]];
f[[2]]=F[[posx+1,posb]];
f[[3]]=F[[posx+1,posb+1]];
f[[4]]=F[[posx,posb+1]];



If[posx==1 && posb==1,

f1[[1]]=polderivative1 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb]],F[[posx+1,posb]],F[[posx+2,posb]]];f1[[2]]=polderivative2 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb]],F[[posx+1,posb]],F[[posx+2,posb]]];
f1[[3]]=polderivative2 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb+1]],F[[posx+1,posb+1]],F[[posx+2,posb+1]]];
f1[[4]]=polderivative1 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb+1]],F[[posx+1,posb+1]],F[[posx+2,posb+1]]];

f2[[1]]=polderivative1 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx,posb]],F[[posx,posb+1]],F[[posx,posb+2]]];
f2[[2]]=polderivative1 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx+1,posb]],F[[posx+1,posb+1]],F[[posx+1,posb+2]]];f2[[3]]=polderivative2 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx+1,posb]],F[[posx+1,posb+1]],F[[posx+1,posb+2]]];f2[[4]]=polderivative2 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx,posb]],F[[posx,posb+1]],F[[posx,posb+2]]];

f12[[1]]=polderivative1 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],f2[[1]],f2[[2]],polderivative1[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx+2,posb]],F[[posx+2,posb+1]],F[[posx+2,posb+2]]]];
f12[[2]]=polderivative2 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],f2[[1]],f2[[2]],polderivative1[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx+2,posb]],F[[posx+2,posb+1]],F[[posx+2,posb+2]]]];f12[[3]]=polderivative2 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],f2[[4]],f2[[3]],polderivative2[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx+2,posb]],F[[posx+2,posb+1]],F[[posx+2,posb+2]]]];f12[[4]]=polderivative1 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],f2[[4]],f2[[3]],polderivative2[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx+2,posb]],F[[posx+2,posb+1]],F[[posx+2,posb+2]]]];

f21[[1]]=polderivative1 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],f1[[1]],f1[[4]],polderivative1[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb+2]],F[[posx+1,posb+2]],F[[posx+2,posb+2]]]];f21[[2]]=polderivative1 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],f1[[2]],f1[[3]],polderivative2[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb+2]],F[[posx+1,posb+2]],F[[posx+2,posb+2]]]];
f21[[3]]=polderivative2 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],f1[[2]],f1[[3]],polderivative2[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb+2]],F[[posx+1,posb+2]],F[[posx+2,posb+2]]]];f21[[4]]=polderivative2 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],f1[[1]],f1[[4]],polderivative1[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb+2]],F[[posx+1,posb+2]],F[[posx+2,posb+2]]]];



f12=0.5*(f12+f21);


];


If[posx==FFnx-1 && posb==FFnb-1,

f1[[1]]=polderivative2 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb]],F[[posx,posb]],F[[posx+1,posb]]];f1[[2]]=polderivative3 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb]],F[[posx,posb]],F[[posx+1,posb]]];f1[[3]]=polderivative3 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb+1]],F[[posx,posb+1]],F[[posx+1,posb+1]]];f1[[4]]=polderivative2 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb+1]],F[[posx,posb+1]],F[[posx+1,posb+1]]];

f2[[1]]=polderivative2 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx,posb-1]],F[[posx,posb]],F[[posx,posb+1]]];f2[[2]]=polderivative2 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx+1,posb-1]],F[[posx+1,posb]],F[[posx+1,posb+1]]];f2[[3]]=polderivative3 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx+1,posb-1]],F[[posx+1,posb]],F[[posx+1,posb+1]]];f2[[4]]=polderivative3 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx,posb-1]],F[[posx,posb]],F[[posx,posb+1]]];

f12[[1]]=polderivative2 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],polderivative2[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx-1,posb-1]],F[[posx-1,posb]],F[[posx-1,posb+1]]],f2[[1]],f2[[2]]];f12[[2]]=polderivative3 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],polderivative2[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx-1,posb-1]],F[[posx-1,posb]],F[[posx-1,posb+1]]],f2[[1]],f2[[2]]];f12[[3]]=polderivative3 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],polderivative3[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx-1,posb-1]],F[[posx-1,posb]],F[[posx-1,posb+1]]],f2[[4]],f2[[3]]];f12[[4]]=polderivative2 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],polderivative3[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx-1,posb-1]],F[[posx-1,posb]],F[[posx-1,posb+1]]],f2[[4]],f2[[3]]];

f21[[1]]=polderivative2 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],polderivative2[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb-1]],F[[posx,posb-1]],F[[posx+1,posb-1]]],f1[[1]],f1[[4]]];f21[[2]]=polderivative2 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],polderivative3[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb-1]],F[[posx,posb-1]],F[[posx+1,posb-1]]],f1[[2]],f1[[3]]];f21[[3]]=polderivative3 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],polderivative3[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb-1]],F[[posx,posb-1]],F[[posx+1,posb-1]]],f1[[2]],f1[[3]]];f21[[4]]=polderivative3 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],polderivative2[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb-1]],F[[posx,posb-1]],F[[posx+1,posb-1]]],f1[[1]],f1[[4]]];

f12 = 0.5*(f12+f21);
];


If[posx==1 && posb==FFnb-1,

f1[[1]]=polderivative1 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb]],F[[posx+1,posb]],F[[posx+2,posb]]];f1[[2]]=polderivative2 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb]],F[[posx+1,posb]],F[[posx+2,posb]]];f1[[3]]=polderivative2 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb+1]],F[[posx+1,posb+1]],F[[posx+2,posb+1]]];f1[[4]]=polderivative1 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb+1]],F[[posx+1,posb+1]],F[[posx+2,posb+1]]];

f2[[1]]=polderivative2 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx,posb-1]],F[[posx,posb]],F[[posx,posb+1]]];f2[[2]]=polderivative2 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx+1,posb-1]],F[[posx+1,posb]],F[[posx+1,posb+1]]];f2[[3]]=polderivative3 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx+1,posb-1]],F[[posx+1,posb]],F[[posx+1,posb+1]]];f2[[4]]=polderivative3 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx,posb-1]],F[[posx,posb]],F[[posx,posb+1]]];

f12[[1]]=polderivative1 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],f2[[1]],f2[[2]],polderivative2[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx+2,posb-1]],F[[posx+2,posb]],F[[posx+2,posb+1]]]];f12[[2]]=polderivative2 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],f2[[1]],f2[[2]],polderivative2[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx+2,posb-1]],F[[posx+2,posb]],F[[posx+2,posb+1]]]];f12[[3]]=polderivative2 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],f2[[4]],f2[[3]],polderivative3[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx+2,posb-1]],F[[posx+2,posb]],F[[posx+2,posb+1]]]];f12[[4]]=polderivative1 [FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],f2[[4]],f2[[3]],polderivative3[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx+2,posb-1]],F[[posx+2,posb]],F[[posx+2,posb+1]]]];

f21[[1]]=polderivative2 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],polderivative1[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb-1]],F[[posx+1,posb-1]],F[[posx+2,posb-1]]],f1[[1]],f1[[4]]];f21[[2]]=polderivative2 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],polderivative2[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb-1]],F[[posx+1,posb-1]],F[[posx+2,posb-1]]],f1[[2]],f1[[3]]];f21[[3]]=polderivative3 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],polderivative2[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb-1]],F[[posx+1,posb-1]],F[[posx+2,posb-1]]],f1[[2]],f1[[3]]];f21[[4]]=polderivative3 [FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],polderivative1[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb-1]],F[[posx+1,posb-1]],F[[posx+2,posb-1]]],f1[[1]],f1[[4]]];

f12=0.5*(f12+f21);


];


If[posx==FFnx-1 && posb==1,

f1[[1]]=polderivative2 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb]],F[[posx,posb]],F[[posx+1,posb]]];f1[[2]]=polderivative3 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb]],F[[posx,posb]],F[[posx+1,posb]]];f1[[3]]=polderivative3 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb+1]],F[[posx,posb+1]],F[[posx+1,posb+1]]];f1[[4]]=polderivative2 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb+1]],F[[posx,posb+1]],F[[posx+1,posb+1]]];

f2[[1]]=polderivative1 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx,posb]],F[[posx,posb+1]],F[[posx,posb+2]]];f2[[2]]=polderivative1 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx+1,posb]],F[[posx+1,posb+1]],F[[posx+1,posb+2]]];f2[[3]]=polderivative2 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx+1,posb]],F[[posx+1,posb+1]],F[[posx+1,posb+2]]];f2[[4]]=polderivative2 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx,posb]],F[[posx,posb+1]],F[[posx,posb+2]]];

f12[[1]]=polderivative2 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],polderivative1[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx-1,posb]],F[[posx-1,posb+1]],F[[posx-1,posb+2]]],f2[[1]],f2[[2]]];f12[[2]]=polderivative3 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],polderivative1[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx-1,posb]],F[[posx-1,posb+1]],F[[posx-1,posb+2]]],f2[[1]],f2[[2]]];f12[[3]]=polderivative3 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],polderivative2[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx-1,posb]],F[[posx-1,posb+1]],F[[posx-1,posb+2]]],f2[[4]],f2[[3]]];f12[[4]]=polderivative2 [FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],polderivative2[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx-1,posb]],F[[posx-1,posb+1]],F[[posx-1,posb+2]]],f2[[4]],f2[[3]]];

f21[[1]]=polderivative1 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],f1[[1]],f1[[4]],polderivative2[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb+2]],F[[posx,posb+2]],F[[posx+1,posb+2]]]];f21[[2]]=polderivative1 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],f1[[2]],f1[[3]],polderivative3[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb+2]],F[[posx,posb+2]],F[[posx+1,posb+2]]]];f21[[3]]=polderivative2 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],f1[[2]],f1[[3]],polderivative3[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb+2]],F[[posx,posb+2]],F[[posx+1,posb+2]]]];f21[[4]]=polderivative2 [FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],f1[[1]],f1[[4]],polderivative2[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb+2]],F[[posx,posb+2]],F[[posx+1,posb+2]]]];

f12 = 0.5*(f12+f21);


];

If[posx!= 1 && posx!= FFnx-1 && posb != 1 && posb != FFnb-1,

f1[[1]]=polderivative2[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb]],F[[posx,posb]],F[[posx+1,posb]]];f1[[2]]=polderivative2[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb]],F[[posx+1,posb]],F[[posx+2,posb]]];f1[[3]]=polderivative2[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb+1]],F[[posx+1,posb+1]],F[[posx+2,posb+1]]];f1[[4]]=polderivative2[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb+1]],F[[posx,posb+1]],F[[posx+1,posb]]];

f2[[1]]=polderivative2[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx,posb-1]],F[[posx,posb]],F[[posx,posb+1]]];f2[[2]]=polderivative2[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx+1,posb-1]],F[[posx+1,posb]],F[[posx+1,posb+1]]];f2[[3]]=polderivative2[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx+1,posb]],F[[posx+1,posb+1]],F[[posx+1,posb+2]]];f2[[4]]=polderivative2[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx,posb]],F[[posx,posb+1]],F[[posx,posb+2]]];

f12[[1]]=polderivative2[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],polderivative2[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx-1,posb-1]],F[[posx-1,posb]],F[[posx-1,posb+1]]],f2[[1]],f2[[2]]];f12[[2]]=polderivative2[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],f2[[1]],f2[[2]],polderivative2[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],F[[posx+2,posb-1]],F[[posx+2,posb]],F[[posx+2,posb+1]]]];f12[[3]]=polderivative2[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],f2[[4]],f2[[3]],polderivative2[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx+2,posb]],F[[posx+2,posb+1]],F[[posx+2,posb+2]]]];f12[[4]]=polderivative2[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],polderivative2[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],F[[posx-1,posb]],F[[posx-1,posb+1]],F[[posx-1,posb+2]]],f2[[4]],f2[[3]]];

f21[[1]]=polderivative2[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],polderivative2[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb-1]],F[[posx,posb-1]],F[[posx+1,posb-1]]],f1[[1]],f1[[4]]];f21[[2]]=polderivative2[FFbGrid[[posb-1]],FFbGrid[[posb]],FFbGrid[[posb+1]],polderivative2[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb-1]],F[[posx+1,posb-1]],F[[posx+2,posb-1]]],f1[[2]],f1[[3]]];f21[[3]]=polderivative2[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],f1[[2]],f1[[3]],polderivative2[FFxGrid[[posx]],FFxGrid[[posx+1]],FFxGrid[[posx+2]],F[[posx,posb+2]],F[[posx+1,posb+2]],F[[posx+2,posb+2]]]];f21[[4]]=polderivative2[FFbGrid[[posb]],FFbGrid[[posb+1]],FFbGrid[[posb+2]],f1[[1]],f1[[4]],polderivative2[FFxGrid[[posx-1]],FFxGrid[[posx]],FFxGrid[[posx+1]],F[[posx-1,posb+2]],F[[posx,posb+2]],F[[posx+1,posb+2]]]];

f12= 0.5*(f12+f21);

];


bicubint[f,f1,f2,f12,FFxGrid[[posx]],FFxGrid[[posx+1]],FFbGrid[[posb]],FFbGrid[[posb+1]],xA,bT]
]




End[]
EndPackage[]
