Bremsstrahlung Beam Calculations

from math import *



eps0 = 8.854187817e-12 #F/m
c    = 299792458.0 #m/s
u    = 1.660539040e-27 #Da (unified atomic mass unit)

def lorentz(vel):
    beta = vel / c
    gamma = (1.0-beta*beta) ** -0.5
    return gamma

def kinetic_energy(kg,vel):
    gamma = lorentz(vel)
    return kg * c*c * (gamma-1.0)
def kinetic_energy_invmass(energy,vel):
    gamma = lorentz(vel)
    return energy / ( c*c * ( gamma - 1.0 ) )

##def parallel_decel_1q( vel,decel ):
##    #https://en.wikipedia.org/wiki/Bremsstrahlung#Total_radiated_power
##    a = decel
##    q = 1.6021766208e-19
##    gamma = lorentz(vel)
##    power = q*q * a*a * (gamma**6) / ( 6.0*pi * eps0 * c*c*c )
##    return power

def parallel_decel_1q( vel0, accel ):
    #https://en.wikipedia.org/wiki/Bremsstrahlung#Total_radiated_power
    a = accel
    v_0 = vel0
    q = 1.6021766208e-19

##    def integrated(t):
##        u = (v_0 + a*t)/c
##        A = ( 10.0*u - 6*u*u*u ) / ( (u*u-1.0)**2 )
##        B = -3.0 * log(1.0-u)
##        C =  3.0 * log(1.0+u)
##        return (1.0/16.0)*( A + B + C )
##    energy = q*q * a / ( 6.0*pi * eps0 * c*c*c )
##    energy *= integrated(vel0/-a) - integrated(0.0)

    def get_power(vel):
        gamma = lorentz(vel)
        power = q*q * a*a * (gamma**6) / ( 6.0*pi * eps0 * c*c*c )
        return power
    def get_vel(t):
        return v_0 - a*t
    stopping_time = particle_speed / -accel
    N = 1000000
    dt = stopping_time / float(N)
    energy = 0.0
    for i in range(N):
        t = (i+0.5)*dt
        energy += get_power(get_vel(t))*dt

    return energy

def parallel_decel( kg, vel,accel ):
    num_atoms = kg / (12.011 * u) #Carbon
    num_protons = 6 * num_atoms
    num_electrons = num_protons
    return (num_protons+num_electrons)*parallel_decel_1q( vel,accel )
def parallel_stop( kg, vel,stop_dist ):
    #d = v v / (2 a), so:
    #2 a d = v*v
    #a = v*v / (2 d)
    accel = -vel*vel / (2.0*stop_dist)
    return parallel_decel( kg, vel,accel )



power_beam     = 1.0e9
particle_speed = 0.3 * c
particle_mass  = 1.0e-9 * 1.0e-3
stopping_dist  = 0.001



accel = -particle_speed*particle_speed / (2.0*stopping_dist)
mass_flow_rate = kinetic_energy_invmass(power_beam,particle_speed)
particles_rate = mass_flow_rate / particle_mass
kinetic_energy_per_particle = kinetic_energy(particle_mass,particle_speed)
bremsstrahlung_per_particle = parallel_stop(particle_mass, particle_speed,stopping_dist)
bremsstrahlung_beam_avg     = bremsstrahlung_per_particle * particles_rate

def fmt(value):
    s = "%g" % value
    if "e" in s:
        s = "%.4g" % value
        mantissa,exponent = s.split("e")
        s = mantissa + "⨯10"
        if exponent[0]=="-": s+="⁻"
        exponent = exponent[1:]
        started = False
        for c in exponent:
            c = "⁰¹²³⁴⁵⁶⁷⁸⁹"[int(c)]
            if c != "⁰":
                started = True
                s += c
            elif started:
                s += c
    return s
print("""
Power (beam):                        %s W
Particle speed:                      %s c
Particle mass:                       %s kg
Stopping distance:                   %s m

Mass flow rate:                      %s kg/s
Particles rate:                      %s particles/s
Particle stopping time:              %s s
Particle acceleration:               %s m/s²
Kinetic energy (per particle):       %s J
Bremsstrahlung (per particle):       %s J (= %s %% of KE)
Bremsstrahlung (beam temporal avg.): %s W"""%(
    fmt( power_beam ),
    fmt( particle_speed / c ),
    fmt( particle_mass ),
    fmt( stopping_dist ),

    fmt( mass_flow_rate ),
    fmt( particles_rate ),
    fmt( particle_speed / -accel ),
    fmt( accel ),
    fmt( kinetic_energy_per_particle ),
    fmt( bremsstrahlung_per_particle ), fmt( 100.0*bremsstrahlung_per_particle/kinetic_energy_per_particle ),
    fmt( bremsstrahlung_beam_avg )
))