What is it that's wrong with the following input in Maple?
Note that all variables are real. The sample worksheet is attached. 

Download Integration-Vec-Example.mw

I am trying to get Maple to understand, and evaluate the following integral analytically AND numerically.  
The original integrand to be evaluated is given by 

(w²  + m²)^(1/4) (p² + mu² )^(1/4) / ( (w² +a²)² (p² +a²) (q² + b2 ²)² ((p - q + w)² + b₁²) )

which is to be integrated simultaneously with respect to dw dp dq (where w, p, q are 3-d vectors).  
Limits of integration are (in cartesian co-ordinates) from -Infinity to + Infinity

All variables contained in the integrand are Real and >0.
a, m, mu, b₁, b₂ are constants and >0.  
dp = d³p, dq = d³q, dw = d³w
p²=p², w²=w², q²=q²

It is of course possible to help Maple along by noting that 
(p - q + w)² = p² + q² + w² - (2 p q Cos(\theta_pq) ) + (2 p w Cos(\theta_pw) ) - (2 q w Cos(\theta_wq) ) + b₁² 
where \theta_pq is the angle between vectors p and q, theta_pw is the angle between vectors p and w, etc. 

However, using the expanded form of (p - q + w)² would imply that the integral would now need to be written out in tems of 
dw_radial, dq_radial, dp_radial, d\theta_w, d\theta_q, d\theta_p, and d\phi_w, d\phi_q, d\phi_p.
Now, the integral in spherical co-ordinates would be given by the pseudocode below (note that now there are no longer any vector symbols, and the limits of integration have changed)

int ( Sin(\theta_w) d\theta_w Limits {0, pi}) int ( w² dw_radial Limits:{0, Infinity} )
int ( Sin(\theta_p) d\theta_p Limits {0, pi}) int ( p² dp_radial Limits:{0, Infinity} )
int ( Sin(\theta_q) d\theta_q Limits {0, pi}) int ( q² dq_radial Limits:{0, Infinity} )

where we have omitted the integration over \phi_w, \phi_p, \phi_q, assuming that there is no dependence on the \phi angles, resulting in all the \phi angles being integrated over, with the result being a constant multiplied to the rest of the remaining integrals. No additional information on the vectors is available, e.g. the angle between p and q, p and w, q and w, are not specified, the magnitude of the vectors is not specified. 

Any ideas if Maple can understand and execute the integral with p, q, w are input as vectors in the integrand directly, as shown earlier?
If it's easier for Maple to evaluate the integral in spherical co-ordinates, any suggestions on setting up the input, i.e. set up that the variables and constants are real, setup the limits of the theta and radial components etc. would also be appreciated. 

Is it possible to get the following integral involving vectors done in Maple 2018, but handle the singularity when p1=p3?

Note that the output of p3 and p1 are slightly different, probably because the earlier version of the worksheet was in Maple 2017? 

(Note: Earlier version - was due to help received from Dr. Edgardo Cheb-Terrab, using Maple 2017). 

with(Physics[Vectors]);
r_ := _i*x+_j*y+_k*z;
                    r_ := _i x + _j y + _k z
p__1_ := _i*`p__1x `+_j*`p__1y `+_k*`p__1z `;
           p__1_ := _i p__1x  + _j p__1y  + _k p__1z 
p__3_ := _i*`p__3x `+_j*`p__3y `+_k*`p__3z `;
           p__3_ := _i p__3x  + _j p__3y  + _k p__3z 
Expression to integrate wrt p3
'exp(i*(p__1_ . r_))/((p__3_ - p__1_)^2)'
             exp(I Physics:-Vectors:-.(p__1_, r_))
             -------------------------------------
                                                2 
              Physics:-Vectors:-+(p__3_, -p__1_)  
Delaying the evaluation above, if evaluation is allowed, the integrand is
exp(I*(p__1_ . r_))/(p__3_-p__1_)^2;
(exp(I (p__1x  x + p__1y  y + p__1z  z)))/((_i (-p__1x  + p__3x )

   + _j (-p__1y  + p__3y ) + _k (-p__1z  + p__3z ))^2)
