(*Remark: Not all the presented options and parameters used for the calculation are necessary. There are probably more efficient methods to implement such derivation in WM14. Therefore, we present here what worked for us, without aiming to optimize the derivations.*)
Clear["Global`*"];

(***************************)
(***************************)
(***************************)
(*We want to approximate point-(infinite) cylinder (P-C) interaction, where the total interaction is described by integrating point-disk (P-D) interaction. The sketch is given in the paper. We model the P-C interaction as an interaction between a point and an infinite number of elliptical plane figures. These plane elliptical figures are sections between a cylinder and planes that pass through a point. The planes are perpendicular to the e1e3-plane. Finally, these elliptical figures are approximated with inscribed disks.*)
(***************************)
(***************************)
(***************************)
SetSystemOptions["SimplificationOptions"->{"AssumptionsMaxNonlinearVariables"->25,"AssumptionsMaxVariables"->25}];
$Assumptions=Rx>0&&Rxk>0&&Ry>0&&p>0&&(*t>=Rx+q2&&*)tk>Rxk&&q2>0&&q1>=0&&p>0&&d>0&&d1>0&&d2>0&&d3>0&&d4>0&&d5>0&&cosa>0;


(*Integration is over a volume. We use modified Langbein coordinates and find the Jacobian.*)

x2[p_,φ_,alfa_]:=t;
y2[p_,φ_,alfa_]:=0;


x1[p_,φ_,alfa_]:=t+p*Cos[φ]*Cos[alfa];
y1[p_,φ_,alfa_]:=p*Sin[φ];
z1[p_,φ_,alfa_]:=p*Cos[φ]*Sin[alfa];

Dg=Grad[{x1[p,φ,alfa],y1[p,φ,alfa],z1[p,φ,alfa]},{p,φ,alfa}];

(*Jacobian of coordinate transformation*)
jac=PowerExpand[FullSimplify[Det[Dg]]]
(*P-P distance*)
Simplify[√((((x2[p,φ,alfa]-x1[p,φ,alfa]))^2)+(((y2[p,φ,alfa]-y1[p,φ,alfa]))^2)+(((0-z1[p,φ,alfa]))^2))]

(*Some geometrical considerations. Let us inscribe a circle inside an ellipse*)
rep={a->2,b->1.5};

(*parametric equation of an ellipse, in 3D and 2D*)
r[a_,b_,t_]={a*Cos[t],b*Sin[t],0};
r2[a_,b_,t_]={a*Cos[t],b*Sin[t]};
(*curvature*)
kk[t_]=((Norm[Cross[D[r[a,b,t],t],D[r[a,b,t],t,t]]])/((Norm[D[r[a,b,t],t]])^3));

ParametricPlot[r2[a,b,t]/.rep,{t,0,2*π}]

(*The radius of inscribed circle is 1/kk[0] and equals the radius of ellipse for t=0. Then, the offset of the inscribed circle is*)
Rxi=1/kk[0];
of=a-Rxi;
ParametricPlot[{r2[a,b,t]/.rep,{of+Rxi*Cos[t],Rxi*Sin[t]}/.rep},{t,0,2*π}]

(*The radius of an ellipse at t=0*)
1/kk[0]

(*The radius of an inscribed circle*)
FullSimplify[Rxi/.{b->Rx,a->Rx/cosa}]

(*Radius Rxi of inscribed disk is Rxi=((Rx^3)/((Rx/cosa)*Rx))=Rx cosa=Rxk. tk coordinate in the kappa-plane is related to t coordinate (distance between cylinder axis and a point) as tk= ((t-Rx+Rx*((Cos[alfa])^2))/(Cos[alfa])). It is derived and sketched in the paper.*)

(*We will derive an approximate P-C law for every even exponent m, m>3. For odd m, we already have nice analytical expressions that follow from applying RIIPI to P-(Subscript[D, IP]) laws, see Notebook 2. P-C interaction requires m>3, due to the 3D integration.*)

