Godot script for simple coulomb (Rutherford) scattering

scattering.gd

extends Control

var ec = 1.602176634/1e19 #elementary charge [C]
var amu = 1.66053906660/1e27 #atomic mass unit [kg]
var vprmtvty = 8.8541878128/1e12 #Vacuum permittivity [F/m]

var scale_m = 1e15 #scale multiplier
var time_m = 1e21 #time divider

var m0 = 1 * amu #Proton
var m1 = 197 * amu #Gold atom nucleus
var r0 = 1/scale_m #pixel size
var r1 = 7.56/1e15 #size approximation
var q0 = 1 * ec
var q1 = 79 * ec
var p0 = Vector2(-1000,300) / scale_m
var p1 = Vector2(800,300) / scale_m
var v0 = Vector2(1,0)
var v1 = Vector2(0,0)

var s = []

func energy_to_velocity(E, m):
	return sqrt(2*E*ec/m)

func _ready():
	v0 *= energy_to_velocity(13e6,m0)
	print("E'/E(proton): " + str((m1-m0)/(m0+m1)))

func _draw():
	draw_circle(p0 * scale_m, r0 * scale_m, Color.red)
	draw_circle(p1 * scale_m, r1 * scale_m, Color.white)
	draw_circle((m0*p0 + m1*p1) / (m0+m1) * scale_m, 2, Color.green)

func get_acceleration(m,p_1,p_2):
	var k = (q0*q1) / (4*PI*vprmtvty)
	var r = p_1 - p_2
	return k/m * r/pow(r.length(),3)

func _process(delta):
	delta /= time_m
	
	p0 += v0 * delta 
	p1 += v1 * delta
	
	update()
	
	v0 += get_acceleration(m0,p0,p1) * delta
	v1 += get_acceleration(m1,p1,p0) * delta
	
	s.append((p1-p0).length())
	print("s/r:" + str(s.min()/r1))
	print(str(m0*v0.length_squared()/2/ec/1e6)+" MeV")