In Goldstein problem 3.5, We have to calculate the scattering angle, $\Phi(s)$, given as a function of input parameter, using Gauss-Legendre Quadrature. The full question is as follows;

Compute numerically $\Phi(s)$ and the differential cross section of $\sigma(\Phi)$ for the repulsive potential $$V=\frac{V_0}{1+r}$$ and for total energy $E=1.2V_0$. It is suggested that the 16-point Gauss-Legendre quadrature will give adequate accuracy. Does the scattering exhibit a rainbow?

So firstly starting with $\Phi(s)$ given as: $$\Phi(s)=\pi-2\displaystyle\int_{r_m}^\infty\frac{sdr}{r\sqrt{r^2 \left(1-\frac{V(r)}{E}\right)^2 - s^2}} \tag{1}$$

Where $s$, is the impact parameter. $r_m$ is the radius at closest approach and $r$ is the radius from the centre of the force. Defining some function $\rho(r)=\sqrt{1-\frac{r_m}{r}}$ and inserting into $(1)$ we get the scattering angle as

$$\Phi(s) = \pi-4s \displaystyle\int _0^1 \displaystyle\dfrac{\rho \ d \rho}{\sqrt{r_m^2(1-\frac{V}{E})^2-s^2(1-\rho^2)^2}} \tag{2}$$

the formulation of the scattering angle in equation $(2)$ allows us to use Gauss-Legendre quadrature mainly because of the definite limits of integration and the fact that the integrand is non singular in between those limits.

So I created a Gauss-Legendre quadrature algorithm in python

import sympy
from sympy import legendre, real_roots, nroots, diff
import numpy as np
import matplotlib.pyplot as plt

def weights_and_roots(p):
    calculates the roots of a P+1 degree Legendre polynomial then uses them to calculate the weights
    c = []
    L = legendre((p+1),x)
    dL = L.diff(x)
    nR = nroots(L)
    for rt in nR:
    return c, nR

def GLI(c,nR,f,lower_limit=-1,upper_limit=1):
    c = list of all the weights calculated in the previous function
    nR = this is a list of all the roots again calculated in the previous function
    f = the function being numerically integrated
    lower_limit = the lower limit of integration
    Upper_limit = the upper limit of integration
    The following sums the weights with the function (f) evaluated at the respective legendre root 
    idx = len(nR)
    S = 0
    for i in range(idx):
        S+= ((upper_limit-lower_limit)/2)*c[i]*f.subs(x,(((upper_limit-lower_limit)/2)*nR[i]+((upper_limit+lower_limit)/2)))
    return S

I tested it on some known functions and it worked well so I tried applying it to the integral in equation $(2)$ with $V$ as given in the question, I also replaced the $r_m^2$ with $r(1-\rho^2)^2$ thus the integral in equation $(2)$ is given as

