Theory detail : please visit my answer on physics.stackexchange.com (i used exactly that in this simulator):
http://physics.stackexchange.com/questi ... 996#305996
To compile use: [path]\fbc\win64-1.05.0\fbc -s console -t 128000 -w 3 -O 2 -v "Relativity_instantfield.bas"
Choose your destiny :) and wait for the trajectory to be computed step by step.
It takes a lot of time - it's relativistic so lot's of Lorentz Boosts in arbitrary directions = lots of calculus.
I've simplified some parts to make it go faster, so some parts of the code will not be easy to understand.
Keys ? Need help ? Just press key "h" during simulation;
ESC key exits simulation
Code: Select all
' Einstein's Special Relativity
' Relativistic Electrodynamics Simulator (free version, but instant field used instead of genius positions)
' Warning: Don't use with medical applications !
' Theory detail : please visit my answer on physics.stackexchange.com :
' http://physics.stackexchange.com/questions/304901/relativistic-2-body-problem-with-instantaneous-colomb-force/305996#305996
' My website http://alphax.info (Warning: opened from 8:00 GMT until 20:00 GMT, Monday to Sunday, closed in the other hours)
' (C) Mihai Barboi 2017 - open source.
' Tested with FreeBASIC Compiler - Version 1.05.0 (01-31-2016), built for win64 (64bit)
' (works with the 32bit fbc compiler too)
' To compile use: [path]\fbc\win64-1.05.0\fbc -s console -t 128000 -w 3 -O 2 -v "Relativity_instantfield.bas"
#include "fbgfx.bi"
Type cube4x4x4
As Double e(0 To 3, 0 To 3, 0 To 3)
End Type
Type matrix4x4
As Double e(0 To 3, 0 To 3)
End Type
Type vector4
As Double e(0 To 3)
End Type
Type quanta
As Double x, y, z
As Double vx, vy, vz
As Double q, m, t
As Double Ex, Ey, Ez
As Double Bx, By, Bz
As Double a, v
End Type
Declare Sub writeF_mu_nu(m As matrix4x4 Ptr , ex1 As Double, ey1 As Double, ez1 As Double, bx1 As Double, by1 As Double, bz1 As Double, c As Double)
Declare Sub writeBoostMatrix_covariant(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Declare Sub writeBoostMatrix_contravariant(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Declare Sub computeElectromagneticTensorLorentzBoost(m1 As matrix4x4 Ptr, m2 As matrix4x4 Ptr, F As matrix4x4 Ptr, res As matrix4x4 Ptr)
Declare Sub writeBoostMatrix_covariant_reverse(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Declare Sub writeBoostMatrix_contravariant_reverse(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Declare Sub electrodynamics_instantfield
Declare Sub show_help
'helper
Declare sub pause1 overload(a As string,x as integer)
Declare sub pause1 overload(a as string)
Declare sub pause1 overload()
Declare sub pause1 overload(x as integer)
Dim Shared As Double pi = 3.14159265359
Dim Shared As Integer maxx, maxy
#Define to_rads(x) x * pi / 180
#Define logit(x) CInt((x * 10 ^ Abs((Log(rref) / Log(10)))))
Screen 19, 32
ScreenInfo maxx, maxy
electrodynamics_instantfield()
End
sub pause1 overload(a As string,x as integer)
print a:while inkey="":sleep x:wend
end sub
sub pause1 overload(a as string)
print a:while inkey="":sleep 333:wend
end sub
sub pause1 overload()
while inkey="":sleep 333:wend
end sub
sub pause1 overload(x as integer)
while inkey="":sleep x:wend
end Sub
'F[a]{b}
Sub writeF_mu_nu(m As matrix4x4 Ptr , ex1 As Double, ey1 As Double, ez1 As Double, bx1 As Double, by1 As Double, bz1 As Double, c As Double)
' F[a]{b}=| 0 -Ex/c -Ey/c -Ez/c |
' | -Ex/c 0 -Bz By |
' | -Ey/c Bz 0 -Bx |
' | -Ez/c -By Bx 0 |
m->e(0,0) = 0
m->e(0,1) = -ex1/c
m->e(0,2) = -ey1/c
m->e(0,3) = -ez1/c
m->e(1,0) = -ex1/c
m->e(1,1) = 0
m->e(1,2) = -bz1
m->e(1,3) = by1
m->e(2,0) = -ey1/c
m->e(2,1) = bz1
m->e(2,2) = 0
m->e(2,3) = -bx1
m->e(3,0) = -ez1/c
m->e(3,1) = -by1
m->e(3,2) = bx1
m->e(3,3) = 0
End Sub
Sub writeBoostMatrix_covariant(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Dim As Double bx, by, bz, b2, g
g = 1 / from->a
bx = from->vx / c
by = from->vy / c
bz = from->vz / c
b2 = (from->v*from->v) / (c*c)
m->e(0,0) = g
m->e(0,1) = bx*g
m->e(0,2) = by*g
m->e(0,3) = bz*g
m->e(1,0) = bx*g
m->e(1,1) = 1 + (g-1)*bx*bx/b2
m->e(1,2) = (g-1)*bx*by/b2
m->e(1,3) = (g-1)*bx*bz/b2
m->e(2,0) = by*g
m->e(2,1) = (g-1)*by*bx/b2
m->e(2,2) = 1 + (g-1)*by*by/b2
m->e(2,3) = (g-1)*by*bz/b2
m->e(3,0) = bz*g
m->e(3,1) = (g-1)*bz*bx/b2
m->e(3,2) = (g-1)*bz*by/b2
m->e(3,3) = 1 + (g-1)*bz*bz/b2
End Sub
Sub writeBoostMatrix_contravariant(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Dim As Double bx, by, bz, b2, g
g = 1 / from->a
bx = from->vx / c
by = from->vy / c
bz = from->vz / c
b2 = (from->v*from->v) / (c*c)
m->e(0,0) = g
m->e(0,1) = -bx*g
m->e(0,2) = -by*g
m->e(0,3) = -bz*g
m->e(1,0) = -bx*g
m->e(1,1) = 1 + (g-1)*bx*bx/b2
m->e(1,2) = (g-1)*bx*by/b2
m->e(1,3) = (g-1)*bx*bz/b2
m->e(2,0) = -by*g
m->e(2,1) = (g-1)*by*bx/b2
m->e(2,2) = 1 + (g-1)*by*by/b2
m->e(2,3) = (g-1)*by*bz/b2
m->e(3,0) = -bz*g
m->e(3,1) = (g-1)*bz*bx/b2
m->e(3,2) = (g-1)*bz*by/b2
m->e(3,3) = 1 + (g-1)*bz*bz/b2
End Sub
Sub computeElectromagneticTensorLorentzBoost(m1 As matrix4x4 Ptr, m2 As matrix4x4 Ptr, F As matrix4x4 Ptr, res As matrix4x4 Ptr)
' F[m]{n} = L{m}[a]L[n]{b}F[a]{b}
' A{nm} : n x m
' B{mp} : m x p
' (AB){ij} = sum{k=1}[m] A{ik} B{kj}
Dim As Integer m,n,a,b
For m = 0 To 3
For n = 0 To 3
res->e(m,n) = 0
For a = 0 To 3
For b = 0 To 3
res->e(m,n) += m1->e(m,a) * m2->e(n,b) * F->e(a,b)
Next b
Next a
next n
next m
End Sub
Sub writeBoostMatrix_covariant_reverse(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Dim As Double bx, by, bz, b2, g
' (AB){ij} = sum{k=1}[m] A{ik} B{kj} ' <= good one
' L{m}[n] = ita{rn} L[r]{s} ita[ms] = ita{rn} ita[ms] L[r]{s}
' L{m}[n} = ita[mr] ita{ns} L[r]{s} ' <= good one
g = 1 / from->a
bx = -from->vx / c
by = -from->vy / c
bz = -from->vz / c
b2 = (from->v*from->v) / (c*c)
m->e(0,0) = g
m->e(0,1) = bx*g
m->e(0,2) = by*g
m->e(0,3) = bz*g
m->e(1,0) = bx*g
m->e(1,1) = 1 + (g-1)*bx*bx/b2
m->e(1,2) = (g-1)*bx*by/b2
m->e(1,3) = (g-1)*bx*bz/b2
m->e(2,0) = by*g
m->e(2,1) = (g-1)*by*bx/b2
m->e(2,2) = 1 + (g-1)*by*by/b2
m->e(2,3) = (g-1)*by*bz/b2
m->e(3,0) = bz*g
m->e(3,1) = (g-1)*bz*bx/b2
m->e(3,2) = (g-1)*bz*by/b2
m->e(3,3) = 1 + (g-1)*bz*bz/b2
End Sub
Sub writeBoostMatrix_contravariant_reverse(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Dim As Double bx, by, bz, b2, g
g = 1 / from->a
bx = -from->vx / c
by = -from->vy / c
bz = -from->vz / c
b2 = (from->v*from->v) / (c*c)
m->e(0,0) = g
m->e(0,1) = -bx*g
m->e(0,2) = -by*g
m->e(0,3) = -bz*g
m->e(1,0) = -bx*g
m->e(1,1) = 1 + (g-1)*bx*bx/b2
m->e(1,2) = (g-1)*bx*by/b2
m->e(1,3) = (g-1)*bx*bz/b2
m->e(2,0) = -by*g
m->e(2,1) = (g-1)*by*bx/b2
m->e(2,2) = 1 + (g-1)*by*by/b2
m->e(2,3) = (g-1)*by*bz/b2
m->e(3,0) = -bz*g
m->e(3,1) = (g-1)*bz*bx/b2
m->e(3,2) = (g-1)*bz*by/b2
m->e(3,3) = 1 + (g-1)*bz*bz/b2
End Sub
Sub electrodynamics_instantfield
Dim As Double c = 299792458 ' speed of light
Dim As Double c2 = c*c ' speed of light to power 2
Dim As Double c4 = c2*c2 ' speed of ligth to power 4
Dim As Double m0 = 1e-27
Dim As Double h = 6.6e-36 ' Plank constant
Dim As Double e = 1.6e-19 ' electric charge of electron
Dim As Double m_elec = 9.1e-31 ' mass of electron
Dim As Double m_prot = 1.627e-27 ' mass of proton
Dim As Double miu = 4*pi*1.0e-7 ' magnetic permeability
Dim As Double epsi = 8.854e-12 ' electric permitivity
Dim As Double grav = 6.673e-11 ' gravitational constant of Newton
Dim As Double pi2 = pi * 2 ' 2 pi
Dim As Double miudiv4pi = miu / (4*pi) ' constant
Dim As Double epsi_ = 4*pi*epsi ' constat
Dim As Double epsi_c2 = (epsi_) * c2 ' constant
Dim As Double var_theta, var_phi, sin1, sin2, x, y, z, phi2, theta2, real_theta, real_phi
Dim As Double angle = 25, angle2 = 0, rref, minz, maxz, tip,cos_opt, sin_opt, logit_z, logit_y, logit_x
Dim As Integer co, xx, yy, xx2, yy2, wx = 0, wy = 0,x_opt,y_opt,tip_base
Dim As double sinangle, cosangle, sinangle2, cosangle2,ttp1,ttp2,maxxshr,maxyshr ' 3d graphics translation
Dim As UInteger color1 = RGB(155,255,255)
Dim As ULongInt num = 0, steps = 0, skip = 0
Dim As Double u
Dim As Double q
Dim As Double r
Dim As Double phi, theta
Dim As Double t
Dim As Double e_calc
Dim As Double timeQuanta
Dim As Integer cont = 0
Dim As String a
Dim As Double startTime, universalTime, rmin
Dim As Integer currentTime, np, maxTime, do_graph = 0,skipper = 5, skippy = 0
Dim As matrix4x4 ma1, ma2, Fma, mres
Dim As vector4 vect1, vect2, vect3
Dim As cube4x4x4 cube1, cube2, cube3
Dim As quanta Ptr quants
dim As quanta Ptr s
dim As quanta Ptr en
dim As quanta Ptr from
Dim As Integer howMany, size, i, k, j, ii, kk, jj
Dim As Double bx, by, bz
Dim As Double tx1, ty1, tz1, flx, fly, flz, d, f, fx, fy, fz, ff
Dim As Double rx1, ry1, rz1, Rmax, ovx, ovy, ovz, oa
Dim As Double ga, b2, gx, gy, gz, gb, sx, sy, sz, ex1, ey1, ez1, bx1, by1, bz1
Dim As Double AA, BB, CC, DD
Dim As Double dtau, temp, dt, temp2
Dim As Double r2, r3
angle = 180
angle2 = 180
rref = 0.5e-14 ' meters
sinangle=Sin(to_rads(90-angle))
cosangle=cos(to_rads(90-angle))
sinangle2=sin(to_rads(90-angle2))
cosangle2=cos(to_rads(90-angle2))
x_opt = CInt(maxx Shr 1)+wx
y_opt = CInt(maxy Shr 1)+wy
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
howMany = 1000 ' 30000000
quants = CPtr(quanta Ptr, Allocate(howMany * SizeOf(quanta)))
For i = 0 To howMany-1
quants[i].x = 0
quants[i].y = 0
quants[i].z = 0
quants[i].vx = 0
quants[i].vy = 0
quants[i].vz = 0
quants[i].v = Sqr(quants[i].vx*quants[i].vx + quants[i].vy*quants[i].vy + quants[i].vz + quants[i].vz)
quants[i].q = 0
quants[i].m = 0
quants[i].t = 0
quants[i].Ex = 0
quants[i].Ey = 0
quants[i].Ez = 0
quants[i].Bx = 0
quants[i].By = 0
quants[i].Bz = 0
quants[i].a = Sqr(1 - quants[i].v*quants[i].v/c2)
Next i
size = howMany * SizeOf(quanta) / (1024*1024)
Dim As Integer iu, num_particles
show_help()
Dim As Integer choice = 0
Print
Print
Print "Select which predefined system of particles do you want to simulate :"
Print "choice = 0 ' electron[0] , positron[1]"
Print "choice = 1 ' electron[0] , proton[1]"
Print "choice = 2 ' electron[0] , proton[1], electron[2] , proton[3]"
Print "choice = 3 ' electron[0] , 2*[q, mass] proton[1], electron[2]"
Print "choice = 4 ' random electron, random proton"
Print "choice = 5 ' collision of electron[0], electron[1]"
Input "Config choice (0..5):",choice
If choice = 0 Then
rref = 0.15e-13
angle = 206.6
angle2 = 180
wx = +200
wy = +200
sinangle = Sin(to_rads(90-angle))
cosangle = Cos(to_rads(90-angle))
sinangle2 = Sin(to_rads(90-angle2))
cosangle2 = Cos(to_rads(90-angle2))
x_opt = CInt(maxx Shr 1)+wx
y_opt = CInt(maxy Shr 1)+wy
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
num_particles = 2
quants[0].q = -e
quants[0].m = m_elec
quants[0].t = 0
quants[0].x = 1e-12 ' good value 1e-12
quants[0].z = 0.5e-12
quants[0].y = 0 ' good value 1e-12
quants[0].vx = 0
quants[0].vy = 0 '7e6
quants[0].vz = 7e6 ' good value 7e6
quants[0].v = Sqr(quants[0].vx*quants[0].vx + quants[0].vy*quants[0].vy + quants[0].vz*quants[0].vz)
quants[0].a = Sqr(1-quants[0].v*quants[0].v/c2)
quants[1].q = e
quants[1].m = m_elec
quants[1].vx = 5e6
quants[1].vy = 0
quants[1].v = Sqr(quants[1].vx*quants[1].vx + quants[1].vy*quants[1].vy + quants[1].vz*quants[1].vz)
quants[1].a = Sqr(1-quants[1].v*quants[1].v/c2)
ElseIf choice = 1 Then
angle = 206.6
angle2 = 180
sinangle = Sin(to_rads(90-angle))
cosangle = Cos(to_rads(90-angle))
sinangle2 = Sin(to_rads(90-angle2))
cosangle2 = Cos(to_rads(90-angle2))
x_opt = CInt(maxx Shr 1)+wx
y_opt = CInt(maxy Shr 1)+wy
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
num_particles = 2
quants[0].q = -e
quants[0].m = m_elec
quants[0].t = 0
quants[0].x = 1e-12 ' good value 1e-12
quants[0].z = 0.5e-12
quants[0].y = 0 ' good value 1e-12
quants[0].vx = 0
quants[0].vy = 0 '7e6
quants[0].vz = 7e6 ' good value 7e6
quants[0].v = Sqr(quants[0].vx*quants[0].vx + quants[0].vy*quants[0].vy + quants[0].vz*quants[0].vz)
quants[0].a = Sqr(1-quants[0].v*quants[0].v/c2)
quants[1].q = e
quants[1].m = m_prot
quants[1].vx = 0
quants[1].vy = 0
quants[1].v = Sqr(quants[1].vx*quants[1].vx + quants[1].vy*quants[1].vy + quants[1].vz*quants[1].vz)
quants[1].a = Sqr(1-quants[1].v*quants[1].v/c2)
ElseIf choice = 2 Then
angle = 206.6
angle2 = 180
sinangle = Sin(to_rads(90-angle))
cosangle = Cos(to_rads(90-angle))
sinangle2 = Sin(to_rads(90-angle2))
cosangle2 = Cos(to_rads(90-angle2))
x_opt = CInt(maxx Shr 1)+wx
y_opt = CInt(maxy Shr 1)+wy
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
num_particles = 4
quants[0].q = -e
quants[0].m = m_elec
quants[0].t = 0
quants[0].x = 1e-12 ' good value 1e-12
quants[0].z = 0.5e-12
quants[0].y = 0 ' good value 1e-12
quants[0].vx = 0
quants[0].vy = 0 '7e6
quants[0].vz = 7e6 ' good value 7e6
quants[0].v = Sqr(quants[0].vx*quants[0].vx + quants[0].vy*quants[0].vy + quants[0].vz*quants[0].vz)
quants[0].a = Sqr(1-quants[0].v*quants[0].v/c2)
quants[1].q = e
quants[1].m = m_prot
quants[1].vx = 0
quants[1].vy = 0
quants[1].v = Sqr(quants[1].vx*quants[1].vx + quants[1].vy*quants[1].vy + quants[1].vz*quants[1].vz)
quants[1].a = Sqr(1-quants[1].v*quants[1].v/c2)
quants[2].q = -e
quants[2].m = m_elec
quants[2].t = 0
quants[2].x = 1e-12 ' good value 1e-12
quants[2].z = -0.5e-12
quants[2].vx = 0
quants[2].vy = 0 '7e6
quants[2].vz = 7e6 ' good value 7e6
quants[2].v = Sqr(quants[2].vx*quants[2].vx + quants[2].vy*quants[2].vy + quants[2].vz*quants[2].vz)
quants[2].a = Sqr(1-quants[2].v*quants[2].v/c2)
quants[3].z = -0.5e-12
quants[3].q = e
quants[3].m = m_prot
quants[3].vx = 0
quants[3].vy = 0
quants[3].v = Sqr(quants[3].vx*quants[3].vx + quants[3].vy*quants[3].vy + quants[3].vz*quants[3].vz)
quants[3].a = Sqr(1-quants[3].v*quants[3].v/c2)
ElseIf choice = 3 Then
angle = 206.6
angle2 = 180
sinangle = Sin(to_rads(90-angle))
cosangle = Cos(to_rads(90-angle))
sinangle2 = Sin(to_rads(90-angle2))
cosangle2 = Cos(to_rads(90-angle2))
x_opt = CInt(maxx Shr 1)+wx
y_opt = CInt(maxy Shr 1)+wy
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
num_particles = 3
quants[0].q = -e
quants[0].m = m_elec
quants[0].t = 0
quants[0].x = 1e-12 ' good value 1e-12
quants[0].z = 0.5e-12
quants[0].y = 0 ' good value 1e-12
quants[0].vx = 0
quants[0].vy = 0 '7e6
quants[0].vz = 0.1e6 ' good value 7e6
quants[0].v = Sqr(quants[0].vx*quants[0].vx + quants[0].vy*quants[0].vy + quants[0].vz*quants[0].vz)
quants[0].a = Sqr(1-quants[0].v*quants[0].v/c2)
quants[1].q = 2*e
quants[1].m = 2*m_prot + 2*m_prot '+2neutrons
quants[1].vx = 0
quants[1].vy = 0
quants[1].v = Sqr(quants[1].vx*quants[1].vx + quants[1].vy*quants[1].vy + quants[1].vz*quants[1].vz)
quants[1].a = Sqr(1-quants[1].v*quants[1].v/c2)
quants[2].q = -e
quants[2].m = m_elec
quants[2].t = 0
quants[2].x = 1e-12 ' good value 1e-12
quants[2].z = -0.5e-12
quants[2].vx = 0
quants[2].vy = 0 '7e6
quants[2].vz = 1e6 ' good value 7e6
quants[2].v = Sqr(quants[2].vx*quants[2].vx + quants[2].vy*quants[2].vy + quants[2].vz*quants[2].vz)
quants[2].a = Sqr(1-quants[2].v*quants[2].v/c2)
ElseIf choice = 4 Then
Randomize Timer
num_particles = 2
quants[0].q = -e
quants[0].m = m_elec
quants[0].t = 0
quants[0].x = 1e-14 + 1.0e-12*rnd ' good value 1e-12
quants[0].z = 1e-14 + 1.0e-12*Rnd
quants[0].y = 1e-14 + 1.0e-12*Rnd
quants[0].vx = 3e6*Rnd
quants[0].vy = 3e6*Rnd '7e6
quants[0].vz = 3e6*rnd ' good value 7e6
quants[0].v = Sqr(quants[0].vx*quants[0].vx + quants[0].vy*quants[0].vy + quants[0].vz*quants[0].vz)
quants[0].a = Sqr(1-quants[0].v*quants[0].v/c2)
quants[1].q = e
quants[1].m = m_prot
quants[1].vx = 3e4*Rnd
quants[1].vy = 3e4*Rnd
quants[1].vz = 3e4*Rnd
quants[1].v = Sqr(quants[1].vx*quants[1].vx + quants[1].vy*quants[1].vy + quants[1].vz*quants[1].vz)
quants[1].a = Sqr(1-quants[1].v*quants[1].v/c2)
ElseIf choice = 5 Then
rref = 0.15e-13
angle = 206.6
angle2 = 180
wx = 0'-200
wy = 0'+200
sinangle = Sin(to_rads(90-angle))
cosangle = Cos(to_rads(90-angle))
sinangle2 = Sin(to_rads(90-angle2))
cosangle2 = Cos(to_rads(90-angle2))
x_opt = CInt(maxx Shr 1)+wx
y_opt = CInt(maxy Shr 1)+wy
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
num_particles = 2
quants[0].q = e
quants[0].m = m_elec+m_prot
quants[0].t = 0
quants[0].x = -1e-12
quants[0].z = 4e-12
quants[0].y = 0
quants[0].vx = 3e7
quants[0].vy = 0
quants[0].vz = -1e8
quants[0].v = Sqr(quants[0].vx*quants[0].vx + quants[0].vy*quants[0].vy + quants[0].vz*quants[0].vz)
quants[0].a = Sqr(1-quants[0].v*quants[0].v/c2)
quants[1].q = e
quants[1].m = m_elec
quants[1].x = -1e-12
quants[1].z = -4e-12
quants[1].vx = 3e7
quants[1].vy = 0
quants[1].vz = 1e8
quants[1].v = Sqr(quants[1].vx*quants[1].vx + quants[1].vy*quants[1].vy + quants[1].vz*quants[1].vz)
quants[1].a = Sqr(1-quants[1].v*quants[1].v/c2)
Else
Print "undefined choice ";choice
pause1
End
End If
Cls
Dim As Double oi = 0
timeQuanta = 2.0e-6/c
maxTime = 100000
Do
' computations
If (skip = 0) Then
currentTime = 0
Do
i = 0
rmin = 10000
i = 0
While i < num_particles
s = @quants[i]
s->Ex = 0 : s->Ey = 0 : s->Ez = 0
s->Bx = 0 : s->By = 0 : s->Bz = 0
ii = 0
While ii < num_particles
If (ii <> i) Then
from = @quants[ii]
tx1 = s->x - from->x
ty1 = s->y - from->y
tz1 = s->z - from->z
r2 = (tx1 * tx1) + (ty1 * ty1) + (tz1 * tz1)
d = Sqr(r2)
Rmax = Abs(s->q * from->q) / (epsi_c2 * ((s->m/s->a) + (from->m/from->a)))
If d > Rmax Then
If rmin > d Then rmin = d
r3 = r2 * d
f = from->q / (epsi_ * r3) ' r^3 = r2 * d = d * d * d
If (oi=0) Then oi = Abs( s->q * f * r2)
ex1 = f * tx1
ey1 = f * ty1
ez1 = f * tz1
f = miudiv4pi * from->q / r3
bx1 = f * (from->vy * tz1 - from->vz * ty1)
by1 = f * (from->vz * tx1 - from->vx * tz1) '+ 0.2e9
bz1 = f * (from->vx * ty1 - from->vy * tx1)
If (from->v > 0) Then
writeF_mu_nu(@Fma, ex1, ey1, ez1, bx1, by1, bz1, c)
writeBoostMatrix_covariant(from, @ma1, c)
writeBoostMatrix_contravariant(from, @ma2, c)
computeElectromagneticTensorLorentzBoost(@ma1, @ma2, @Fma, @mres)
ex1 = -mres.e(0,1)*c
ey1 = -mres.e(0,2)*c
ez1 = -mres.e(0,3)*c
bx1 = mres.e(3,2)
by1 = mres.e(1,3)
bz1 = mres.e(2,1)
EndIf
/' If (s->v > 0) Then
writeF_mu_nu(@Fma, ex1, ey1, ez1, bx1, by1, bz1, c)
writeBoostMatrix_covariant_reverse(s, @ma1, c)
writeBoostMatrix_contravariant_reverse(s, @ma2, c)
computeElectromagneticTensorLorentzBoost(@ma1, @ma2, @Fma, @mres)
ex1 = -mres.e(0,1)*c
ey1 = -mres.e(0,2)*c
ez1 = -mres.e(0,3)*c
bx1 = mres.e(3,2)
by1 = mres.e(1,3)
bz1 = mres.e(2,1)
EndIf '/
s->Ex += ex1
s->Ey += ey1
s->Ez += ez1
s->Bx += bx1
s->By += by1
s->Bz += bz1
Else
rmin = rmax
End If
End If
ii+=1
Wend
i += 1
Wend
If (rmin <= 0) Then
Locate 1,1
Print "ERROR ";rmin; " "
pause1
EndIf
dt = rmin * timeQuanta
i = 0
While i < num_particles
s = @quants[i]
'dtau = dt * s->a ' dtau = tau1 - tau0 = t1*a1 - t0*a0 = dt
temp = (s->q) * dt ' /s->a * dtau
temp2 = s->m / s->a
DD = temp/c * (s->Ex * s->vx + s->Ey * s->vy + s->Ez * s->vz) + temp2 * c
AA = temp * (s->Ex + s->Bz * s->vy - s->By * s->vz) + temp2 * s->vx
BB = temp * (s->Ey - s->Bz * s->vx + s->Bx * s->vz) + temp2 * s->vy
CC = temp * (s->Ez + s->By * s->vx - s->Bx * s->vy) + temp2 * s->vz
oa = s->a
ovx = s->vx
ovy = s->vy
ovz = s->vz
s->t += dt
temp = c / DD
s->vx = AA * temp
s->vy = BB * temp
s->vz = CC * temp
's->a = s->m * c / DD
s->v = Sqr(s->vx * s->vx + s->vy * s->vy + s->vz * s->vz)
temp = s->v / c
s->a = Sqr(1 - temp * temp)
/'
temp = dt'*s->a
s->x = s->x + ( s->vx/s->a + ovx/oa )/2 * temp
s->y = s->y + ( s->vy/s->a + ovy/oa )/2 * temp
s->z = s->z + ( s->vz/s->a + ovz/oa )/2 * temp
'/
' compute new position relative to our "observer"
temp = dt
s->x = s->x + ( s->vx + ovx )/2 * temp
s->y = s->y + ( s->vy + ovy )/2 * temp
s->z = s->z + ( s->vz + ovz )/2 * temp
i+=1
Wend
currentTime+=1
Loop Until currentTime >= maxTime
steps += currentTime
Else
Sleep 100
End If
' -- user input ---
a = InKey
Select Case a
Case "h"
show_help
Print "press any key to continue."
pause1
Case Chr(9)
skip = IIf(skip = 0, 1, 0)
Case "c"
Cls
Case "i"
Locate 2,1
Input "i=",i
If (i < 0) Then i =0
If (i >= num_particles) Then i = num_particles - 1
Print "m ";quants[i].m;" q ";quants[i].q; " steps ";steps
Print "x ";quants[i].x;" ";quants[i].y;" ";quants[i].z;" t ";quants[i].t;" "
Print "a ";quants[i].a;" V ";quants[i].v; " ct ";c*quants[i].t;" "
Print "v ";quants[i].vx;" ";quants[i].vy;" ";quants[i].vz;" "
Print "E ";quants[i].Ex;" ";quants[i].Ey;" ";quants[i].Ez;" "
Print "B ";quants[i].Bx;" ";quants[i].By;" ";quants[i].Bz;" "
Print "K ";quants[i].m*c2*((1/quants[i].a) - 1)/h;" Hz => " ;Log(quants[i].m*c2*((1/quants[i].a) - 1)/h)/Log(10);" (log)Hz"
Print
pause1
Case "v"
do_graph = IIf(do_graph, 0, 1)
Case "r"
Print "resolution ";maxTime;" (=num computations before drawing)"
Print "default = 100000 ; decrease it for better graphics; increase it for succint graphics "
Input "new resolution :",maxTime
If (maxTime < 1) Then maxTime = 1
If (maxTime > 10000000) Then maxTime = 10000000
Case "1"
angle-=2
sinangle=Sin(to_rads(90-angle))
cosangle=cos(to_rads(90-angle))
sinangle2=sin(to_rads(90-angle2))
cosangle2=cos(to_rads(90-angle2))
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
Locate 1,1
Print angle;" ";angle2; " ";
Case "2"
angle+=2
sinangle=Sin(to_rads(90-angle))
cosangle=cos(to_rads(90-angle))
sinangle2=sin(to_rads(90-angle2))
cosangle2=cos(to_rads(90-angle2))
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
Locate 1,1
Print angle;" ";angle2; " ";
Case "3"
angle2-=2
sinangle=Sin(to_rads(90-angle))
cosangle=cos(to_rads(90-angle))
sinangle2=sin(to_rads(90-angle2))
cosangle2=cos(to_rads(90-angle2))
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
Locate 1,1
Print angle;" ";angle2; " ";
Case "4"
angle2+=2
sinangle=Sin(to_rads(90-angle))
cosangle=cos(to_rads(90-angle))
sinangle2=sin(to_rads(90-angle2))
cosangle2=cos(to_rads(90-angle2))
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
Locate 1,1
Print angle;" ";angle2; " ";
Case Chr(255)+"H"'up
wy-=30
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
Case Chr(255)+"P"'down
wy+=30
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
Case Chr(255)+"K"'left
wx-=30
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
Case Chr(255)+"M"'right
wx+=30
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
Case "-"
rref=rref*2.222222
'wx=wx\10
'wy=wy\10
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
Locate 1,1
Print "rref ";rref;" ";
Case "+"
rref=rref/2.222222
ttp1=cosangle*cosangle2
ttp2=sinangle*cosangle2
maxxshr=maxx Shr 1
maxyshr=maxy Shr 1
maxxshr+=wx
maxyshr+=wy
Locate 1,1
Print "rref ";rref;" ";
Case Chr(27)
Exit do
End Select
If do_graph Then
skippy +=1
If (skippy >= skipper) Then
Scope
Dim As Double x,y,z,bx,by,bz,ex,ey,ez, kk,min, ee2, bb, minEE = 1e30, minBB = 1e30, maxEE = 0, maxBB = 0
Dim As UInteger color2
For k = 0 To 480
For i = 0 To 700
x = rref * (i-350)
z = rref * (k-240)
y = 0
ii = 0
ex = 0: ey = 0: ez = 0: bx = 0: by = 0: bz = 0
s = @quants[kk]
While ii < num_particles
from = @quants[ii]
tx1 = x - from->x
ty1 = y - from->y
tz1 = z - from->z
r2 = (tx1 * tx1) + (ty1 * ty1) + (tz1 * tz1)
d = Sqr(r2)
If (d > 1e-18) Then
r3 = r2 * d
f = from->q / (epsi_ * r3) ' r^3 = r2 * d = d * d * d
ff = miudiv4pi * from->q / r3
'' ---
ex1 = f * tx1
ey1 = f * ty1
ez1 = f * tz1
f = miudiv4pi * from->q / r3
bx1 = f * (from->vy * tz1 - from->vz * ty1)
by1 = f * (from->vz * tx1 - from->vx * tz1) '+ 0.2e9
bz1 = f * (from->vx * ty1 - from->vy * tx1)
If (from->v > 0) Then
writeF_mu_nu(@Fma, ex1, ey1, ez1, bx1, by1, bz1, c)
writeBoostMatrix_covariant(from, @ma1, c)
writeBoostMatrix_contravariant(from, @ma2, c)
computeElectromagneticTensorLorentzBoost(@ma1, @ma2, @Fma, @mres)
ex1 = -mres.e(0,1)*c
ey1 = -mres.e(0,2)*c
ez1 = -mres.e(0,3)*c
bx1 = mres.e(3,2)
by1 = mres.e(1,3)
bz1 = mres.e(2,1)
EndIf
Ex += ex1
Ey += ey1
Ez += ez1
Bx += bx1
By += by1
Bz += bz1
EndIf
ii += 1
Wend
' -min 0-255 r g b
ee2 = Log(Sqr(ex*ex + ey*ey + ez*ez))
bb = Log(Sqr(bx*bx + by*by + bz*bz))
If (ee2 < 0) Then ee2 = 0
If (bb < 0) Then bb = 0
If (ee2 < minEE) Then minEE = ee2
If (ee2 > maxEE) Then maxEE = ee2
If (bb < minBB) Then minBB = bb
If (bb > maxBB) Then maxBB = bb
Next i
Next k
If (minBB < 0) Then minBB = 0
If (minEE < 0) Then minEE = 0
'For j = 0 To 10
For k = 0 To 480
For i = 0 To 700
x = rref * (i-350)
z = rref * (k-240)
y = 0 'rref * (j-5)*10
ii = 0
ex = 0: ey = 0: ez = 0: bx = 0: by = 0: bz = 0
s = @quants[kk]
While ii < num_particles
from = @quants[ii]
tx1 = x - from->x
ty1 = y - from->y
tz1 = z - from->z
r2 = (tx1 * tx1) + (ty1 * ty1) + (tz1 * tz1)
d = Sqr(r2)
If (d > 1e-18) Then
r3 = r2 * d
f = from->q / (epsi_ * r3) ' r^3 = r2 * d = d * d * d
ff = miudiv4pi * from->q / r3
'' ---
ex1 = f * tx1
ey1 = f * ty1
ez1 = f * tz1
f = miudiv4pi * from->q / r3
bx1 = f * (from->vy * tz1 - from->vz * ty1)
by1 = f * (from->vz * tx1 - from->vx * tz1) '+ 0.2e9
bz1 = f * (from->vx * ty1 - from->vy * tx1)
If (from->v > 0) Then
writeF_mu_nu(@Fma, ex1, ey1, ez1, bx1, by1, bz1, c)
writeBoostMatrix_covariant(from, @ma1, c)
writeBoostMatrix_contravariant(from, @ma2, c)
computeElectromagneticTensorLorentzBoost(@ma1, @ma2, @Fma, @mres)
ex1 = -mres.e(0,1)*c
ey1 = -mres.e(0,2)*c
ez1 = -mres.e(0,3)*c
bx1 = mres.e(3,2)
by1 = mres.e(1,3)
bz1 = mres.e(2,1)
EndIf
Ex += ex1
Ey += ey1
Ez += ez1
Bx += bx1
By += by1
Bz += bz1
EndIf
ii += 1
Wend
ee2 = Log(Sqr(ex*ex + ey*ey + ez*ez))
bb = Log(Sqr(bx*bx + by*by + bz*bz))
If (ee2 < 0) Then ee2 = 0
If (bb < 0) Then bb = 0
color2 = RGB(0,((ee2-minEE)*3550/(maxEE-minEE)) And 255, ((bb-minBB)*3550/(maxBB-minBB))And 255)
'color2 = RGB(0,0, ((bb-minBB)*5550/(maxBB-minBB))And 255)
'color2 = RGB(0,((ee2-minEE)*5550/(maxEE-minEE)) And 255, 0)
'color2 = RGB(0,0, ((bb-minBB)*255/(maxBB-minBB))And 255)
'color2 = RGB(0,((ee2-minEE)*255/(maxEE-minEE)) And 255, 0)
'color2 = RGB(0,((ee2-minEE)*255/(maxEE-minEE)) And 255, ((bb-minBB)*72/(maxBB-minBB))And 255)
logit_y = logit(y)
logit_x = logit(x)
xx = -logit_x * sinangle + logit_y * cosangle + maxxshr
yy = -logit_x*ttp1 - logit_y*ttp2 + logit(z) * sinangle2 + maxyshr
PSet (xx,yy), color2
Next i
Next k
'BSave "relativity\sim_"+Str(steps)+".bmp",0
'Next j
End Scope
skippy = 0
End If
End If
' graphic output
If do_graph <> 1 Then
i = 0
While i < num_particles
'rotated by 2 angles
logit_y = logit(quants[i].y)
logit_x = logit(quants[i].x)
xx = -logit_x * sinangle + logit_y * cosangle + maxxshr
yy = -logit_x*ttp1 - logit_y*ttp2 + logit(quants[i].z) * sinangle2 + maxyshr
PSet(xx,yy), color1
i+=1
Wend
End If
Loop While 1 = 1
DeAllocate quants
End Sub
Sub show_help
Print "Usable keys during simulation:"
Print "h = This help"
Print "c = CLS "
Print "i = info about particle (x,y,z, vx, vy, vz, etc) "
Print "v = toggle FIELD Graph/particle trajectory "
Print "1 = angle -= 2 "
Print "2 = angle += 2 "
Print "3 = angle2 -= 2 "
Print "4 = angle2 += 2 "
Print "KEY UP = move UP viewport "
Print "KEY DOWN = move DOWN viewport "
Print "KEY LEFT = move LEFT viewport "
Print "KEY RIGHT = move RIGHT viewport "
Print "- = zoom out (rref *= 2.222222) "
Print "+ = zoom in (rref /= 2.222222) "
Print "TAB = toggle simulation run (on/off"
End Sub