T i is the imaginary unit
interface(imaginaryunit = I);
                               I
Because (4), has the value of i before being the imaginary unit, input the integrand again
'exp(i*(p__1_ . r_))/((p__3_ - p__1_)^2)'
             exp(I Physics:-Vectors:-.(p__1_, r_))
             -------------------------------------
                                                2 
              Physics:-Vectors:-+(p__3_, -p__1_)  


Now on the integration
Int(exp(I*(p__1_ . r_))/(p__3_-p__1_)^2, [`p__3x ` = -infinity .. infinity, `p__3y ` = -infinity .. infinity, `p__3y ` = -infinity .. infinity]);
  /infinity   /infinity   /infinity                            
 |           |           |                                     
 |           |           |          (exp(I (p__1x  x + p__1y  y
 |           |           |                                     
/-infinity  /-infinity  / -infinity                            

   + p__1z  z)))/((_i (-p__1x  + p__3x ) + _j (-p__1y  + p__3y )

   + _k (-p__1z  + p__3z ))^2) dp__3x  dp__3y  dp__3y 
value(Int(exp(I*(`p__1x `*x+`p__1y `*y+`p__1z `*z))/(_i*(-`p__1x `+`p__3x `)+_j*(-`p__1y `+`p__3y `)+_k*(-`p__1z `+`p__3z `))^2, [`p__3x ` = -infinity .. infinity, `p__3y ` = -infinity .. infinity, `p__3y ` = -infinity .. infinity]));
                                         /         /[          
 exp(I (p__1x  x + p__1y  y + p__1z  z)) |PIECEWISE|[infinity, 
                                         \         \[          

   Im(_i p__1x  + _j p__1y  + _k p__1z  - _k p__3z )    ]  
   ------------------------------------------------- = 0], 
                          _i                            ]  

                 \\         
   [0, otherwise]|| infinity
                 //         
In the above, there is the product of three Dirac delta functions, that can be represented as a single 3D Dirac delta
combine(exp(I*(`p__1x `*x+`p__1y `*y+`p__1z `*z))*piecewise(Im(_i*`p__1x `+_j*`p__1y `+_k*`p__1z `-_k*`p__3z `)/_i = 0, infinity, 0)*infinity);
                                         /         /[          
 exp(I (p__1x  x + p__1y  y + p__1z  z)) |PIECEWISE|[infinity, 
                                         \         \[          

   Im(_i p__1x  + _j p__1y  + _k p__1z  - _k p__3z )    ]  
   ------------------------------------------------- = 0], 
                          _i                            ]  

                 \\         
   [0, otherwise]|| infinity
                 //         

Is there a way in Maple to evaluate the following integral - where the "w" are set to go to zero to avoid singularities? B.t.w. Maple just returns the integral unevaluated with or without the "w" with the simple "int" command. 

int(ln(((p1-p3)^2+w^2)/(p1+p3)^(w^2+2))*ln(((p1-p4)^2+w^2)/((p1+p4)^2+w^2))*sin(p1)/(p1*p3*p4), p1 = 0 .. infinity);

or even

int(ln(((p1-p3)^2+w^2)/(p1+p3)^(w^2+2))*ln(((p1-p4)^2+w^2)/((p1+p4)^2+w^2))*sin(p1)/(p1*p3*p4), p1 = 0 .. 100, numeric);

How does one declare 3d vector in Maple, say p1, p3, and x, and integrate the following function from -Infinity to q, and q to + Infinity? Note that there is an iota in the exponential, and p1.x is the dot product of p1 and x. 

Integrate from -Infinity to q, q to infinity w.r.t. p1 and p3 

e^( - i p1 . x ) * 1/((p3 -p1)^2)

The integration should preferably be in spherical co-ordinates, but a demonstration in Cartesian co-ordinates would also be helpful.