(*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&&Rx+2*Ry+q2>=t>=Rx+q2&&2*R1+2*R2+q2>=p>=q2&&q2>0&&(*Element[q1,Reals]*)q1>0&&d>Rx+Ry&&d>Rx+q2&&d>Ry+q2&&d>0&&d<Rx+q2+Ry+Rx&&d<Rx+q2+Ry+Ry&&d<Rx+q2+Ry+q2&&t>0&&p>0&&π/2>th>0;

(*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]]]

(**********************************)
(**********************************)
(**********************************)
(*We derive a vdW disk-disk law with an offset (q2-gap, q1-offset) using RPCS1. The integral to solve is 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))]*((Sqrt[(p^2)+(q1^2)])^(-m)) p t,p]*)
(**********************************)
(**********************************)

m=6;
(*First, let us find a P-D law*)
integrandP=2 ArcCos[(p^2+t^2-Rx^2)/(2*p*t)]*((Sqrt[p^2+q1^2])^(-m))*p;
(*WM 14 cannot integrate this directly*)
Integrate[integrandP,{p,t-Rx,t+Rx} ]

(*Therefore, we look for the indefinite integral*)
indefP6=Simplify[Integrate[integrandP,p ]]

(*Plot of the primitive function shows that its real part has a discontinuity*)
r={t->1,Rx->0.5,q1->1};
ReImPlot[indefP6/.r,{p,t-Rx,t+Rx}/.r,PlotRange->All]

(*WM14 can find a discontinuity point*)
FunctionDiscontinuities[indefP6,p,Complexes];
discP=Simplify[SolveValues[%,p][[1,1]]]
discP/.r

(*When we know the discontinuity point, then the definite integral follows from the fundamental theorem of calculus. It turns out that if we use the assumption q1∈Reals, a wrong results follows. If we use assumption q1>0, the result is correct. This is the consequence of the complex-valued integral, and we can approach it in several ways, but the simplest one is to assume q1>0, which does not influence subsequent derivation*)

i61Lower=Simplify[Limit[indefP6,p->t-Rx,Direction->"FromAbove"]]
i61Upper=Simplify[Limit[indefP6,p->discP,Direction->"FromBelow"]]
i62Lower=Simplify[Limit[indefP6,p->discP,Direction->"FromAbove"]]
i62Upper=Simplify[Limit[indefP6,p->t+Rx,Direction->"FromBelow"]]
idef6Final=FullSimplify[i61Upper-i61Lower+i62Upper-i62Lower]

(*The result*)
idef6Final=(π (-(q1^6)+(q1^4) ((Rx^2)-3 (t^2)+(Sqrt[((((q1^2)+(Rx^2)))^2)+2 (q1-Rx) (q1+Rx) (t^2)+(t^4)]))+(((Rx-t))^2) (((Rx+t))^2) ((Rx^2)-(t^2)+(Sqrt[((((q1^2)+(Rx^2)))^2)+2 (q1-Rx) (q1+Rx) (t^2)+(t^4)]))+(q1^2) ((Rx^2)+(t^2)) (3 (Rx^2)-3 (t^2)+2 (Sqrt[((((q1^2)+(Rx^2)))^2)+2 (q1-Rx) (q1+Rx) (t^2)+(t^4)]))))/(4 (q1^4) (((((((q1^2)+(Rx^2)))^2)+2 (q1-Rx) (q1+Rx) (t^2)+(t^4)))^(3/2)));

(*This expression is correct, but has a numerical issue for q1->0 which requires arbitrary precision*)
rep={Rx->1,Ry->2,q1->1/100,t->2};
NIntegrate[2*ArcCos[(((p^2)+(t^2)-(Rx^2))/(2*p*t))]*((((Sqrt[((p^2)+(q1^2))])))^(-6))*p/.rep,{p,t-Rx,t+Rx}/.rep]
1.*idef6Final/.rep
N[idef6Final/.rep,20]

(*We can solve this numerical issue. Let us "take out" q2 from the square roots of the numerator. In essence, find a-(Sqrt[b]), multiply by a+(Sqrt[b]), and get (a^2)-b. Here, oppNum is the numerator but with an opposite sign in front of each square root*)
num=Numerator[idef6Final]
den=Denominator[idef6Final]
oppNum=num/.{(Sqrt[x_])->-(Sqrt[x])};
idefP6Reg=Simplify[((Simplify[Expand[num*oppNum]])/(den*oppNum))]

(*The expression can be more simplified, but this form is fine as well, and does not influence further integration. 
Numerical test of accuracy*)
m=6;
rep={Rx->1,Ry->2,q1->1/10000.,t->2};
NIntegrate[2*ArcCos[(((p^2)+(t^2)-(Rx^2))/(2*p*t))]*((((Sqrt[((p^2)+(q1^2))])))^(-m))*p/.rep,{p,t-Rx,t+Rx}/.rep]
idefP6Reg/.rep
idef6Final/.rep

