```python
from sympy.interactive import init_printing
init_printing()
from sympy import sqrt, sin, cos, pi, var, integrate, Symbol, S, Integral, symbols, solve
def rpoisson(f, r):
    return -(r**2*f.diff(r)).diff(r)/r**2/(4*pi)
def rpoisson_int(nr, r):
    # Boundary conditions:
    # r_c^2*v'(r_c) = -1       (1)
    # v(r_c) = 1/r_c           (2)
    first = ((-4*pi*nr)*r**2).integrate(r)
    first = first - first.subs(r, r_c) + (-1)     # Impose (1)
    second = (first/r**2).integrate(r)
    second = second - second.subs(r, r_c) + 1/r_c # Impose (2)
    return second
def check_gr_vr(gr, vr, r, r_c):
    # Some checks:
    # density is zero at r=r_c
    assert gr.subs(r, r_c).simplify() == 0
    # density is normalized to 1
    assert integrate(4*pi*r**2*gr, (r, 0, r_c)).simplify() == 1
    # vr is the potential corresponding to the density gr
    assert (gr-rpoisson(vr, r)).simplify() == 0
    # gr is the density corresponding to vr
    assert (vr-rpoisson_int(gr, r)).simplify() == 0
```


```python
var("r r_c")
gr = -21*(r-r_c)**3*(6*r**2+3*r*r_c+r_c**2)/(5*pi*r_c**8)
vr = (9*r**7 - 30*r**6*r_c + 28*r**5*r_c**2 - 14*r**2*r_c**5 + 12*r_c**7)/(5*r_c**8)
```


```python
check_gr_vr(gr, vr, r, r_c)
```


```python
Ig = integrate(4*pi*r**2*gr*(1/r-vr), (r, 0, r_c))
assert Ig == 10976/(17875*r_c)
Ig
```




    10976/(17875*r_c)




```python
Isph = integrate(4*pi*r**2*(1/r-vr), (r, 0, r_c))
assert Isph == 14*pi*r_c**2/75
Isph
```




    14*pi*r_c**2/75



# Example 1


```python
#gr = cos(pi*r/r_c/2) /( 8*r_c**3*(-8 + pi**2)/pi**2 )
gr = 1-(r/r_c)**2
gr = gr / integrate(4*pi*r**2*gr, (r, 0, r_c))
gr = gr.simplify()
gr
```




    15*(-r**2 + r_c**2)/(8*pi*r_c**5)




```python
vr = rpoisson_int(gr, r).simplify()
vr
```




    (3*r**4 - 10*r**2*r_c**2 + 15*r_c**4)/(8*r_c**5)




```python
check_gr_vr(gr, vr, r, r_c)
```


```python
Ig = integrate(4*pi*r**2*gr*(1/r-vr), (r, 0, r_c))
Ig
```




    25/(56*r_c)




```python
Isph = integrate(4*pi*r**2*(1/r-vr), (r, 0, r_c))
Isph
```




    2*pi*r_c**2/7




```python
print "Ig =", Ig
print "Isph =", Isph
print "v0 =", vr.limit(r, 0)
print
print "gr =", gr
```

    Ig = 25/(56*r_c)
    Isph = 2*pi*r_c**2/7
    v0 = 15/(8*r_c)
    
    gr = 15*(-r**2 + r_c**2)/(8*pi*r_c**5)


# Example 2


```python
n = 10
p = r**(n-1)
syms = symbols("a0:%d" % (n-1))
for i in range(n-1):
    p += syms[i]*r**i
f = (1-r)**n * p
eqs = [f.diff(r, i).subs(r, 0) for i in range(1, n)]
d = solve(eqs, syms)
f = f.subs(d)
gr = f.subs(r, r/r_c)
gr = gr / integrate(4*pi*r**2*gr, (r, 0, r_c))
gr = gr.simplify()
gr
```




    21*(r - r_c)**10*(48620*r**9 + 24310*r**8*r_c + 11440*r**7*r_c**2 + 5005*r**6*r_c**3 + 2002*r**5*r_c**4 + 715*r**4*r_c**5 + 220*r**3*r_c**6 + 55*r**2*r_c**7 + 10*r*r_c**8 + r_c**9)/(4*pi*r_c**22)




```python
vr = rpoisson_int(gr, r).simplify()
vr
```




    (-8840*r**21 + 92378*r**20*r_c - 432432*r**19*r_c**2 + 1191190*r**18*r_c**3 - 2130128*r**17*r_c**4 + 2567565*r**16*r_c**5 - 2089164*r**15*r_c**6 + 1108536*r**14*r_c**7 - 348840*r**13*r_c**8 + 49742*r**12*r_c**9 - 14*r**2*r_c**19 + 11*r_c**21)/(4*r_c**22)




```python
check_gr_vr(gr, vr, r, r_c)
```


