(*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.*)

(*allow more assumptions*)
SetSystemOptions["SimplificationOptions"->{"AssumptionsMaxNonlinearVariables"->25,"AssumptionsMaxVariables"->25}];
(*Expressions for in-plane interaction of two rectangular sections in parallel orientation. Dimensions of rectangles are 2*bx/2 hx and 2 by/ 2 hy. *)
$Assumptions=bx>0 && by>0&&hx>0 && hy>0&& q2>0&&  q1 ∈ Reals && bx+q2+2*by>=y1>=bx+q2&&bx+q2>x1 &&Element[m,Integers]&&Element[v,Integers]&&m>3&&v>1&&p>0 && Element[x1,Reals]&& Element[x2,Reals]&& Element[y1,Reals]&& Element[y2,Reals]&&y1>x1&&a>0&&b>0;

(*A function that helps integration*)
integrateD[x_,v_,l_,u_]:=Module[{int=x,var=v,low=l,up=u},
indef=FullSimplify[Integrate[int,var],TimeConstraint->1];
upper=Simplify[Limit[indef,var->up,Direction->"FromBelow"],TimeConstraint->1];
lower=Simplify[Limit[indef,var->low,Direction->"FromAbove"],TimeConstraint->1];
res=FullSimplify[upper-lower,TimeConstraint->1];
Return[res]]

(**************************************)
(******************rectangle-rectangle in-plane for arbitrary m********************)
(**************************************)
(*It is possible to derive R-(Subscript[R, IP]) laws for specific value of m. However, it turns out that we can find a general expression for arbitrary integer m. For this, let us first derive an infinite half-strip - infinite half-strip in-plane law, and then we can apply RIHBI. The integration in WM14 does not work for an arbitrary integer m. On the other hand, it works if we distinguish between odd and even m. First, we are considering the odd values, m=2*v+1. We move the origin of the coordinate system to (bx,0,0), that is, at the edge of the rectangle x.*)
$Assumptions=bx>0 && by>0&&hx>0 && hy>0&& q2>0&&  q1 ∈ Reals && y1>=q2&&0>=x1 &&Element[m,Integers]&&Element[v,Integers]&&m>3&&v>1&&p>0 (*&& Element[x1,Reals]&& Element[x2,Reals]&& Element[y1,Reals]&& Element[y2,Reals]&&y1>x1*)&&a>0&&b>0;

defX2odd=integrateD[1/((y1-x1)^2+(y2-x2)^2)^((2*v+1)/2),x2,-hx,hx]
defX2X1odd=integrateD[defX2odd,x1,-Infinity,0]
defX2X1Y2odd=integrateD[defX2X1odd,y2,-hy,hy]
defX2X1Y2Y1odd=integrateD[defX2X1Y2odd,y1,q2,Infinity]

(*Let us replace v=((m-1)/2), so that we have in-plane infinite half-strip-infinite half-strip for m=5,7,9,...*)
finalOdd=defX2X1Y2Y1odd/.v->((m-1)/2)

m

(*The result*)
finalOdd=-(1/(2 (((2-m))^2)))(((2 (-2+m) ((((((hx-hy))^2)+(q2^2)))^((3/2)+((1-m)/2))))/((-1+(1/2) (-1+m)) (-4+m)))-((2 (-2+m) ((((((hx+hy))^2)+(q2^2)))^((3/2)+((1-m)/2))))/((-1+(1/2) (-1+m)) (-4+m)))+4 (((hx-hy))^2) (q2^(2-m)) Hypergeometric2F1[(1/2),-(1/2)+(1/2) (-1+m),(3/2),-((((hx-hy))^2)/(q2^2))]-4 (((hx+hy))^2) (q2^(2-m)) Hypergeometric2F1[(1/2),-(1/2)+(1/2) (-1+m),(3/2),-((((hx+hy))^2)/(q2^2))]+(((hx-hy))^2) (q2^(2-m)) Gamma[-1+(1/2) (-1+m)] Hypergeometric2F1Regularized[-(1/2)+(1/2) (-1+m),(1/2) (-1+m),1+(1/2) (-1+m),-((((hx-hy))^2)/(q2^2))]-(((hx+hy))^2) (q2^(2-m)) Gamma[-1+(1/2) (-1+m)] Hypergeometric2F1Regularized[-(1/2)+(1/2) (-1+m),(1/2) (-1+m),1+(1/2) (-1+m),-((((hx+hy))^2)/(q2^2))]);