(*Define this integral explicitly, for further evaluation. It is the P-D law for m=6.*)
pointDisk6=(π (Rx^2) (2 (q1^6)+(Rx^6)-3 (Rx^2) (t^4)+2 (t^6)+(q1^4) (5 (Rx^2)+6 (t^2))+(q1^2) (4 (Rx^4)+6 (Rx^2) (t^2)+6 (t^4))))/(((((q1^4)+((((Rx^2)-(t^2)))^2)+2 (q1^2) ((Rx^2)+(t^2))))^(3/2)) ((q1^6)+((((Rx^2)-(t^2)))^2) (-(Rx^2)+(t^2)+(Sqrt[(q1^4)+((((Rx^2)-(t^2)))^2)+2 (q1^2) ((Rx^2)+(t^2))]))+(q1^4) (-(Rx^2)+3 (t^2)+(Sqrt[(q1^4)+((((Rx^2)-(t^2)))^2)+2 (q1^2) ((Rx^2)+(t^2))]))+(q1^2) ((Rx^2)+(t^2)) (-3 (Rx^2)+3 (t^2)+2 (Sqrt[(q1^4)+((((Rx^2)-(t^2)))^2)+2 (q1^2) ((Rx^2)+(t^2))]))));

(*The limit value for q1->0 is well-defined and gives the expression that we already have for q1=0, P-(Subscript[D, IP]) vdW law.*)
Simplify[pointDisk6/.q1->0]

(*The pointDisk6 expression is the result of integrating P-P law over the area of disk X. It is a point-disk law with an offset. The distance of a point Y from the projection of the center of disk X is t. For example, if we integrate this P-D law over q1, we obtain the point-cylinder law.*)
(*The point-disk law, a numerical test*)
m=6;
Rr1=1;
Rr2=1;
d=q2+Rx+Ry;
d1=q2+Ry;
rep={Rx->Rr1,Ry->Rr2,t->1.1*(Rr1),q2->1/10,q1->0};
rep1={Rx->Rr1,Ry->Rr2,q2->1/10,q1->0};

(*two different approaches, P-D w.r.t. disk X and disk Y*)
2*NIntegrate[ArcCos[(t^2+p^2-Rx^2)/(2*t*p)]*(Sqrt[p^2+q1^2])^-m*p/. rep,{p,t-Rx,t+Rx}/. rep]
2*NIntegrate[ArcCos[(t^2+d1^2-Ry^2)/(2*t*d1)]*(Sqrt[t^2+q1^2])^-m*t/. rep1,{t,q2,q2+2*Ry}/. rep1]
pointDisk6/. rep

(*To derive the D-D law, it is left to integrate w.r.t. the area of disk Y. Since the integration w.r.t. the angle is already made, it is left to integrate over coordinate t. If we use variable d=(q2+Rx+Ry), a simpler primitive function is obtained. Unfortunately, the limits are then more complicated. In what follows, one limit is very simple, without the elliptic integrals. For the approach with d parameter, both limits have elliptic integrals. Let us note that the primitive function is very complicated and it is quite interesting that it exists*)

indef6T=Integrate[2*pointDisk6*ArcCos[((-(Ry^2)+(((q2+Rx+Ry))^2)+(t^2))/(2 (q2+Rx+Ry) t))]*t,t]

(*There is no discontinuity of the primitive function, and the imaginary part is constant*)
rep={Rx->1.5,Ry->0.1,q1->0.01,q2->0.01};
Plot[Re[indef6T/.rep],{t,Rx+q2,Rx+q2+2*Ry}/.rep]
Plot[Im[indef6T]/.rep,{t,Rx+q2,Rx+q2+2*Ry}/.rep]

(*Extracting the real part is slow, and it seems that we get no benefit from it*)
reIndef6T=ComplexExpand[Re[indef6T]]
i6TLower=Limit[reIndef6T,t->Rx+q2,Direction->"FromAbove"]
i6TUpper=Limit[reIndef6T,t->Rx+q2+2*Ry,Direction->"FromBelow"]

(*Therefore, we search for limits of the complete primitive function. The lower limit follows easily, but the upper requires some additional effort*)
i6TLower=Limit[indef6T,t->Rx+q2,Direction->"FromAbove"]
i6TUpper=Limit[indef6T,t->Rx+q2+2*Ry,Direction->"FromBelow"]

(*The result for the lower limit*)
i6TLower=((π (Ry^2) (π-I Log[4 Ry (q2+Rx+Ry)]))/(4 (q1^4)));

