(*** Testing 2D scalar conformal blocks ***)
(*** Define projective lightcone variables for the four operators ***)
Do[X[i]=Table[x[i,j],{j,1,4}],{i,1,4}];

(*** Lorentzian inner product on the auxiliary (R^(3,1)) space ***)
Inn[Y_,Z_]:=-Y[[1]] Z[[1]]+Sum[Y[[i]] Z[[i]],{i,2,4}]  

(*** Note that at this stage we haven't yet restricted to the lightcone ***)

(*** The scalar 4-point conformal block, expressed in terms of G[u,v] ***)
FourPointFunction=Inn[X[1],X[2]]^(-(Δ1+Δ2)/2) Inn[X[3],X[4]]^(-(Δ3+Δ4)/2) (Inn[X[2],X[4]]/Inn[X[1],X[4]])^((Δ1-Δ2)/2)(Inn[X[1],X[4]]/Inn[X[1],X[3]])^((Δ3-Δ4)/2) G[Inn[X[1],X[2]] Inn[X[3],X[4]]/(Inn[X[1],X[3]] Inn[X[2],X[4]]),Inn[X[1],X[4]] Inn[X[2],X[3]]/(Inn[X[1],X[3]] Inn[X[2],X[4]])];

(*** Define the conformal generators L[i,j] in projective lightcone representation, acting on the first two fields in the 4-point function, with coordinates X[1] and X[2] ***)
L[i_,j_][f_]:=If[i==1,-1,1] (x[1,i] D[f,x[1,j]]+x[2,i] D[f,x[2,j]])-If[j==1,-1,1] (x[1,j] D[f,x[1,i]]+x[2,j] D[f,x[2,i]])

(*** L2[f] is minus the Casimir operator, (L^2)=(1/2)(L^AB)(Subscript[L, AB]), acting on the function f ***)
L2[f_]:=1/2 Sum[If[i==1,-1,1] If[j==1,-1,1] L[i,j][L[i,j][f]],{i,1,4},{j,1,4}]

(*** We will set X[1]={1,1,0,0}, X[2]={((1+(a^2)+(b^2))/2),((1-(a^2)-(b^2))/2),a,b}, X[3]={1,0,1,0}, X[4]={((1+(ϵ^(-2)))/2),((1-(ϵ^(-2)))/2),0,(1/ϵ)} ***)
CoordinateSubstitution={x[1,1]->1,x[1,2]->1,x[1,3]->0,x[1,4]->0,x[2,1]->(1+a^2+b^2)/2,x[2,2]->(1-a^2-b^2)/2,x[2,3]->a,x[2,4]->b,x[3,1]->1,x[3,2]->0,x[3,3]->1,x[3,4]->0,x[4,1]->(1+1/ϵ^2)/2,x[4,2]->(1-1/ϵ^2)/2,x[4,3]->0,x[4,4]->1/ϵ};
(*** After taking ϵ to zero limit, (a,b) will be related to z by z=a+ib ***)

(*** This expression is supposed to be equal to be the Casimir eigenvalue when acting on the conformal block ***)
Casimir=Assuming[a>0 && b>0 && ϵ>0 && Δ1>0 && Δ2>0 && Δ3>0 && Δ4>0,(-L2[FourPointFunction]/FourPointFunction/.CoordinateSubstitution//Simplify)/.ϵ->0]

(*** Changing variables, from G[u,v] to F[z,zbar] ***)
G[u_,v_]:=F[(1/2) (1+u+(Sqrt[-4 u+(((1+u-v))^2)])-v),(1/2) (1+u-(Sqrt[-4 u+(((1+u-v))^2)])-v)]
(*** verifying that the change of coordinate is correct ***)
Assuming[b>0,G[(a^2)+(b^2),1-2 a+(a^2)+(b^2)]//Simplify]

(*** Casimir eigenvalue, now expressed in terms of F[z,zbar] ***)
CasimirZ=Assuming[b>0,(Casimir//Simplify)/.{a->(z+zbar)/2,b->(z-zbar)/(2I)}//Simplify]

(*** Checking the leading term in the conformal block from OPE ***)
H[z_,zbar_]:=z^((Δ+s)/2) zbar^((Δ-s)/2)

((CasimirZ/.F->H)//Simplify)/.{z->0,zbar->0}//Simplify

(*** The 2D scalar conformal block ***)
Block2D[z_,zbar_]:=z^((Δ+s)/2) zbar^((Δ-s)/2) Hypergeometric2F1[(Δ+s)/2-(Δ1-Δ2)/2,(Δ+s)/2+(Δ3-Δ4)/2,Δ+s,z] Hypergeometric2F1[(Δ-s)/2-(Δ1-Δ2)/2,(Δ-s)/2+(Δ3-Δ4)/2,Δ-s,zbar]

(*** Checking that it has the correct Casimir by Taylor expansion in z, zbar ***)
(CasimirZ-((s^2)+(-2+Δ) Δ)/.F->Block2D/.{z->ϵ z,zbar->ϵ zbar}/.{Hypergeometric2F1[a_,b_,c_,ϵ x_]:>(Hypergeometric2F1[a,b,c,ϵ x]+O[ϵ]^4)})//Simplify

(*** Checking that it has the correct Casimir by plugging in random numerical values ***)
CasimirZ-((s^2)+(-2+Δ) Δ)/.F->Block2D/.{z->.3,zbar->.4,Δ->1.3,s->1,Δ1->1.7,Δ2->1.8,Δ3->1.22,Δ4->1.25}

(*** testing numerical values of the conformal block itself ***)
Plot[(Block2D[x,x])/.{Δ1->1/8,Δ2->1/8,Δ3->1/8,Δ4->1/8,Δ->1.5,s->0},{x,0,1}]

(*** After setting Δ1=Δ2=Δ3=Δ4, contribution to the crossing equation, (v^(Δ1))F[u,v]-(u^(Δ1))F[v,u] ***)
CrossingBlock[z_,zbar_]=(1-z)^Δ1(1-zbar)^Δ1 Block2D[z,zbar]-z^Δ1 zbar^Δ1 Block2D[1-z,1-zbar]/.{Δ2->Δ1,Δ3->Δ1,Δ4->Δ1}//Simplify

(*** Extracting the x and x^3 terms in the z=zbar=1/2+x expansion of CrossingBlock ***)
P2D:=(series=CrossingBlock[1/2+x,1/2+x]+O[x]^4//N;{SeriesCoefficient[series,1],SeriesCoefficient[series,3]})

(*** Some numerical values of P2D, with Δ1=2, at different values of Δ, spin s=0 ***)
Table[{weight,P2D/.{Δ1->1/8,Δ->weight,s->0}},{weight,0,6,1}]

maxweight=30;weightstep=0.05;
PlotP2D[weight1_,spin_,gap_,color_]:=ListPlot[Table[P2D/1.1^(spin)/2^weight/.{Δ1->weight1,Δ->weight,s->spin},{weight,Max[spin,gap],maxweight,weightstep}],PlotRange->{2{-1,1},5{-1,1}},AspectRatio->1,Axes->True,Frame->True,PlotStyle->{PointSize[.01],color}]
PlotIdentity[weight1_]:=ListPlot[{P2D/.{Δ1->weight1,Δ->0,s->0}},PlotRange->{5{-1,1},50{-1,1}},AspectRatio->1,Axes->True,Frame->True,PlotStyle->{PointSize[.015],Black}]

GapPlot[weight1_,gap_]:=Show[PlotP2D[weight1,0,gap,Red],PlotP2D[weight1,2,0,Brown],PlotP2D[weight1,4,0,Blue],PlotIdentity[weight1]]

(*** Plotting the order x and x^3 components of CrossingBlock (suitably rescaled for visualization) for external weight Δ1=1/8, a range of internal weight Δ, spin s=0,2,4 (in red, brown, blue) with scalar gap 0, 1, and 2 respectively ***)
GapPlot[1/8,0]
GapPlot[1/8,1]
GapPlot[1/8,2]