fork download
copy 
  1. import numpy as np
  2. from scipy.integrate import quad
  3. import matplotlib.pyplot as plt
  4.  
  5. # Constants
  6. a = 1.0
  7. v = 1.0
  8. E = 0.5 * v**2
  9.  
  10. def V(r, a=1.0):
  11.     if r < a:
  12.         return 1/r + r**2 / (2*a**2) - 3/(2*a)
  13.     else:
  14.         return 0.0
  15.  
  16. def integrand(r, b, E):
  17.     U = V(r)
  18.     inside = 1 - b**2 / r**2 - 2 * U / E
  19.     if inside < 0:
  20.         return 0
  21.     return b / (r**2 * np.sqrt(inside))
  22.  
  23. def scattering_angle(b):
  24.     r_min = find_r_min(b)
  25.     if r_min is None:
  26.         return 0
  27.     integral, _ = quad(integrand, r_min, 10*a, args=(b, E))
  28.     return np.pi - 2 * integral
  29.  
  30. def find_r_min(b):
  31.     from scipy.optimize import brentq
  32.     def equation(r):
  33.         return 0.5 * b**2 / r**2 + V(r) - E
  34.     try:
  35.         return brentq(equation, 1e-5, a)
  36.     except ValueError:
  37.         return None
  38.  
  39. # Generate data
  40. b_vals = np.linspace(0.01, a, 100)
  41. theta_vals = [scattering_angle(b) for b in b_vals]
  42.  
  43. # Plot θ vs b
  44. plt.plot(np.degrees(theta_vals), b_vals)
  45. plt.xlabel("Scattering angle θ (degrees)")
  46. plt.ylabel("Impact parameter b")
  47. plt.title("Scattering angle vs impact parameter")
  48. plt.grid(True)
  49. plt.show()
  50.  
Success #stdin #stdout #stderr 1.79s 70180KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
prog:27: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
https://ideone.com/g7K9FV
language:
Python (cpython 2.7.16)
created:
10 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!