(*Now we introduce functions and make a some manual simplifications.*)
k[s_]=hx+s*hy;
f[a_,s_]=s (2/((((m-2))^2) (m-1)(m-3))) *((k[s])^2) (a^(2-m)) Hypergeometric2F1[((m-2)/2),((m-1)/2) ,((m+1)/2) ,-(((k[s])^2)/(a^2))];
g[a_,s_]=s*(2/(((m-2))^2)) ((k[s])^2) (a^(2-m)) Hypergeometric2F1[(1/2),((m-2)/2) ,(3/2),-(((k[s])^2)/(a^2))];
u[a_,s_]=s*(2/((m-4) (m-3) (m-2)))(((((k[s])^2)+(a^2)))^(2-(m/2)));

(*test*)
FullSimplify[f[q2,1]+f[q2,-1]+g[q2,1]+g[q2,-1]+u[q2,1]+u[q2,-1]-finalOdd]

(*The final form of the in-plane infinite half-strip-infinite half-strip law for odd m*)
infHSTInfHSTOdd[q2_]=f[q2,1]+f[q2,-1]+g[q2,1]+g[q2,-1]+u[q2,1]+u[q2,-1];

(*Now, let us consider even values m=2*v*)
defX2even=integrateD[1/((y1-x1)^2+(y2-x2)^2)^((2*v)/2),x2,-hx,hx]
defX2X1even=integrateD[defX2even,x1,-Infinity,0]
defX2X1Y2even=integrateD[defX2X1even,y2,-hy,hy]
defX2X1Y2Y1even=integrateD[defX2X1Y2even,y1,q2,Infinity]

(*Let us replace v=(m/2), so that we have the infinite half-strip-infinite half-strip law for m=4,6,8,...*)
finalEven=defX2X1Y2Y1even/.v->(m/2)

(*The result*)
finalEven=(1/(2 (3-m) (((-1+(m/2)))^2)))((1/(-2+(m/2)))(-1+(m/2)) ((((hx^4)+2 (hx^2) (-(hy^2)+(q2^2))+((((hy^2)+(q2^2)))^2)))^(-m/2)) (-((((((hx-hy))^2)+(q2^2)))^(m/2)) ((((((hx+hy))^2)+(q2^2)))^2)+((((((hx-hy))^2)+(q2^2)))^2) ((((((hx+hy))^2)+(q2^2)))^(m/2)))+(1/(-1+m))(q2^(2-m)) ((((hx-hy))^2) (3+2 (-2+(m/2)) m) Hypergeometric2F1[(1/2),-1+(m/2),(3/2),-((((hx-hy))^2)/(q2^2))]-(((hx+hy))^2) (3+2 (-2+(m/2)) m) Hypergeometric2F1[(1/2),-1+(m/2),(3/2),-((((hx+hy))^2)/(q2^2))]+(((hx-hy))^2) Hypergeometric2F1[-1+(m/2),-(1/2)+(m/2),(1/2)+(m/2),-((((hx-hy))^2)/(q2^2))]-(((hx+hy))^2) Hypergeometric2F1[-1+(m/2),-(1/2)+(m/2),(1/2)+(m/2),-((((hx+hy))^2)/(q2^2))]));

(*We would prefer to write the first term in brackets in a different form*)
term=finalEven[[4,1]]
n=Numerator[term]
d=Denominator[term]
termNew=FullSimplify[((n[[1]]*n[[2,1]])/d)]+FullSimplify[((n[[1]]*n[[2,2]])/d)]
FullSimplify[term-termNew]