tab={};
For[m=4,m<13,m=m+2,

(*first, we integrate of angle φ*)
defFi=Integrate[(p^(-m))*(p^2) Cos[φ],{φ,-ArcCos[(((p^2)+(tk^2)-(Rxk^2))/(2*tk*p))],ArcCos[(((p^2)+(tk^2)-(Rxk^2))/(2*tk*p))]}];
(*Then we integrate over p*)
indefP=Integrate[defFi,p];

u=Limit[indefP,p->tk+Rxk,Direction->"FromBelow"];
l=Limit[indefP,p->tk-Rxk,Direction->"FromAbove"];
defP=FullSimplify[u-l];

(*Finally, we integrate over the angle alpha and the approximate P-C laws follow:*)

indefAlfa=Integrate[(defP/. {Rxk->Rx*Cos[alfa],tk->(t-Rx+Rx*Cos[alfa]^2)/Cos[alfa]}),alfa];
u=Limit[indefAlfa,alfa->π/2,Direction->"FromBelow"];
l=Limit[indefAlfa,alfa->0,Direction->"FromAbove"];
defAlfa=2*FullSimplify[u-l];

AppendTo[tab,defAlfa];
Print[m,"  ",defAlfa]
]

(*The results: the approximate P-(Subscript[C, app]) laws for m=4,6,8,10,12. Let us compare it with the exact ones. The exact P-C laws follow from the P-(Subscript[D, IP]) laws. P-(Subscript[D, IP]) laws are collected in the table pointSecEx for m=3,4,5,6,...11, see Notebook 2. We need to multiply P-(Subscript[D, IP]) law with a factor fm to get the P-C laws for m=4-12*)
pointCylApp={(1/2) π ((2/(-Rx+t))-((4 ArcTanh[(1/(Sqrt[(t/Rx)]))])/(Sqrt[Rx t]))+(((Sqrt[2]) ArcTanh[(((Sqrt[2]) Rx)/(Sqrt[Rx (Rx+t)]))])/(Sqrt[Rx (Rx+t)]))),(π (-166 (Rx^(7/2)) (Sqrt[Rx+t])-86 (Rx^(5/2)) t (Sqrt[Rx+t])-34 (Rx^(3/2)) (t^2) (Sqrt[Rx+t])+30 (t^3) (Sqrt[Rx (Rx+t)])+3 (Sqrt[2]) (((Rx-t))^3) (11 Rx+5 t) ArcTanh[(((Sqrt[2]) Rx)/(Sqrt[Rx (Rx+t)]))]))/(384 (Rx^(3/2)) (((Rx-t))^3) (((Rx+t))^(5/2))),(1/(122880 (Rx^3) ((((Rx^2)-(t^2)))^5)))(-2 π Rx (Rx+t) (12031 (Rx^6)+16474 (Rx^5) t+40941 (Rx^4) (t^2)-1956 (Rx^3) (t^3)-5399 (Rx^2) (t^4)+1450 Rx (t^5)+1995 (t^6))+15 (Sqrt[2]) π (235 (Rx^(15/2)) (Sqrt[Rx+t])-1469 (Rx^(13/2)) t (Sqrt[Rx+t])+2151 (Rx^(11/2)) (t^2) (Sqrt[Rx+t])-529 (Rx^(9/2)) (t^3) (Sqrt[Rx+t])-287 (Rx^(7/2)) (t^4) (Sqrt[Rx+t])-375 (Rx^(5/2)) (t^5) (Sqrt[Rx+t])+141 (Rx^(3/2)) (t^6) (Sqrt[Rx+t])+133 (t^7) (Sqrt[Rx (Rx+t)])) ArcTanh[(((Sqrt[2]) Rx)/(Sqrt[Rx (Rx+t)]))]),(1/(22020096 (Rx^4) (((Rx-t))^7) (((Rx+t))^7)))π (-2 Rx (Rx+t) (1305141 (Rx^9)+2936943 (Rx^8) t+12255444 (Rx^7) (t^2)+3274924 (Rx^6) (t^3)+4143158 (Rx^5) (t^4)-220302 (Rx^4) (t^5)+1459204 (Rx^3) (t^6)+577052 (Rx^2) (t^7)-370755 Rx (t^8)-194985 (t^9))+105 (Sqrt[2]) (((Rx-t))^4) (Sqrt[Rx (Rx+t)]) (3231 (Rx^6)-19642 (Rx^5) t-29327 (Rx^4) (t^2)-41708 (Rx^3) (t^3)-30191 (Rx^2) (t^4)-11578 Rx (t^5)-1857 (t^6)) ArcTanh[(((Sqrt[2]) Rx)/(Sqrt[Rx (Rx+t)]))]),-(1/(12079595520 (Rx^(9/2)) (((Rx-t))^9) (((Rx+t))^(17/2))))π (4 (246027571 (Rx^(25/2)) (Sqrt[Rx+t])+771616208 (Rx^(23/2)) t (Sqrt[Rx+t])+4483631626 (Rx^(21/2)) (t^2) (Sqrt[Rx+t])+2601988760 (Rx^(19/2)) (t^3) (Sqrt[Rx+t])+5958326329 (Rx^(17/2)) (t^4) (Sqrt[Rx+t])+613476416 (Rx^(15/2)) (t^5) (Sqrt[Rx+t])+2242836604 (Rx^(13/2)) (t^6) (Sqrt[Rx+t])+801195920 (Rx^(11/2)) (t^7) (Sqrt[Rx+t])-160039979 (Rx^(9/2)) (t^8) (Sqrt[Rx+t])-449773392 (Rx^(7/2)) (t^9) (Sqrt[Rx+t])-63121254 (Rx^(5/2)) (t^10) (Sqrt[Rx+t])+100083480 (Rx^(3/2)) (t^11) (Sqrt[Rx+t])+33620895 (t^12) (Sqrt[Rx (Rx+t)]))+45 (Sqrt[2]) (((Rx-t))^5) (1314405 (Rx^8)-12139080 (Rx^7) t-23275284 (Rx^6) (t^2)-53318584 (Rx^5) (t^3)-62126850 (Rx^4) (t^4)-47345528 (Rx^3) (t^5)-22645588 (Rx^2) (t^6)-6208776 Rx (t^7)-747131 (t^8)) (Log[-2 (Sqrt[Rx])+(Sqrt[2]) (Sqrt[Rx+t])]-Log[2 (Sqrt[Rx])+(Sqrt[2]) (Sqrt[Rx+t])]))};
pointDiskEx={2 (((EllipticE[((4 Rx t)/(((Rx+t))^2))])/(-Rx+t))-((EllipticK[((4 Rx t)/(((Rx+t))^2))])/(Rx+t))),((π (Rx^2))/((((Rx^2)-(t^2)))^2)),((2 (-((7 (Rx^2)+(t^2)) EllipticE[((4 Rx t)/(((Rx+t))^2))])+(((Rx-t))^2) EllipticK[((4 Rx t)/(((Rx+t))^2))]))/(9 (((Rx-t))^3) (((Rx+t))^2))),((π (Rx^2) ((Rx^2)+2 (t^2)))/(2 ((((Rx^2)-(t^2)))^4))),((2 (-((43 (Rx^4)+82 (Rx^2) (t^2)+3 (t^4)) EllipticE[((4 Rx t)/(((Rx+t))^2))])+(((Rx-t))^2) (13 (Rx^2)+3 (t^2)) EllipticK[((4 Rx t)/(((Rx+t))^2))]))/(75 (((Rx-t))^5) (((Rx+t))^4))),((π (Rx^2) ((Rx^4)+6 (Rx^2) (t^2)+3 (t^4)))/(3 ((((Rx^2)-(t^2)))^6))),-(1/(735 (((Rx-t))^7) (((Rx+t))^6)))2 ((337 (Rx^6)+1773 (Rx^4) (t^2)+947 (Rx^2) (t^4)+15 (t^6)) EllipticE[((4 Rx t)/(((Rx+t))^2))]-(((Rx-t))^2) (127 (Rx^4)+242 (Rx^2) (t^2)+15 (t^4)) EllipticK[((4 Rx t)/(((Rx+t))^2))]),((π (Rx^2) ((Rx^6)+12 (Rx^4) (t^2)+18 (Rx^2) (t^4)+4 (t^6)))/(4 ((((Rx^2)-(t^2)))^8))),-(1/(2835 (((Rx-t))^9) (((Rx+t))^8)))2 ((1091 (Rx^8)+11156 (Rx^6) (t^2)+16434 (Rx^4) (t^4)+4052 (Rx^2) (t^6)+35 (t^8)) EllipticE[((4 Rx t)/(((Rx+t))^2))]-(((Rx-t))^2) (461 (Rx^6)+2313 (Rx^4) (t^2)+1287 (Rx^2) (t^4)+35 (t^6)) EllipticK[((4 Rx t)/(((Rx+t))^2))])};
pointCylEx=Flatten[Table[(Sqrt[π])*((Gamma[((i+2)/2)])/(Gamma[((i+3)/2)]))*pointDiskEx[[i]],{i,1,Length[pointDiskEx]}]];

