#Calculation of flux in x direction (Qx) because of force at node j
#r = (x, y, z)
#r' is taken to be zero
#Variable declaration
var('x, y, z, fx, fy, fz, ux, pi, mu, dj')
#uxx is the ux because of force is x direction - fx
uxx=3/2/pi/mu*dj *x*z*x/(x^2+y^2+z^2)^(5/2) * fx
#uxy is the ux because of force is y direction - fy
uxy=3/2/pi/mu*dj *x*z*y/(x^2+y^2+z^2)^(5/2) * fy
#some assumptions needed by sage
assume(y^2 > 0, x^2 > 0, z^2 > 0 )
#Qxx is the contribtion to Qx because of force is x direction - fx
Qxx = integrate(uxx,z,0,infinity) # integrarion with respect to z
Qxx = integrate(Qxx,y,-infinity,infinity) # integrarion with respect to y
#Qxy is the contribtion to Qx because of force is y direction - fy
Qxy = integrate(uxy,z,0,infinity) # integrarion with respect to z
Qxy = integrate(Qxy,y,-infinity,infinity) # integrarion with respect to y
#summing up
Qx = Qxx+Qxy
#displaying Qx - the flux in x direction (Qx) because of force at node j
Qx