(*Functions and simplifications*)
k[s_]=hx+s*hy;
f[a_,s_]=s*((k[s])^2) (2/(((-2+m))^2)) (a^(2-m)) Hypergeometric2F1[(1/2),-1+(m/2),(3/2),-(((k[s])^2)/(a^2))];
g[a_,s_]=s*(2/((-3+m) (((-2+m))^2) (-1+m))) ((k[s])^2) (a^(2-m)) Hypergeometric2F1[-1+(m/2),-(1/2)+(m/2),(1/2)+(m/2),-(((k[s])^2)/(a^2))];
u[a_,s_]=s*((-2+m)/(2 (m-3) (((-1+(m/2)))^2) (-4+m)))(((((k[s])^2)+(a^2)))^(2-(m/2)));
(*test*)
FullSimplify[f[q2,1]+f[q2,-1]+g[q2,1]+g[q2,-1]+w[q2,1]+w[q2,-1]-finalEven]

(*The final form of the in-plane infinite half-strip-infinite half-strip law*)
infHSTInfHSTEven[q2_]=f[q2,1]+f[q2,-1]+g[q2,1]+g[q2,-1]+u[q2,1]+u[q2,-1];

(*However, it turns out that both expressions, for odd and even exponent m are the same!*)
FullSimplify[infHSTInfHSTEven[q2]-infHSTInfHSTOdd[q2]]
FullSimplify[finalEven-finalOdd]

(*It is interesting to observe that the derived expression is well defined only for m>4. For m=4, it has a problem*)
infHSTInfHSTOdd[q2]/.m->4

(*However, the limit m->4 exists, so we can find it and define our final function as a piecewise expression*)
specm4=Evaluate[Limit[infHSTInfHSTOdd[q2],m->4]]

(*Therefore, the final expression for the in-plane infinite half-strip-infinite half-strip interaction potential*)
infHSTHSTIP[q2_,hx_,hy_,m_]=Piecewise[{{specm4,m==4},{infHSTInfHSTOdd[q2],m>4}}]

(*Numerical test. Everything is fine. However, there is a singularity for hx->hy*)
r={q2->2.,hx->1,hy->1+1/10000,m->4};
infHSTHSTIP@@r[[All,2]]
NIntegrate[1/((y1-x1)^2+(y2-x2)^2)^(m/2)/.r,{x1,-1000,0}/.r,{x2,-hx,hx}/.r,{y1,q2,1000}/.r,{y2,-hy,hy}/.r]
r={q2->2.,hx->1,hy->1+1/10000,m->5};
infHSTHSTIP@@r[[All,2]]
NIntegrate[1/((y1-x1)^2+(y2-x2)^2)^(m/2)/.r,{x1,-1000,0}/.r,{x2,-hx,hx}/.r,{y1,q2,1000}/.r,{y2,-hy,hy}/.r]
r={q2->2.,hx->1,hy->1+1/10000,m->6};
infHSTHSTIP@@r[[All,2]]
NIntegrate[1/((y1-x1)^2+(y2-x2)^2)^(m/2)/.r,{x1,-1000,0}/.r,{x2,-hx,hx}/.r,{y1,q2,1000}/.r,{y2,-hy,hy}/.r]

(*Yet, the limit for hx->hy is well-defined*)
infHSTHSTIPSpec[q2_,hy_,m_]=Piecewise[{{Limit[infHSTHSTIP[q2,hx,hy,m],hx->hy,Assumptions->m>4],m>4},{Limit[infHSTHSTIP[q2,hx,hy,4],hx->hy],m==4}}]

(*The expressions for the special case hx=hy and integer m are quite simple*)
Table[FullSimplify[infHSTHSTIPSpec[q2,hx,i]],{i,4,8}]

(*Now, we apply RIHBI and obtain R-(Subscript[R, IP]) law for m>=4 as a piecewise function*)
recRecIP[q2_,bx_,by_,hx_,hy_,m_]=infHSTHSTIP[q2,hx,hy,m]+infHSTHSTIP[q2+2bx+2by,hx,hy,m]-infHSTHSTIP[q2+2bx,hx,hy,m]-infHSTHSTIP[q2+2by,hx,hy,m]