(*So, WM14 cannot find the upper limit. Let us find the limit in an alternative way. First, we decompose indef6T into 16 additive parts - dec1add16. It turns out that the limit of the 13th part is not well-defined - it consists of elliptic functions, and we will deal with it separately*)
decIadd16=Level[Apart[indef6T],1];
res1={};
decIadd16//Length
For[i=1,i<13,i++,
apaL=Limit[decIadd16[[i]],t->q2+Rx+2 Ry,Direction->"FromBelow"];
Print[i,": ",apaL];
AppendTo[res1,apaL]]

For[i=14,i<17,i++,
apaL=Limit[decIadd16[[i]],t->q2+Rx+2 Ry,Direction->"FromBelow"];
Print[i,": ",apaL];
AppendTo[res1,apaL]]

(*So, we need to find the limit of the 13th part of decIadd16.
Then, the upper limit is: Total[res1]+limDecIadd16[[13]]. Let us designate i13=limDecIadd16[[13]]*)
(*i6TUpper=Total[res1]+i13limit;*)

i13=Simplify[decIadd16[[13]]]

(*Now, we decompose i13=decIadd16[[13]] into two multiplicative parts. The first part has a limit, and the second is problematic for WM14, but we can show that the correct result is obtained by a straightforward replacement of t with a limit value. As a result, we we obtain the upper limit of the primitive function*)

le=Length[i13]
i6TUpper=Simplify[Total[res1]+Simplify[Limit[i13[[1;;(le-1)]],t->Rx+q2+2*Ry,Direction->"FromBelow"]]*Simplify[i13[[le]]/.t->Rx+q2+2*Ry]]

(*The D-D interaction potential for m=6*)
i6Tdef=Simplify[Expand[i6TUpper-i6TLower]];

(*The final result*)
i6Tdef=(π (-(1/(Sqrt[(((q1^2)+2 I q1 Rx+q2 (q2+2 Rx))/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))]))(q1+I q2) (q1+I Rx) (q2+Rx) (q1-I (q2+2 Rx)) ((q1^2)+2 I q1 Rx+q2 (q2+2 Rx)) (Sqrt[(((q1-I q2) (q1+I (q2+2 Rx)) ((((q1^2)-2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))^2))/((q1^2)-2 I q1 Rx+q2 (q2+2 Rx)))]) EllipticE[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Rx)) (q1-I (q2+2 Ry)) (q1+I (q2+2 (Rx+Ry)))))]+(q2+Rx) (Sqrt[(((q1-I q2) (q1+I (q2+2 Rx)) ((((q1^2)-2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))^2))/((q1^2)-2 I q1 Rx+q2 (q2+2 Rx)))]) (Sqrt[(((q1^2)+2 I q1 Rx+q2 (q2+2 Rx))/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))]) ((q1^5)+I (q1^4) Rx+I q2 Rx (q2+2 Rx) (q2+2 Ry) (q2+2 (Rx+Ry))+2 I (q1^2) Rx ((q2^2)-2 (Rx^2)+2 Rx Ry+6 (Ry^2)+2 q2 (Rx+Ry))+2 (q1^3) ((q2^2)+2 q2 (Rx+Ry)+2 Ry (Rx+Ry))+q1 ((q2^4)+4 (q2^3) (Rx+Ry)-8 (Rx^2) Ry (Rx+Ry)+4 (q2^2) Ry (3 Rx+Ry)-8 q2 ((Rx^3)-Rx (Ry^2)))) EllipticK[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Rx)) (q1-I (q2+2 Ry)) (q1+I (q2+2 (Rx+Ry)))))]+2 ((q1^2)+2 I q1 Rx+q2 (q2+2 Rx)) ((q1-I Rx) (Rx^2) (q2+Rx+2 Ry) √(((q1-I q2) (q1+I (q2+2 Rx)) ((((q1^2)-2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))^2))/((q1^2)-2 I q1 Rx+q2 (q2+2 Rx))) (Sqrt[(((q1^2)+2 I q1 Rx+q2 (q2+2 Rx))/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))]) EllipticPi[-((4 (((q1+I Rx))^2) Ry (q2+Rx+Ry))/((((q2+Rx))^2) (I q1+q2+2 Ry) (-I q1+q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Rx)) (q1-I (q2+2 Ry)) (q1+I (q2+2 (Rx+Ry)))))]+(q1+I Rx) (q2+Rx) (π (Ry^2) ((q1^2)-2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry)))+((Rx^2)-2 (Ry^2)) √(((q1-I q2) (q1+I (q2+2 Rx)) ((((q1^2)-2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))^2))/((q1^2)-2 I q1 Rx+q2 (q2+2 Rx))) (Sqrt[(((q1^2)+2 I q1 Rx+q2 (q2+2 Rx))/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))]) EllipticPi[((4 Ry (q2+Rx+Ry))/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Rx)) (q1-I (q2+2 Ry)) (q1+I (q2+2 (Rx+Ry)))))]))))/(8 (q1^4) (q1-I q2) (q1+I Rx) (q2+Rx) (q1+I (q2+2 Rx)) (q1+I (q2+2 Ry)) (q1-I (q2+2 (Rx+Ry))));

