u12=(SQRT(b1*x+a1)*SQRT(b2*x+a2)*SQRT(b3*y+a3)+SQRT(b1*y+a1)*SQRT(b2*y+a2)*SQ~
RT(b3*x+a3))/(x-y)

u13=(SQRT(b1*x+a1)*SQRT(b2*y+a2)*SQRT(b3*x+a3)+SQRT(b1*y+a1)*SQRT(b2*x+a2)*SQ~
RT(b3*y+a3))/(x-y)

u23=(SQRT(b1*x+a1)*SQRT(b2*y+a2)*SQRT(b3*y+a3)+SQRT(b1*y+a1)*SQRT(b2*x+a2)*SQ~
RT(b3*x+a3))/(x-y)

IHAT(0)=4*R_F((u12+u13)*(u12+u23),(u12+u13)*(u13+u23),(u12+u23)*(u13+u23))

IHAT(-e1)=2*b1*(2*b1*b2*(SQRT(b1*x+a1)+SQRT(b1*y+a1))*(a1*b3-a3*b1)*(a2*b3-a3~
*b2)*SQRT(b1*x+a1)*SQRT(b1*y+a1)*R_D((u12+u13)*(u13+u23),(u12+u23)*(u13+u23),~
(u12+u13)*(u12+u23))+3*(2*b2*(a1*b3-a3*b1)*SQRT(b1*y+a1)*(SQRT(b1*x+a1)*SQRT(~
b1*y+a1)+b1*x+a1)*R_F((u12+u13)*(u12+u23),(u12+u13)*(u13+u23),(u12+u23)*(u13+~
u23))+b1*(SQRT(b1*x+a1)*SQRT(b1*y+a1)*(SQRT(b2*x+a2)*SQRT(b3*y+a3)-SQRT(b2*y+~
a2)*SQRT(b3*x+a3))-(b1*y+a1)*SQRT(b2*x+a2)*SQRT(b3*x+a3)+b1*x*SQRT(b2*y+a2)*S~
QRT(b3*y+a3)+a1*SQRT(b2*y+a2)*SQRT(b3*y+a3))))/(3*(SQRT(b1*x+a1)+SQRT(b1*y+a1~
))*(a1*b2-a2*b1)*(a1*b3-a3*b1)*SQRT(b1*x+a1)*SQRT(b1*y+a1))

u14=SQRT((a1*(a2*b1*b3*b4-a4*b1*b2*b3+b4*u12^2)-b1*(a2*a3*b1*b4-a4*(a3*b1*b2-~
u12^2)))/(a1*b4-a4*b1))

q14=u14*SQRT(b4*x+a4)*SQRT(b4*y+a4)/(SQRT(b1*x+a1)*SQRT(b1*y+a1))

s14=(SQRT(b1*x+a1)*SQRT(b2*y+a2)*SQRT(b3*y+a3)*(b4*x+a4)+SQRT(b1*y+a1)*SQRT(b~
2*x+a2)*SQRT(b3*x+a3)*(b4*y+a4))/((x-y)*SQRT(b1*x+a1)*SQRT(b1*y+a1))

b14=(u12+u14)*(u13+u14)*(u14+u23)*(a4*b1-a1*b4)/(b1*(a1*b2-a2*b1)*(a1*b3-a3*b~
1))

IHAT(-e4)=4*b4*(3*b4*(b14+q14+s14)*(a1*b4-a4*b1)^2*R_C(1,b4*(b14+q14+s14)^2*(~
a1*b4-a4*b1)*(a2*b4-a4*b2)*(a4*b3-a3*b4)/(a1*b14*b4*(q14+s14)+a2*b4^2*(a3*b4-~
a4*b3)-a4*(a3*b2*b4^2-a4*b2*b3*b4+b1*b14*(q14+s14)))^2+1)-b1*(3*(a1*b4-a4*b1)~
*R_F((u12+u13)*(u12+u23),(u12+u13)*(u13+u23),(u12+u23)*(u13+u23))+b4*(a1*b2-a~
2*b1)*(a1*b3-a3*b1)*R_J((u12+u13)*(u12+u23),(u12+u13)*(u13+u23),(u12+u23)*(u1~
3+u23),u12*(u13+u23)+u13*u23+u14^2))*(a1*b14*b4*(q14+s14)+a2*b4^2*(a3*b4-a4*b~
3)-a4*(a3*b2*b4^2-a4*b2*b3*b4+b1*b14*(q14+s14))))/(3*(a1*b14*b4*(q14+s14)+a2*~
b4^2*(a3*b4-a4*b3)-a4*(a3*b2*b4^2-a4*b2*b3*b4+b1*b14*(q14+s14)))*(a1*b4-a4*b1~
)^2)