$$\displaystyle\int _0^1 \displaystyle\dfrac{\rho \ d \rho}{\sqrt{(r(1-\rho^2)^2(1-\frac{1}{1.2(1+r)})^2-s^2(1-\rho^2)^2}} \tag{3}$$

My reasoning for replacing $r_m^2$ is because we don't know this value before hand so just sticking a value in felt wrong. I then ran equation $(3)$ through my Gauss-Legendre quadrature code. I set $r$ as $20$ and I set $V_0$ as $10$ as well just to give them values. The code is as follows

r = sympy.Symbol("r") #radius
s = sympy.Symbol("s") # impact parameter
x = sympy.Symbol("x") # I use x instead of rho 

a = 20 #setting the radius

f = x/((((r*(1-x**2))**2)*(1-(1/1.2*(1+r)))**2)-(s*(1-x**2))**2)**(1/2)
p =16 # number of points

c,nR = weights_and_roots(p)
S = GLI(c,nR,f,0,1)

phi = np.pi - 4*s*S #Calculate Phi
phi = phi.subs(r,a) # substitute a value for the radius r given by a

Plotting the potential as a function of radius gave the following Graph

enter image description here

However, when I plotted $\Phi(s)$ calculated using the script shown above as a function of impact parameter, $s$. I got the following graph

enter image description here

However, this doesn't look right to me at all! I would assume the scattering angle to fall away to zero as the impact parameter increases. However I am not sure where I have gone wrong? Furthermore, the question askes if the scattering exhibits a rainbow but I am not sure what sort of scattering profile is evidence of a rainbow scattering effect.

  • $\begingroup$ An easy way to troubleshoot an algorithm applied to a problem with an unknown answer is to apply a simpler, though perhaps less efficient, algorithm which you can be more confident in. Shouldn't take long at all to program a simple Riemann sum. It might not be as good an algorithm, but there are very few ways it can go wrong, and it should be plenty for just brute forcing an answer to see what you get here. $\endgroup$ May 11, 2021 at 17:26
  • $\begingroup$ Two things to mention: $1)$ Your equation $(1)$ seems wrong when comparing it with the corresponding formula in Goldstein's book (equation $(3.96)$ in my edition), with no square power for the $\left(1-\frac{V(r)}{E}\right)$ term. $2)$ It may be important to take into account that $r_m$ can be found from the equation for the energy with $\dot{r} = 0$. $r_m$ then becomes a function of $s$ as well as $V_0$ and $E$, although we know that $E = 1.2 V_0$. $\endgroup$
    – secavara
    May 11, 2021 at 19:10

1 Answer 1


To keep track of units, let me write the potential and the energy in the exercise as \begin{eqnarray} V &=& V_0 \frac{r_0}{r_0 + r} \\ E &=& \kappa V_0 \, , \end{eqnarray} where $r_0 = 1$ and $\kappa = 1.2$. By definition, the $z$ component of the angular momentum is $l = m r^2 \dot{\phi}$, but we may use the initial conditions, that is, the impact parameter $s$ and initial velocity $v_0$, to write \begin{equation} l = m s v_0 = m s \sqrt{\frac{2 E}{m}} = s \sqrt{2 E m} = s \sqrt{2 \kappa V_0 m} \, . \end{equation} The total energy is hence given by \begin{eqnarray} \kappa V_0 &=& E \\ &=& K + V \\ &=& \frac{1}{2} m \dot{r}^2 + \frac{1}{2} m r^2 \dot{\phi}^2 + V_0 \frac{r_0}{r_0 + r} \\ &=& \frac{1}{2} m \dot{r}^2 + \frac{1}{2} \frac{l^2}{m r^2} + V_0 \frac{r_0}{r_0 + r} \\ &=& \frac{1}{2} m \dot{r}^2 + \kappa V_0 \frac{s^2}{r^2} + V_0 \frac{r_0}{r_0 + r} \, . \end{eqnarray} When $r = r_m$ we have $\dot{r} = 0$. From the previous equation we then find an equation for $r_m$, \begin{equation} 0 = r_m^2 \left[ 1 - \frac{r_0}{\kappa \left( r_0 + r_m \right)} \right] - s^2 \, . \end{equation} With symbolic solving software you can find the solutions for this equation and single out the real and semi-positive root, for the values of $\kappa$, $r_0$ and $s$ of interest. We end up with an expression for the scattering angle $\Theta$, \begin{equation} \Theta = \pi - 2 \int_{r_m}^{\infty} \frac{s \, dr}{r \sqrt{r^2 \left[ 1 - \frac{V}{E} \right] - s^2}} \, , \end{equation} of the form \begin{equation} \Theta \left(s, \kappa, r_0 \right) = \pi - 2 \int_{r_m\left(s, \kappa, r_0 \right)}^{\infty} \frac{s \, dr}{r \sqrt{r^2 \left[ 1 - \frac{r_0}{\kappa \left( r_0 + r \right)} \right] - s^2}} \, . \end{equation} Notice that the denominator in the integrand vanishes when $r = r_m$. Mathematica seems to handle this well:

enter image description here

  • $\begingroup$ Awesome thanks! I didn't think about getting rm from the total energy. I'll try applying it to what I've done and see if I can match your graph output. $\endgroup$
    – seraphimk
    May 12, 2021 at 19:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.