recRecIP[q2_,bx_,by_,hx_,hy_,m_]=(({{\[Piecewise], ({{(1/2) ((((-hx+hy+((q2^2)/(hx-hy))) ArcTan[((hx-hy)/q2)])/q2)+((((hx^2)+2 hx hy+(hy^2)-(q2^2)) ArcTan[((hx+hy)/q2)])/((hx+hy) q2))+Log[(((hx^2)-2 hx hy+(hy^2)+(q2^2))/((hx^2)+2 hx hy+(hy^2)+(q2^2)))]), m==4}, {-((2 ((((((hx-hy))^2)+(q2^2)))^(2-(m/2))))/((-4+m) (-3+m) (-2+m)))+((2 ((((((hx+hy))^2)+(q2^2)))^(2-(m/2))))/((-4+m) (-3+m) (-2+m)))-((2 (((hx-hy))^2) (q2^(2-m)) Hypergeometric2F1[(1/2),(1/2) (-2+m),(3/2),-((((hx-hy))^2)/(q2^2))])/(((-2+m))^2))+((2 (((hx+hy))^2) (q2^(2-m)) Hypergeometric2F1[(1/2),(1/2) (-2+m),(3/2),-((((hx+hy))^2)/(q2^2))])/(((-2+m))^2))-((2 (((hx-hy))^2) (q2^(2-m)) Hypergeometric2F1[(1/2) (-2+m),(1/2) (-1+m),((1+m)/2),-((((hx-hy))^2)/(q2^2))])/((-3+m) (((-2+m))^2) (-1+m)))+((2 (((hx+hy))^2) (q2^(2-m)) Hypergeometric2F1[(1/2) (-2+m),(1/2) (-1+m),((1+m)/2),-((((hx+hy))^2)/(q2^2))])/((-3+m) (((-2+m))^2) (-1+m))), m>4}, {0, True}})}}))-(({{\[Piecewise], ({{(1/2) ((((-hx+hy+((((2 bx+q2))^2)/(hx-hy))) ArcTan[((hx-hy)/(2 bx+q2))])/(2 bx+q2))+((((hx^2)+2 hx hy+(hy^2)-(((2 bx+q2))^2)) ArcTan[((hx+hy)/(2 bx+q2))])/((hx+hy) (2 bx+q2)))+Log[(((hx^2)-2 hx hy+(hy^2)+(((2 bx+q2))^2))/((hx^2)+2 hx hy+(hy^2)+(((2 bx+q2))^2)))]), m==4}, {-((2 ((((((hx-hy))^2)+(((2 bx+q2))^2)))^(2-(m/2))))/((-4+m) (-3+m) (-2+m)))+((2 ((((((hx+hy))^2)+(((2 bx+q2))^2)))^(2-(m/2))))/((-4+m) (-3+m) (-2+m)))-((2 (((hx-hy))^2) (((2 bx+q2))^(2-m)) Hypergeometric2F1[(1/2),(1/2) (-2+m),(3/2),-((((hx-hy))^2)/(((2 bx+q2))^2))])/(((-2+m))^2))+((2 (((hx+hy))^2) (((2 bx+q2))^(2-m)) Hypergeometric2F1[(1/2),(1/2) (-2+m),(3/2),-((((hx+hy))^2)/(((2 bx+q2))^2))])/(((-2+m))^2))-((2 (((hx-hy))^2) (((2 bx+q2))^(2-m)) Hypergeometric2F1[(1/2) (-2+m),(1/2) (-1+m),((1+m)/2),-((((hx-hy))^2)/(((2 bx+q2))^2))])/((-3+m) (((-2+m))^2) (-1+m)))+((2 (((hx+hy))^2) (((2 bx+q2))^(2-m)) Hypergeometric2F1[(1/2) (-2+m),(1/2) (-1+m),((1+m)/2),-((((hx+hy))^2)/(((2 bx+q2))^2))])/((-3+m) (((-2+m))^2) (-1+m))), m>4}, {0, True}})}}))-(({{\[Piecewise], ({{(1/2) ((((-hx+hy+((((2 by+q2))^2)/(hx-hy))) ArcTan[((hx-hy)/(2 by+q2))])/(2 by+q2))+((((hx^2)+2 hx hy+(hy^2)-(((2 by+q2))^2)) ArcTan[((hx+hy)/(2 by+q2))])/((hx+hy) (2 by+q2)))+Log[(((hx^2)-2 hx hy+(hy^2)+(((2 by+q2))^2))/((hx^2)+2 hx hy+(hy^2)+(((2 by+q2))^2)))]), m==4}, {-((2 ((((((hx-hy))^2)+(((2 by+q2))^2)))^(2-(m/2))))/((-4+m) (-3+m) (-2+m)))+((2 ((((((hx+hy))^2)+(((2 by+q2))^2)))^(2-(m/2))))/((-4+m) (-3+m) (-2+m)))-((2 (((hx-hy))^2) (((2 by+q2))^(2-m)) Hypergeometric2F1[(1/2),(1/2) (-2+m),(3/2),-((((hx-hy))^2)/(((2 by+q2))^2))])/(((-2+m))^2))+((2 (((hx+hy))^2) (((2 by+q2))^(2-m)) Hypergeometric2F1[(1/2),(1/2) (-2+m),(3/2),-((((hx+hy))^2)/(((2 by+q2))^2))])/(((-2+m))^2))-((2 (((hx-hy))^2) (((2 by+q2))^(2-m)) Hypergeometric2F1[(1/2) (-2+m),(1/2) (-1+m),((1+m)/2),-((((hx-hy))^2)/(((2 by+q2))^2))])/((-3+m) (((-2+m))^2) (-1+m)))+((2 (((hx+hy))^2) (((2 by+q2))^(2-m)) Hypergeometric2F1[(1/2) (-2+m),(1/2) (-1+m),((1+m)/2),-((((hx+hy))^2)/(((2 by+q2))^2))])/((-3+m) (((-2+m))^2) (-1+m))), m>4}, {0, True}})}}))+(({{\[Piecewise], ({{(1/2) ((((-hx+hy+((((2 bx+2 by+q2))^2)/(hx-hy))) ArcTan[((hx-hy)/(2 bx+2 by+q2))])/(2 bx+2 by+q2))+((((hx^2)+2 hx hy+(hy^2)-(((2 bx+2 by+q2))^2)) ArcTan[((hx+hy)/(2 bx+2 by+q2))])/((hx+hy) (2 bx+2 by+q2)))+Log[(((hx^2)-2 hx hy+(hy^2)+(((2 bx+2 by+q2))^2))/((hx^2)+2 hx hy+(hy^2)+(((2 bx+2 by+q2))^2)))]), m==4}, {-((2 ((((((hx-hy))^2)+(((2 bx+2 by+q2))^2)))^(2-(m/2))))/((-4+m) (-3+m) (-2+m)))+((2 ((((((hx+hy))^2)+(((2 bx+2 by+q2))^2)))^(2-(m/2))))/((-4+m) (-3+m) (-2+m)))-((2 (((hx-hy))^2) (((2 bx+2 by+q2))^(2-m)) Hypergeometric2F1[(1/2),(1/2) (-2+m),(3/2),-((((hx-hy))^2)/(((2 bx+2 by+q2))^2))])/(((-2+m))^2))+((2 (((hx+hy))^2) (((2 bx+2 by+q2))^(2-m)) Hypergeometric2F1[(1/2),(1/2) (-2+m),(3/2),-((((hx+hy))^2)/(((2 bx+2 by+q2))^2))])/(((-2+m))^2))-((2 (((hx-hy))^2) (((2 bx+2 by+q2))^(2-m)) Hypergeometric2F1[(1/2) (-2+m),(1/2) (-1+m),((1+m)/2),-((((hx-hy))^2)/(((2 bx+2 by+q2))^2))])/((-3+m) (((-2+m))^2) (-1+m)))+((2 (((hx+hy))^2) (((2 bx+2 by+q2))^(2-m)) Hypergeometric2F1[(1/2) (-2+m),(1/2) (-1+m),((1+m)/2),-((((hx+hy))^2)/(((2 bx+2 by+q2))^2))])/((-3+m) (((-2+m))^2) (-1+m))), m>4}, {0, True}})}}));

