fork download 
  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. # your code goes here
Success #stdin #stdout #stderr 1.87s 70176KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
stderr 
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/6GB1Ae
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!