# gnuplot commands
# suffisso r => riferimento in cui q e' a riposo
  c = 1.

# v tale che: gamma = 2
 v = sqrt(3.)/2. 

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

# v/c = .99 (gamma = 10)
# v = .999

  eps0 = 1.
  q = 1.

  beta = v/c
  gamma = 1./sqrt(1-beta**2)
  tr(x,y,z,t) = gamma*(t-beta*x)
  xr(x,y,z,t) = gamma*(x-v*t)
  yr(x,y,z,t) = y
  zr(x,y,z,t) = z

  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)

  min(x,y) = (x < y) ? x : y
# trucchetto: i campi vanno all'infinito, quindi li limito a un valore plottabile

  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) = gamma*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) = gamma*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) = min (100., 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) = gamma*(beta*Ez(x,y,z,t) )    
  Bz(x,y,z,t) = gamma*(-beta*Ey(x,y,z,t) )  
  scala = beta/.01
# riporto tutti i valori di B a quelli della "velocita' non relativistica"
  B(x,y,z,t) = min (1.,sqrt (Bx(x,y,z,t)**2 + By(x,y,z,t)**2 + Bz(x,y,z,t)**2)/scala )

 set view map
 set aspect 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 E(x,y,0,0), B(x,y,0,0)