(*Numerical test. Everything is fine.*)
r={q2->1/1.,bx->1,by->2,hx->1,hy->2,m->6};
recRecIP@@r[[All,2]]
NIntegrate[1/((y1-x1)^2+(y2-x2)^2)^(m/2)/.r,{x1,-bx,bx}/.r,{x2,-hx,hx}/.r,{y1,bx+q2,bx+q2+2 by}/.r,{y2,-hy,hy}/.r]
r={q2->1/1.,bx->1,by->2,hx->1,hy->2,m->12};
recRecIP@@r[[All,2]]
NIntegrate[1/((y1-x1)^2+(y2-x2)^2)^(m/2)/.r,{x1,-bx,bx}/.r,{x2,-hx,hx}/.r,{y1,bx+q2,bx+q2+2 by}/.r,{y2,-hy,hy}/.r]

(*The expression is valid for arbitrary m>=4. We could make it work also for m<4, but due to the asymptotic scaling, this definition is well-suited.*)
r={q2->1/10.,bx->1,by->2,hx->1,hy->2,m->4.5};
recRecIP@@r[[All,2]]
NIntegrate[1/((y1-x1)^2+(y2-x2)^2)^(m/2)/.r,{x1,-bx,bx}/.r,{x2,-hx,hx}/.r,{y1,bx+q2,bx+q2+2 by}/.r,{y2,-hy,hy}/.r]

