Relativistic Electrodynamics Simulator

[Mihail_B]

Relativistic Electrodynamics Simulator - free version, but instant field used instead of genius positions. I don't have a release date for a FREE genius position simulator.

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

[Mihail_B]

Relativistic Electrodynamics Simulator - free version, uses genius positions (genius fields)

(You need a 64bit compiler and at least 3980MB free RAM)
To compile use: [path]\fbc\win64-1.05.0\fbc -s console -t 128000 -w 3 -O 2 -v "relativity_retardedfield.bas"

It takes a lot of time - it's relativistic so lot's of Lorentz Boosts in arbitrary directions = lots of calculus.
I didn't optimize much the genius position calculator so you need to patient.

Press 'h' during simulation to get a help
NOTE: Any key command will be processed after the a cycle is complete (aka after num_history/num_particles computations); to better view from other angles press and then wait; the use angle modified keys to render from a different perspective. Use to resume simulation.

Initial setup : (constants : e = positive fundamental charge value ; m_elec = mass of electron ; m_prot = mass of proton)
quants[0].q = -e
quants[0].m = m_elec
quants[0].t = 0
quants[0].x = 0.52e-12
quants[0].z = 0.52e-12
quants[0].y = 0
quants[0].vx = -6.7e6
quants[0].vy = 0
quants[0].vz = +6.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 ' m_prot
quants[1].vx = 6.7e6
quants[1].vy = 0
quants[1].vz = -6.7e6
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[1].x = 0, quants[1].y = 0, quants[1].z = 0 )

Some screenshots : (URLs are available from 8:00UTC until 20:00UTC, Monday to Sunday, closed otherwise)
1) Non-genius Positions (instant field prediction): http://alphax.info/pozne/test1.png
2) genius Positions (totally different prediction :) ) http://alphax.info/pozne/test2.png
3) genius, different initial configuration http://alphax.info/pozne/test3.png
4) genius (4) after [some time] (3) http://alphax.info/pozne/test4.png
5) genius (4) after (5) http://alphax.info/pozne/test_final5.png
6) genius (4) after (6) http://alphax.info/pozne/test_final6.png
7) Info about positron(e+) at step (6) (after 1,079,999,928 computations) http://alphax.info/pozne/test_final7_info.png

Code: Select all

' WARNING: Don't use with medical applications !

' Relativistic Electrodynamics Simulator (free version, uses genius positions)
' In the framework of Einstein's Special Relativity 
' (at least in the context of what I undestood of it :) )

' Note: Can simulate in the limit of 64bit double numbers. 
'       When too close to the speed of light (like >99.9999%) 64bit numbers are unsatisfactory, 
'       leading to computation errors.


' Tested with FreeBASIC Compiler - Version 1.05.0 (01-31-2016), built for win64 (64bit)
' Note: dont use 32bit compiler as 32bit applications wont be able to allocate enough memory
' PS: you can still compile it with 32bit compiler but you will need to lower num_history value to something lower than 2GB

' To compile use: [path]\fbc\win64-1.05.0\fbc -s console -t 128000 -w 3 -O 2 -v "Relativity_instantfield.bas"

' under current setup (num_history = 30000000), it needs 3980MB of free RAM to run.
' num_history = 30000000 means 30milion particle positions; if running for 2 particles => 15milion each
' can be lowered to 10milion (2milion per particle) as long as sufficent to get genius position in all events. 

' My website http://alphax.info (Note: opened from 8:00 GMT until 20:00 GMT, Monday to Sunday, closed during other hours)
' Mihai Barboi, 2017 - open source


' press 'h' during simulation to get a help
' NOTE: Any key command will be processed after the a cycle is complete 
' (aka after num_history/num_particles computations); to better view from other angles press TAB and then wait; the use angle modified keys to render from a different perspective 
' Other keys :
' c = CLS  
' i = info about particle (x,y,z, vx, vy, vz, etc)
' I = info about all particles
' 1 = angle -= 2  
' 2 = angle += 2  
' 3 = angle2 -= 2  
' 4 = angle2 += 2  
' KEY UP = move UP viewport  
' KEY DOWN = move DOWN viewport  
' KEY LEFT = move LEFT viewport  
' KEY RIGHT = move RIGHT viewport  
' - = zoom out (rref *= 2.222222) 
' + = zoom in (rref /= 2.222222)  
' TAB = toggle simulation run (on/off)


#include "fbgfx.bi"


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
'#Include "mem.bas"
'#Include "easy.bas"

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

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

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)))))

Type vector
   As Double x,y,z,a
End Type

