Solve integral in an annular domain in Python

Question

I am trying to solve a function in an annular domain that has a change of phase with respect to the angular direction of the annulus.

My attempt to solve it is the following:

import numpy as np
from scipy import integrate

def f(x0, y0):
    
    r = np.sqrt(x0**2 + y0**2)
    
    if r >= rIn and r <= rOut:
        theta = np.arctan(y0 / x0)
        R = np.sqrt((x - x0)**2 + (y - y0)**2 + z**2)
        integrand = (np.exp(-1j * (k*R + theta))) / R
        return integrand
    else:
        return 0


# Test

rIn = 0.5
rOut = 1.5
x = 1
y = 1
z = 1
k = 3.66

I = integrate.dblquad(f, -rOut, rOut, lambda x0: -rOut, lambda x0: rOut)

My problem is that I don't know how to get rid of the division by zero occuring when I evaluate theta.

Any help will be more than appreciated!

score 1 · Accepted Answer · answered May 17 '22 at 16:30

Use numpy.arctan2 instead, it will have problems only if both x and y are zero, in which case the angle is undetermined.

Also I see you that your integrand is complex, in this case you will probably have to handle real and imaginary part separately, as done here.

Solve integral in an annular domain in Python

1 Answers1