(*scaling*)
tab1={};
tab2={};
For[i=4,i<13,i++,

f = -D[Log[recRecIP[q2,bx,by,hx,hy,i]/.{q2->Exp[q2]}],q2]/.{hx->1,hy->1/2,bx->1/2,by->2/1};
AppendTo[tab1,Re[f]];
AppendTo[tab2,i]];

Plot[tab1,{q2,-10,10},AxesLabel->{"Log[q2]","Scaling factor of (Subscript[Π, R - (Subscript[R, IP])])"},PlotLegends->tab2,WorkingPrecision->100]

(*We can compare our result for the vdW attraction, m=6, with the result of De Rocco and Hoover 1960. They consider rectangles of equal size and m=6 exclusively. *)

recRecIP6spec=FullSimplify[FunctionExpand[infHRHRIPSpec[q2,hy,6]+infHRHRIPSpec[q2+2 by+2 by,hy,6]-infHRHRIPSpec[q2+2 by,hy,6]-infHRHRIPSpec[q2+2 by,hy,6]]]

(*The final obtained expression is the same as in De Rocco and Hoover 1960, although in a different form. Let us emphasize that the expression derived here is valid for arbitrary dimensions of rectangles*)
recRecIP6deRocco=(1/48) (-(1/(q2^2))+(2/(((a+q2))^2))-(1/(((2 a+q2))^2))+((3 q2 ArcCot[(b/q2)])/(b^3))-((6 (a+q2) ArcCot[(b/(a+q2))])/(b^3))+((3 (2 a+q2) ArcCot[(b/(2 a+q2))])/(b^3))+((3 b ArcTan[(b/q2)])/(q2^3))-((6 b ArcTan[(b/(a+q2))])/(((a+q2))^3))+((3 b ArcTan[(b/(2 a+q2))])/(((2 a+q2))^3)));
FullSimplify[recRecIP6deRocco-recRecIP6spec/.{bx->(a/2),by->(a/2),hy->b/2}]

(*Now, it is straightforward to apply the RIIPI, and obtain R-RP laws for m>=5*)
recRP[q2_,bx_,by_,hx_,hy_,m_]=(Sqrt[π])*((Gamma[((m-1)/2)])/(Gamma[(m/2)]))*recRecIP[q2,bx,by,hx,hy,m-1];