Declare Sub show_help
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 matrixMatrixMult(m1 As matrix4x4 Ptr, m2  As matrix4x4 Ptr, res As matrix4x4 Ptr)
Declare Sub contravariantVelocityLorentzBoost(m1 As matrix4x4 Ptr, c As double, v As vector4 ptr , result As vector4 Ptr)
Declare Sub computeElectromagneticTensorLorentzBoost(m1 As matrix4x4 Ptr, m2  As matrix4x4 Ptr, F As matrix4x4 Ptr, res As matrix4x4 Ptr)
Declare Sub writeBoostMatrix_contravariant(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Declare Sub writeBoostMatrix_covariant(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Declare Sub writeBoostMatrixFromVect_contravariant(from As vector4 Ptr, m As matrix4x4 Ptr, c As Double)
Declare Sub velocityAdition(from As quanta ptr, s As quanta Ptr, vect1 As vector4 Ptr, c2 As Double)
Declare Sub writeBoostMatrixFromVect_covariant(from As vector4 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 info_quanta(q as quanta Ptr, steps As ULongInt, c As Double, h As Double)
Declare Sub coordTrans(s As quanta Ptr, c As Double, dt As Double)
Declare Sub lorentzBoostEB(s As quanta ptr, from As quanta Ptr, c As Double)
Declare Sub computeXU(s As quanta Ptr, en As quanta Ptr, dt As Double, c As Double) 
Declare function computeEB(s As quanta Ptr, from As quanta Ptr, epsi_ As Double, epsi_c2 As Double, miudiv4pi As Double, c As Double, ma1 As matrix4x4 Ptr, ma2 As matrix4x4 Ptr, fma As matrix4x4 Ptr, mres As matrix4x4 ptr) As Double
Declare Sub CopyQuanta(s As quanta Ptr, p As quanta Ptr)
Declare function distance(s As quanta Ptr, p As quanta Ptr)  As Double
Declare Sub electrodynamics()

/'
Dim As Double c,v,m0,e,de,vf,c2,h,l10,e0,nf,freq
Dim As Integer i,k,j,x,y,x0,y0
Dim As String a
Dim As Double relativeEpsi, chargeElec, massElec, massProt, epsi, grav, c4
c = 3e8
c2 = c*c
c4=c2*c2
m0=1e-27
h=6.6e-36
relativeEpsi = 1.216e-17
chargeElec = 1.6e-19
massElec  = 9.1e-31
massProt = 1.627e-27
epsi = 8.854e-12
grav = 6.673e-11
'/
Screen 19,32
ScreenInfo maxx,maxy

electrodynamics()

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 "I = info about all particles"
   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

'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 matrixMatrixMult(m1 As matrix4x4 Ptr, m2  As matrix4x4 Ptr, res As matrix4x4 Ptr)
   ' A{nm} : n x m
   ' B{mp} : m x p
   ' (AB){ij} = sum{k=1}[m] A{ik} B{kj}
   Dim As Integer i,k,j
   For i = 0 To 3
      For j = 0 To 3
         res->e(i,j) = 0
         For k = 0 To 3
            res->e(i,j) += m1->e(i, k) * m2->e(k, j)
         Next k
      Next j
   Next i
End Sub
Sub contravariantVelocityLorentzBoost(m1 As matrix4x4 Ptr, c As double, v As vector4 ptr , result As vector4 Ptr)
   Dim As Double m,n
   For m = 0 To 3
      result->e(m) = 0
      For n = 0 To 3
         result->e(m) += m1->e(m,n) * v->e(n)
      Next n
   Next m
   result->e(0) /= c
   result->e(1) /= result->e(0)
   result->e(2) /= result->e(0)
   result->e(3) /= result->e(0)
End Sub
Sub convariantVelocityLorentzBoost(m1 As matrix4x4 Ptr, c As double, v As vector4 ptr , result As vector4 Ptr)
   Dim As Double m,n
   For m = 0 To 3
      result->e(m) = 0
      For n = 0 To 3
         result->e(m) += m1->e(m,n) * v->e(n)
      Next n
   Next m
   result->e(0) /= c
   result->e(1) /= result->e(0)
   result->e(2) /= result->e(0)
   result->e(3) /= result->e(0)
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_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 writeBoostMatrix_covariant(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 writeBoostMatrixFromVect_contravariant(from As vector4 Ptr, m As matrix4x4 Ptr, c As Double)
   Dim As Double bx, by, bz, b2, g
   
   g = from->e(0)
   bx = from->e(1) / c 
   by = from->e(2) / c 
   bz = from->e(2) / c 
   b2 = (bx * bx + by * by + bz * bz)

   
   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 velocityAdition(from As quanta ptr, s As quanta Ptr, vect1 As vector4 Ptr, c2 As Double)
   ' w2 = ((u-v).(u-v) - (u x v)2/c2 )/(1- (u.v)/c2)2
   ' u.v = uxvx+uyvy+uzvz ; (u x v)2 = (u.u)(v.v) - (u.v)2 ; u-v=(ux-vx, uy-vy, uz-vz)
   Dim As Double a1, a2, a3, w2, tx1, ty1, tz1, r, a4
   a1 = (s->vx - from->vx)*(s->vx - from->vx) + (s->vy - from->vy)*(s->vy - from->vy) + (s->vz - from->vz)*(s->vz - from->vz)
   a2 = (from->vx*from->vx + from->vy*from->vy+from->vz*from->vz)*(s->vx*s->vx + s->vy*s->vy+s->vz*s->vz)
   a3 = (from->vx*s->vx + from->vy*s->vy + from->vz*s->vz)
   a4 = (1-a3/c2)
   w2 = (a1 - (a2 - a3*a3)/c2)/(a4*a4)
   
   tx1 = s->x - from->x
	ty1 = s->y - from->y
	tz1 = s->z - from->z
	r = Sqr((tx1 * tx1) + (ty1 * ty1) + (tz1 * tz1))
	
	w2 = w2/r
	vect1->e(1) = w2 * tx1
	vect1->e(2) = w2 * ty1
	vect1->e(3) = w2 * tz1 

   
End Sub
Sub writeBoostMatrixFromVect_covariant(from As vector4 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 = from->e(0)
   bx = from->e(1) / c 
   by = from->e(2) / c 
   bz = from->e(2) / c 
   b2 = (bx * bx + by * by + bz * bz)
   
   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 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 info_quanta(q as quanta Ptr, steps As ULongInt, c As Double, h As Double)
   Print "m ";q->m;" q ";q->q; " steps ";steps;"  "
   Print "x ";q->x;" ";q->y;" ";q->z;" t ";q->t;"  "
   Print "a ";q->a;" V ";q->v; " ct ";c*q->t;"  "
   Print "v ";q->vx;" ";q->vy;" ";q->vz;"  "
   Print "E ";q->Ex;" ";q->Ey;" ";q->Ez;"  "
   Print "B ";q->Bx;" ";q->By;" ";q->Bz;"  "
   Print "K ";q->m*c*c*((1/q->a) - 1)/h;" Hz => " ;Log(q->m*c*c*((1/q->a) - 1)/h)/Log(10);" (log)Hz"
   Print 
End Sub
Sub coordTrans(s As quanta Ptr, c As Double, dt As Double) 
   ' L{ik} = |  g      -Bxg                 -Byg                -Bzg               | |ct
   '         | -Bxg    1 + (g - 1)Bx2/B2    (g - 1)BxBy/B2      (g - 1)BxBz/B2     | | x
   '         | -Byg    (g - 1)ByBx/B2       1 + (g - 1)By2/B2   (g - 1)ByBz/B2     | | y
   '         | -Bzg    (g - 1)BzBx/B2       (g - 1)BzBy/B2       1 + (g - 1)Bz2/B2 | | z
   Dim As Double g, bx, by, bz, b, b2
   Dim As Double x, y, z, t, ct
   g = 1 / s->a
   bx = s->vx / c
   by = s->vy / c
   bz = s->vz / c
   b = s->v / c
   b2 = b * b
   If (b = 0) Then return
   ct = s->t * c
   t = g*ct - bx*g*s->x - by*g*s->y - bz*g*s->z
   x = -Bx*g*ct + (1+(g-1)*bx*bx/b2)*s->x + (g-1)*bx*by/b2*s->y + (g-1)*bx*bz/b2*s->z
   y = -by*g*ct + (g-1)*by*bz/b2*s->x + ((g-1)*by*bz/b2)*s->y + (g-1)*by*bz/b2*s->z
   z = -bz*g*ct + (g-1)*bz*bx/b2*s->x + (g-1)*bz*by/b2*s->y + (g-1)*bz*bz/b2*s->z
   
   't = ( 1 - Bx  - by  - bz) * g * ct
   'x = (-bx * g +  1+(g-1)*bx*bx/b2 + (g-1)*bx*by/b2 + (g-1)*bx*bz/b2) * s->x
   'y = (-by * g +  (g-1)*by*bz/b2 + 1 + (g-1)*by*by/b2 + (g-1)*by*bz/b2) * s->y
   'z = (-bz * g +  (g-1)*bz*bx/b2 + (g-1)*bz*by/b2 + 1 + (g-1)*bz*bz/b2) * s->z 
   
   
   
   Locate 1,1
   Print s->x;" ";s->y;" ";s->z;" ";s->t;"   "
   Print x; " ";y; " "; z; " ";t/c; " "
   's->x = x
   's->y = y
   's->z = z
   's->t += dt 't/c
   
   pause1
   
End Sub

Sub lorentzBoostEB(s As quanta ptr, from As quanta Ptr, c As Double)
   Dim As Double ex, ey, ez, bx, by, bz
   Dim As Double matrix(0 To 3,0 To 3)
   ' F{ab} = |  0     Ex/c  Ey/c  Ez/c |  
   '         | -Ex/c     0   -Bz    By | 
   '         | -Ey/c    Bz     0   -Bx | 
   '         | -Ez/c   -By    Bx     0 | 
   '
   ' g = 1/a = 1/sqr(1-v2/c2)
   ' B{k} = v{k}/c
   ' F'{mn} = Lmi Fik Lnk* ?
   ' F'{mn} = L{m}[i] L{n}[k] F{ik}
   
   ' L{ik} = |  g      -Bxg                 -Byg                -Bzg               |
   '         | -Bxg    1 + (g - 1)Bx2/B2    (g - 1)BxBy/B2      (g - 1)BxBz/B2     |
   '         | -Byg    (g - 1)ByBz/B2       1 + (g - 1)By2/B2   (g - 1)ByBz/B2     |
   '         | -Bzg    (g - 1)BzBx/B2       (g - 1)BzBy/B2       1 + (g - 1)Bz2/B2 |
   
   ' L{ik} = |  1/a    -vx/ca                 -vy/ca                -vz/ca             |
   '         | -vx/ca  1 + (1/a - 1)vx2/v2    (1/a - 1)vxvy/v2      (1/a - 1)vxvz/v2   |
   '         | -vy/ca  (1/a - 1)vyvz/v2       1 + (1/a - 1)vy2/v2   (1/a - 1)vyvz/v2   |
   '         | -vz/ca  (1/a - 1)vzvx/v2       (1/a - 1)vzvy/v2      1 + (1/a -1)vz2/v2 |
   
   ' F'{ik} = L[i]{a} L[k]{b} F{ab}
   ' F{00} = L[0]{a} L[0]{b} F{ab} =   
   
   ' [vx, vy, vz] : velocity relative to  frame 2 from frame 1; speed composition
   ' but we use "vacuum" and we are considering field in general; no need velocity transformation
   
   ' F'{mn} = deron(x[a])/deron(x[m]) deron(x[b])/deron(x[n]) F{ab} 
   ' ^^^^^^^^^^ much better
   ' X{1} = dX[1]/dY[0] Y{0} + dX[1]/dY[1] Y{1} + dX[1]/dY[2] Y{2} + dX[1]/dY[3] Y{3}
   ' x = vx/a * T + 1/a X + dx/dy Y + dx/dz Z 
   ' x = (vx*T + X)/a + vx Y/vya + vx Z/vza
   ' x = (vx*T + X + vx*Y/vy + vx*z/vz)/a
   ' 
   ' matrix(0,0) = 
   
End Sub
' computes new X4 (position) and U4 (velocity)
' en = (ending) = final
' s = (source) = initial
Sub computeXU(s As quanta Ptr, en As quanta Ptr, dt As Double, c As Double) 
   'all up
   ' (0)  m * c /a' = e / (a*c) * (Ex*vx+ Ey*vy + Ez*vz) * dtau + m*c/a = (D)
   ' (1)  m * vx' / a' = e/a * (Ex + Bz*vy - By*vz) * dtau + m*vx/a = (A)
   ' (2)  m * vy' / a' = e/a * (Ey - Bz*vx + Bx*vz) * dtau + m*vy/a = (B)
   ' (3)  m * vz' / a' = e/a * (Ez + By*vx - Bx*vy) * dtau + m*vz/a = (D)
   
   ' v^2 = vx^2 + vy^2 + vz^2   
   ' A/vx' = B/vy' = C/vz' = D/c = m/a' 
   ' =>
   ' A/vx' = D/c   <=> vx' = A*c/D 
   ' B/vy' = D/c   <=> vy' = B*c/D
   ' C/vz' = D/c   <=> vz' = C*c/D
   ' D/c = m/a'    <=> a' = m*c/D
   ' a' = sqr(1 - v^2/c^2) = sqr(1 - (vx'^2 + vy'^2 + vz'^2)/c^2)
   
   Dim As Double AA, BB, CC, DD
   Dim As Double dtau, temp, temp2
   
   /'
   dtau = dt * s->a ' dtau = tau1 - tau0 = t1*a1 - t0*a0 = dt
     
   DD = (s->q / (s->a * c)) * ( s->Ex * s->vx + s->Ey * s->vy + s->Ez * s->vz) * dtau + s->m * c / s->a
   
   AA = (s->q / s->a) * (s->Ex + s->Bz * s->vy - s->By * s->vz) * dtau + s->m * s->vx / s->a
   
   BB = (s->q / s->a) * (s->Ey - s->Bz * s->vx + s->Bx * s->vz) * dtau + s->m * s->vy / s->a
   
   CC = (s->q / s->a) * (s->Ez + s->By * s->vx - s->Bx * s->vy) * dtau + s->m * s->vz / s->a
   '/

   temp2 = s->q * dt
   temp = s->m / s->a
     
   DD = (temp2/c) * ( s->Ex * s->vx + s->Ey * s->vy + s->Ez * s->vz) + c * temp
   
   AA = temp2 * (s->Ex + s->Bz * s->vy - s->By * s->vz) + s->vx * temp
   
   BB = temp2 * (s->Ey - s->Bz * s->vx + s->Bx * s->vz) + s->vy * temp
   
   CC = temp2 * (s->Ez + s->By * s->vx - s->Bx * s->vy) + s->vz * temp
   
   en->m = s->m
   en->q = s->q
   en->t = s->t + dt
   
   temp = c / DD
   en->vx = AA * temp ' = AA * c / DD
   en->vy = BB * temp ' = BB * c / DD
   en->vz = CC * temp ' = CC * c / DD
   'en->a = s->m * temp ' = s->m * c / DD
   
   en->v = Sqr(en->vx * en->vx + en->vy * en->vy + en->vz * en->vz)  
   
   temp = en->v / c
   en->a = Sqr(1 - temp * temp)
   
   'coordTrans(en, c, dt)
   
   ' compute new position relative to our "observer"
   
   
   'Return
   ' temp = dt/en->a
   
   en->x = s->x + (s->vx + en->vx)/2 * dt ' temp
   en->y = s->y + (s->vy + en->vy)/2 * dt ' temp
   en->z = s->z + (s->vz + en->vz)/2 * dt ' temp
   
   /'
   If (en->a > 1) Then
      Locate 1,1
      Print "D";dd; " A";aa; " B"; bb; " C";cc;" a";en->a; " q";en->q
      pause1
   End If
   '/
    
   /'
   en->Bx = 0
   en->By = 0
   en->Bz = 0
   en->Ex = 0
   en->Ey = 0
   en->Ez = 0
   '/
End Sub

' from any particle in radius r = - ct (t = light time, field change propagation speed)
' s = final
' from = initial
function computeEB(s As quanta Ptr, from As quanta Ptr, epsi_ As Double, epsi_c2 As Double, miudiv4pi As Double, c As Double, ma1 As matrix4x4 Ptr, ma2 As matrix4x4 Ptr, fma As matrix4x4 Ptr, mres As matrix4x4 ptr) As double
   Dim As Double tx1, ty1, tz1, bx, by, bz, fx, fy, fz
   Dim As Double d, f, Rmax, r2
   
   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
	   f = from->q / (epsi_* r2 * d) ' r^3
	   fx = f * tx1 
	   fy = f * ty1 
	   fz = f * tz1 
	   
      f = miudiv4pi * from->q / (r2 * d)
		bx = f * (from->vy * tz1 - from->vz * ty1)  
		by = f * (from->vz * tx1 - from->vx * tz1)
		bz = f * (from->vx * ty1 - from->vy * tx1)
		
		If (from->v > 0) Then
	      writeF_mu_nu(Fma, fx, fy, fz, bx, by, bz, c)
         writeBoostMatrix_covariant(from, ma1, c)
         writeBoostMatrix_contravariant(from, ma2, c)
	      
	      computeElectromagneticTensorLorentzBoost(ma1, ma2, Fma, mres)
	      
	      fx = -mres->e(0,1)*c
	      fy = -mres->e(0,2)*c
	      fz = -mres->e(0,3)*c
	      bx = mres->e(3,2)
	      by = mres->e(1,3)
	      bz = mres->e(2,1)
		End If
		
		/'
		If (1 = 0 And s->v > 0) Then
	      writeF_mu_nu(Fma, fx, fy, fz, bx, by, bz, c)
         writeBoostMatrix_covariant_reverse(s, ma1, c) '_reverse
         writeBoostMatrix_contravariant_reverse(s, ma2, c) '_reverse
	      
	      computeElectromagneticTensorLorentzBoost(ma1, ma2, Fma, mres)
	      
	      fx = -mres->e(0,1)*c
	      fy = -mres->e(0,2)*c
	      fz = -mres->e(0,3)*c
	      bx = mres->e(3,2)
	      by = mres->e(1,3)
	      bz = mres->e(2,1)
		EndIf
		'/
		s->Ex += fx
      s->Ey += fy
      s->Ez += fz 
      
		s->Bx += bx
      s->By += by
      s->Bz += bz
      
      ' we have to apply a general boost transformation to (fx,fy,fz, bx,by,bz) for (vx,vy,vz) at new coords 
      
	Else  ' review this
	   Return rmax 
	End If
   
    
   '             | i  j  k  |			| i  j  k  |
	' a x b = det | a1 a2 a3 | = det | x1 y1 z1 |
	'				  | b1 b2 b3 |			| x2 y2 z2 |
	' a x b = (a2*b3 - a3*b2)i + (a3*b1 - a1*b3)j + (a1*b2 - a2*b1)k
   Return d
End Function

Sub CopyQuanta(s As quanta Ptr, p As quanta Ptr)
   s->x = p->x
   s->y = p->y
   s->z = p->z
   s->vx = p->vx
   s->vy = p->vy
   s->vz = p->vz
   s->q = p->q
   s->m = p->m
   s->t = p->t
   s->Ex = p->Ex
   s->Ex = p->Ey
   s->Ez = p->Ez
   s->Bx = p->Bx
   s->By = p->By
   s->Bz = p->Bz
   s->a = p->a
   s->v = p->v
End Sub
function distance(s As quanta Ptr, p As quanta Ptr)  As double
   Dim As Double dr, temp
   dr = Abs(s->x - p->x) 
   dr = dr * dr
   temp = Abs(s->y - p->y)
   temp = temp*temp 
   dr += temp
   temp = abs(s->z - p->z)
   temp = temp*temp
   dr += temp
   dr = Sqr(dr)
   return dr
End Function
Sub electrodynamics
   ' From covariant formulation of electromagnetism we have
   ' p{a}/dtau = e * F{ab} * u[b]
   ' p[a]/dtau = e * F[a]{b} *u[b]
   ' {..} = covariant (lower index} ; [..] = contravariant (upper index) 
   ' a = sqr(1 - v^2/c^2),
   ' gamma = 1/a , 
   ' m = m0 (rest mass), 
   ' e = charge,  
   ' c = speed of light 3e8 m/s
   ' dtau = dt * a, (proper time)
   ' v = dx / dt
   ' v^2 = vx^2 + vy^2 + vz^2
   ' Ex, Ey, Ez = electric field 
   ' Bx, By, Bz = magnetic field
   ' E = mc^2 (energy, in Joules)
   ' ^ = power (takes highest precedence here)
   
   ' under Minkowski metric :
   ' (proper time dtau^2 = dt^2 - (1/(c^2))*(dx^2 + dy^2 + dz^2) 
   ' (proper distance dsigma^2 = -(c*t)^2 + dx^2 + dy^2 + dz^2)
   
   ' x[a] = distance/position etc in [m];
   ' u[a] = gamma(c, u)  ; u[a] velocity in [m/s]
   ' p{a} = m0 * u[a] ; p = impulse
   ' F{ab} = electromagnetic tensor
   ' ita = [-1, 1, 1, 1]
   
   ' x[0] = c*t , x[1] = x , x[2] = y , x[3] = z
   ' u[0] = c/a ; u[1] = vx/a ; u[1] = vy/a ; u[2] = vz/a
   ' p[0] =  m*c/a ; p[1] = m*vx/a ; p[2] = m*vy/a ; p[3] = m*vz/a
   ' p{0} = -m*c/a ; p{1} = m*vx/a ; p{2} = m*vy/a ; p{3} = m*vz/a
   ' F[ab] = |  0    -Ex/c -Ey/c -Ez/c |  
   '         |  Ex/c     0   -Bz    By | 
   '         |  Ey/c    Bz     0   -Bx | 
   '         |  Ez/c   -By    Bx     0 |
   ' F{ab} = |  0     Ex/c  Ey/c  Ez/c |  
   '         | -Ex/c     0   -Bz    By | 
   '         | -Ey/c    Bz     0   -Bx | 
   '         | -Ez/c   -By    Bx     0 |
   ' 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 | 
   
   ' {ab} -> a=row pos, b=col pos (like in LOCATE 10,10)
   ' P{a} = g{ab}P[b]
   ' http://www.feynmanlectures.caltech.edu/II_26.html Feynman uses F{ab} like [0, -Ex/c, -Ey/c, -Ez/c] for time component.
   
   ' for [0,0 to 3]
   ' dp[0]/dtau = -e * F[0]{0} * u[0] =>  m*c/(a*dtau) = -e * 0 * -c/a
   ' dp[0]/dtau = -e * F[0]{1} * u[1] =>  m*c/(a*dtau) = -e * (-Ex/c) * vx/a
   ' dp[0]/dtau = -e * F[0]{2} * u[2] =>  m*c/(a*dtau) = -e * (-Ey/c) * vy/a
   ' dp[0]/dtau = -e * F[0]{3} * u[3] =>  m*c/(a*dtau) = -e * (-Ez/c) * vz/a
   ' ----- sum ----- (add all rows)
   ' d m*c/a / dtau = e/a *(-(-Ex/c)*vx -(-Ey/c*vy) -(- Ez/c*vz))
   ' d m*c/a =  e / (a*c) * (Ex*vx+ Ey*vy + Ez*vz) * dtau
   ' m * c /a' = e / (a*c) * (Ex*vx+ Ey*vy + Ez*vz) * dtau + m*c/a
   
   ' for [1, 0 to 3]
   ' p[1]/dtau = e * F[1]{0} * u[0] => m*vx/(a*dtau) = -e * (-Ex/c) * c/a
   ' p[1]/dtau = e * F[1]{1} * u[1] => m*vx/(a*dtau) = -e * (0) * vx/a
   ' p[1]/dtau = e * F[1]{2} * u[2] => m*vx/(a*dtau) = -e * (-Bz) * vy/a
   ' p[1]/dtau = e * F[1]{3} * u[3] => m*vx/(a*dtau) = -e * (By) * vz/a
   ' ------ sum ----- (add all rows)
   ' d m*vx/(a*dtau) = -e/a * (-Ex*c/c - (Bz*vy) + (By*vz))
   ' d m*vx/a  = e/a *(Ex + Bz*vy - By*vz) * dtau
   ' m * vx' / a' = e/a *(Ex + Bz*vy - By*vz) * dtau + m*vx/a
   
   ' for [2, 0 to 3]
   ' p[2]/dtau = e * F[2]{0} * u[0] => m*vy/(a*dtau) = -e * (-Ey/c) * c/a
   ' p[2]/dtau = e * F[2]{1} * u[1] => m*vy/(a*dtau) = -e * (Bz) * vx/a
   ' p[2]/dtau = e * F[2]{2} * u[2] => m*vy/(a*dtau) = -e * (0) * vy/a
   ' p[2]/dtau = e * F[2]{3} * u[3] => m*vy/(a*dtau) = -e * (-Bx) * vz/a
   ' ------ sum ------ (add all rows)
   ' d m*vy/(a*dtau) = -e/a * (-Ey*c/c + Bz*vx + 0 - (Bx*vz))
   ' d m*vy/a = e/a * (Ey - Bz*vx + Bx*vz)*dtau
   ' m * vy' /a' = e/a * (Ey - Bz*vx + Bx*vz)*dtau + m*vy/a
   
   ' for [3, 0..3]
   ' dp[3]/dtau = e * F{3]{0} * u[0] => m*vz/(a*dtau) = -e * (-Ez/c) * c/a
   ' dp[3]/dtau = e * F{3]{1} * u[1] => m*vz/(a*dtau) = -e * (-By) * vx/a
   ' dp[3]/dtau = e * F{3]{2} * u[2] => m*vz/(a*dtau) = -e * (Bx) * vy/a
   ' dp[3]/dtau = e * F{3]{3} * u[3] => m*vz/(a*dtau) = -e * (0) * vz/a
   ' ------ sum ----- (add all rows)
   ' d m*vz/(a*dtau) = - e/a * (-Ez - By*vx + Bx*vy)
   ' d m*vz/a = e/a * (Ez + By*vx - Bx*vy)*dtau
   ' m * vz' / a' = e/a * (Ez + By*vx - Bx*vy)*dtau + m*vz/a
   
   'all up
   ' (0)  m * c /a' = e / (a*c) * (Ex*vx+ Ey*vy + Ez*vz) * dtau + m*c/a = (D)
   ' (1)  m * vx' / a' = e/a * (Ex + Bz*vy - By*vz) * dtau + m*vx/a = (A)
   ' (2)  m * vy' / a' = e/a * (Ey - Bz*vx + Bx*vz) * dtau + m*vy/a = (B)
   ' (3)  m * vz' / a' = e/a * (Ez + By*vx - Bx*vy) * dtau + m*vz/a = (D)
   
   ' v^2 = vx^2 + vy^2 + vz^2   
   ' A/vx' = B/vy' = C/vz' = D/c = m/a' 
   ' =>
   ' A/vx' = D/c   <=> vx' = A*c/D 
   ' B/vy' = D/c   <=> vy' = B*c/D
   ' C/vz' = D/c   <=> vz' = C*c/D
   ' D/c = m/a'    <=> a' = m*c/D
   ' a' = sqr(1 - v^2/c^2) = sqr(1 - (vx'^2 + vy'^2 + vz'^2)/c^2)
   
   ' @(x,y,z,t) [vx, vy, vz, Ex(t), Ey(t), Ez(t), Bx(t), By(t), Bz(t)]
   ' stuff : (x,y,z,t), vx, vy, vz, m, q
   ' traveling with c, from p(i)(x,y,z,t) to p(me)(x', y', z', t');  Ex(t), Ey(t), Ez(t), Bx(t), By(t), Bz(t)
   
   Dim As Double c = 299792458
   Dim As Double c2 = c*c
   Dim As Double c4 = c2*c2
   Dim As Double m0 = 1e-27
   Dim As Double h = 6.6e-36
   Dim As Double e = 1.6e-19
   Dim As Double m_elec  = 9.1e-31
   Dim As Double m_prot = 1.627e-27
   Dim As Double m_up_quark = 4.08888e-30 'up quark +2/3e ('2.3 1e6  1.6e-19 1/9 1e-16) 29 ' radius(u-u) Rmax= 0.90097 e-16 
   dim as double m_down_quark = 8.53333e-30 'down quark -1/3e  (u+u+d = 16.71109e-30 kg)
   Dim As Double miu = 4*pi*1.0e-7
   Dim As Double epsi = 8.854e-12
   Dim As Double grav = 6.673e-11
   Dim As Double pi2 = pi * 2
   Dim As Double miudiv4pi = miu / (4*pi)
   Dim As Double epsi_ = 4*pi*epsi
   Dim As Double epsi_c2 = (epsi_) * c2
   
   Dim As Double var_theta, var_phi, sin1, sin2, x, y, z, phi2, theta2, real_theta, real_phi
   Dim As Double unghi = 25, unghi2 = 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 sinunghi, cosunghi, sinunghi2, cosunghi2, ttp1, ttp2, maxxshr, maxyshr ' 3d graphics translation
   Dim As UInteger color1 = RGB(155,255,255)
   Dim As UInteger color2 = RGB(255,055,055)
   Dim As UInteger color3 = RGB(055,255,055)
   Dim As ULongInt num = 0, steps = 0, skip = 0
   
   Dim As Double s
   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, spaceQuanta
   Dim As Integer cont = 0
   Dim As String a
   Dim As Double startTime, universalTime, rmin
   Dim As Integer currentTime, np, maxTime
   Dim As Double tempt, tempx, tempy, tempz
   Dim As Integer tempi, tempk, maxtempi, tempdelta
   
   Dim As matrix4x4 ma1, ma2, fma, mres
   Dim As quanta Ptr quants
   Dim As Integer num_history, size, i, k, j, ii, kk, jj
   
   
   unghi = 206.6
   unghi2 = 180
   rref = 0.5e-14
   sinunghi=Sin(to_rads(90-unghi))
	cosunghi=cos(to_rads(90-unghi))
	sinunghi2=sin(to_rads(90-unghi2))
	cosunghi2=cos(to_rads(90-unghi2))
   x_opt = CInt(maxx Shr 1)+wx
   y_opt = CInt(maxy Shr 1)+wy  
    
   ttp1=cosunghi*cosunghi2
	ttp2=sinunghi*cosunghi2
	maxxshr=maxx Shr 1
	maxyshr=maxy Shr 1
	maxxshr+=wx
	maxyshr+=wy
   
   num_history =  30000000 ' 10000000 ' 30000000
   quants = CPtr(quanta Ptr, Allocate(num_history * SizeOf(quanta)))
   For i = 0 To num_history-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 = num_history * SizeOf(quanta) / (1024*1024)
   Print "allocated ";size;"MB ";IIf(quants <> 0, "yes", "no")
   
   Dim As Integer iu, num_particles
   
   spaceQuanta = 0.9e-18
   timeQuanta = timeQuanta/c 'sec ct = 1e-19 't = 0.3333333e-27
   num_particles = 2 '4
   np = 1'2
   cls
   quants[0].q = -e
   quants[0].m = m_elec
   quants[0].t = 0
   quants[0].x = 0.52e-12     ' good value 1e-12
   quants[0].z = 0.52e-12
   quants[0].y = 0     ' good value 1e-12     
   quants[0].vx = -6.7e6
   quants[0].vy = 0 '7e6
   quants[0].vz = +6.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 'm_elec ' m_prot
   quants[1].vx = 6.7e6
   quants[1].vy = 0
   quants[1].vz = -6.7e6
   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)
   '/
   
   
   ' time{start} = dr/c 
   startTime = distance(@quants[0], @quants[1])/c
   universalTime = startTime
   maxtempi = 0
   
   Do
      ' computations
      If (maxtempi > 0) Then 
         currentTime = ((num_history-1-maxtempi) Shr np)
         ' Locate 2,1: Print "currentTime ";currentTime;" ";maxtempi;"      ";
         maxtempi = 0
      Else
         currentTime = 0
         maxtempi = 0
      EndIf
      maxTime = (num_history Shr np) - np
      If (skip = 0) then
         Do
            'If (currentTime And &hfffff) = 0 Then Locate 2,1: Print currentTime;"   ";
            k = currentTime Shl np
            j = k + num_particles
            i = k
            r = 10000         
            While i < j
               ii = k
               jj = j
               While ii < jj
                  If (ii <> i) Then
                     ' go back in time and compute to quants[i] from quants[ii] 
                     ' where ii is the appropiate moment in time when c * dt = dr 
                     ' (dr distance between the 2 charges) ; dr = rmin/c
                     ' choose ii such that (quants[i].t - quants[ii].t)*c = sqrt((quants[i].x - quants[ii].x)^2+(quants[i].y - quants[ii].y)^2+(quants[i].z - quants[ii].z)^2)
                     tempdelta = 17000
                     If (ii < tempdelta Shl np) Then 
                        rmin = computeEB(@quants[i], @quants[ii], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                        If (r > rmin) Then r = rmin
                     Else
                        tempi = ii
                        dummy_search:
                        tempx = (quants[i].x - quants[tempi].x): tempx *= tempx
                        tempy = (quants[i].y - quants[tempi].y): tempy *= tempy
                        tempz = (quants[i].z - quants[tempi].z): tempz *= tempz
                        tempx = tempx+tempy+tempz
                        tempt = c*(quants[i].t - quants[tempi].t)
                        tempt *= tempt
                        If (tempt >= tempx) Then
                           'If (currentTime And &hfffff) = 0 Then Locate 1,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                           While (tempt > tempx)
                              tempi+=1000
                              tempx = (quants[i].x - quants[tempi].x): tempx *= tempx
                              tempy = (quants[i].y - quants[tempi].y): tempy *= tempy
                              tempz = (quants[i].z - quants[tempi].z): tempz *= tempz
                              tempx = tempx+tempy+tempz
                              tempt = c*(quants[i].t - quants[tempi].t)
                              tempt *= tempt
                           Wend
                           'If (currentTime And &hfffff) = 0 Then Locate 2,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                           While (tempt < tempx)
                              tempi-=100
                              tempx = (quants[i].x - quants[tempi].x): tempx *= tempx
                              tempy = (quants[i].y - quants[tempi].y): tempy *= tempy
                              tempz = (quants[i].z - quants[tempi].z): tempz *= tempz
                              tempx = tempx+tempy+tempz
                              tempt = c*(quants[i].t - quants[tempi].t)
                              tempt *= tempt
                           Wend
                           'If (currentTime And &hfffff) = 0 Then Locate 3,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                           While (tempt > tempx)
                              tempi+=10
                              tempx = (quants[i].x - quants[tempi].x): tempx *= tempx
                              tempy = (quants[i].y - quants[tempi].y): tempy *= tempy
                              tempz = (quants[i].z - quants[tempi].z): tempz *= tempz
                              tempx = tempx+tempy+tempz
                              tempt = c*(quants[i].t - quants[tempi].t)
                              tempt *= tempt
                           Wend
                           
                           'If (currentTime And &hfffff) = 0 Then Locate 4,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                           If maxtempi < tempi Then maxtempi = tempi
                           rmin = computeEB(@quants[i], @quants[tempi], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                           If (r > rmin) Then r = rmin
                           If (currentTime And &h3ffff) = 0 Then
                              logit_y = logit(quants[ii].y)
                           	logit_x = logit(quants[ii].x) 
                     	      xx = -logit_x * Sinunghi + logit_y * cosunghi + maxxshr
                     			yy = -logit_x*ttp1 - logit_y*ttp2 + logit(quants[ii].z) * sinunghi2 + maxyshr
                     			PSet(xx,yy), color1
                              logit_y = logit(quants[tempi].y)
                           	logit_x = logit(quants[tempi].x) 
                     	      xx = -logit_x * Sinunghi + logit_y * cosunghi + maxxshr
                     			yy = -logit_x*ttp1 - logit_y*ttp2 + logit(quants[tempi].z) * sinunghi2 + maxyshr
                     			PSet(xx,yy), color2
                  			endif
                        Else
                           tempi -= num_particles*tempdelta
                           If tempi < 0 Then
                              rmin = computeEB(@quants[i], @quants[ii], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                              If (r > rmin) Then r = rmin
                           Else
                              GoTo dummy_search
                           EndIf
                        EndIf
                     EndIf
                  End If
                  ii+=1
               Wend
               i += 1
            Wend
            If (r <= 0) Then
               Locate 1,1
               Print "ERROR ";timeQuanta
               pause1
            EndIf
            timeQuanta = (r * 2.0e-6 /(c+quants[i].v)) 
            i = k
            While i < j
               ' spaceQuanta/(c+quants[i].v)
               computeXU(@quants[i], @quants[i+num_particles], timeQuanta, c)
               i+=1
            Wend
            /'
            If (currentTime And &hfffff) = 0 Then 
               i = k               
               Locate 3,1:
      			info_quanta(@quants[i], steps, c,h)
   			EndIf
   			'/
            currentTime+=1
         Loop Until currentTime >= maxTime
         steps += currentTime
      Else 
         Sleep 100
      End If
      
      maxtempi = num_history Shr 1
      ' -- user input ---
      
      a = InKey
      
      Select Case a
         Case "h"
            show_help
            Print "press any key to continue."
            pause1
         Case "s"
            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
            If skip Then
               i = ((0) Shl np) + i
            Else
               i = ((currentTime -1) Shl np) + i
            EndIf
            
            info_quanta(@quants[i], steps, c,h)
            pause1
         Case "I"
            Locate 2,1
            For k = 0 To num_particles - 1 
               i = ((currentTime -1) Shl np) + k
               Print "particle=";k
               info_quanta(@quants[i], steps, c,h)
            Next k
            pause1
         Case "1"
   			unghi-=2
   			sinunghi=Sin(to_rads(90-unghi))
   			cosunghi=cos(to_rads(90-unghi))
   			sinunghi2=sin(to_rads(90-unghi2))
   			cosunghi2=cos(to_rads(90-unghi2))
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print unghi;" ";unghi2; "   ";
   		Case "2"
   			unghi+=2
   			sinunghi=Sin(to_rads(90-unghi))
   			cosunghi=cos(to_rads(90-unghi))
   			sinunghi2=sin(to_rads(90-unghi2))
   			cosunghi2=cos(to_rads(90-unghi2))
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print unghi;" ";unghi2; "   ";
   	   Case "3"
   	      unghi2-=2
   			sinunghi=Sin(to_rads(90-unghi))
   			cosunghi=cos(to_rads(90-unghi))
   			sinunghi2=sin(to_rads(90-unghi2))
   			cosunghi2=cos(to_rads(90-unghi2))
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print unghi;" ";unghi2; "   ";
   	   Case "4"
   	      unghi2+=2
   			sinunghi=Sin(to_rads(90-unghi))
   			cosunghi=cos(to_rads(90-unghi))
   			sinunghi2=sin(to_rads(90-unghi2))
   			cosunghi2=cos(to_rads(90-unghi2))
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print unghi;" ";unghi2; "   ";
   		Case Chr(255)+"H"'up
   			wy-=30
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print wx;",";wy; "      ";
   		Case Chr(255)+"P"'down
   			wy+=30
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print wx;",";wy; "      ";
   		Case Chr(255)+"K"'left
   			wx-=30
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print wx;",";wy; "      ";
   		Case Chr(255)+"M"'right
   			wx+=30
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print wx;",";wy; "      ";
         Case "-"
   			rref=rref*1.1111111
   			'wx=wx\2
   			'wy=wy\2
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print rref; "      ";
   		Case "+"
   			rref=rref/1.1111111
   			'wx=wx*2
   			'wy=wy*2
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print rref; "      ";
   		Case Chr(27)
   			Exit do
   	End Select
   	
   	
   	
   	' graphic output
   
   	
      currentTime = 0
      Do
         k = currentTime Shl np
         j = k + num_particles
         i = k
         i = k
         While i < j
            ' in 2D plane
            'PSet(logit(quants[i].x) + x_opt,y_opt - logit(quants[i].z)), color1
            ' rotated by 1 angle
            'logit_y = logit(quants[i].y) 
            'xx = -logit_y * sin_opt + logit(quants[i].x) + x_opt
         	'yy = -logit_y * cos_opt - logit(quants[i].z) + y_opt
         	'PSet(xx,yy), color1
         	'rotated bt 2 angles
         	logit_y = logit(quants[i].y)
         	logit_x = logit(quants[i].x) 
   	      xx = -logit_x * Sinunghi + logit_y * cosunghi + maxxshr
   			yy = -logit_x*ttp1 - logit_y*ttp2 + logit(quants[i].z) * sinunghi2 + maxyshr
   			PSet(xx,yy), color1
         	i+=1
         Wend 
      	currentTime += 50000
      Loop Until currentTime >= maxtime
      /'
      currentTime = maxtempi Shr np
      Do
         k = currentTime Shl np
         j = k + num_particles
         i = k
         i = k
         While i < j
            ' in 2D plane
            'PSet(logit(quants[i].x) + x_opt,y_opt - logit(quants[i].z)), color1
            ' rotated by 1 angle
            'logit_y = logit(quants[i].y) 
            'xx = -logit_y * sin_opt + logit(quants[i].x) + x_opt
         	'yy = -logit_y * cos_opt - logit(quants[i].z) + y_opt
         	'PSet(xx,yy), color1
         	'rotated bt 2 angles
         	logit_y = logit(quants[i].y)
         	logit_x = logit(quants[i].x) 
   	      xx = -logit_x * Sinunghi + logit_y * cosunghi + maxxshr
   			yy = -logit_x*ttp1 - logit_y*ttp2 + logit(quants[i].z) * sinunghi2 + maxyshr
   			PSet(xx,yy), color3
         	i+=1
         Wend 
      	currentTime += 50000
      Loop Until currentTime >= maxtime
      '/    
      ' prepare for next cycle
      'currentTime = maxTime - 1
      'i = ((currentTime) Shl np ) 
      i = maxtempi 
      For k = 0 To (num_history-1-maxtempi)
         'coordTrans(@quants[i+k], c, 0)
   	   CopyQuanta(@quants[k], @quants[i+k])
      Next k
   
      ' clear for next cycle'(num_history-maxtempi-num_particles) ' (num_history-maxtempi-num_particles-2) 
      For i = 0 To num_history-1 '0 to num_history-1
         quants[i].Ex = 0
         quants[i].Ey = 0
         quants[i].Ez = 0
         quants[i].Bx = 0
         quants[i].By = 0
         quants[i].Bz = 0
      Next i
      

   Loop While 1 = 1 
   DeAllocate quants
 
End Sub

[Mihail_B]

Relativistic Electrodynamics Simulator - genius positions (genius fields), field propagation graphics (drawn after history elements are fully occupied)

(You need a 64bit compiler and at least 3980MB free RAM)
To compile use: [path]\fbc\win64-1.05.0\fbc -s console -t 128000 -w 3 -O 2 -v "relativity_retardedfield.bas"
It takes a lot of time - it's relativistic so lot's of Lorentz Boosts in arbitrary directions = lots of calculus.
I didn't optimize much the genius position calculator so you need to be patient.

This analysis program sample will use the history of each particle to display electric and magnetic field propagation in time from the charges that generated it. We use the framework of special relativity and we treat electric charges as field singularities. We take the upper limit of the potential energy as being at most max(E) = E = (m1+m2) * c^2 (where m1, m2 are the relativistic masses, ex: m1 = restMass(m1) / sqrt(1 - (v^2)/(c^2)). Collisions are treated by them selves by the equations of motions but up to a maximum energy of max(E), otherwise i choose to let the particle freely pass through the other particle since I can not treat here particle annihilation/creation properly.

At approx. line 1093 we define the number of maximum history positions held in memory. Feel free to increase this up to whatever limit of RAM you have, as explained below.
For a 8GB RAM PC (Win 10/7 with no other applications opened) : num_history = 31000000
For a 16GB RAM PC (Win 7 with no other applications opened) : num_history = 95000000

So for 2 particles and a total of 31000000 history elements we get 15,500,000 (15million) positions per particle.
This is barely enough to show full genius field propagation for a screen of 800x600 pixels where the radius of reference is 0.5e-13meters between 2pixels.

I used 95000000 elements with a machine that has 16GB of RAM and a screen resolution of 1280x1024 but still this is not enough to cover all the screen with the appropriate field values..

The field at each pixel is drawn as :

Code: Select all

colorf = RGB(Log(Sqr(temp_quanta.Ex*temp_quanta.Ex+temp_quanta.Ey*temp_quanta.Ey+temp_quanta.Ez*temp_quanta.Ez))*255, Log(Sqr(temp_quanta.Bx*temp_quanta.Bx+temp_quanta.By*temp_quanta.By+temp_quanta.Bz*temp_quanta.Bz))*255,10)

where .Ex, .Ey, .Ez are the components of the electric field vector, and .Bx, .By, .Bz are the components of the magnetic induction field vector (the components are taken from the relativistic electromagnetic tensor after appropriate field transformation takes place), and we display the modulus of the vector (scalar value at that point) otherwise we would have needed too many colors which makes harder to understand the field values.

Every time we complete a cycle of num_history elements we display the graphics and then we cut half of the upper positions and copy them on the lower half and in the next cycle we start from the half of the array. It takes about 10-20 seconds on a 3.2GHz cpu to complete a cycle.

To start the simulation with a specific zoom/scale modify variable rref at approx line 1077

Code: Select all

rref = 0.5e-13

( value is in meters and represents the distance between 2pixels when the angles of view are : angle1 = 206.6, angle2 = 180 )


16GB Ram machine, screen 21,
first cycle:


second cycle

Code: Select all

' WARNING: Don't use with medical applications !

' Relativistic Electrodynamics Simulator (free version, uses genius positions)
' In the framework of Einstein's Special Relativity 
' (at least in the context of what I undestood of it :) )

' Note: Can simulate in the limit of 64bit double numbers. 
'       When too close to the speed of light (like >99.9999%) 64bit numbers are unsatisfactory, 
'       leading to computation errors.


' Tested with FreeBASIC Compiler - Version 1.05.0 (01-31-2016), built for win64 (64bit)
' Note: dont use 32bit compiler as 32bit applications wont be able to allocate enough memory
' PS: you can still compile it with 32bit compiler but you will need to lower num_history value to something lower than 2GB

' To compile use: [path]\fbc\win64-1.05.0\fbc -s console -t 128000 -w 3 -O 2 -v "Relativity_instantfield.bas"
' My website http://alphax.info
' Mihai Barboi, 2017 - open source

' press 'h' during simulation to get a help
' NOTE: Any key command will be processed after the a cycle is complete 
' (aka after num_history/num_particles computations); to better view from other angles press TAB and then wait; the use angle modified keys to render from a different perspective 
' Other keys :
' c = CLS  
' i = info about particle (x,y,z, vx, vy, vz, etc)
' I = info about all particles
' 1 = angle -= 2  
' 2 = angle += 2  
' 3 = angle2 -= 2  
' 4 = angle2 += 2  
' KEY UP = move UP viewport  
' KEY DOWN = move DOWN viewport  
' KEY LEFT = move LEFT viewport  
' KEY RIGHT = move RIGHT viewport  
' - = zoom out (rref *= 2.222222) 
' + = zoom in (rref /= 2.222222)  
' s = toggle simulation run (on/off)


#include "fbgfx.bi"


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
'#Include "mem.bas"
'#Include "easy.bas"

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

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

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)))))

Type vector
   As Double x,y,z,a
End Type

Declare Sub show_help
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 matrixMatrixMult(m1 As matrix4x4 Ptr, m2  As matrix4x4 Ptr, res As matrix4x4 Ptr)
Declare Sub contravariantVelocityLorentzBoost(m1 As matrix4x4 Ptr, c As double, v As vector4 ptr , result As vector4 Ptr)
Declare Sub computeElectromagneticTensorLorentzBoost(m1 As matrix4x4 Ptr, m2  As matrix4x4 Ptr, F As matrix4x4 Ptr, res As matrix4x4 Ptr)
Declare Sub writeBoostMatrix_contravariant(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Declare Sub writeBoostMatrix_covariant(from As quanta Ptr, m As matrix4x4 Ptr, c As Double)
Declare Sub writeBoostMatrixFromVect_contravariant(from As vector4 Ptr, m As matrix4x4 Ptr, c As Double)
Declare Sub velocityAdition(from As quanta ptr, s As quanta Ptr, vect1 As vector4 Ptr, c2 As Double)
Declare Sub writeBoostMatrixFromVect_covariant(from As vector4 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 info_quanta(q as quanta Ptr, steps As ULongInt, c As Double, h As Double)
Declare Sub coordTrans(s As quanta Ptr, c As Double, dt As Double)
Declare Sub lorentzBoostEB(s As quanta ptr, from As quanta Ptr, c As Double)
Declare Sub computeXU(s As quanta Ptr, en As quanta Ptr, dt As Double, c As Double) 
Declare function computeEB(s As quanta Ptr, from As quanta Ptr, epsi_ As Double, epsi_c2 As Double, miudiv4pi As Double, c As Double, ma1 As matrix4x4 Ptr, ma2 As matrix4x4 Ptr, fma As matrix4x4 Ptr, mres As matrix4x4 ptr) As Double
Declare Sub CopyQuanta(s As quanta Ptr, p As quanta Ptr)
Declare function distance(s As quanta Ptr, p As quanta Ptr)  As Double
Declare Sub electrodynamics()

/'
Dim As Double c,v,m0,e,de,vf,c2,h,l10,e0,nf,freq
Dim As Integer i,k,j,x,y,x0,y0
Dim As String a
Dim As Double relativeEpsi, chargeElec, massElec, massProt, epsi, grav, c4
c = 3e8
c2 = c*c
c4=c2*c2
m0=1e-27
h=6.6e-36
relativeEpsi = 1.216e-17
chargeElec = 1.6e-19
massElec  = 9.1e-31
massProt = 1.627e-27
epsi = 8.854e-12
grav = 6.673e-11
'/
Screen 19,32
ScreenInfo maxx,maxy

electrodynamics()

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 "I = info about all particles"
   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

'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 matrixMatrixMult(m1 As matrix4x4 Ptr, m2  As matrix4x4 Ptr, res As matrix4x4 Ptr)
   ' A{nm} : n x m
   ' B{mp} : m x p
   ' (AB){ij} = sum{k=1}[m] A{ik} B{kj}
   Dim As Integer i,k,j
   For i = 0 To 3
      For j = 0 To 3
         res->e(i,j) = 0
         For k = 0 To 3
            res->e(i,j) += m1->e(i, k) * m2->e(k, j)
         Next k
      Next j
   Next i
End Sub
Sub contravariantVelocityLorentzBoost(m1 As matrix4x4 Ptr, c As double, v As vector4 ptr , result As vector4 Ptr)
   Dim As Double m,n
   For m = 0 To 3
      result->e(m) = 0
      For n = 0 To 3
         result->e(m) += m1->e(m,n) * v->e(n)
      Next n
   Next m
   result->e(0) /= c
   result->e(1) /= result->e(0)
   result->e(2) /= result->e(0)
   result->e(3) /= result->e(0)
End Sub
Sub convariantVelocityLorentzBoost(m1 As matrix4x4 Ptr, c As double, v As vector4 ptr , result As vector4 Ptr)
   Dim As Double m,n
   For m = 0 To 3
      result->e(m) = 0
      For n = 0 To 3
         result->e(m) += m1->e(m,n) * v->e(n)
      Next n
   Next m
   result->e(0) /= c
   result->e(1) /= result->e(0)
   result->e(2) /= result->e(0)
   result->e(3) /= result->e(0)
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_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 writeBoostMatrix_covariant(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 writeBoostMatrixFromVect_contravariant(from As vector4 Ptr, m As matrix4x4 Ptr, c As Double)
   Dim As Double bx, by, bz, b2, g
   
   g = from->e(0)
   bx = from->e(1) / c 
   by = from->e(2) / c 
   bz = from->e(2) / c 
   b2 = (bx * bx + by * by + bz * bz)

   
   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 velocityAdition(from As quanta ptr, s As quanta Ptr, vect1 As vector4 Ptr, c2 As Double)
   ' w2 = ((u-v).(u-v) - (u x v)2/c2 )/(1- (u.v)/c2)2
   ' u.v = uxvx+uyvy+uzvz ; (u x v)2 = (u.u)(v.v) - (u.v)2 ; u-v=(ux-vx, uy-vy, uz-vz)
   Dim As Double a1, a2, a3, w2, tx1, ty1, tz1, r, a4
   a1 = (s->vx - from->vx)*(s->vx - from->vx) + (s->vy - from->vy)*(s->vy - from->vy) + (s->vz - from->vz)*(s->vz - from->vz)
   a2 = (from->vx*from->vx + from->vy*from->vy+from->vz*from->vz)*(s->vx*s->vx + s->vy*s->vy+s->vz*s->vz)
   a3 = (from->vx*s->vx + from->vy*s->vy + from->vz*s->vz)
   a4 = (1-a3/c2)
   w2 = (a1 - (a2 - a3*a3)/c2)/(a4*a4)
   
   tx1 = s->x - from->x
	ty1 = s->y - from->y
	tz1 = s->z - from->z
	r = Sqr((tx1 * tx1) + (ty1 * ty1) + (tz1 * tz1))
	
	w2 = w2/r
	vect1->e(1) = w2 * tx1
	vect1->e(2) = w2 * ty1
	vect1->e(3) = w2 * tz1 

   
End Sub
Sub writeBoostMatrixFromVect_covariant(from As vector4 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 = from->e(0)
   bx = from->e(1) / c 
   by = from->e(2) / c 
   bz = from->e(2) / c 
   b2 = (bx * bx + by * by + bz * bz)
   
   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 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 info_quanta(q as quanta Ptr, steps As ULongInt, c As Double, h As Double)
   Print "m ";q->m;" q ";q->q; " steps ";steps;"  "
   Print "x ";q->x;" ";q->y;" ";q->z;" t ";q->t;"  "
   Print "a ";q->a;" V ";q->v; " ct ";c*q->t;"  "
   Print "v ";q->vx;" ";q->vy;" ";q->vz;"  "
   Print "E ";q->Ex;" ";q->Ey;" ";q->Ez;"  "
   Print "B ";q->Bx;" ";q->By;" ";q->Bz;"  "
   Print "K ";q->m*c*c*((1/q->a) - 1)/h;" Hz => " ;Log(q->m*c*c*((1/q->a) - 1)/h)/Log(10);" (log)Hz"
   Print 
End Sub
Sub coordTrans(s As quanta Ptr, c As Double, dt As Double) 
   ' L{ik} = |  g      -Bxg                 -Byg                -Bzg               | |ct
   '         | -Bxg    1 + (g - 1)Bx2/B2    (g - 1)BxBy/B2      (g - 1)BxBz/B2     | | x
   '         | -Byg    (g - 1)ByBx/B2       1 + (g - 1)By2/B2   (g - 1)ByBz/B2     | | y
   '         | -Bzg    (g - 1)BzBx/B2       (g - 1)BzBy/B2       1 + (g - 1)Bz2/B2 | | z
   Dim As Double g, bx, by, bz, b, b2
   Dim As Double x, y, z, t, ct
   g = 1 / s->a
   bx = s->vx / c
   by = s->vy / c
   bz = s->vz / c
   b = s->v / c
   b2 = b * b
   If (b = 0) Then return
   ct = s->t * c
   t = g*ct - bx*g*s->x - by*g*s->y - bz*g*s->z
   x = -Bx*g*ct + (1+(g-1)*bx*bx/b2)*s->x + (g-1)*bx*by/b2*s->y + (g-1)*bx*bz/b2*s->z
   y = -by*g*ct + (g-1)*by*bz/b2*s->x + ((g-1)*by*bz/b2)*s->y + (g-1)*by*bz/b2*s->z
   z = -bz*g*ct + (g-1)*bz*bx/b2*s->x + (g-1)*bz*by/b2*s->y + (g-1)*bz*bz/b2*s->z
   
   't = ( 1 - Bx  - by  - bz) * g * ct
   'x = (-bx * g +  1+(g-1)*bx*bx/b2 + (g-1)*bx*by/b2 + (g-1)*bx*bz/b2) * s->x
   'y = (-by * g +  (g-1)*by*bz/b2 + 1 + (g-1)*by*by/b2 + (g-1)*by*bz/b2) * s->y
   'z = (-bz * g +  (g-1)*bz*bx/b2 + (g-1)*bz*by/b2 + 1 + (g-1)*bz*bz/b2) * s->z 
   
   
   
   Locate 1,1
   Print s->x;" ";s->y;" ";s->z;" ";s->t;"   "
   Print x; " ";y; " "; z; " ";t/c; " "
   's->x = x
   's->y = y
   's->z = z
   's->t += dt 't/c
   
   pause1
   
End Sub

Sub lorentzBoostEB(s As quanta ptr, from As quanta Ptr, c As Double)
   Dim As Double ex, ey, ez, bx, by, bz
   Dim As Double matrix(0 To 3,0 To 3)
   ' F{ab} = |  0     Ex/c  Ey/c  Ez/c |  
   '         | -Ex/c     0   -Bz    By | 
   '         | -Ey/c    Bz     0   -Bx | 
   '         | -Ez/c   -By    Bx     0 | 
   '
   ' g = 1/a = 1/sqr(1-v2/c2)
   ' B{k} = v{k}/c
   ' F'{mn} = Lmi Fik Lnk* ?
   ' F'{mn} = L{m}[i] L{n}[k] F{ik}
   
   ' L{ik} = |  g      -Bxg                 -Byg                -Bzg               |
   '         | -Bxg    1 + (g - 1)Bx2/B2    (g - 1)BxBy/B2      (g - 1)BxBz/B2     |
   '         | -Byg    (g - 1)ByBz/B2       1 + (g - 1)By2/B2   (g - 1)ByBz/B2     |
   '         | -Bzg    (g - 1)BzBx/B2       (g - 1)BzBy/B2       1 + (g - 1)Bz2/B2 |
   
   ' L{ik} = |  1/a    -vx/ca                 -vy/ca                -vz/ca             |
   '         | -vx/ca  1 + (1/a - 1)vx2/v2    (1/a - 1)vxvy/v2      (1/a - 1)vxvz/v2   |
   '         | -vy/ca  (1/a - 1)vyvz/v2       1 + (1/a - 1)vy2/v2   (1/a - 1)vyvz/v2   |
   '         | -vz/ca  (1/a - 1)vzvx/v2       (1/a - 1)vzvy/v2      1 + (1/a -1)vz2/v2 |
   
   ' F'{ik} = L[i]{a} L[k]{b} F{ab}
   ' F{00} = L[0]{a} L[0]{b} F{ab} =   
   
   ' [vx, vy, vz] : velocity relative to  frame 2 from frame 1; speed composition
   ' but we use "vacuum" and we are considering field in general; no need velocity transformation
   
   ' F'{mn} = deron(x[a])/deron(x[m]) deron(x[b])/deron(x[n]) F{ab} 
   ' ^^^^^^^^^^ much better
   ' X{1} = dX[1]/dY[0] Y{0} + dX[1]/dY[1] Y{1} + dX[1]/dY[2] Y{2} + dX[1]/dY[3] Y{3}
   ' x = vx/a * T + 1/a X + dx/dy Y + dx/dz Z 
   ' x = (vx*T + X)/a + vx Y/vya + vx Z/vza
   ' x = (vx*T + X + vx*Y/vy + vx*z/vz)/a
   ' 
   ' matrix(0,0) = 
   
End Sub
' computes new X4 (position) and U4 (velocity)
' en = (ending) = final
' s = (source) = initial
Sub computeXU(s As quanta Ptr, en As quanta Ptr, dt As Double, c As Double) 
   'all up
   ' (0)  m * c /a' = e / (a*c) * (Ex*vx+ Ey*vy + Ez*vz) * dtau + m*c/a = (D)
   ' (1)  m * vx' / a' = e/a * (Ex + Bz*vy - By*vz) * dtau + m*vx/a = (A)
   ' (2)  m * vy' / a' = e/a * (Ey - Bz*vx + Bx*vz) * dtau + m*vy/a = (B)
   ' (3)  m * vz' / a' = e/a * (Ez + By*vx - Bx*vy) * dtau + m*vz/a = (D)
   
   ' v^2 = vx^2 + vy^2 + vz^2   
   ' A/vx' = B/vy' = C/vz' = D/c = m/a' 
   ' =>
   ' A/vx' = D/c   <=> vx' = A*c/D 
   ' B/vy' = D/c   <=> vy' = B*c/D
   ' C/vz' = D/c   <=> vz' = C*c/D
   ' D/c = m/a'    <=> a' = m*c/D
   ' a' = sqr(1 - v^2/c^2) = sqr(1 - (vx'^2 + vy'^2 + vz'^2)/c^2)
   
   Dim As Double AA, BB, CC, DD
   Dim As Double dtau, temp, temp2
   
   /'
   dtau = dt * s->a ' dtau = tau1 - tau0 = t1*a1 - t0*a0 = dt
     
   DD = (s->q / (s->a * c)) * ( s->Ex * s->vx + s->Ey * s->vy + s->Ez * s->vz) * dtau + s->m * c / s->a
   
   AA = (s->q / s->a) * (s->Ex + s->Bz * s->vy - s->By * s->vz) * dtau + s->m * s->vx / s->a
   
   BB = (s->q / s->a) * (s->Ey - s->Bz * s->vx + s->Bx * s->vz) * dtau + s->m * s->vy / s->a
   
   CC = (s->q / s->a) * (s->Ez + s->By * s->vx - s->Bx * s->vy) * dtau + s->m * s->vz / s->a
   '/

   temp2 = s->q * dt
   temp = s->m / s->a
     
   DD = (temp2/c) * ( s->Ex * s->vx + s->Ey * s->vy + s->Ez * s->vz) + c * temp
   
   AA = temp2 * (s->Ex + s->Bz * s->vy - s->By * s->vz) + s->vx * temp
   
   BB = temp2 * (s->Ey - s->Bz * s->vx + s->Bx * s->vz) + s->vy * temp
   
   CC = temp2 * (s->Ez + s->By * s->vx - s->Bx * s->vy) + s->vz * temp
   
   en->m = s->m
   en->q = s->q
   en->t = s->t + dt
   
   temp = c / DD
   en->vx = AA * temp ' = AA * c / DD
   en->vy = BB * temp ' = BB * c / DD
   en->vz = CC * temp ' = CC * c / DD
   'en->a = s->m * temp ' = s->m * c / DD
   
   en->v = Sqr(en->vx * en->vx + en->vy * en->vy + en->vz * en->vz)  
   
   temp = en->v / c
   en->a = Sqr(1 - temp * temp)
   
   'coordTrans(en, c, dt)
   
   ' compute new position relative to our "observer"
   
   
   'Return
   ' temp = dt/en->a
   
   en->x = s->x + (s->vx + en->vx)/2 * dt ' temp
   en->y = s->y + (s->vy + en->vy)/2 * dt ' temp
   en->z = s->z + (s->vz + en->vz)/2 * dt ' temp
   
   /'
   If (en->a > 1) Then
      Locate 1,1
      Print "D";dd; " A";aa; " B"; bb; " C";cc;" a";en->a; " q";en->q
      pause1
   End If
   '/
    
   /'
   en->Bx = 0
   en->By = 0
   en->Bz = 0
   en->Ex = 0
   en->Ey = 0
   en->Ez = 0
   '/
End Sub

' from any particle in radius r = - ct (t = light time, field change propagation speed)
' s = final
' from = initial
function computeEB(s As quanta Ptr, from As quanta Ptr, epsi_ As Double, epsi_c2 As Double, miudiv4pi As Double, c As Double, ma1 As matrix4x4 Ptr, ma2 As matrix4x4 Ptr, fma As matrix4x4 Ptr, mres As matrix4x4 ptr) As double
   Dim As Double tx1, ty1, tz1, bx, by, bz, fx, fy, fz
   Dim As Double d, f, Rmax, r2
   
   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
	   f = from->q / (epsi_* r2 * d) ' r^3
	   fx = f * tx1 
	   fy = f * ty1 
	   fz = f * tz1 
	   
      f = miudiv4pi * from->q / (r2 * d)
		bx = f * (from->vy * tz1 - from->vz * ty1)  
		by = f * (from->vz * tx1 - from->vx * tz1)
		bz = f * (from->vx * ty1 - from->vy * tx1)
		
		If (from->v > 0) Then
	      writeF_mu_nu(Fma, fx, fy, fz, bx, by, bz, c)
         writeBoostMatrix_covariant(from, ma1, c)
         writeBoostMatrix_contravariant(from, ma2, c)
	      
	      computeElectromagneticTensorLorentzBoost(ma1, ma2, Fma, mres)
	      
	      fx = -mres->e(0,1)*c
	      fy = -mres->e(0,2)*c
	      fz = -mres->e(0,3)*c
	      bx = mres->e(3,2)
	      by = mres->e(1,3)
	      bz = mres->e(2,1)
		End If
		
		/'
		If (1 = 0 And s->v > 0) Then
	      writeF_mu_nu(Fma, fx, fy, fz, bx, by, bz, c)
         writeBoostMatrix_covariant_reverse(s, ma1, c) '_reverse
         writeBoostMatrix_contravariant_reverse(s, ma2, c) '_reverse
	      
	      computeElectromagneticTensorLorentzBoost(ma1, ma2, Fma, mres)
	      
	      fx = -mres->e(0,1)*c
	      fy = -mres->e(0,2)*c
	      fz = -mres->e(0,3)*c
	      bx = mres->e(3,2)
	      by = mres->e(1,3)
	      bz = mres->e(2,1)
		EndIf
		'/
		s->Ex += fx
      s->Ey += fy
      s->Ez += fz 
      
		s->Bx += bx
      s->By += by
      s->Bz += bz
      
      ' we have to apply a general boost transformation to (fx,fy,fz, bx,by,bz) for (vx,vy,vz) at new coords 
      
	Else  ' review this
	   Return rmax 
	End If
   
    
   '             | i  j  k  |			| i  j  k  |
	' a x b = det | a1 a2 a3 | = det | x1 y1 z1 |
	'				  | b1 b2 b3 |			| x2 y2 z2 |
	' a x b = (a2*b3 - a3*b2)i + (a3*b1 - a1*b3)j + (a1*b2 - a2*b1)k
   Return d
End Function

Sub CopyQuanta(s As quanta Ptr, p As quanta Ptr)
   s->x = p->x
   s->y = p->y
   s->z = p->z
   s->vx = p->vx
   s->vy = p->vy
   s->vz = p->vz
   s->q = p->q
   s->m = p->m
   s->t = p->t
   s->Ex = p->Ex
   s->Ex = p->Ey
   s->Ez = p->Ez
   s->Bx = p->Bx
   s->By = p->By
   s->Bz = p->Bz
   s->a = p->a
   s->v = p->v
End Sub
function distance(s As quanta Ptr, p As quanta Ptr)  As double
   Dim As Double dr, temp
   dr = Abs(s->x - p->x) 
   dr = dr * dr
   temp = Abs(s->y - p->y)
   temp = temp*temp 
   dr += temp
   temp = abs(s->z - p->z)
   temp = temp*temp
   dr += temp
   dr = Sqr(dr)
   return dr
End Function
Sub fields_graphics(rref As Double, Sinunghi As Double, cosunghi As Double, ttp1 As Double, ttp2 As Double, _
   maxxshr As double, maxyshr As double, e as double, m_prot as double, quants As quanta Ptr, num_history as ULongInt, _
   sinunghi2 as double, num_particles as integer, tempdelta as integer, epsi_ as double,  tempi as integer, tempx as double, tempy as double, tempz as double, _
   np as integer, epsi_c2 as double, tempt as Double, miudiv4pi as double, c as double, maxTime as Integer, _
   byref ma1 As matrix4x4, byref ma2 As matrix4x4, byref fma As matrix4x4, byref mres As matrix4x4)
   Dim As Integer scr_x, scr_y, scr_z
   Dim As Integer i_p, k_p
   Dim As UInteger colorf
   Dim As Double x5,y5,z5
   Dim As Double logit_const, logit_z2, const_stuff_x, const_stuff_y, const_1, const_2
   Dim As ULongInt ok = 0, bad = 0, over = 0, same = 0, Step1 = 0, Step2 = 0, Step3 = 0
   Dim As quanta temp_quanta
  
               logit_const = (10 ^ Abs((Log(rref)/Log(10))))
               
               scr_z = 0
               z5 = scr_z / logit_const   
               logit_z2 = logit(z5)
               
               const_stuff_x = ((Sinunghi * ttp2 + ttp1*cosunghi) * logit_const)
               const_stuff_y = ((cosunghi * ttp1 + ttp2*sinunghi) * logit_const)
               const_1 = maxyshr - logit_z2 * sinunghi2
               
               
               temp_quanta.q = e
               temp_quanta.m = m_prot
               temp_quanta.a = 0.5
               temp_quanta.t = quants[num_history - num_particles].t
               For scr_y = 0 To maxy-1
                  For scr_x = 0 To maxx - 1
            			x5 = -(ttp2*(scr_x - maxxshr) + cosunghi*(scr_y - const_1)) / const_stuff_x
            			y5 =  (ttp1*(scr_x - maxxshr) - sinunghi*(scr_y - const_1)) / const_stuff_y
            			' field at (x,y,z)
            			   temp_quanta.x = x5
            			   temp_quanta.y = z5
            			   temp_quanta.z = -y5
            			   temp_quanta.Ex = 0
            			   temp_quanta.Ey = 0
            			   temp_quanta.Ez = 0
            			   temp_quanta.Bx = 0
            			   temp_quanta.by = 0
            			   temp_quanta.Bz = 0
            			   
            		      k_p = num_history - num_particles' should start from last particle state ;
            		      While k_p < num_history
                           /'If (k_p < (tempdelta Shl np)) Then 
                              computeEB(@temp_quanta, @quants[k_p And np], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                              direct +=1
                           Else'/
                              tempi = k_p
                              dummy_search2:
                              tempx = (temp_quanta.x - quants[tempi].x): tempx *= tempx
                              tempy = (temp_quanta.y - quants[tempi].y): tempy *= tempy
                              tempz = (temp_quanta.z - quants[tempi].z): tempz *= tempz
                              tempx = tempx+tempy+tempz
                              tempt = c*(temp_quanta.t - quants[tempi].t)
                              tempt *= tempt
                              If (tempt >= tempx) Then
                                 'If (currentTime And &hfffff) = 0 Then Locate 1,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                                 While (tempt > tempx And tempi < num_history)
                                    tempi+=1000
                                    tempx = (temp_quanta.x - quants[tempi].x): tempx *= tempx
                                    tempy = (temp_quanta.y - quants[tempi].y): tempy *= tempy
                                    tempz = (temp_quanta.z - quants[tempi].z): tempz *= tempz
                                    tempx = tempx+tempy+tempz
                                    tempt = c*(temp_quanta.t - quants[tempi].t)
                                    tempt *= tempt
                                    step1+=1
                                 Wend
                                 'If (currentTime And &hfffff) = 0 Then Locate 2,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                                 While (tempt < tempx And tempi > 0)
                                    tempi-=100
                                    tempx = (temp_quanta.x - quants[tempi].x): tempx *= tempx
                                    tempy = (temp_quanta.y - quants[tempi].y): tempy *= tempy
                                    tempz = (temp_quanta.z - quants[tempi].z): tempz *= tempz
                                    tempx = tempx+tempy+tempz
                                    tempt = c*(temp_quanta.t - quants[tempi].t)
                                    tempt *= tempt
                                    step2+=1
                                 Wend
                                 'If (currentTime And &hfffff) = 0 Then Locate 3,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                                 While (tempt > tempx And tempi < num_history)
                                    tempi+=10
                                    tempx = (temp_quanta.x - quants[tempi].x): tempx *= tempx
                                    tempy = (temp_quanta.y - quants[tempi].y): tempy *= tempy
                                    tempz = (temp_quanta.z - quants[tempi].z): tempz *= tempz
                                    tempx = tempx+tempy+tempz
                                    tempt = c*(temp_quanta.t - quants[tempi].t)
                                    tempt *= tempt
                                    step3+=1
                                 Wend
                                 If tempi < 0 Then
                                    computeEB(@temp_quanta, @quants[k_p And np], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                                    bad+=1
                                 ElseIf tempi >= num_history Then
                                    computeEB(@temp_quanta, @quants[num_history - num_particles + (k_p And np)], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                                    bad+=1
                                 Else 
                                    computeEB(@temp_quanta, @quants[tempi], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                                    If tempi <> k_p Then
                                       ok+=1
                                    Else
                                       same+=1
                                    EndIf
                                 EndIf
                              Else
                                 tempi -= num_particles * tempdelta
                                 If tempi < 0 Then
                                    computeEB(@temp_quanta, @quants[k_p And np], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                                    over+=1
                                 Else
                                    GoTo dummy_search2
                                 EndIf
                              EndIf
                           'EndIf
            		         k_p += 1
            		      Wend
            		   
            			colorf = RGB(Log(Sqr(temp_quanta.ex*temp_quanta.ex+temp_quanta.ey*temp_quanta.ey+temp_quanta.ez*temp_quanta.ez))*255, _
            			             Log(Sqr(temp_quanta.bx*temp_quanta.bx+temp_quanta.by*temp_quanta.by+temp_quanta.bz*temp_quanta.bz))*255,10)
            			             
   		             ' colorf = RGB(10, 10, Log(Sqr(temp_quanta.bx*temp_quanta.bx+temp_quanta.by*temp_quanta.by+temp_quanta.bz*temp_quanta.bz))*255)
            			PSet (scr_x,scr_y), colorf
                  Next scr_x
               Next scr_y
               
   			
   			Locate 1,1:
   			Print "press ok=";ok;" bad=";bad;" over=";over;" same=";same;" step1=";step1;" step2=";step2;" step3=";step3;
   			
   			pause1
 
   			   
End Sub

Sub electrodynamics
' From covariant formulation of electromagnetism we have
   ' p{a}/dtau = e * F{ab} * u[b]
   ' p[a]/dtau = e * F[a]{b} *u[b]
   ' {..} = covariant (lower index} ; [..] = contravariant (upper index) 
   ' a = sqr(1 - v^2/c^2),
   ' gamma = 1/a , 
   ' m = m0 (rest mass), 
   ' e = charge,  
   ' c = speed of light 3e8 m/s
   ' dtau = dt * a, (proper time)
   ' v = dx / dt
   ' v^2 = vx^2 + vy^2 + vz^2
   ' Ex, Ey, Ez = electric field 
   ' Bx, By, Bz = magnetic field
   ' E = mc^2 (energy, in Joules)
   ' ^ = power (takes highest precedence here)
   
   ' under Minkowski metric :
   ' (proper time dtau^2 = dt^2 - (1/(c^2))*(dx^2 + dy^2 + dz^2) 
   ' (proper distance dsigma^2 = -(c*t)^2 + dx^2 + dy^2 + dz^2)
   
   ' x[a] = distance/position etc in [m];
   ' u[a] = gamma(c, u)  ; u[a] velocity in [m/s]
   ' p{a} = m0 * u[a] ; p = impulse
   ' F{ab} = electromagnetic tensor
   ' ita = [-1, 1, 1, 1]
   
   ' x[0] = c*t , x[1] = x , x[2] = y , x[3] = z
   ' u[0] = c/a ; u[1] = vx/a ; u[1] = vy/a ; u[2] = vz/a
   ' p[0] =  m*c/a ; p[1] = m*vx/a ; p[2] = m*vy/a ; p[3] = m*vz/a
   ' p{0} = -m*c/a ; p{1} = m*vx/a ; p{2} = m*vy/a ; p{3} = m*vz/a
   ' F[ab] = |  0    -Ex/c -Ey/c -Ez/c |  
   '         |  Ex/c     0   -Bz    By | 
   '         |  Ey/c    Bz     0   -Bx | 
   '         |  Ez/c   -By    Bx     0 |
   ' F{ab} = |  0     Ex/c  Ey/c  Ez/c |  
   '         | -Ex/c     0   -Bz    By | 
   '         | -Ey/c    Bz     0   -Bx | 
   '         | -Ez/c   -By    Bx     0 |
   ' 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 | 
   
   ' @(x,y,z,t) [vx, vy, vz, Ex(t), Ey(t), Ez(t), Bx(t), By(t), Bz(t)]
   ' stuff : (x,y,z,t), vx, vy, vz, m, q
   ' traveling with c, from p(i)(x,y,z,t) to p(me)(x', y', z', t');  Ex(t), Ey(t), Ez(t), Bx(t), By(t), Bz(t)


      
   Dim As Double c = 299792458
   Dim As Double c2 = c*c
   Dim As Double c4 = c2*c2
   Dim As Double m0 = 1e-27
   Dim As Double h = 6.6e-36
   Dim As Double e = 1.6e-19
   Dim As Double m_elec  = 9.1e-31
   Dim As Double m_prot = 1.627e-27
   Dim As Double m_up_quark = 4.08888e-30 'up quark +2/3e ('2.3 1e6  1.6e-19 1/9 1e-16) 29 ' radius(u-u) Rmax= 0.90097 e-16 
   dim as double m_down_quark = 8.53333e-30 'down quark -1/3e  (u+u+d = 16.71109e-30 kg)
   Dim As Double miu = 4*pi*1.0e-7
   Dim As Double epsi = 8.854e-12
   Dim As Double grav = 6.673e-11
   Dim As Double pi2 = pi * 2
   Dim As Double miudiv4pi = miu / (4*pi)
   Dim As Double epsi_ = 4*pi*epsi
   Dim As Double epsi_c2 = (epsi_) * c2
   
   Dim As Double var_theta, var_phi, sin1, sin2, x, y, z, phi2, theta2, real_theta, real_phi
   Dim As Double unghi = 25, unghi2 = 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 sinunghi, cosunghi, sinunghi2, cosunghi2, ttp1, ttp2, maxxshr, maxyshr ' 3d graphics translation
   Dim As UInteger color1 = RGB(155,255,255)
   Dim As UInteger color2 = RGB(255,055,055)
   Dim As UInteger color3 = RGB(055,255,055)
   Dim As ULongInt num = 0, steps = 0, skip = 0
   
   Dim As Double s
   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, spaceQuanta
   Dim As Integer cont = 0
   Dim As String a
   Dim As Double startTime, universalTime, rmin
   Dim As Integer currentTime, np, maxTime
   Dim As Double tempt, tempx, tempy, tempz
   Dim As Integer tempi, tempk, maxtempi, tempdelta
   
   Dim As matrix4x4 ma1, ma2, fma, mres
   Dim As quanta Ptr quants
   Dim As Integer num_history, size, i, k, j, ii, kk, jj
   
   
   
   unghi = 206.6
   unghi2 = 180
   rref = 0.5e-13'0.5e-14
   sinunghi=Sin(to_rads(90-unghi))
	cosunghi=cos(to_rads(90-unghi))
	sinunghi2=sin(to_rads(90-unghi2))
	cosunghi2=cos(to_rads(90-unghi2))
   x_opt = CInt(maxx Shr 1)+wx
   y_opt = CInt(maxy Shr 1)+wy  
    
   ttp1=cosunghi*cosunghi2
	ttp2=sinunghi*cosunghi2
	maxxshr=maxx Shr 1
	maxyshr=maxy Shr 1
	maxxshr+=wx
	maxyshr+=wy
   
   tempdelta = 17000
   num_history = 31000000 ' 10000000 ' 30000000 
   '  8GB RAM, Win 10/7 : num_history = 31000000 ' (with no other applications opened)
   ' 16GB RAM, Win 7: num_history = 95000000 ' (with no other applications opened)
   quants = CPtr(quanta Ptr, Allocate(num_history * SizeOf(quanta)))
   For i = 0 To num_history-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 = num_history * SizeOf(quanta) / (1024*1024)
   Print "allocated ";size;"MB ";IIf(quants <> 0, "yes", "no")
   
   Dim As Integer iu, num_particles
   
   spaceQuanta = 0.9e-18
   timeQuanta = timeQuanta/c 'sec ct = 1e-19 't = 0.3333333e-27
   Cls
   
   ' e+ , e-
   num_particles = 2 '4
   np = 1'2
   quants[0].q = -e
   quants[0].m = m_elec
   quants[0].t = 0
   quants[0].x = 0.52e-12     ' good value 1e-12
   quants[0].z = 0.52e-12
   quants[0].y = 0     ' good value 1e-12     
   quants[0].vx = -6.7e6
   quants[0].vy = 0 '7e6
   quants[0].vz = +6.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 'm_elec ' m_prot
   quants[1].vx = 6.7e6
   quants[1].vy = 0
   quants[1].vz = -6.7e6
   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)
 
   
   
   ' time{start} = dr/c 
   startTime = distance(@quants[0], @quants[1])/c
   universalTime = startTime
   maxtempi = 0
   
   Do
      ' computations
      If (maxtempi > 0) Then 
         currentTime = ((num_history-1-maxtempi) Shr np)
         ' Locate 2,1: Print "currentTime ";currentTime;" ";maxtempi;"      ";
         maxtempi = 0
      Else
         currentTime = 0
         maxtempi = 0
      EndIf
      maxTime = (num_history Shr np) - np
      If (skip = 0) then
         Do
            'If (currentTime And &hfffff) = 0 Then Locate 2,1: Print currentTime;"   ";
            k = currentTime Shl np
            j = k + num_particles
            i = k
            r = 10000         
            While i < j
               ii = k
               jj = j
               While ii < jj
                  If (ii <> i) Then
                     ' go back in time and compute to quants[i] from quants[ii] 
                     ' where ii is the appropiate moment in time when c * dt = dr 
                     ' (dr distance between the 2 charges) ; dr = rmin/c
                     ' choose ii such that (quants[i].t - quants[ii].t)*c = sqrt((quants[i].x - quants[ii].x)^2+(quants[i].y - quants[ii].y)^2+(quants[i].z - quants[ii].z)^2)
                     
                     /' If (ii < tempdelta Shl np) Then 
                        rmin = computeEB(@quants[i], @quants[ii And np], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                        If (r > rmin) Then r = rmin
                     Else'/
                        tempi = ii
                        dummy_search:
                        tempx = (quants[i].x - quants[tempi].x): tempx *= tempx
                        tempy = (quants[i].y - quants[tempi].y): tempy *= tempy
                        tempz = (quants[i].z - quants[tempi].z): tempz *= tempz
                        tempx = tempx+tempy+tempz
                        tempt = c*(quants[i].t - quants[tempi].t)
                        tempt *= tempt
                        If (tempt >= tempx) Then
                           'If (currentTime And &hfffff) = 0 Then Locate 1,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                           While (tempt > tempx And tempi < num_history)
                              tempi+=1000
                              tempx = (quants[i].x - quants[tempi].x): tempx *= tempx
                              tempy = (quants[i].y - quants[tempi].y): tempy *= tempy
                              tempz = (quants[i].z - quants[tempi].z): tempz *= tempz
                              tempx = tempx+tempy+tempz
                              tempt = c*(quants[i].t - quants[tempi].t)
                              tempt *= tempt
                           Wend
                           'If (currentTime And &hfffff) = 0 Then Locate 2,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                           While (tempt < tempx And tempi > 0)
                              tempi-=100
                              tempx = (quants[i].x - quants[tempi].x): tempx *= tempx
                              tempy = (quants[i].y - quants[tempi].y): tempy *= tempy
                              tempz = (quants[i].z - quants[tempi].z): tempz *= tempz
                              tempx = tempx+tempy+tempz
                              tempt = c*(quants[i].t - quants[tempi].t)
                              tempt *= tempt
                           Wend
                           'If (currentTime And &hfffff) = 0 Then Locate 3,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                           While (tempt > tempx And tempi < num_history)
                              tempi+=10
                              tempx = (quants[i].x - quants[tempi].x): tempx *= tempx
                              tempy = (quants[i].y - quants[tempi].y): tempy *= tempy
                              tempz = (quants[i].z - quants[tempi].z): tempz *= tempz
                              tempx = tempx+tempy+tempz
                              tempt = c*(quants[i].t - quants[tempi].t)
                              tempt *= tempt
                           Wend
                           
                           'If (currentTime And &hfffff) = 0 Then Locate 4,1:Print tempt;" >= ";tempx;" ? ";ii-tempi;"       ";
                           If tempi < 0 Then
                              rmin = computeEB(@quants[i], @quants[0 +(ii And np)], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                           ElseIf tempi >= num_history Then
                              rmin = computeEB(@quants[i], @quants[num_history - num_particles+(ii And np)], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                           Else
                              If maxtempi < tempi Then maxtempi = tempi
                              rmin = computeEB(@quants[i], @quants[tempi], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                              If (r > rmin) Then r = rmin
                              If (currentTime And &h3ffff) = 0 Then
                                 logit_y = logit(quants[ii].y)
                              	logit_x = logit(quants[ii].x) 
                        	      xx = -logit_x * Sinunghi + logit_y * cosunghi + maxxshr
                        			yy = -logit_x*ttp1 - logit_y*ttp2 + logit(quants[ii].z) * sinunghi2 + maxyshr
                        			PSet(xx,yy), color1
                                 logit_y = logit(quants[tempi].y)
                              	logit_x = logit(quants[tempi].x) 
                        	      xx = -logit_x * Sinunghi + logit_y * cosunghi + maxxshr
                        			yy = -logit_x*ttp1 - logit_y*ttp2 + logit(quants[tempi].z) * sinunghi2 + maxyshr
                        			PSet(xx,yy), color2
                     			EndIf
                           EndIf
                           
                        Else
                           tempi -= num_particles*tempdelta
                           If tempi < 0 Then
                              rmin = computeEB(@quants[i], @quants[ii And np], epsi_, epsi_c2, miudiv4pi, c, @ma1, @ma2, @fma, @mres)
                              If (r > rmin) Then r = rmin
                           Else
                              GoTo dummy_search
                           EndIf
                        EndIf
                     'EndIf
                  End If
                  ii+=1
               Wend
               i += 1
            Wend
            If (r <= 0) Then
               Locate 1,1
               Print "ERROR ";timeQuanta
               pause1
            EndIf
            timeQuanta = (r * 2.0e-6 /(c+quants[i].v)) 
            i = k
            While i < j
               ' spaceQuanta/(c+quants[i].v)
               computeXU(@quants[i], @quants[i+num_particles], timeQuanta, c)
               i+=1
            Wend
            /'
            If (currentTime And &hfffff) = 0 Then 
               i = k               
               Locate 3,1:
      			info_quanta(@quants[i], steps, c,h)
   			EndIf
   			'/
   			
   			/'
   			If (currentTime > ((num_history-1-maxtempi) Shr np) And ((currentTime And &hfffff) = 0)) Then
      			fields_graphics(rref, Sinunghi, cosunghi, ttp1, ttp2, _
                  maxxshr, maxyshr, e, m_prot, quants, currentTime Shl 1, _
                  sinunghi2, num_particles, tempdelta, epsi_, tempi, tempx, tempy, tempz, _
                  np, epsi_c2, tempt, miudiv4pi, c, currentTime, _
                  ma1, ma2, fma, mres)
   			EndIf
   			'/
            currentTime+=1
         Loop Until currentTime >= maxTime
         steps += currentTime
      Else 
         Sleep 100
      End If
      
      maxtempi = num_history Shr 1
      ' -- user input ---
      
      a = InKey
      
      If a = "" Then a= "f" ' if no key is pressed just show field lines
      
      Select Case a
         Case "s"
            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
            If skip Then
               i = ((0) Shl np) + i
            Else
               i = ((currentTime -1) Shl np) + i
            EndIf
            
            info_quanta(@quants[i], steps, c,h)
            pause1
         Case "I"
            Locate 2,1
            For k = 0 To num_particles - 1 
               i = ((currentTime -1) Shl np) + k
               Print "particle=";k
               info_quanta(@quants[i], steps, c,h)
            Next k
            pause1
         Case "1"
   			unghi-=2
   			sinunghi=Sin(to_rads(90-unghi))
   			cosunghi=cos(to_rads(90-unghi))
   			sinunghi2=sin(to_rads(90-unghi2))
   			cosunghi2=cos(to_rads(90-unghi2))
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print unghi;" ";unghi2; "   ";
   		Case "2"
   			unghi+=2
   			sinunghi=Sin(to_rads(90-unghi))
   			cosunghi=cos(to_rads(90-unghi))
   			sinunghi2=sin(to_rads(90-unghi2))
   			cosunghi2=cos(to_rads(90-unghi2))
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print unghi;" ";unghi2; "   ";
   	   Case "3"
   	      unghi2-=2
   			sinunghi=Sin(to_rads(90-unghi))
   			cosunghi=cos(to_rads(90-unghi))
   			sinunghi2=sin(to_rads(90-unghi2))
   			cosunghi2=cos(to_rads(90-unghi2))
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print unghi;" ";unghi2; "   ";
   	   Case "4"
   	      unghi2+=2
   			sinunghi=Sin(to_rads(90-unghi))
   			cosunghi=cos(to_rads(90-unghi))
   			sinunghi2=sin(to_rads(90-unghi2))
   			cosunghi2=cos(to_rads(90-unghi2))
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print unghi;" ";unghi2; "   ";
   		Case Chr(255)+"H"'up
   			wy-=30
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print wx;",";wy; "      ";
   		Case Chr(255)+"P"'down
   			wy+=30
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print wx;",";wy; "      ";
   		Case Chr(255)+"K"'left
   			wx-=30
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print wx;",";wy; "      ";
   		Case Chr(255)+"M"'right
   			wx+=30
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print wx;",";wy; "      ";
   		Case "-"
   			rref=rref*1.1111111
   			'wx=wx\2
   			'wy=wy\2
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print rref; "      ";
         Case "+"
   			rref=rref/1.1111111
   			'wx=wx*2
   			'wy=wy*2
   			
   			ttp1=cosunghi*cosunghi2
         	ttp2=sinunghi*cosunghi2
         	maxxshr=maxx Shr 1
         	maxyshr=maxy Shr 1
         	maxxshr+=wx
         	maxyshr+=wy
         	Locate 1,1
         	Print rref; "      ";
         Case "f"
            ' show genius fields
            ' for screen x,y = (0,0) to screen x,y = (max_x, max_y) -> tranform positions
            ' use position history to draw genius fields lorentz transformed
            ' we can draw in a buffer the instantaneous field and compare results 
               
               ' logit(fz) = constant for this
               
            fields_graphics(rref, Sinunghi, cosunghi, ttp1, ttp2, _
               maxxshr, maxyshr, e, m_prot, quants, num_history, _
               sinunghi2, num_particles, tempdelta, epsi_, tempi, tempx, tempy, tempz, _
               np, epsi_c2, tempt, miudiv4pi, c, maxTime, _
               ma1, ma2, fma, mres)
         Case Chr(27)
   			Exit do
   	End Select
   	
   	
   	
   	' graphic output
   
   	
      currentTime = 0
      Do
         k = currentTime Shl np
         j = k + num_particles
         i = k
         i = k
         While i < j
            ' in 2D plane
            'PSet(logit(quants[i].x) + x_opt,y_opt - logit(quants[i].z)), color1
            ' rotated by 1 angle
            'logit_y = logit(quants[i].y) 
            'xx = -logit_y * sin_opt + logit(quants[i].x) + x_opt
         	'yy = -logit_y * cos_opt - logit(quants[i].z) + y_opt
         	'PSet(xx,yy), color1
         	'rotated bt 2 angles
         	logit_y = logit(quants[i].y)
         	logit_x = logit(quants[i].x) 
   	      xx = -logit_x * Sinunghi + logit_y * cosunghi + maxxshr
   			yy = -logit_x*ttp1 - logit_y*ttp2 + logit(quants[i].z) * sinunghi2 + maxyshr
   			PSet(xx,yy), color1
         	i+=1
         Wend 
      	currentTime += 50000
      Loop Until currentTime >= maxtime
      /'
      currentTime = maxtempi Shr np
      Do
         k = currentTime Shl np
         j = k + num_particles
         i = k
         i = k
         While i < j
            ' in 2D plane
            'PSet(logit(quants[i].x) + x_opt,y_opt - logit(quants[i].z)), color1
            ' rotated by 1 angle
            'logit_y = logit(quants[i].y) 
            'xx = -logit_y * sin_opt + logit(quants[i].x) + x_opt
         	'yy = -logit_y * cos_opt - logit(quants[i].z) + y_opt
         	'PSet(xx,yy), color1
         	'rotated bt 2 angles
         	logit_y = logit(quants[i].y)
         	logit_x = logit(quants[i].x) 
   	      xx = -logit_x * Sinunghi + logit_y * cosunghi + maxxshr
   			yy = -logit_x*ttp1 - logit_y*ttp2 + logit(quants[i].z) * sinunghi2 + maxyshr
   			PSet(xx,yy), color3
         	i+=1
         Wend 
      	currentTime += 50000
      Loop Until currentTime >= maxtime
      '/    
      ' prepare for next cycle
      'currentTime = maxTime - 1
      'i = ((currentTime) Shl np ) 
      i = maxtempi 
      For k = 0 To (num_history-1-maxtempi)
         'coordTrans(@quants[i+k], c, 0)
   	   CopyQuanta(@quants[k], @quants[i+k])
      Next k
   
      ' clear for next cycle'(num_history-maxtempi-num_particles) ' (num_history-maxtempi-num_particles-2) 
      For i = 0 To num_history-1 '0 to num_history-1
         quants[i].Ex = 0
         quants[i].Ey = 0
         quants[i].Ez = 0
         quants[i].Bx = 0
         quants[i].By = 0
         quants[i].Bz = 0
      Next i
      

   Loop While 1 = 1 
   DeAllocate quants
 
End Sub

[dodicat]

Thanks Mihail_B
I had to reduce line 985
num_history = 3100000
to get it running
win 10, but only 2Gb memory.
Looks good.
It's been a long time since I studied Physics.
B fields/E fields ~~ tensors, I remember, but I would need some serious revision to appreciate your work.