(*Let us simplify it a bit*)

i6TdefFS=FullSimplify[Expand[i6Tdef]];

(*the result*)
i6TdefFS=(π (2 π (q1-I q2) (q1+I Rx) (q2+Rx) (q1+I (q2+2 Rx)) (Ry^2) ((q1^2)-2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry)))+(Sqrt[(((q1-I q2) (q1+I (q2+2 Rx)) ((((q1^2)-2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))^2))/((q1^2)-2 I q1 Rx+q2 (q2+2 Rx)))]) (Sqrt[(((q1^2)+2 I q1 Rx+q2 (q2+2 Rx))/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))]) (-((q1+I q2) (q1+I Rx) (q2+Rx) (q1-I (q2+2 Rx)) ((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))) EllipticE[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))])+(q2+Rx) ((q1-I q2) (q1 ((q1^2)+(q2^2))-I (((q1+I q2))^2) Rx-2 (q1-I q2) (Rx^2)) (q1+I (q2+2 Rx))+4 (q2+Rx) ((q1^3)+I (q1^2) Rx+I q2 Rx (q2+2 Rx)+q1 ((q2^2)+2 q2 Rx-2 (Rx^2))) Ry+4 (q1+I Rx) ((q1^2)+2 I q1 Rx+q2 (q2+2 Rx)) (Ry^2)) EllipticK[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]+2 (q1-I q2) (-I q1+q2+2 Rx) ((Rx^2) (I q1+Rx) (q2+Rx+2 Ry) EllipticPi[-((4 (((q1+I Rx))^2) Ry (q2+Rx+Ry))/((((q2+Rx))^2) (I q1+q2+2 Ry) (-I q1+q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]+I (q1+I Rx) (q2+Rx) ((Rx^2)-2 (Ry^2)) EllipticPi[((4 Ry (q2+Rx+Ry))/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]))))/(8 (q1^4) (q1-I q2) (q1+I Rx) (q2+Rx) (q1+I (q2+2 Rx)) (q1+I (q2+2 Ry)) (q1-I (q2+2 (Rx+Ry))));
(*We are combining square roots together and the result is*)
i6TdefFS1=(π (2 π  (q1+I Rx) (q2+Rx)  (Ry^2) ((q1^2)-2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry)))+(Sqrt[(1/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]) *  ((q1^2)-2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry)))(-((q1+I q2) (q1+I Rx) (q2+Rx) (q1-I (q2+2 Rx)) ((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))) EllipticE[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))])+(q2+Rx) ((q1-I q2) (q1 ((q1^2)+(q2^2))-I (((q1+I q2))^2) Rx-2 (q1-I q2) (Rx^2)) (q1+I (q2+2 Rx))+4 (q2+Rx) ((q1^3)+I (q1^2) Rx+I q2 Rx (q2+2 Rx)+q1 ((q2^2)+2 q2 Rx-2 (Rx^2))) Ry+4 (q1+I Rx) ((q1^2)+2 I q1 Rx+q2 (q2+2 Rx)) (Ry^2)) EllipticK[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]+2 (q1-I q2) (-I q1+q2+2 Rx) ((Rx^2) (I q1+Rx) (q2+Rx+2 Ry) EllipticPi[-((4 (((q1+I Rx))^2) Ry (q2+Rx+Ry))/((((q2+Rx))^2) (I q1+q2+2 Ry) (-I q1+q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]+I (q1+I Rx) (q2+Rx) ((Rx^2)-2 (Ry^2)) EllipticPi[((4 Ry (q2+Rx+Ry))/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]))))/(8 (q1^4)  (q1+I Rx) (q2+Rx)  (q1+I (q2+2 Ry)) (q1-I (q2+2 (Rx+Ry))));

(*another FullSimplify*)
i6TdefFS2=FullSimplify[i6TdefFS1];

(*the result*)
i6TdefFS2=(π ((q1^2)-2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))) (2 π (q1+I Rx) (q2+Rx) (Ry^2)+(Sqrt[(1/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]) (-((q1+I q2) (q1+I Rx) (q2+Rx) (q1-I (q2+2 Rx)) ((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))) EllipticE[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))])+(q2+Rx) ((q1-I q2) (q1 ((q1^2)+(q2^2))-I (((q1+I q2))^2) Rx-2 (q1-I q2) (Rx^2)) (q1+I (q2+2 Rx))+4 (q2+Rx) ((q1^3)+I (q1^2) Rx+I q2 Rx (q2+2 Rx)+q1 ((q2^2)+2 q2 Rx-2 (Rx^2))) Ry+4 (q1+I Rx) ((q1^2)+2 I q1 Rx+q2 (q2+2 Rx)) (Ry^2)) EllipticK[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]+2 (q1-I q2) (-I q1+q2+2 Rx) ((Rx^2) (I q1+Rx) (q2+Rx+2 Ry) EllipticPi[-((4 (((q1+I Rx))^2) Ry (q2+Rx+Ry))/((((q2+Rx))^2) (I q1+q2+2 Ry) (-I q1+q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]+I (q1+I Rx) (q2+Rx) ((Rx^2)-2 (Ry^2)) EllipticPi[((4 Ry (q2+Rx+Ry))/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]))))/(8 (q1^4) (q1+I Rx) (q2+Rx) (q1+I (q2+2 Ry)) (q1-I (q2+2 (Rx+Ry))));

(*Canceling some terms manually. This is the simplest form. We can choose to introduce some substitutions here*)
i6TdefFS3=(1/(8 (q1^4)))π (2 π (Ry^2)+(Sqrt[(1/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]) (-((q1+I q2) (q1-I (q2+2 Rx)) ((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))) EllipticE[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))])+(1/(q1+I Rx))((q1-I q2) (q1 ((q1^2)+(q2^2))-I (((q1+I q2))^2) Rx-2 (q1-I q2) (Rx^2)) (q1+I (q2+2 Rx))+4 (q2+Rx) ((q1^3)+I (q1^2) Rx+I q2 Rx (q2+2 Rx)+q1 ((q2^2)+2 q2 Rx-2 (Rx^2))) Ry+4 (q1+I Rx) ((q1^2)+2 I q1 Rx+q2 (q2+2 Rx)) (Ry^2)) EllipticK[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]+2 (q1-I q2) (-I q1+q2+2 Rx) ((1/((q1+I Rx) (q2+Rx)))(Rx^2) (I q1+Rx) (q2+Rx+2 Ry) EllipticPi[-((4 (((q1+I Rx))^2) Ry (q2+Rx+Ry))/((((q2+Rx))^2) (I q1+q2+2 Ry) (-I q1+q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]+I ((Rx^2)-2 (Ry^2)) EllipticPi[((4 Ry (q2+Rx+Ry))/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))])));
(*replacements*)
drep={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};
drepRev={d->q2+Rx+Ry,d1->q2+Rx+2 Ry,d2->q2+2 (Rx+Ry),d3->q2+Rx,d4->q2+2 Rx,d5->q2+2 Ry};

(*this is the expression for the paper*)
i6TdefFSsubs=i6TdefFS3/.{1/((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry))))->1/a1}/.{-((q1+I q2) (q1-I (q2+2 Rx)) ((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))) EllipticE[x_])->-a1 EllipticE[x]}/.{(1/((((q2+Rx))^2) (I q1+q2+2 Ry) (-I q1+q2+2 (Rx+Ry))))->(1/a2),(1/((q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))))->(1/a3)}/.{Ry (q2+Rx+Ry)->a4}/.{-16 I a4 q1 Rx->a5,-4 a4 (((q1+I Rx))^2)->a6}/.{(1/(q1+I Rx))((q1-I q2) (q1 ((q1^2)+(q2^2))-I (((q1+I q2))^2) Rx-2 (q1-I q2) (Rx^2)) (q1+I (q2+2 Rx))+4 (q2+Rx) ((q1^3)+I (q1^2) Rx+I q2 Rx (q2+2 Rx)+q1 ((q2^2)+2 q2 Rx-2 (Rx^2))) Ry+4 (q1+I Rx) ((q1^2)+2 I q1 Rx+q2 (q2+2 Rx)) (Ry^2))->a1-(1/(q1+I Rx))4 q1 Rx ((q1^2) Rx+Rx (q2+2 Ry) (q2+2 (Rx+Ry))+2 I q1 (-q2 Ry+(Rx-2 Ry) (Rx+Ry)))}/.{-(1/(q1+I Rx))4 q1 Rx ((q1^2) Rx+Rx (q2+2 Ry) (q2+2 (Rx+Ry))+2 I q1 (-q2 Ry+(Rx-2 Ry) (Rx+Ry)))->a7}/.drep;