(*The relative error of the approximate P-C laws for m=4,6,8,10,12. For small separations, the relative error is small, for large separations, it is bounded to approximately 30%.*)
tab1={};
tab2={};
For[i=1,i<6,i++,

error=((pointCylEx[[2*i-1]]-pointCylApp[[i]])/(pointCylEx[[2*i-1]]))*100/.{t->Rx+q2}/.{Rx->1};
AppendTo[tab1,error];
AppendTo[tab2,2+2*i]]
LogLogPlot[tab1,{q2,(10^(-4)),(10^3)},PlotRange->All,AxesLabel->{"q2","Error [%]"},PlotLegends->tab2,WorkingPrecision->30,PlotLabel->"Error of the P-(Subscript[C, app]) law for even m"]

(*Scaling of the approximate P-(Subscript[C, app]) laws for m=4,6,8,10,12. For small separations,scaling is m-3, for large, it is m-1.*)
tab1={};
tab2={};
For[i=1,i<6,i++,f=-D[Log[pointCylApp[[i]]/.t->Rx+q2/. {q2->Exp[q2]}],q2]/. {Rx->1};
AppendTo[tab1,f];
AppendTo[tab2,2*i+2]];

Plot[tab1,{q2,-10,10},AxesLabel->{"Log[q2]","Scaling factor"},PlotLegends->tab2,Evaluated->True,WorkingPrecision->80]