```python
Ig = integrate(4*pi*r**2*gr*(1/r-vr), (r, 0, r_c))
Ig
```




    517399110137/(809268832200*r_c)




```python
Isph = integrate(4*pi*r**2*(1/r-vr), (r, 0, r_c))
Isph
```




    91*pi*r_c**2/690




```python
print "Ig =", Ig
print "Isph =", Isph
print "v0 =", vr.limit(r, 0)
print
print "gr =", gr
```

    Ig = 517399110137/(809268832200*r_c)
    Isph = 91*pi*r_c**2/690
    v0 = 11/(4*r_c)
    
    gr = 21*(r - r_c)**10*(48620*r**9 + 24310*r**8*r_c + 11440*r**7*r_c**2 + 5005*r**6*r_c**3 + 2002*r**5*r_c**4 + 715*r**4*r_c**5 + 220*r**3*r_c**6 + 55*r**2*r_c**7 + 10*r*r_c**8 + r_c**9)/(4*pi*r_c**22)



```python
print "Ig =", str(Ig).replace("r_c", "rc")
print "Isph =", str(Isph).replace("r_c", "rc")
print "v0 =", str(vr.limit(r, 0)).replace("r_c", "rc")
print
print "gr =", str(gr).replace("r_c", "rc")
```

    Ig = 517399110137/(809268832200*rc)
    Isph = 91*pi*rc**2/690
    v0 = 11/(4*rc)
    
    gr = 21*(r - rc)**10*(48620*r**9 + 24310*r**8*rc + 11440*r**7*rc**2 + 5005*r**6*rc**3 + 2002*r**5*rc**4 + 715*r**4*rc**5 + 220*r**3*rc**6 + 55*r**2*rc**7 + 10*r*rc**8 + rc**9)/(4*pi*rc**22)


# Example 3


```python
from scipy.integrate import cumtrapz
from scipy.special import erf
from numpy import linspace, trapz, pi, sin, empty, size, exp, sqrt

def integrate(R, f):
    return cumtrapz(f, R, initial=0)
    
def rpoisson_int(gr, R):
    rc = 1
    first = integrate(R, (-4*pi*gr)*R**2)
    first = first - first[-1] + (-1)
    second = integrate(R, first/R**2)
    second = second - second[-1] + 1/rc
    return second
```


```python
alpha = 37
C = 1/(pi**(3./2) * erf(sqrt(alpha)) / alpha**(3./2) - 2*pi*exp(-alpha)/alpha)

for i in range(3, 7):
    R = linspace(1e-5, 1, 10**i)
    gr = C * exp(-37*R**2)

    print "gr norm:", 1/trapz(4*pi*gr*R**2, R)

    vr = rpoisson_int(gr, R)

    Ig = trapz(4*pi*R**2*gr*(1/R-vr), R)
    Isph = trapz(4*pi*R**2*(1/R-vr), R)
    v0 = vr[0]

    print "Ig = %.14f_dp / rc" % Ig
    print "Isph = %.14f_dp * rc**2" % Isph
    print "v0 = %.14f_dp / rc" % v0
    print
```

    gr norm: 1.00000000085
    Ig = 2.01027837029584_dp / rc
    Isph = 0.08490520777895_dp * rc**2
    v0 = 6.86830681054508_dp / rc
    
    gr norm: 1.00000000001
    Ig = 2.01031976486788_dp / rc
    Isph = 0.08490788196539_dp * rc**2
    v0 = 6.86367256211640_dp / rc
    
    gr norm: 1.0
    Ig = 2.01032017805483_dp / rc
    Isph = 0.08490790865855_dp * rc**2
    v0 = 6.86366259004854_dp / rc
    
    gr norm: 1.0
    Ig = 2.01032018218137_dp / rc
    Isph = 0.08490790892400_dp * rc**2
    v0 = 6.86366256490661_dp / rc
    


491d5c7229b8355c8e7888b6641a45f18730bd85
ipynb
Jupyter Notebook
src/tests/nb/Radial Neutralizing Density.ipynb
certik/hfsolver
b4c50c1979fb7e468b1852b144ba756f5a51788d
BSD-2-Clause
2015-03-24T13:06:39.000Z
2022-03-29T00:14:02.000Z
src/tests/nb/Radial Neutralizing Density.ipynb
certik/hfsolver
b4c50c1979fb7e468b1852b144ba756f5a51788d
BSD-2-Clause
2015-03-25T04:59:43.000Z
2017-06-06T23:00:09.000Z
src/tests/nb/Radial Neutralizing Density.ipynb
certik/hfsolver
b4c50c1979fb7e468b1852b144ba756f5a51788d
BSD-2-Clause
2016-01-20T13:38:22.000Z
2020-11-24T15:35:43.000Z
Qwen/Qwen-72B
1. YES
2. YES