a1s=((q1+I q2) (q1-I (q2+2 Ry)) ((q1^2)+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry))))/.drep;
a2s=(((q2+Rx))^2) (I q1+q2+2 Ry) (-I q1+q2+2 (Rx+Ry))/.drep;
a3s=(q1^2)+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))/.drep;
a4s=Ry (q2+Rx+Ry)/.drep;
a5s=-16 I a4s q1 Rx/.drep;
a6s=-4 a4s (((q1+I Rx))^2)/.drep;
a7s=-(1/(q1+I Rx))4 q1 Rx ((q1^2) Rx+Rx (q2+2 Ry) (q2+2 (Rx+Ry))+2 I q1 (-q2 Ry+(Rx-2 Ry) (Rx+Ry)))/.drep;

(*Check if the replacements are correct*)
FullSimplify[(i6TdefFSsubs/.{a1->a1s,a2->a2s,a3->a3s,a4->a4s,a5->a5s,a6->a6s,a7->a7s}/.drepRev)-i6TdefFS1]

(*Final numerical check, All the expressions are OK. There is a numerical issue for q1->0.*)
rep={Rx->1,Ry->2,q2->1/10,q1->1/10};
r2=a1->(q1+I  q2)  (q1-I  (q2+2  Rx))  (q1^2+2  I  q1  Rx+(q2+2  Ry)  (q2+2  (Rx+Ry))) ;
4*NIntegrate[ArcCos[(((t^2)+(((q2+Rx+Ry))^2)-(Ry^2))/(2*t*(q2+Rx+Ry)))]*ArcCos[(((t^2)+(p^2)-(Rx^2))/(2*t*p))]*((((Sqrt[(p^2)+(q1^2)])))^(-m))*p*t/.rep,{t,q2+Rx,q2+Rx+2*Ry}/.rep,{p,t-Rx,t+Rx}/.rep]
N[i6Tdef/.rep,16]
1.*i6Tdef/.rep
1.*i6TdefFS/.rep
1.*i6TdefFS1/.rep
N[i6TdefFS1/.rep,16]
1.*i6TdefFS2/.rep
N[i6TdefFS2/.rep,16]
1.*i6TdefFS3/.r2/.rep
N[i6TdefFS3/.r2/.rep,16]