(***************************)
(***************************)
(***************************)
(*Now, we want to derive an approximate disk-(infinite) cylinder law (or unit-length cylinder-infinite cylinder), D-C, law for m=6. We can do it for other exponents as well, but they require different approximations*)
(*We are using RCPS1. In essence, we approximate the D-C integral, Integrate[(Subscript[Π, P-C])*2*ArcCos[((-(Ry^2)+(((q2+Rx+Ry))^2)+(t^2))/(2 (q2+Rx+Ry) t))]*t,{t,Rx+q2,q2+Rx+2*Ry}]. The first approximation is to replace (Subscript[Π, P-C]) with our approximate P-C interaction potential. However, additional approximations are required. *)
(***************************)
(***************************)
(*After careful consideration, it turns out that, for Rx>=Ry, the angle psi is relatively small. We can show it both graphically and analytically.
*)

(*extremum of the psi angle is at the minimum of Cos[psi]*)
sol=Solve[D[((-(Ry^2)+(((q2+Rx+Ry))^2)+(t^2))/(2 (q2+Rx+Ry) t)),t]==0,t]
(*insert the value of t where Cos[psi] is maximum, and plot the max value of Cos[psi] for different ratios s=Rx/Ry and q2*)
maxCosPsi=FullSimplify[((-(Ry^2)+(((q2+Rx+Ry))^2)+(t^2))/(2 (q2+Rx+Ry) t))/.t->(Sqrt[(q2+Rx) (q2+Rx+2 Ry)])]
Plot3D[maxCosPsi/.Ry->s*Rx/.Rx->1,{q2,1/10000,10},{s,1/100,1},AxesLabel->{"q2","s"},PlotRange->All]

