(*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`*"];

(*set global assumptions*)
SetSystemOptions["SimplificationOptions"->{"AssumptionsMaxNonlinearVariables"->25,"AssumptionsMaxVariables"->25}];
$Assumptions=Rx>0&&Ry>0&&t+Rx>=p≥t-Rx&&Rx+2*Ry+q2>=t>=Rx+q2&&q2>0&&Element[q1,Reals]&&d>0&&m≥3(*&&Element[m,Integers]*)&&n≥0&&Element[n,Integers]&&k≥0&&Element[k,Integers]&&t>0&&rx>0&&ry>0&&π>fx>0&&π>fy>0;

(*A function that helps integration. We restrict the time allowed for the simplification to speed-up the process at the expense of additional transformations.*)
integrateD[x_,v_,l_,u_]:=Module[{int=x,var=v,low=l,up=u},
indef=FullSimplify[Integrate[int,var],TimeConstraint->2];
upper=Limit[indef,var->up,Direction->"FromBelow"];
lower=Limit[indef,var->low,Direction->"FromAbove"];
res=FullSimplify[upper-lower,TimeConstraint->2];
Return[res]
]

(*Coordinate transformation used in Langbein 1972. Two relative polar coordinate systems, RPCS1*)

x2[p_,t_,φ_,ψ_]=t*Cos[ψ];
y2[p_,t_,φ_,ψ_]=t*Sin[ψ];

x1[p_,t_,φ_,ψ_]=x2[p,t,φ,ψ]-p*Cos[-φ+ψ];
y1[p_,t_,φ_,ψ_]=y2[p,t,φ,ψ]-p*Sin[-φ+ψ];

Dg=Grad[{x1[p,t,φ,ψ],y1[p,t,φ,ψ],x2[p,t,φ,ψ],y2[p,t,φ,ψ]},{p,t,φ,ψ}];

(*The Jacobian of the coordinate transformation*)
jac=PowerExpand[FullSimplify[Det[Dg]]]

(*Two polar coordinate systems, for the reference*)

x2[p_,t_,φ_,ψ_]=t*Cos[ψ];
y2[p_,t_,φ_,ψ_]=t*Sin[ψ];

x1[p_,t_,φ_,ψ_]=Rx+q2+Ry-p*Cos[φ];
y1[p_,t_,φ_,ψ_]=p*Sin[φ];

Dg=Grad[{x1[p,t,φ,ψ],y1[p,t,φ,ψ],x2[p,t,φ,ψ],y2[p,t,φ,ψ]},{p,t,φ,ψ}];

(*The Jacobian of the coordinate transformation*)
jac=PowerExpand[FullSimplify[Det[Dg]]]

(*Coordinate transformation used in Langbein 1972. Two relative polar coordinate systems, RPCS2. q2 is the gap*)

dd=Rx+Ry+q2;
tt=(Sqrt[(rx^2)+(d^2)-2*rx*d*Cos[fx]])/.d->dd;
cosa=(((d^2)+(t^2)-(rx^2))/(2*t*d))/.t->tt/.d->dd;
sina=(rx/t)*Sin[fx]/.t->tt;
x2[rx_,ry_,fx_,fy_]=d+ry*Cos[fy]*cosa+ry*Sin[fy]*sina/.d->dd;
y2[rx_,ry_,fx_,fy_]=ry*Cos[fy]*sina-ry*Sin[fy]*cosa;

x1[rx_,ry_,fx_,fy_]=rx*Cos[fx];
y1[rx_,ry_,fx_,fy_]=rx*Sin[fx];

Dg=Grad[{x1[rx,ry,fx,fy],y1[rx,ry,fx,fy],x2[rx,ry,fx,fy],y2[rx,ry,fx,fy]},{rx,ry,fx,fy}];

(*The Jacobian of the coordinate transformation*)
jac=PowerExpand[FullSimplify[Det[Dg]]]

(*plotting two in-plane disks to check the coordinate transformations*)
rep={Rx->0.8,Ry->1.2,q2->0.1};

d1=8;
d2=8;
aPlot=Flatten[Table[{x1[rx,0,fx,0],y1[rx,0,fx,0]},{rx,Table[y,{y,0,Rx,Rx/d1}]/.rep},{fx,Table[i,{i,0,2*π,2*π/d1}]}],1];
bPlot=Flatten[Table[{x2[rx,ry,fx,fy],y2[rx,ry,fx,fy]}/.rep,{fy,Table[y,{y,0,π*2,π*2/d2}]},{ry,Table[y,{y,0,Ry,Ry/d2}]/.rep},{rx,Table[y,{y,0,Rx,Rx/d1}]/.rep},{fx,Table[i,{i,0,2*π,2*π/d1}]}],3];

ListPlot[{aPlot,bPlot},PlotStyle->{PointSize[0.01]},PlotRange->All]

(*Using RPCS1, the in-plane disk-disk potential of (m^th) order is obtained as Integrate[4*ArcCos[(((t^2)+(((Rx+Ry+q2))^2)-(Ry^2))/(2*t*(Rx+Ry+q2)))]ArcCos[(((p^2)+(t^2)-(Rx^2))/(2*t*p))]*(p^(-m)) p t,{p,t-Rx,t+Rx},{t,q2+Rx,q2+Rx+2*Ry}]. To simplify, we are first integrating over the p coordinate, and the result is an in-plane point-disk law. Therefore, we pursue integration w.r.t. p for m∈[3,12]. We apply limits on additively decomposed primitive functions, since it turns out to be the most efficient approach*)
$Assumptions=t>0&&t>Rx&&Rx>0&&p>0;
indef=FullSimplify[Table[Integrate[2*ArcCos[(p^2+t^2-Rx^2)/(2*t*p)]*p^(-m+1),p],{m,3,12}],TimeConstraint->1];
lev=Level[#,1]&/@Apart/@indef;
limL=Limit[#,p->t-Rx,Direction->"FromAbove"]&/@lev;
limU=Limit[#,p->t+Rx,Direction->"FromBelow"]&/@lev;
tabIntP=FullSimplify[Total/@(limU-limL)];

(*This is the result, P-(Subscript[D, IP]) laws for m=3,4,5,...12*)
tabIntP={((2 EllipticE[((4 Rx t)/(((Rx+t))^2))])/(-Rx+t))-((2 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))]+2 (((Rx-t))^2) EllipticK[((4 Rx t)/(((Rx+t))^2))])/(9 (((Rx-t))^3) (((Rx+t))^2))),((π ((Rx^4)+2 (Rx^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))]+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))),((-2 (337 (Rx^6)+1773 (Rx^4) (t^2)+947 (Rx^2) (t^4)+15 (t^6)) EllipticE[((4 Rx t)/(((Rx+t))^2))]+2 (((Rx-t))^2) (127 (Rx^4)+242 (Rx^2) (t^2)+15 (t^4)) EllipticK[((4 Rx t)/(((Rx+t))^2))])/(735 (((Rx-t))^7) (((Rx+t))^6))),((π (Rx^2) ((Rx^6)+12 (Rx^4) (t^2)+18 (Rx^2) (t^4)+4 (t^6)))/(4 ((((Rx^2)-(t^2)))^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))]+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))])/(2835 (((Rx-t))^9) (((Rx+t))^8))),((π (Rx^2) ((Rx^8)+20 (Rx^6) (t^2)+60 (Rx^4) (t^4)+40 (Rx^2) (t^6)+5 (t^8)))/(5 ((((Rx^2)-(t^2)))^10)))};

(*Let us check the scaling of the P-(Subscript[D, IP]) laws of order m. For small separations, the asymptotic scaling goes with m-2, and for large as m, because then it approaches the point-point interaction.*)
tab1={};
tab2={};
For[i=1,i<Length[tabIntP]+1,i++,

f = D[Log[tabIntP[[i]]/.{t->Rx+q2}/.{q2->Exp[q2]}],q2]/.{Rx->1};
AppendTo[tab1,f];
AppendTo[tab2,i+2]];

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

(*Now, we have the P-(Subscript[D, IP]) law. If we multiply it by a factor fm,it becomes a P-C law of (((m+1))^th) order (RIIPI). By integrating P-(Subscript[D, IP]) w.r.t. the other disk, we obtain D-(Subscript[D, IP]) laws. It turns out that the integration is not possible for functions involving elliptic integrals, that is for P-(Subscript[D, IP]) with odd m. We will pursue an infinite series representation later. For now, we only integrate P-(Subscript[D, IP]) laws for even m,and obtain D-(Subscript[D, IP]) laws. *)

tabIntPT=Table[integrateD[2 ArcCos[(((t^2)+(((Rx+Ry+q2))^2)-(Ry^2))/(2*t*(Rx+Ry+q2)))] tabIntP[[i]] t,t,Rx+q2,Rx+q2+2*Ry],{i,2,Length[tabIntP],2}]

(*The result (with some additional simplification): D-(Subscript[D, IP]) laws for m=4,6,8,10,12*)
tabIntPT=({{(((π^2) ((q2^2)+2 Rx Ry+2 q2 (Rx+Ry)-(Sqrt[q2 (q2+2 Rx) (q2+2 Ry) (q2+2 (Rx+Ry))])))/(2 (Sqrt[q2 (q2+2 Rx) (q2+2 Ry) (q2+2 (Rx+Ry))])))}, {(((π^2) (Rx^2) (Ry^2) (2 (q2^4)+8 (q2^3) (Rx+Ry)+6 Rx Ry (((Rx+Ry))^2)+2 q2 (Rx+Ry) (3 (Rx^2)+8 Rx Ry+3 (Ry^2))+(q2^2) (11 (Rx^2)+24 Rx Ry+11 (Ry^2))))/(2 (((q2 (q2+2 Rx) (q2+2 Ry) (q2+2 (Rx+Ry))))^(5/2))))}, {(((π^2) (Rx^2) (Ry^2) (3 (q2^10)+30 (q2^9) (Rx+Ry)+140 (Rx^3) (Ry^3) (((Rx+Ry))^4)+60 q2 (Rx^2) (Ry^2) (((Rx+Ry))^3) (5 (Rx^2)+14 Rx Ry+5 (Ry^2))+48 (q2^7) (Rx+Ry) (7 (Rx^2)+15 Rx Ry+7 (Ry^2))+12 (q2^5) (Rx+Ry) (5 (Rx^2)+11 Rx Ry+5 (Ry^2)) (9 (Rx^2)+25 Rx Ry+9 (Ry^2))+6 (q2^8) (22 (Rx^2)+45 Rx Ry+22 (Ry^2))+6 (q2^2) Rx Ry (((Rx+Ry))^2) (40 (Rx^4)+273 (Rx^3) Ry+494 (Rx^2) (Ry^2)+273 Rx (Ry^3)+40 (Ry^4))+(q2^6) (538 (Rx^4)+2352 (Rx^3) Ry+3642 (Rx^2) (Ry^2)+2352 Rx (Ry^3)+538 (Ry^4))+16 (q2^3) (Rx+Ry) (5 (Rx^6)+73 (Rx^5) Ry+287 (Rx^4) (Ry^2)+441 (Rx^3) (Ry^3)+287 (Rx^2) (Ry^4)+73 Rx (Ry^5)+5 (Ry^6))+12 (q2^4) (26 (Rx^6)+225 (Rx^5) Ry+691 (Rx^4) (Ry^2)+985 (Rx^3) (Ry^3)+691 (Rx^2) (Ry^4)+225 Rx (Ry^5)+26 (Ry^6))))/(3 (((q2 (q2+2 Rx) (q2+2 Ry) (q2+2 (Rx+Ry))))^(9/2))))}, {(((π^2) (Rx^2) (Ry^2) (4 (q2^16)+64 (q2^15) (Rx+Ry)+4620 (Rx^5) (Ry^5) (((Rx+Ry))^6)+2100 q2 (Rx^4) (Ry^4) (((Rx+Ry))^5) (7 (Rx^2)+20 Rx Ry+7 (Ry^2))+28 (q2^13) (Rx+Ry) (79 (Rx^2)+160 Rx Ry+79 (Ry^2))+(q2^14) (478 (Rx^2)+960 Rx Ry+478 (Ry^2))+70 (q2^2) (Rx^3) (Ry^3) (((Rx+Ry))^4) (280 (Rx^4)+1801 (Rx^3) Ry+3252 (Rx^2) (Ry^2)+1801 Rx (Ry^3)+280 (Ry^4))+16 (q2^11) (Rx+Ry) (1018 (Rx^4)+4277 (Rx^3) Ry+6551 (Rx^2) (Ry^2)+4277 Rx (Ry^3)+1018 (Ry^4))+4 (q2^12) (1765 (Rx^4)+7189 (Rx^3) Ry+10859 (Rx^2) (Ry^2)+7189 Rx (Ry^3)+1765 (Ry^4))+80 (q2^3) (Rx^2) (Ry^2) (((Rx+Ry))^3) (175 (Rx^6)+1973 (Rx^5) Ry+7160 (Rx^4) (Ry^2)+10839 (Rx^3) (Ry^3)+7160 (Rx^2) (Ry^4)+1973 Rx (Ry^5)+175 (Ry^6))+2 (q2^9) (Rx+Ry) (17103 (Rx^6)+120912 (Rx^5) Ry+337516 (Rx^4) (Ry^2)+467984 (Rx^3) (Ry^3)+337516 (Rx^2) (Ry^4)+120912 Rx (Ry^5)+17103 (Ry^6))+(q2^10) (27603 (Rx^6)+179168 (Rx^5) Ry+471740 (Rx^4) (Ry^2)+640464 (Rx^3) (Ry^3)+471740 (Rx^2) (Ry^4)+179168 Rx (Ry^5)+27603 (Ry^6))+16 (q2^7) (Rx+Ry) (1135 (Rx^8)+13997 (Rx^7) Ry+67315 (Rx^6) (Ry^2)+163172 (Rx^5) (Ry^3)+217520 (Rx^4) (Ry^4)+163172 (Rx^3) (Ry^5)+67315 (Rx^2) (Ry^6)+13997 Rx (Ry^7)+1135 (Ry^8))+4 (q2^4) Rx Ry (((Rx+Ry))^2) (1400 (Rx^8)+26290 (Rx^7) Ry+160545 (Rx^6) (Ry^2)+441316 (Rx^5) (Ry^3)+612147 (Rx^4) (Ry^4)+441316 (Rx^3) (Ry^5)+160545 (Rx^2) (Ry^6)+26290 Rx (Ry^7)+1400 (Ry^8))+2 (q2^8) (15132 (Rx^8)+153927 (Rx^7) Ry+644532 (Rx^6) (Ry^2)+1438332 (Rx^5) (Ry^3)+1865272 (Rx^4) (Ry^4)+1438332 (Rx^3) (Ry^5)+644532 (Rx^2) (Ry^6)+153927 Rx (Ry^7)+15132 (Ry^8))+4 (q2^5) (Rx+Ry) (280 (Rx^10)+9680 (Rx^9) Ry+96785 (Rx^8) (Ry^2)+430672 (Rx^7) (Ry^3)+1005900 (Rx^6) (Ry^4)+1325376 (Rx^5) (Ry^5)+1005900 (Rx^4) (Ry^6)+430672 (Rx^3) (Ry^7)+96785 (Rx^2) (Ry^8)+9680 Rx (Ry^9)+280 (Ry^10))+2 (q2^6) (3320 (Rx^10)+63560 (Rx^9) Ry+458283 (Rx^8) (Ry^2)+1680168 (Rx^7) (Ry^3)+3533748 (Rx^6) (Ry^4)+4503296 (Rx^5) (Ry^5)+3533748 (Rx^4) (Ry^6)+1680168 (Rx^3) (Ry^7)+458283 (Rx^2) (Ry^8)+63560 Rx (Ry^9)+3320 (Ry^10))))/(4 (((q2 (q2+2 Rx) (q2+2 Ry) (q2+2 (Rx+Ry))))^(13/2))))}, {(((π^2) (Rx^2) (Ry^2) (5 (q2^22)+110 (q2^21) (Rx+Ry)+180180 (Rx^7) (Ry^7) (((Rx+Ry))^8)+84084 q2 (Rx^6) (Ry^6) (((Rx+Ry))^7) (9 (Rx^2)+26 Rx Ry+9 (Ry^2))+200 (q2^19) (Rx+Ry) (39 (Rx^2)+77 Rx Ry+39 (Ry^2))+10 (q2^20) (116 (Rx^2)+231 Rx Ry+116 (Ry^2))+15 (q2^18) (2494 (Rx^4)+9880 (Rx^3) Ry+14779 (Rx^2) (Ry^2)+9880 Rx (Ry^3)+2494 (Ry^4))+462 (q2^2) (Rx^5) (Ry^5) (((Rx+Ry))^6) (3024 (Rx^4)+19043 (Rx^3) Ry+34446 (Rx^2) (Ry^2)+19043 Rx (Ry^3)+3024 (Ry^4))+30 (q2^17) (Rx+Ry) (4510 (Rx^4)+17936 (Rx^3) Ry+26915 (Rx^2) (Ry^2)+17936 Rx (Ry^3)+4510 (Ry^4))+672 (q2^3) (Rx^4) (Ry^4) (((Rx+Ry))^5) (2205 (Rx^6)+23109 (Rx^5) Ry+81679 (Rx^4) (Ry^2)+123125 (Rx^3) (Ry^3)+81679 (Rx^2) (Ry^4)+23109 Rx (Ry^5)+2205 (Ry^6))+160 (q2^15) (Rx+Ry) (5254 (Rx^6)+32742 (Rx^5) Ry+84411 (Rx^4) (Ry^2)+113883 (Rx^3) (Ry^3)+84411 (Rx^2) (Ry^4)+32742 Rx (Ry^5)+5254 (Ry^6))+10 (q2^16) (37996 (Rx^6)+230010 (Rx^5) Ry+579876 (Rx^4) (Ry^2)+775761 (Rx^3) (Ry^3)+579876 (Rx^2) (Ry^4)+230010 Rx (Ry^5)+37996 (Ry^6))+56 (q2^4) (Rx^3) (Ry^3) (((Rx+Ry))^4) (17640 (Rx^8)+280490 (Rx^7) Ry+1585497 (Rx^6) (Ry^2)+4221567 (Rx^5) (Ry^3)+5809145 (Rx^4) (Ry^4)+4221567 (Rx^3) (Ry^5)+1585497 (Rx^2) (Ry^6)+280490 Rx (Ry^7)+17640 (Ry^8))+4 (q2^13) (Rx+Ry) (509081 (Rx^8)+4644480 (Rx^7) Ry+18166580 (Rx^6) (Ry^2)+39097760 (Rx^5) (Ry^3)+50136910 (Rx^4) (Ry^4)+39097760 (Rx^3) (Ry^5)+18166580 (Rx^2) (Ry^6)+4644480 Rx (Ry^7)+509081 (Ry^8))+2 (q2^14) (736223 (Rx^8)+6304800 (Rx^7) Ry+23379020 (Rx^6) (Ry^2)+48579600 (Rx^5) (Ry^3)+61538850 (Rx^4) (Ry^4)+48579600 (Rx^3) (Ry^5)+23379020 (Rx^2) (Ry^6)+6304800 Rx (Ry^7)+736223 (Ry^8))+40 (q2^5) (Rx^2) (Ry^2) (((Rx+Ry))^3) (10584 (Rx^10)+248864 (Rx^9) Ry+2090319 (Rx^8) (Ry^2)+8552402 (Rx^7) (Ry^3)+19236866 (Rx^6) (Ry^4)+25079082 (Rx^5) (Ry^5)+19236866 (Rx^4) (Ry^6)+8552402 (Rx^3) (Ry^7)+2090319 (Rx^2) (Ry^8)+248864 Rx (Ry^9)+10584 (Ry^10))+16 (q2^11) (Rx+Ry) (113885 (Rx^10)+1536613 (Rx^9) Ry+8943355 (Rx^8) (Ry^2)+28940340 (Rx^7) (Ry^3)+57056850 (Rx^6) (Ry^4)+71275750 (Rx^5) (Ry^5)+57056850 (Rx^4) (Ry^6)+28940340 (Rx^3) (Ry^7)+8943355 (Rx^2) (Ry^8)+1536613 Rx (Ry^9)+113885 (Ry^10))+4 (q2^12) (550166 (Rx^10)+6618053 (Rx^9) Ry+34959996 (Rx^8) (Ry^2)+105298180 (Rx^7) (Ry^3)+198719380 (Rx^6) (Ry^4)+244627110 (Rx^5) (Ry^5)+198719380 (Rx^4) (Ry^6)+105298180 (Rx^3) (Ry^7)+34959996 (Rx^2) (Ry^8)+6618053 Rx (Ry^9)+550166 (Ry^10))+4 (q2^6) Rx Ry (((Rx+Ry))^2) (28224 (Rx^12)+1009176 (Rx^11) Ry+12380984 (Rx^10) (Ry^2)+73703881 (Rx^9) (Ry^3)+245974094 (Rx^8) (Ry^4)+495036590 (Rx^7) (Ry^5)+622747390 (Rx^6) (Ry^6)+495036590 (Rx^5) (Ry^7)+245974094 (Rx^4) (Ry^8)+73703881 (Rx^3) (Ry^9)+12380984 (Rx^2) (Ry^10)+1009176 Rx (Ry^11)+28224 (Ry^12))+10 (q2^10) (111760 (Rx^12)+2004376 (Rx^11) Ry+15549259 (Rx^10) (Ry^2)+67969704 (Rx^9) (Ry^3)+185623382 (Rx^8) (Ry^4)+333487088 (Rx^7) (Ry^5)+404353743 (Rx^6) (Ry^6)+333487088 (Rx^5) (Ry^7)+185623382 (Rx^4) (Ry^8)+67969704 (Rx^3) (Ry^9)+15549259 (Rx^2) (Ry^10)+2004376 Rx (Ry^11)+111760 (Ry^12))+4 (q2^9) (Rx+Ry) (119728 (Rx^12)+2674272 (Rx^11) Ry+24616083 (Rx^10) (Ry^2)+122042352 (Rx^9) (Ry^3)+363270198 (Rx^8) (Ry^4)+686263680 (Rx^7) (Ry^5)+845953815 (Rx^6) (Ry^6)+686263680 (Rx^5) (Ry^7)+363270198 (Rx^4) (Ry^8)+122042352 (Rx^3) (Ry^9)+24616083 (Rx^2) (Ry^10)+2674272 Rx (Ry^11)+119728 (Ry^12))+32 (q2^7) (Rx+Ry) (504 (Rx^14)+31528 (Rx^13) Ry+585291 (Rx^12) (Ry^2)+5029009 (Rx^11) (Ry^3)+23956541 (Rx^10) (Ry^4)+69493498 (Rx^9) (Ry^5)+129331219 (Rx^8) (Ry^6)+158639285 (Rx^7) (Ry^7)+129331219 (Rx^6) (Ry^8)+69493498 (Rx^5) (Ry^9)+23956541 (Rx^4) (Ry^10)+5029009 (Rx^3) (Ry^11)+585291 (Rx^2) (Ry^12)+31528 Rx (Ry^13)+504 (Ry^14))+4 (q2^8) (32032 (Rx^14)+1077552 (Rx^13) Ry+13996300 (Rx^12) (Ry^2)+95284995 (Rx^11) (Ry^3)+388181520 (Rx^10) (Ry^4)+1017251142 (Rx^9) (Ry^5)+1787285202 (Rx^8) (Ry^6)+2151762975 (Rx^7) (Ry^7)+1787285202 (Rx^6) (Ry^8)+1017251142 (Rx^5) (Ry^9)+388181520 (Rx^4) (Ry^10)+95284995 (Rx^3) (Ry^11)+13996300 (Rx^2) (Ry^12)+1077552 Rx (Ry^13)+32032 (Ry^14))))/(5 (((q2 (q2+2 Rx) (q2+2 Ry) (q2+2 (Rx+Ry))))^(17/2))))}});

(*Let us check the scaling. For small separations, the asymptotic scaling is m-7/2. So, although we did the integration over the area of the second disk, the asymptotic scaling factor has reduced only for a factor of 3/2. For large separations, the asymptotic scaling factor is m (resembles the P-P interaction).*)

tab1={};
tab2={};
For[i=1,i<6,i++,

f = D[Log[tabIntPT[[i,1]]/.{q2->Exp[q2]}],q2]/.{q2->q2f,Rx->1,Ry->1};
AppendTo[tab1,f];
AppendTo[tab2,2*i+2]];

Plot[tab1,{q2f,-8,8},AxesLabel->{"gap Log[q2]","Scaling factor"},PlotLegends->tab2,WorkingPrecision->50]

(*For example, D-(Subscript[D, IP]) for m=6 can be found in Dube 1939. The expressions are the same, just in a different form*)
dube=(((π^2) (a^2) (b^2))/2)*(((c^2) (2 (c^2)-(a^2)-(b^2))-((((b^2)-(a^2)))^2))/(((((((c+b))^2)-(a^2)))^(5/2))*((((((c-b))^2)-(a^2)))^(5/2))))/.{a->Rx,b->Ry,c->Rx+Ry+q2}

(*some additional replacements and simplification*)
tabIntPTdFS=FullSimplify[PowerExpand[tabIntPT]/. q2+2 (Rx+Ry)->d+Rx+Ry/. q2+2 Ry->d-Rx+Ry/. q2+2 Rx->d-Ry+Rx/. q2->d-Rx-Ry]

(*The result after the simplification*)
tabIntPTdFS=({{-(((π^2) (-(d^2)+(Rx^2)+(Ry^2)+(Sqrt[d-Rx-Ry]) (Sqrt[(d+Rx-Ry) (d-Rx+Ry) (d+Rx+Ry)])))/(2 (Sqrt[d-Rx-Ry]) (Sqrt[(d+Rx-Ry) (d-Rx+Ry) (d+Rx+Ry)])))}, {-(((π^2) (Rx^2) (Ry^2) (-2 (d^4)+((((Rx^2)-(Ry^2)))^2)+(d^2) ((Rx^2)+(Ry^2))))/(2 (((d-Rx-Ry))^(5/2)) ((((d+Rx-Ry) (d-Rx+Ry) (d+Rx+Ry)))^(5/2))))}, {-(((π^2) (Rx^2) (Ry^2) (-3 (d^10)+3 (d^8) ((Rx^2)+(Ry^2))+((((Rx^2)-(Ry^2)))^4) ((Rx^2)+(Ry^2))+3 (d^2) ((((Rx^2)-(Ry^2)))^2) ((Rx^4)+6 (Rx^2) (Ry^2)+(Ry^4))+(d^6) (8 (Rx^4)-30 (Rx^2) (Ry^2)+8 (Ry^4))-6 (d^4) (2 (Rx^6)-3 (Rx^4) (Ry^2)-3 (Rx^2) (Ry^4)+2 (Ry^6))))/(3 (((d-Rx-Ry))^(9/2)) ((((d+Rx-Ry) (d-Rx+Ry) (d+Rx+Ry)))^(9/2))))}, {(((π^2) (Rx^2) (Ry^2) (4 (d^16)-2 (d^14) ((Rx^2)+(Ry^2))-((((Rx^2)-(Ry^2)))^6) ((Rx^4)+3 (Rx^2) (Ry^2)+(Ry^4))-(d^2) (((Rx-Ry))^4) (((Rx+Ry))^4) ((Rx^2)+(Ry^2)) (7 (Rx^4)+61 (Rx^2) (Ry^2)+7 (Ry^4))-2 (d^12) (19 (Rx^4)-60 (Rx^2) (Ry^2)+19 (Ry^4))+3 (d^10) ((Rx^2)+(Ry^2)) (27 (Rx^4)-73 (Rx^2) (Ry^2)+27 (Ry^4))-2 (d^6) ((Rx^2)+(Ry^2)) (12 (Rx^8)-222 (Rx^6) (Ry^2)+445 (Rx^4) (Ry^4)-222 (Rx^2) (Ry^6)+12 (Ry^8))+2 (d^4) ((((Rx^2)-(Ry^2)))^2) (16 (Rx^8)-22 (Rx^6) (Ry^2)-213 (Rx^4) (Ry^4)-22 (Rx^2) (Ry^6)+16 (Ry^8))-(d^8) (45 (Rx^8)+255 (Rx^6) (Ry^2)-764 (Rx^4) (Ry^4)+255 (Rx^2) (Ry^6)+45 (Ry^8))))/(4 (((d-Rx-Ry))^(13/2)) ((((d+Rx-Ry) (d-Rx+Ry) (d+Rx+Ry)))^(13/2))))}, {(((π^2) (Rx^2) (Ry^2) (5 (d^22)+5 (d^20) ((Rx^2)+(Ry^2))-(((Rx-Ry))^8) (((Rx+Ry))^8) ((Rx^2)+(Ry^2)) ((Rx^4)+5 (Rx^2) (Ry^2)+(Ry^4))-5 (d^18) (23 (Rx^4)-67 (Rx^2) (Ry^2)+23 (Ry^4))+5 (d^16) ((Rx^2)+(Ry^2)) (53 (Rx^4)-143 (Rx^2) (Ry^2)+53 (Ry^4))-(d^2) ((((Rx^2)-(Ry^2)))^6) (13 (Rx^8)+181 (Rx^6) (Ry^2)+396 (Rx^4) (Ry^4)+181 (Rx^2) (Ry^6)+13 (Ry^8))-4 (d^14) (26 (Rx^8)+515 (Rx^6) (Ry^2)-1350 (Rx^4) (Ry^4)+515 (Rx^2) (Ry^6)+26 (Ry^8))+(d^4) (((Rx-Ry))^4) (((Rx+Ry))^4) ((Rx^2)+(Ry^2)) (59 (Rx^8)-228 (Rx^6) (Ry^2)-2602 (Rx^4) (Ry^4)-228 (Rx^2) (Ry^6)+59 (Ry^8))-4 (d^12) ((Rx^2)+(Ry^2)) (98 (Rx^8)-1365 (Rx^6) (Ry^2)+2780 (Rx^4) (Ry^4)-1365 (Rx^2) (Ry^6)+98 (Ry^8))-5 (d^6) ((((Rx^2)-(Ry^2)))^2) ((Rx^12)-438 (Rx^10) (Ry^2)+3 (Rx^8) (Ry^4)+2436 (Rx^6) (Ry^6)+3 (Rx^4) (Ry^8)-438 (Rx^2) (Ry^10)+(Ry^12))-10 (d^8) ((Rx^2)+(Ry^2)) (32 (Rx^12)+150 (Rx^10) (Ry^2)-1698 (Rx^8) (Ry^4)+3081 (Rx^6) (Ry^6)-1698 (Rx^4) (Ry^8)+150 (Rx^2) (Ry^10)+32 (Ry^12))+(d^10) (616 (Rx^12)-2772 (Rx^10) (Ry^2)-8964 (Rx^8) (Ry^4)+24330 (Rx^6) (Ry^6)-8964 (Rx^4) (Ry^8)-2772 (Rx^2) (Ry^10)+616 (Ry^12))))/(5 (((d-Rx-Ry))^(17/2)) ((((d+Rx-Ry) (d-Rx+Ry) (d+Rx+Ry)))^(17/2))))}});

(***********************************)
(*Infinite series approach for the P-(Subscript[D, IP]) interaction using the table integral and RPCS2 coordinates.*)
(***********************************)
(*The expression is general and the results are equal to the expressions obtained by the direct integration shown above.*)
PDdouble[m_,Rx_,t_]=FullSimplify[π*Sum[((Pochhammer[m/2,k]*Pochhammer[m/2,k])/(Pochhammer[1,k]*(k+1)!))*((Rx^(2 k+2))/(t^(2 k+m))),{k,0,Infinity}]]

(*The expression is valid for arbitrary m*)
rep={m->0.5,Rx->1,t->2};
NIntegrate[2*ArcCos[(p^2+t^2-Rx^2)/(2*t*p)]*p^(-m+1)/.rep,{p,t-Rx,t+Rx}/.rep]
PDdouble@@rep[[All,2]]

(*However, since the asymptotic scaling factor for small separations is m-2, the law is valid only for m>2*)
tab1={};
tab2={};

For[mm=3,mm<13,mm++,

f = D[Log[PDdouble[mm,Rx,Rx+q2]/.{q2->Exp[q2]}],q2]/.{Rx->1};

AppendTo[tab1,Table[{q2,N[Re[f],16]},{q2,-10,10,1}]];

AppendTo[tab2,mm]];


ListLinePlot[tab1,AxesLabel->{"Log[q2]","Scaling factor of (Subscript[Π, D - (Subscript[D, IP])])"},PlotLegends->tab2]

(*WM14 cannot recognize that these P-(Subscript[D, IP]) laws are the same for m=3,5,7,...*)
Table[FullSimplify[FunctionExpand[PDdouble[m,Rx,t]]-tabIntP[[m-2]]],{m,3,12}]

(*But numerically, the results are the same*)
Table[FullSimplify[N[(PDdouble[m,Rx,t]/.{Rx->1,t->110/10})-(tabIntP[[m-2]]/.{Rx->1,t->110/10}),16]],{m,3,12}]

(*Therefore, for even m, our general P-(Subscript[D, IP]) in the form of a hypergeometric series returns rational functions. For odd m, we get relations between elliptic and hypergeometric functions that are numerically equivalent but WM14 does not recognize them analytically equal.*)
(*******************************)
(*The next step is to find the general D-(Subscript[D, IP]) law. It follows from another integration using the table integral, as done in Langbein 1972. We could not do this integration with Mathematica.*)
(*This is the final expression to be used, with a finite limit of the outer sum*)
DDdouble[m_,Rx_,Ry_,d_,lim_]=FullSimplify[((π^2)/(d^(m-4)))*Sum[Sum[FunctionExpand[(((Pochhammer[n+m/2,k])^2)/(Pochhammer[1,k]))*(((Pochhammer[m/2,n])^2)/(Pochhammer[1,n]))*(1/((k+1)!))*(1/((n+1)!))]*((((Rx/d)))^(2k+2))*((((Ry/d)))^(2n+2)),{k,0,Infinity}],{n,0,lim}],Assumptions->k∈Integers&&n∈Integers&&m∈Integers&&k>0&&n>0&&v>4&&Rx>0&&Ry>0&&d>0]

(*If we want a general expression, for the paper, we insert infinity as the upper limit*)
FullSimplify[((π^2)/(d^(m-4)))*Sum[Sum[FunctionExpand[(((Pochhammer[n+m/2,k])^2)/(Pochhammer[1,k]))*(((Pochhammer[m/2,n])^2)/(Pochhammer[1,n]))*(1/((k+1)!))*(1/((n+1)!))]*((((Rx/d)))^(2k+2))*((((Ry/d)))^(2n+2)),{k,0,Infinity}],{n,0,Infinity}],Assumptions->k∈Integers&&n∈Integers&&m∈Integers&&k>0&&n>0&&v>4&&Rx>0&&Ry>0&&d>0]

(*The expression is valid for arbitrary m. The sum converges slowly for small gaps*)
rep={m->2.2,Rx->1,Ry->1.5,d->2.51,lim->100};
NIntegrate[4*ArcCos[(((t^2)+(d^2)-(Ry^2))/(2*t*d))]ArcCos[(((p^2)+(t^2)-(Rx^2))/(2*t*p))]*(p^(-m)) p t/.rep,{t,d-Ry,d+Ry}/.rep,{p,t-Rx,t+Rx}/.rep]
DDdouble@@rep[[All,2]]

(*WM14 does not recognize that these D-(Subscript[D, IP]) laws for even m are the same*)
Table[FullSimplify[FunctionExpand[DDdouble[m,Rx,Ry,d,Infinity]]-tabIntPT[[m-3]]],{m,4,12,2}]

(*It turns out that this expression in the form of a double series cannot return rational functions for even m. On the other hand, for the special case R=Rx=Ry, we obtain an expression in the form of a generalized hypergeometric function. For this expression, it is necessary to write down the double series in WM14 inside one Sum function*)
DDdoubleRR[m_,R_,d_]=FullSimplify[((π^2)/(d^(m-4)))*Sum[FunctionExpand[(((Pochhammer[n+m/2,k])^2)/(Pochhammer[1,k]))*(((Pochhammer[m/2,n])^2)/(Pochhammer[1,n]))*(1/((k+1)!))*(1/((n+1)!))]*((((Rx/d)))^(2k+2))*((((Ry/d)))^(2n+2))/.{Rx->R,Ry->R},{k,0,Infinity},{n,0,Infinity}],Assumptions->k∈Integers&&n∈Integers&&m∈Integers&&k>0&&n>0&&v>4&&Rx>0&&Ry>0&&d>0]

(*Again, the expression is valid for arbitrary m*)
rep={m->0.2,Rx->1,d->2.01};
NIntegrate[4*ArcCos[(((t^2)+(d^2)-(Rx^2))/(2*t*d))]ArcCos[(((p^2)+(t^2)-(Rx^2))/(2*t*p))]*(p^(-m)) p t/.rep,{t,d-Rx,d+Rx}/.rep,{p,t-Rx,t+Rx}/.rep]
DDdoubleRR@@rep[[All,2]]

(*However, since the asymptotic scaling factor for small separations is m-3.5, the law is valid only for m>3.5*)
tab1={};
tab2={};
For[mm=4,mm<13,mm++,
f = D[Log[DDdoubleRR[mm,Rx,Rx+Rx+q2]/.{q2->Exp[q2]}],q2]/.{Rx->1};
AppendTo[tab1,Table[{q2,N[Re[f],16]},{q2,-10,10,1}]];
AppendTo[tab2,mm]];

ListLinePlot[tab1,AxesLabel->{"Log[q2]","Scaling factor of (Subscript[Π, D - (Subscript[D, IP])])"},PlotLegends->tab2]

(*For this special case, WM14 recognize that even terms are the same*)
Table[FullSimplify[DDdoubleRR[m*2+2,R,2*R+q2]-(tabIntPT[[m,1]]/.{Rx->R,Ry->R,d->2*R+q2}),Assumptions->R>0&&q2>0],{m,1,5}]

(*For the reference, let us write here the D-(Subscript[D, IP]) law for small separations, as derived in Langbein 1972.*)
langShortApp[m_,Rx_,Ry_,q2_]=(2^(5/2-m))*(π^(3/2))*(Sqrt[Rx*(Ry/(Rx+Ry))])*((Gamma[m-7/2])/((Gamma[m/2])^2))*(q2^(-m+7/2));
(*Final steric (m=12) and vdW (m=6) D-(Subscript[D, IP]) laws. We will use these expressions for the approximate D-D law (with offset)*)
exDD6[Rx_,Ry_,q2_]=-(((π^2) (Rx^2) (Ry^2) (-2 (d^4)+((((Rx^2)-(Ry^2)))^2)+(d^2) ((Rx^2)+(Ry^2))))/(2 (((d-Rx-Ry))^(5/2)) ((((d+Rx-Ry) (d-Rx+Ry) (d+Rx+Ry)))^(5/2))))/.d->q2+Rx+Ry;
exDD12[Rx_,Ry_,q2_]=((π^2) (Rx^2) (Ry^2) (5 (d^22)+5 (d^20) ((Rx^2)+(Ry^2))-(((Rx-Ry))^8) (((Rx+Ry))^8) ((Rx^2)+(Ry^2)) ((Rx^4)+5 (Rx^2) (Ry^2)+(Ry^4))-5 (d^18) (23 (Rx^4)-67 (Rx^2) (Ry^2)+23 (Ry^4))+5 (d^16) ((Rx^2)+(Ry^2)) (53 (Rx^4)-143 (Rx^2) (Ry^2)+53 (Ry^4))-(d^2) ((((Rx^2)-(Ry^2)))^6) (13 (Rx^8)+181 (Rx^6) (Ry^2)+396 (Rx^4) (Ry^4)+181 (Rx^2) (Ry^6)+13 (Ry^8))-4 (d^14) (26 (Rx^8)+515 (Rx^6) (Ry^2)-1350 (Rx^4) (Ry^4)+515 (Rx^2) (Ry^6)+26 (Ry^8))+(d^4) (((Rx-Ry))^4) (((Rx+Ry))^4) ((Rx^2)+(Ry^2)) (59 (Rx^8)-228 (Rx^6) (Ry^2)-2602 (Rx^4) (Ry^4)-228 (Rx^2) (Ry^6)+59 (Ry^8))-4 (d^12) ((Rx^2)+(Ry^2)) (98 (Rx^8)-1365 (Rx^6) (Ry^2)+2780 (Rx^4) (Ry^4)-1365 (Rx^2) (Ry^6)+98 (Ry^8))-5 (d^6) ((((Rx^2)-(Ry^2)))^2) ((Rx^12)-438 (Rx^10) (Ry^2)+3 (Rx^8) (Ry^4)+2436 (Rx^6) (Ry^6)+3 (Rx^4) (Ry^8)-438 (Rx^2) (Ry^10)+(Ry^12))-10 (d^8) ((Rx^2)+(Ry^2)) (32 (Rx^12)+150 (Rx^10) (Ry^2)-1698 (Rx^8) (Ry^4)+3081 (Rx^6) (Ry^6)-1698 (Rx^4) (Ry^8)+150 (Rx^2) (Ry^10)+32 (Ry^12))+(d^10) (616 (Rx^12)-2772 (Rx^10) (Ry^2)-8964 (Rx^8) (Ry^4)+24330 (Rx^6) (Ry^6)-8964 (Rx^4) (Ry^8)-2772 (Rx^2) (Ry^10)+616 (Ry^12))))/(5 (((d-Rx-Ry))^(17/2)) ((((d+Rx-Ry) (d-Rx+Ry) (d+Rx+Ry)))^(17/2)))/.d->q2+Rx+Ry;