# Animated 3D plot in spherical coordinates with gnuplot

Published:

I am working on a physics experiment that involves muons that are stored in a magnetic field and are traveling along circular orbits. Their motion inside the magnetic field is quite complex and involves oscillations with different periods – essentially they are like a spinning top, but instead of precessing around only one axis, they can precess around three.

What I am interested in is to track the spin of the muons through time and in order to visualize things better I wanted to have an animated 3D plot of the vector of the spin as the muons experience the changing magnetic field. My favorite plotting tool for 3D graphs is `gnuplot`. The code to generate this plot is a bit large and messy, but it may be useful to someone so I’m posting it here:

``````# Initial spin vector
Sx = -0.0
Sy = -1.0
Sz = -0.0

# Transform to spherical coordinates
phi_0 = atan2(Sx, Sy) # rad
theta_0 = atan2(sqrt(Sx**2+Sy**2), Sz)-pi/2 # rad

# Spin motion equations
Eex = 0.05

f(x) = (191.49*(Eex) + 178.05*0.06/27.94/omega)*(cos(omega*x))
l(x) = -(((sin(omega*x+phi_0)))*81.75e-3/omega)
p(x) =  phi_0 -om_z*x
if (om_z != 0) {
e(x) = -e_z/om_z*sin(p(x))
} else {
e(x) = e_z*x*cos(phi_0)
}
c =    (f(0))*sin(p(0)) + l(0)*cos(p(0)) + e(0)
g(x) = (f(x))*sin(p(x)) + l(x)*cos(p(x)) + e(x) + theta_0 * 1e6 - c
# End of spin motion equations

# Plotting
set key samplen 1 height 5 bottom
set view 70, 340, 1.6, 1.6
set view equal xyz
set parametric
unset xtics
unset ytics
unset ztics
set border 0
set xlabel 'p'
set ylabel 'z'
set isosamples 12
set xyplane -0.5
set grid x y
set xrange [-1:1]
set yrange [-1:1]
set zrange [-1:1]
set urange [0:2*pi]
set vrange [-pi:pi]

phi(x) = (g(x) - theta_0*1e6 + c) / 20 + theta_0 - c / 1e6
the(x) = p(x) + f(x) * sin(phi(x)) / 20
av_the(x) = p(x)
av_phi(x) = e(x) / 20

av_x(u, v) = v*cos(av_the(u))*sin(av_phi(u))
av_y(u, v) = v*cos(av_the(u))*cos(av_phi(u))
av_z(u, v) = v*sin(av_the(u))

fx(u,v) = v*cos(the(u))*sin(phi(u))
fy(u,v) = v*cos(the(u))*cos(phi(u))
fz(u,v) = v*sin(the(u))

# Parametric functions for the sphere
r = 1.0
sx(u,v) = r*cos(u)*sin(v)
sy(u,v) = r*cos(u)*cos(v)
sz(u,v) = r*sin(u)

# Plot axis vectors
set arrow 2 from 0,0,0 to 0,-1.2,0 lc rgb 'gray40' lw 2 back
set arrow 5 from 0,0,0 to -1.2,0,0 lc rgb 'gray40' lw 2 back
set arrow 6 from 0,0,0 to 0,0,-1.2 lc rgb 'gray40' lw 2 back

# Plot magnetic fields in the top right corner
o = 0.3
o0 = 1.5
set arrow 7 from o0,o0,o to o0,o0,(o+sin(omega*t)/10)
set arrow 8 from o0,o0,o to 1.2,o0,o
set arrow 9 from o0,o0,o to o0,(o0+sin(omega*t)/8),o
set label 1 'B_r' at o0+0.05, o0+0.08, o+0.05
set label 2 'B_z' at 1.2,o0,o-0.05
set label 3 'B_p' at 1.6, 1.6, o-0.15

set term pngcairo size 800,800 enhanced
outfile = sprintf('animation/3d_vector_plot%05.0f.png', t)
newtitle = sprintf('t = %i ns', t)
set title newtitle
set output outfile

set object 99 circle at 0,0,0 radius 0.025 front fc rgb "red" fs solid
set arrow 1 from fx(t,-1),fy(t,-1),fz(t,-1) to fx(t, 1), fy(t, 1), fz(t, 1) back lc rgb 'black' lw 2

set multiplot
# Plot the sphere, the equator and zero meridian
splot sx(u,v),sy(u,v),sz(u,0) lc 'black' dt 3 notitle,\
sx(0, v), sy(0,v), sz(0, 0) lc rgb '#8888dd' notitle,\
0, sy(pi,u), sz(u,0) lc rgb '#dd8888' notitle

# Plot an arrow pointing to the average movement
set arrow 3 from 0,0,0 to av_x(t, 0.5), av_y(t, 0.5), 0 lc rgb 'red' back

# Change the range
set urange [0:t]
splot av_x(u, 0.46), av_y(u, 0.46),av_z(0, 0.46)  lw 2 lc rgb 'red' notitle,\
av_x(u, 1), av_y(u,1), av_z(u,1) lw 2 lc rgb 'red' t 'Average value'

# Plot the spin trajectory, but only a part of it
set key height 10 bottom
start = t - 400 > 0 ? t - 400 : 0
set urange [start:t]
splot fx(u, 1), fy(u, 1), fz(u, 1) t 'Spin position' lc rgb '#57a757'

unset multiplot

print 'Time step: ', t, ' ns'

t = t + step
unset arrow 1
``````\$ gnuplot
Note: be sure that you have a folder like `animation/` to store all the files!