(*This allows us to approximate the arc length (t ArcCos[psi]) with its series expansion at t=Rx+q2. We have compared several approaches and here present only the one that provides reasonable balance between efficiency and accuracy*)

pointCylApp6=(1/(384 (Rx^(3/2)) (((Rx-t))^3) (((Rx+t))^(5/2))))π (30 (t^3) (Sqrt[Rx (Rx+t)])-34 (t^2) (Sqrt[(Rx^3) (Rx+t)])-86 t (Sqrt[(Rx^5) (Rx+t)])-166 (Sqrt[(Rx^7) (Rx+t)])+3 (Sqrt[2]) (((Rx-t))^3) (11 Rx+5 t) ArcTanh[(Sqrt[2]) (Sqrt[(Rx/(Rx+t))])]);
pointCylEx6=((π (-((7 (Rx^2)+(t^2)) EllipticE[((4 Rx t)/(((Rx+t))^2))])+(((Rx-t))^2) EllipticK[((4 Rx t)/(((Rx+t))^2))]))/(12 (((Rx-t))^3) (((Rx+t))^2)));

(*A first order series approximation of the arc length ArcCos[x]*t at t=Rx+q2. *)
appACt=FullSimplify[Normal[Series[ArcCos[((-(Ry^2)+(((q2+Rx+Ry))^2)+(t^2))/(2 (q2+Rx+Ry) t))]*t,{t,Rx+q2,1}]]]


(*Another question is where to approximate the ArcTanh function. If we chose Infinity, we obtain elliptic integrals as solutions. If we choose Rx+q2+2*Ry or Rx+q2+Ry, we end up with a more complicated expressions with Log and ArcTanh functions, but without elliptic integrals. Since WM14 calculates elliptic integrals efficiently, we choose the first option.*)
appATanh=FullSimplify[Normal[Series[ArcTanh[(Sqrt[2]) (Sqrt[(Rx/(Rx+t))])],{t,Infinity,2}]]]

(*Therefore, by using these approximations, the definite integral is*)

indef=Simplify[Total[Integrate[FullSimplify[Level[Collect[pointCylApp6,{ArcTanh[_]}]/.ArcTanh[_]->appATanh,1]*2*appACt],t]]]
l=Limit[indef,t->Rx+q2,Direction->"FromAbove"]
u=Limit[indef,t->Rx+q2+2*Ry,Direction->"FromBelow"]
def=u-l;
FullSimplify[Expand[def]]