1.*(i6TdefFSsubs/.{a1->a1s,a2->a2s,a3->a3s,a4->a4s,a5->a5s,a6->a6s,a7->a7s})/.drepRev/.rep
N[(i6TdefFSsubs/.{a1->a1s,a2->a2s,a3->a3s,a4->a4s,a5->a5s,a6->a6s,a7->a7s})/.drepRev/.rep,16]

(*The derived D-D law is correct and accurate. However, it has an issue with q1->0 which requires arbitrary precision arithmetics. Also, it cannot be integrated along q1, which is somehow expected, since we know that the vdW C-C potential has a double series solution, see Langbein 1972. *)

(*On the other hand, having the exact D-D law, allows us to precisely estimate the error of various approximate D-D laws or even derive new approximate laws. Let us compare the exact D-D law with the two approximate D-D laws*)
ddEx[Rx_,Ry_,q1_,q2_]=1/(8 q1^4) π (2 π Ry^2+Sqrt[1/((q1+I q2) (q1-I (q2+2 Ry)) (q1^2+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry))))] (-((q1+I q2) (q1-I (q2+2 Rx)) (q1^2+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))) EllipticE[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) (q1^2+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))])+1/(q1+I Rx) ((q1-I q2) (q1 (q1^2+q2^2)-I (q1+I q2)^2 Rx-2 (q1-I q2) Rx^2) (q1+I (q2+2 Rx))+4 (q2+Rx) (q1^3+I q1^2 Rx+I q2 Rx (q2+2 Rx)+q1 (q2^2+2 q2 Rx-2 Rx^2)) Ry+4 (q1+I Rx) (q1^2+2 I q1 Rx+q2 (q2+2 Rx)) Ry^2) EllipticK[-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) (q1^2+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]+2 (q1-I q2) (-I q1+q2+2 Rx) (1/((q1+I Rx) (q2+Rx)) Rx^2 (I q1+Rx) (q2+Rx+2 Ry) EllipticPi[-((4 (q1+I Rx)^2 Ry (q2+Rx+Ry))/((q2+Rx)^2 (I q1+q2+2 Ry) (-I q1+q2+2 (Rx+Ry)))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) (q1^2+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))]+I (Rx^2-2 Ry^2) EllipticPi[(4 Ry (q2+Rx+Ry))/(q1^2+2 I q1 Rx+(q2+2 Ry) (q2+2 (Rx+Ry))),-((16 I q1 Rx Ry (q2+Rx+Ry))/((q1+I q2) (q1-I (q2+2 Ry)) (q1^2+2 I q1 Ry+(q2+2 Rx) (q2+2 (Rx+Ry)))))])));

(*The approximate disk-disk law named ISSIP is derived in 2024 Borkovic. It is represented as a product of the approximate D-(Subscript[D, IP]) law and a hypergeometric function. Here we improve the ISSIP law by employing the exact D-(Subscript[D, IP]) law. The new approximate law is designated D-(Subscript[D, app])*)

(*Exact vdW D-(Subscript[D, IP])*)

inPlaneVdWExact=((π^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)));

(*Approximate vdW D-(Subscript[D, IP]), Langbein 1972*)