(*Next, let us check the scaling of the R-RP laws. For small separations, the asymptotic scaling is m-4, and for large it is m-1, because it approaches point-line interaction.*)
tab1={};
tab2={};
For[i=5,i<13,i++,

f = -D[Log[recRP[q2,bx,by,hx,hy,i]/.{q2->Exp[q2]}],q2]/.{hx->1,hy->1/2,bx->1/2,by->2/1};
AppendTo[tab1,Re[f]];
AppendTo[tab2,i]];

Plot[tab1,{q2,-10,10},AxesLabel->{"Log[q2]","Scaling factor of (Subscript[Π, R - RP])"},PlotLegends->tab2,WorkingPrecision->100]

(*If we want to find a simplified version of expressions, it can be simpler to go back to the infHSTHSTIPSpec[q2,hx,hy,m] function. For example, let us find R-RP for m=6. First we simplify HST-(Subscript[HST, IP]) for m=5, then calculate R-(Subscript[R, IP]), multiply it by f6, and simplify again.*)
simHSTHST5[q2_]= FullSimplify[Collect[FullSimplify[infHSTHSTIP[q2,hx,hy,5]],{(Sqrt[(((hx-hy))^2)+(q2^2)]),(Sqrt[(((hx+hy))^2)+(q2^2)])}]]
RRP6=FullSimplify[Expand[(Sqrt[π])*((Gamma[((6-1)/2)])/(Gamma[(6/2)]))(simHSTHST5[q2]+simHSTHST5[q2+2bx+2by]-simHSTHST5[q2+2bx]-simHSTHST5[q2+2by])],TimeConstraint->1]
simRRP6=FullSimplify[Collect[RRP6,{(Sqrt[(((hx-hy))^2)+((_)^2)]),(Sqrt[(((hx+hy))^2)+((_)^2)])}]]

(*The special case, hx=hy*)
simHSTHST5spec[q2_]=FullSimplify[Collect[FullSimplify[infHSTHSTIPSpec[q2,hy,5]],{(Sqrt[(hy^2)+(q2^2)]),(Sqrt[4 (hy^2)+(q2^2)])}]]
RRP6spec=FullSimplify[Expand[(Sqrt[π])*((Gamma[((6-1)/2)])/(Gamma[(6/2)]))(simHSTHST5spec[q2]+simHSTHST5spec[q2+2bx+2by]-simHSTHST5spec[q2+2bx]-simHSTHST5spec[q2+2by])],TimeConstraint->1]
simRRP6spec=FullSimplify[Collect[RRP6spec,{(Sqrt[(hy^2)+((_)^2)]),(Sqrt[4 (hy^2)+((_)^2)])}]]

(*For the paper, we can first find the infinite half-strip-infinite half-plate law, and then apply RIHBI*)
infHSTinfHP[q2_]=FullSimplify[Collect[FullSimplify[((Sqrt[π])*((Gamma[((m-1)/2)])/(Gamma[(m/2)]))/.m->6)*infHSTHSTIPSpec[q2,hy,5]],{(Sqrt[4 (hy^2)+(q2^2)])}]]

(*The term linear w.r.t. the gap q2 cancels out in the RIHBI expression and the final representation is very simple*)

ξ[a_]=((π ((((4 (hy^2)+(a^2)))^(3/2))-6 (hy^2) a))/(48 (hy^2) (a^2)));
recRPSpec6=ξ[q2]+ξ[q2+2bx+2by]-ξ[q2+2bx]-ξ[q2+2by];
(*The result is the same as in the above derivation*)
FullSimplify[simRRP6spec-recRPSpec6]

(*Let us consider the scaling factor of the R-RP interaction for m=6 when we vary by*)
tab1={};
tab2=Table[(10^k),{k,-3,3,1}];
For[i=1,i<Length[tab2]+1,i++,

f = -D[Log[recRP[q2,bx,by,hx,hy,6]/.{q2->Exp[q2]}],q2]/.{hx->1,hy->1/2,bx->1,by->tab2[[i]]};
AppendTo[tab1,Re[f]]];

Plot[tab1,{q2,-10,14},AxesLabel->{"Log[q2]","Scaling factor of (Subscript[Π, R - RP]) for m=6 and different hy/hx ratios"},PlotLegends->tab2,WorkingPrecision->100]