(*The result: the approximate disk-cylinder vdW law, D-(Subscript[C, app])*)
diskCylApp6=(1/1920)π (Sqrt[((q2+Rx)/(q2+Rx+Ry))]) ((160/q2)-(160/(q2+2 Ry))-(77/(Sqrt[(q2+Rx+2 Ry) (q2+2 (Rx+Ry))]))+((11 (((7 (Sqrt[(q2+Rx+2 Ry) (q2+2 (Rx+Ry))]))/(q2+Rx))-((22 (Sqrt[(q2+Rx+2 Ry) (q2+2 (Rx+Ry))]))/(q2+2 Rx))+15 (Sqrt[1-(Rx/(q2+2 (Rx+Ry)))])))/Rx)+(1/(Rx^2))Ry (3024 (Sqrt[1-(Rx/(q2+2 (Rx+Ry)))])+2 Rx (((160 (3 q2-2 Rx+6 Ry))/(((q2+2 Ry))^2))+242 (Sqrt[((q2+Rx+2 Ry)/(((q2+2 (Rx+Ry)))^3))])+(180/(q2+2 (Rx+Ry)))-77 (Sqrt[((q2+2 (Rx+Ry))/(((q2+Rx+2 Ry))^3))])+(291/(Sqrt[(q2+Rx+2 Ry) (q2+2 (Rx+Ry))]))))+(1/((Rx^2) (q2+2 Rx) (((q2 Ry))^(3/2))))20 (Sqrt[2]) (Ry^2) (-2 (q2+2 Rx) (9 (q2^2)+12 q2 Rx-4 (Rx^2)) ArcTan[(Sqrt[2]) (Sqrt[(Ry/q2)])]+9 (3 Rx (Sqrt[(q2^3) (q2+2 Rx)])+2 (Sqrt[(q2^5) (q2+2 Rx)])) ArcTan[(Sqrt[2]) (Sqrt[(Ry/(q2+2 Rx))])])+(1/((Rx^2) (q2+Rx) (Sqrt[(q2+2 Rx) Ry])))I (Sqrt[2]) Ry ((1512 (q2^2)+4371 q2 Rx+2936 (Rx^2)) EllipticE[I ArcSinh[(Sqrt[2]) (Sqrt[(Ry/(q2+Rx))])],((q2+Rx)/(q2+2 Rx))]-11 Rx (96 q2+103 Rx) EllipticF[I ArcSinh[(Sqrt[2]) (Sqrt[(Ry/(q2+Rx))])],((q2+Rx)/(q2+2 Rx))]));

(*Let us calculate the relative error for different ratios of Ry/Rx*)

tab2={1/4,1/2,1,2,4};
tab3={};
tab4={};

exp=Range[-4,5/2,1/32];
For[j=1,j<Length[tab2]+1,j++,
tab1={};
tab5={};
For[i=1,i<Length[exp]+1,i++,
r={Rx->1,Ry->tab2[[j]],q2->(10^(exp[[i]]))};
ex=Re[NIntegrate[pointCylEx6*2*ArcCos[(((t^2)+(((q2+Rx+Ry))^2)-(Ry^2))/(2*t*(q2+Rx+Ry)))]*t/.r,{t,Rx+q2,q2+Rx+2*Ry}/.r]];

app=Re[diskCylApp6/.r];
(*Print[{ex,1.*app}];*)

error={(10^(exp[[i]])),N[Abs[((ex-app)/ex)]*100,160]};
AppendTo[tab5,{error,ex}];

AppendTo[tab1,error]];
AppendTo[tab4,tab5];
AppendTo[tab3,tab1]

]

ListLogLogPlot[tab3,Joined->True,PlotLegends->tab2,PlotLabel->"Relative error of the approximate ((Subscript[Π, D - C])^vdW) for different ratios Ry/Rx",AxesLabel->{"q2","Relative error [%]"}]

(*Scaling is correct, 3/2 for small and 5 for large separations*)

f025 = -D[Log[diskCylApp6/.{q2->Exp[q2]}],q2]/.{Rx->1,Ry->1/4};
f05 = -D[Log[diskCylApp6/.{q2->Exp[q2]}],q2]/.{Rx->1,Ry->1/2};
f1 = -D[Log[diskCylApp6/.{q2->Exp[q2]}],q2]/.{Rx->1,Ry->1};
f2= -D[Log[diskCylApp6/.{q2->Exp[q2]}],q2]/.{Rx->1,Ry->2};
f4= -D[Log[diskCylApp6/.{q2->Exp[q2]}],q2]/.{Rx->1,Ry->4};

