# gnuplot commands

# natural units system:
  c = 1.
  eps0 = 1.

# unitary charge
  q = 1.

# uncomment desired velocity
# or write it in gnuplot command line 
#   before pasting the whole script

# v << c (gamma = 1.00005..):
# v = .01

# v = c/2 (gamma = 1.1547..):
# v = .5

# v such that gamma = 2
# v = sqrt(3.)/2. 

# v/c = .99 (gamma = 7.09..)
# v = .99

# v/c = .999 (gamma = 22.4)
# v = .999

# v/c = .9999 (gamma = 70)
# v = .9999

  beta(v) = v/c
  Lorentz(v) = 1./sqrt(1-beta(v)**2)
# note: gamma() cannot be used, is a gnuplot intrinsic function

# Lorentz transformation of coordinates 
# from laboratory's reference system to 
# reference system where q is at rest
# (suffix r indicates the reference system where q is at Rest)

  tr(x,y,z,t) = Lorentz(v)*(t+beta(v)*x)
  xr(x,y,z,t) = Lorentz(v)*(x+v*t)
  yr(x,y,z,t) = y
  zr(x,y,z,t) = z

# E computed in reference where q is at rest (and B is zero):

  rsquare(xr,yr,zr) = xr**2+yr**2+zr**2
  r(xr,yr,zr) = sqrt(rsquare(xr,yr,zr))
  Er(xr,yr,zr,tr) = q/(4*pi*eps0)/rsquare(xr,yr,zr)
  Erx(xr,yr,zr,tr) = Er(xr,yr,zr,tr)*xr/r(xr,yr,zr)
  Ery(xr,yr,zr,tr) = Er(xr,yr,zr,tr)*yr/r(xr,yr,zr)
  Erz(xr,yr,zr,tr) = Er(xr,yr,zr,tr)*zr/r(xr,yr,zr)

# Lorentz transformation of EM tensor to Laboratory' reference (SI units):
#
#  E'x = Ex      
#  E'y = gamma(Ey - v Bz )  (but Bz=0)
#  E'z = gamma(Ez + v By )  (but By=0)
#
#  B'x = Bx
#  B'y = gamma(By + beta/c Ez )  (but By=0)
#  B'z = gamma(Bz - beta/c Ey )  (but Bz=0)


  Ex(x,y,z,t) =  Erx(xr(x,y,z,t),yr(x,y,z,t),zr(x,y,z,t),tr(x,y,z,t))
  Ey(x,y,z,t) = Lorentz(v)*Ery(xr(x,y,z,t),yr(x,y,z,t),zr(x,y,z,t),tr(x,y,z,t))
  Ez(x,y,z,t) = Lorentz(v)*Erz(xr(x,y,z,t),yr(x,y,z,t),zr(x,y,z,t),tr(x,y,z,t))

  E(x,y,z,t) = sqrt (Ex(x,y,z,t)**2 + Ey(x,y,z,t)**2 + Ez(x,y,z,t)**2)


  Bx(x,y,z,t) = 0.
  By(x,y,z,t) = Lorentz(v)*(-beta(v)/c*Ez(x,y,z,t) )    
  Bz(x,y,z,t) = Lorentz(v)*(+beta(v)/c*Ey(x,y,z,t) )  

  B(x,y,z,t) = sqrt (Bx(x,y,z,t)**2 + By(x,y,z,t)**2 + Bz(x,y,z,t)**2)


# graphic trick: fields goes to infinity approaching origin, 
# let's limit them to a max plottable value

  min(x,y) = (x < y) ? x : y
  Eplot(x,y,z,t) = min (1., E(x,y,z,t) )

# graphic trick: B increases with v even in non-relativistic cases,
# let's divide it by v so as to compare plots at different v.

  Bplot(x,y,z,t) = min (1., B(x,y,z,t)/v )

 set view map
 set size ratio -1
 unset surface
 set xrange [-1:1]
 set yrange [-1:1]
 set samples 500
 set isosamples 500
 set contour base
 set cntrparam levels auto 10 

 splot Eplot(x,y,0,0), Bplot(x,y,0,0)