inPlaneVdWApp=((3 (π^2) (Sqrt[((Rx Ry)/(Rx+Ry))]))/(128 (Sqrt[2]) (q2^(5/2))));

(*Approximate  disk-disk law derived in Borkovic 2024,ISSIP  (improved section-section interaction potential)*)

ISSIP[Rx_,Ry_,q1_,q2_]=inPlaneVdWApp*Hypergeometric2F1[1/4 (-7+2 m),1/4 (-5+2 m),m/2,-(q1/q2)^2]/. m->6;

(*The new approximate disk-disk law (D-(Subscript[D, app]))*)

ddApp[Rx_,Ry_,q1_,q2_]=inPlaneVdWExact*Hypergeometric2F1[(1/4) (-7+2 m),(1/4) (-5+2 m),(m/2),-((((q1/q2)))^2)]/.m->6;

(*We can also find forces along the offset q1 (f1) and the gap q2 (f2)*)
ddF1=FullSimplify[D[ddEx[Rx,Ry,q1,q2],q1]];
ddF2=FullSimplify[D[ddEx[Rx,Ry,q1,q2],q2]];

(*As an additional check, let us find interaction force f2 for two in-plane disks. We compare results from the general D-D law and the special D-(Subscript[D, IP]) case, see Notebook 2. The same expressions are obtained which confirms that the D-D law is general and correct*)
fIP2=FullSimplify[Limit[ddF2,q1->0]]
fIP2spec=FullSimplify[D[inPlaneVdWExact,q2]]
FullSimplify[fIP2-fIP2spec]

(*Test the scaling w.r.t. q2, for different values of parameter c=q1/q2. The exact law shows that the factor c does not influence results for c<0.1 (all graphs overlap). Then, we see that the scaling factor of the D-(Subscript[D, app]) law is insensitive to the ratio q1/q2=c. Therefore, although it is asymptotically correct, it gives wrong results for c>0.1. This confirms that it is usable for short-range interactions when only small parts of disks interact. As the separation increases, the areas of disks that contribute to the interaction increase, and the error of D-(Subscript[D, app]) increases*)

tab={1/100,1/10,1,10,100};
f =- D[Log[ddEx[Rx,Ry,c*q2,q2]/.{q2->Exp[q2]}],q2]/.{Rx->1,Ry->1};
as1=Table[{q2,Re[N[f,100]]},{q2,-8,8,1/4},{c,tab}];
ListLinePlot[Transpose@as1,AxesLabel->{"Log[q2]","Scaling factor of ΠvdW - Exact D-D law"},PlotLegends->{tab}]

tab={1/100,1,100};
f = -D[Log[ddApp[Rx,Ry,c*q2,q2]/.{q2->Exp[q2]}],q2]/.{Rx->1,Ry->1};
as2=Table[{q2,Re[N[f,100]]},{q2,-8,8,1/4},{c,tab}];
ListLinePlot[Transpose@as2,AxesLabel->{"Log[q2]","Scaling factor of ΠvdW - D-(Subscript[D, app])"},PlotLegends->{tab}]

tab={1/100,1,100};
f =- D[Log[ISSIP[Rx,Ry,c*q2,q2]/.{q2->Exp[q2]}],q2]/.{Rx->1,Ry->1};
as3=Table[{q2,Re[N[f,100]]},{q2,-8,8,1/4},{c,tab}];
ListLinePlot[Transpose@as3,AxesLabel->{"Log[q2]","Scaling factor of ΠvdW - ISSIP"},PlotLegends->{tab}]

(*Scaling with respect to the distance of closest points - h. Again, this shows that the scaling factor of the D-(Subscript[D, app]) law is accurate for small values of c=q1/q2*)

a={π/1000,Pi/8,Pi/4,3Pi/8,99Pi/200};

f =- D[Log[ddEx[Rx,Ry,q1,q2]/.{q1->h Sin[alpha],q2->h Cos[alpha],Rx->1,Ry->1}/.{h->Exp[h]}],h];
as=Table[{h,Re[N[f,100]]},{h,-5,10,1},{alpha,a}];
ListLinePlot[Join[Transpose@as],AxesLabel->{"Log[h]","Scaling factor of ΠvdW - Exact D-D"},PlotLegends->a]

f1 = -D[Log[ddApp[Rx,Ry,q1,q2]/.{q1->h Sin[alpha],q2->h Cos[alpha],Rx->1,Ry->1}/.{h->Exp[h]}],h];
as1=Table[{h,Re[N[f1,100]]},{h,-5,10,1},{alpha,a}];
ListLinePlot[Transpose@as1,AxesLabel->{"Log[h]","Scaling factor of ΠvdW - D-(Subscript[D, app])"},PlotLegends->a]