a025=Table[{q2,Re[N[f025,100]]},{q2,-10,10,1/4}];
a05=Table[{q2,Re[N[f05,100]]},{q2,-10,10,1/4}];
a1=Table[{q2,Re[N[f1,100]]},{q2,-10,10,1/4}];
a2=Table[{q2,Re[N[f2,100]]},{q2,-10,10,1/4}];
a4=Table[{q2,Re[N[f4,100]]},{q2,-10,10,1/4}];
tabX={a025,a05,a1,a2,a4};

ListLinePlot[tabX,AxesLabel->{"gap: Log[q2]","Scaling factor of ΠvdW"}]

(*Replacements in the expression for the paper. Few additional cancelations are done in the paper*)
defFSd=diskCylApp6/.{q2+Rx+Ry->d,q2+Rx+2 Ry->d1,q2+2 (Rx+Ry)->d2,q2+Rx->d3,q2+2*Rx->d4,q2+2 Ry->d5}/.{1-(Rx/d2)->(d1/d2)}

(******************************)
(******************************)
(*Let us now derive some special cases of interaction, for general m. First, from point-infinite plane law, we can get P-PT law, where thickness of the plate is 2by. Then, by letting by->Infinity, it is a P-HS law. See also Notebook 1. *)
(******************************)
(******************************)
Clear[m];
$Assumptions=q2>0&&m>3&&by>0;
pointPN=((2 π (((-x1+y1))^(2-m)))/(-2+m));

pointPTa[m_,q2_]=FullSimplify[Integrate[pointPN,{y1,x1+q2,x1+q2+2*by}]]

pointHSa[m_,q2_]=FullSimplify[Limit[pointPTa[m,q2],by->Infinity]]

(*We can do it also vice versa, by first deriving P-HS, and then finding P-PT law.**)

pointPN=((2 π (((-x1+y1))^(2-m)))/(-2+m));

pointHS[m_,q2_]=FullSimplify[Integrate[pointPN,{y1,x1+q2,Infinity}]]

pointPT[m_,q2_]=FullSimplify[pointHS[m,q2]-pointHS[m,q2+2 by]]

(*The results are the same, just in a different form*)
FullSimplify[pointPTa[m,q2]-pointPT[m,q2]]

(*Now, let us derive a general D-HS law. We use Cartesian coordinates (x1,x2) with the origin at the boundary of the HS. The limits are x1 ∈ [q2,q2+2Rx] and x2 ∈ [-(Sqrt[(Rx^2)-(((Rx+q2-x1))^2)]),(Sqrt[(Rx^2)-(((Rx+q2-x1))^2)])]*)
$Assumptions=q2>0&&x1>0&&Rx>0&&m>4&&by>0;
defX2=FullSimplify[Integrate[(pointHS[m,q2])/. q2->x1,{x2,-Sqrt[Rx^2-(Rx+q2-x1)^2],Sqrt[Rx^2-(Rx+q2-x1)^2]}]]
indef=Integrate[defX2,x1]
l=Limit[indef,x1->q2,Direction->"FromAbove"]
u=Limit[indef,x1->q2+2*Rx,Direction->"FromBelow"]
def=FullSimplify[u-l]

diskHS[m_,q2_,Rx_]=(((π^2) (q2^(-m)) (Rx^2) (4 (q2^3) Hypergeometric2F1[(3/2),m,3,-((2 Rx)/q2)]+3 q2 Rx (4 q2 Hypergeometric2F1[(5/2),m,4,-((2 Rx)/q2)]+5 Rx Hypergeometric2F1[(7/2),m,5,-((2 Rx)/q2)])+7 (Rx^3) Hypergeometric2F1[(9/2),m,6,-((2 Rx)/q2)]))/(2 (-3+m) (-2+m)));

