R = 1.
def f(ph):
    return 1./(2+sin(ph))

def integrando_a(ph,n):
    return f(ph)*cos(n*ph)
def a(n):
    return (1./(pi*R**n))*quad(integrando_a, -pi, pi, args=n)[0]
def integrando_b(ph,n):
    return f(ph)*sin(n*ph)
def b(n):
    return (1./(pi*R**n))*quad(integrando_b, -pi, pi, args=n)[0]

def S(nmax,ph,r):
    S = 0
    for n in range(nmax+1):
        S = S + (a(n)*cos(n*ph)+b(n)*sin(n*ph))*(r**n)
    return S

"""
Explanation: Ecuación de Laplace en coordenadas polares:
$$ \psi(r,\varphi)=c_{0}+c_{1}\ln(r)+\sum_{n=1}^{\infty}A_{n}\cos(n\varphi)+B_{n}\sin(n\varphi)$$
Si se buscan soluciones finitas dentro de un círculo $r< R$ y con $\psi(r=R,\varphi)=f(\varphi)$, entonces
$$ \psi(r,\varphi)=\sum_{n=0}^{\infty}A_{n}\cos(n\varphi)r^{n}+\sum_{n=1}^{\infty}B_{n}\sin(n\varphi)r^{n},$$
con
$$A_{n}:=\frac{1}{\pi R^{n}}\int_{-\pi}^{\pi}f(\varphi) \cos(n \varphi) d\varphi,\qquad n=0,1,2,\ldots$$
y
$$B_{n}:=\frac{1}{\pi R^{n}}\int_{-\pi}^{\pi}f(\varphi) \sin(n \varphi) d\varphi,\qquad n=1,2,\ldots$$
Consideremos el caso particular en que $f(\varphi):=\frac{1}{2+\sin(\varphi)}$ y calculemos numéricamente los coeficientes $A_{n}$ y $B_{n}$.
End of explanation
"""


n_ang = 120 #número de divisiones de la coordenada angular
n_radio = 50 #número de divisiones de la coordenada radial
radio = linspace(0.001, R, n_radio) #se define coordenada radial entre 0 y R
ang = linspace(0, 2*pi, n_ang) #se define coordenada angular entre 0 y 2 veces pi

r,p = meshgrid(radio,ang)

X = r*cos(p) #Define dominio en x
Y = r*sin(p) #Define dominio en y
n_max = 20 #número máximo de terminos de la serie truncada
Z = S(n_max, p, r) #valor del potencial

"""
Explanation: Ahora evaluamos algunos parámetros que usaremos para confeccionar el gráfico:
End of explanation
"""


fig = figure()
pcolormesh(X, Y, Z, shading='gouraud')
xlabel('$x$',fontsize=15)
ylabel('$y$',fontsize=15)
title('$\psi(x,y)$',fontsize=15)
colorbar()
axes().set_aspect('equal', 'datalim')

"""
Explanation: A continuación, un gráfico en colores en el dominio.
End of explanation
"""