(*ListLinePlot[Join[Transpose@as,Transpose@as1],AxesLabel->{"gap: Log[d]","Scaling factor of ϕvdW"},PlotLegends->a]*)

(*Relative error of the D-(Subscript[D, app]) law blows up for large separations and large c. However, for small c, the error is reasonably small for large separations.*)
tab=Table[(10^ee),{ee,-2,2,1}];
error=Abs[Re[((ddEx[Rx,Ry,c*q2,q2]-ddApp[Rx,Ry,c*q2,q2])/(ddEx[Rx,Ry,c*q2,q2]))*100]]/.{Rx->1,Ry->1};
as=Table[N[error,100],{c,tab}];

LogLogPlot[as,{q2,(10^(-4)),(10^3)},AxesLabel->{"q2","Relative error of the approximate D-(Subscript[D, app]) law [%]"},PlotLegends->{tab},WorkingPrecision->50]

(*Regarding the absolute error, it blows up for small separations and large c*)
tab=Table[(10^ee),{ee,-4,3,1}];
error=Abs[Re[ddEx[Rx,Ry,c*q2,q2]-ddApp[Rx,Ry,c*q2,q2]]]/.{Rx->1,Ry->1};
as=Table[N[error,100],{c,tab}];

LogLogPlot[as,{q2,(10^(-4)),(10^3)},AxesLabel->{"q2","Absolute error of the D-(Subscript[D, app]) law"},PlotLegends->{tab},WorkingPrecision->50]

(*The relative difference between the ISSIP and D-(Subscript[D, app]) laws grows linearly until, approximately, q2=Rx=1. At that point, the difference is 120%*)
tab=Table[(10^ee),{ee,-2,0,1}];
error=Abs[Re[((ddApp[Rx,Ry,c*q2,q2]-ISSIP[Rx,Ry,c*q2,q2])/(ddApp[Rx,Ry,c*q2,q2]))*100]]/.{Rx->1,Ry->1};
as=Table[N[error,100],{c,tab}];

LogLogPlot[as,{q2,(10^(-4)),(10^3)},AxesLabel->{"q2","Relative difference between D-(Subscript[D, app]) and ISSIP [%]"},PlotLegends->{tab},WorkingPrecision->50]

(*Relative error of the D-(Subscript[D, app]) law and ISSIP. ISSIP blows up for large separations. D-(Subscript[D, app]) is excellent for small separations. As c increases, both laws are similar for small separations, but D-(Subscript[D, app]) is better for large separations.*)
tab=Table[(10^ee),{ee,-2,0,1}];
errorApp=Abs[Re[((ddEx[Rx,Ry,c*q2,q2]-ddApp[Rx,Ry,c*q2,q2])/(ddEx[Rx,Ry,c*q2,q2]))*100]]/.{Rx->1,Ry->1};
errorISSIP=Abs[Re[((ddEx[Rx,Ry,c*q2,q2]-ISSIP[Rx,Ry,c*q2,q2])/(ddEx[Rx,Ry,c*q2,q2]))*100]]/.{Rx->1,Ry->1};
asApp=Table[N[errorApp,100],{c,tab}];
asISSIP=Table[N[errorISSIP,100],{c,tab}];

LogLogPlot[{asApp,asISSIP},{q2,(10^(-4)),(10^3)},AxesLabel->{"q2","Relative error of D-(Subscript[D, app]) and ISSIP[%]"},PlotLegends->{"DDapp: c=1/100","DDapp: c=1/10","DDapp: c=1","ISSIP: c=1/100","ISSIP: c=1/10","ISSIP: c=1"},WorkingPrecision->100]

(*Comparison of the absolute errors of the D-(Subscript[D, app]) law and ISSIP. The error of the ISSIP decreases linearly, while the error of the D-(Subscript[D, app]) drops with increased rate for q2>1*)
tab=Table[(10^ee),{ee,-2,0,1}];
errorAbsApp=Abs[Re[ddEx[Rx,Ry,c*q2,q2]-ddApp[Rx,Ry,c*q2,q2]]]/.{Rx->1,Ry->1};
errorAbsISSIP=Abs[Re[ddEx[Rx,Ry,c*q2,q2]-ISSIP[Rx,Ry,c*q2,q2]]]/.{Rx->1,Ry->1};
asAbsApp=Table[N[errorAbsApp,100],{c,tab}];
asAbsISSIP=Table[N[errorAbsISSIP,100],{c,tab}];

LogLogPlot[{asAbsApp,asAbsISSIP},{q2,(10^(-4)),(10^3)},AxesLabel->{"q2","Abssolute error of D-(Subscript[D, app]) and ISSIP"},PlotLegends->{"DDapp: c=1/100","DDapp: c=1/10","DDapp: c=1","ISSIP: c=1/100","ISSIP: c=1/10","ISSIP: c=1"},WorkingPrecision->100]