(*For special cases, integer m, this law returns simple rational functions*)
FullSimplify[diskHS[#,q2,Rx]&/@Range[5,12]]

(*From the D-HS laws, it is straightforward to find the D-PT laws*)
diskPlate[m_,q2_,Rx_,by_]:=FullSimplify[diskHS[m,q2,Rx]]-FullSimplify[diskHS[m,q2+2 by,Rx]]
FullSimplify[diskPlate[#,q2,Rx,by]&/@Range[5,12]]

(*The expression is valid for arbitrary m, but the law is valid only for m>4.5, due to the asymptotic scaling at small separations*)
rep={m->5.2,q2->1/1000,Rx->1,by->2};
NIntegrate[pointPT[m,q2]/. q2->x1/.rep,{x1,q2,q2+2*Rx}/.rep,{x2,-Sqrt[Rx^2-(Rx+q2-x1)^2],Sqrt[Rx^2-(Rx+q2-x1)^2]}/.rep,AccuracyGoal->20]
1.*diskPlate@@rep[[All,2]]

(*Scaling of the exact disk-HS laws for m=5,6,7,...12. For small separations, the scaling factor is m-4.5, for large, m-3. By multiplying with 2Lx, it is finite cylinder-HS law, where the length of the cylinder is 2Lx*)
tab1={};
tab2={};
For[i=5,i<12+1,i++,
f=-D[Log[diskHS[i,q2,Rx]/. {q2->Exp[q2]}],q2]/. {Rx->1};
AppendTo[tab1,f];
AppendTo[tab2,i]];

Plot[tab1,{q2,-10,10},AxesLabel->{"Log[q2]","Scaling factor"},PlotLegends->tab2,Evaluated->True,WorkingPrecision->40]

(*Scaling of the exact D-PT laws for m=5,6,7,...12. For small separations, the scaling factor is m-4.5 (same as for D-HS), for large separations, it is m-2 (different than for D-HS). By multiplying with 2Lx, it is finite cylinder-(infinite) plate law, where the length of a cylinder is 2Lx*)
tab1={};
tab2={};
For[i=5,i<12+1,i++,
f=-D[Log[diskPlate[i,q2,Rx,by]/. {q2->Exp[q2]}],q2]/. {Rx->1,by->1};
AppendTo[tab1,f];
AppendTo[tab2,i]];

Plot[tab1,{q2,-10,10},AxesLabel->{"Log[q2]","Scaling factor"},PlotLegends->tab2,Evaluated->True,WorkingPrecision->40]

(*Scaling of the exact D-PT laws for m=6 and different values of plate thickness by*)
tab1={};
tab2={};
For[i=1,i<5+1,i++,
f=-D[Log[diskPlate[6,q2,Rx,by]/. {q2->Exp[q2]}],q2]/. {Rx->1,by->(10^(i-2))};
AppendTo[tab1,f];
AppendTo[tab2,(10^(i-2))]];

Plot[tab1,{q2,-10,15},AxesLabel->{"Log[q2]","Scaling factor"},PlotLegends->tab2,Evaluated->True,WorkingPrecision->40]

(********************************)
(********************************)
(*Finally, let us consider the approximate D-HS laws that follow from approximate D-C law for m=6. We observe two variants, one where Rx->Infinity, and the other where Ry->Infinity. The first one is much simpler and in line with our assumption Rx>Ry. However, it turns out that the approximate D-HS law obtained with Ry->Infinity is more accurate for large separations. However, the exact D-HS is much simpler than both approximate ones!*)
(********************************)
(********************************)
diskHS6appRx=Limit[diskCylApp6,Rx->Infinity]/.Ry->Rx
diskHS6appRy=Limit[diskCylApp6,Ry->Infinity]
diskHS6ex=(((π^2) (Rx^2))/(6 (((q2 (q2+2 Rx)))^(3/2))))
LogLogPlot[{(Abs[((diskHS6ex-diskHS6appRx)/diskHS6ex)]*100/.Rx->1),(Abs[((diskHS6ex-diskHS6appRy)/diskHS6ex)]*100/.Rx->1)},{q2,(10^(-4)),(10^3)},PlotLegends->{"Rx->Infinity","Ry->Infinity"}]