# -*- coding: utf-8 -*-
from load import *
from fft import *
from plots import *

print('\nplotting fields\n')
outdir = './fig_fields/'

# Load 2D cut
ncfile = netcdf.netcdf_file(input_dir+runname+'.out.2D.nc'+restart_num, 'r')   

tt_fld  = np.copy(ncfile.variables['tt' ][:]); tt_fld  = np.delete(tt_fld , ignored_points_fld, axis = 0)
xx_fld  = np.copy(ncfile.variables['xx' ][:])
yy_fld  = np.copy(ncfile.variables['yy' ][:])
zz_fld  = np.copy(ncfile.variables['zz' ][:])

rho_r_z0 = np.copy(ncfile.variables['rho_r_z0'][:]); rho_r_z0 = np.delete(rho_r_z0, ignored_points_fld, axis = 0)
rho_r_x0 = np.copy(ncfile.variables['rho_r_x0'][:]); rho_r_x0 = np.delete(rho_r_x0, ignored_points_fld, axis = 0)
rho_r_y0 = np.copy(ncfile.variables['rho_r_y0'][:]); rho_r_y0 = np.delete(rho_r_y0, ignored_points_fld, axis = 0)

mx_r_z0  = np.copy(ncfile.variables['mx_r_z0' ][:]); mx_r_z0  = np.delete(mx_r_z0 , ignored_points_fld, axis = 0)
mx_r_x0  = np.copy(ncfile.variables['mx_r_x0' ][:]); mx_r_x0  = np.delete(mx_r_x0 , ignored_points_fld, axis = 0)
mx_r_y0  = np.copy(ncfile.variables['mx_r_y0' ][:]); mx_r_y0  = np.delete(mx_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
my_r_z0  = np.copy(ncfile.variables['my_r_z0' ][:]); my_r_z0  = np.delete(my_r_z0 , ignored_points_fld, axis = 0)
my_r_x0  = np.copy(ncfile.variables['my_r_x0' ][:]); my_r_x0  = np.delete(my_r_x0 , ignored_points_fld, axis = 0)
my_r_y0  = np.copy(ncfile.variables['my_r_y0' ][:]); my_r_y0  = np.delete(my_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
mz_r_z0  = np.copy(ncfile.variables['mz_r_z0' ][:]); mz_r_z0  = np.delete(mz_r_z0 , ignored_points_fld, axis = 0)
mz_r_x0  = np.copy(ncfile.variables['mz_r_x0' ][:]); mz_r_x0  = np.delete(mz_r_x0 , ignored_points_fld, axis = 0)
mz_r_y0  = np.copy(ncfile.variables['mz_r_y0' ][:]); mz_r_y0  = np.delete(mz_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
wx_r_z0  = np.copy(ncfile.variables['wx_r_z0' ][:]); wx_r_z0  = np.delete(wx_r_z0 , ignored_points_fld, axis = 0)
wx_r_x0  = np.copy(ncfile.variables['wx_r_x0' ][:]); wx_r_x0  = np.delete(wx_r_x0 , ignored_points_fld, axis = 0)
wx_r_y0  = np.copy(ncfile.variables['wx_r_y0' ][:]); wx_r_y0  = np.delete(wx_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
wy_r_z0  = np.copy(ncfile.variables['wy_r_z0' ][:]); wy_r_z0  = np.delete(wy_r_z0 , ignored_points_fld, axis = 0)
wy_r_x0  = np.copy(ncfile.variables['wy_r_x0' ][:]); wy_r_x0  = np.delete(wy_r_x0 , ignored_points_fld, axis = 0)
wy_r_y0  = np.copy(ncfile.variables['wy_r_y0' ][:]); wy_r_y0  = np.delete(wy_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
wz_r_z0  = np.copy(ncfile.variables['wz_r_z0' ][:]); wz_r_z0  = np.delete(wz_r_z0 , ignored_points_fld, axis = 0)
wz_r_x0  = np.copy(ncfile.variables['wz_r_x0' ][:]); wz_r_x0  = np.delete(wz_r_x0 , ignored_points_fld, axis = 0)
wz_r_y0  = np.copy(ncfile.variables['wz_r_y0' ][:]); wz_r_y0  = np.delete(wz_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
bx_r_z0  = np.copy(ncfile.variables['bx_r_z0' ][:]); bx_r_z0  = np.delete(bx_r_z0 , ignored_points_fld, axis = 0)
bx_r_x0  = np.copy(ncfile.variables['bx_r_x0' ][:]); bx_r_x0  = np.delete(bx_r_x0 , ignored_points_fld, axis = 0)
bx_r_y0  = np.copy(ncfile.variables['bx_r_y0' ][:]); bx_r_y0  = np.delete(bx_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
by_r_z0  = np.copy(ncfile.variables['by_r_z0' ][:]); by_r_z0  = np.delete(by_r_z0 , ignored_points_fld, axis = 0)
by_r_x0  = np.copy(ncfile.variables['by_r_x0' ][:]); by_r_x0  = np.delete(by_r_x0 , ignored_points_fld, axis = 0)
by_r_y0  = np.copy(ncfile.variables['by_r_y0' ][:]); by_r_y0  = np.delete(by_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
bz_r_z0  = np.copy(ncfile.variables['bz_r_z0' ][:]); bz_r_z0  = np.delete(bz_r_z0 , ignored_points_fld, axis = 0)
bz_r_x0  = np.copy(ncfile.variables['bz_r_x0' ][:]); bz_r_x0  = np.delete(bz_r_x0 , ignored_points_fld, axis = 0)
bz_r_y0  = np.copy(ncfile.variables['bz_r_y0' ][:]); bz_r_y0  = np.delete(bz_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
jx_r_z0  = np.copy(ncfile.variables['jx_r_z0' ][:]); jx_r_z0  = np.delete(jx_r_z0 , ignored_points_fld, axis = 0)
jx_r_x0  = np.copy(ncfile.variables['jx_r_x0' ][:]); jx_r_x0  = np.delete(jx_r_x0 , ignored_points_fld, axis = 0)
jx_r_y0  = np.copy(ncfile.variables['jx_r_y0' ][:]); jx_r_y0  = np.delete(jx_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
jy_r_z0  = np.copy(ncfile.variables['jy_r_z0' ][:]); jy_r_z0  = np.delete(jy_r_z0 , ignored_points_fld, axis = 0)
jy_r_x0  = np.copy(ncfile.variables['jy_r_x0' ][:]); jy_r_x0  = np.delete(jy_r_x0 , ignored_points_fld, axis = 0)
jy_r_y0  = np.copy(ncfile.variables['jy_r_y0' ][:]); jy_r_y0  = np.delete(jy_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
jz_r_z0  = np.copy(ncfile.variables['jz_r_z0' ][:]); jz_r_z0  = np.delete(jz_r_z0 , ignored_points_fld, axis = 0)
jz_r_x0  = np.copy(ncfile.variables['jz_r_x0' ][:]); jz_r_x0  = np.delete(jz_r_x0 , ignored_points_fld, axis = 0)
jz_r_y0  = np.copy(ncfile.variables['jz_r_y0' ][:]); jz_r_y0  = np.delete(jz_r_y0 , ignored_points_fld, axis = 0)

ux_r_z0  = mx_r_z0/rho_r_z0 
ux_r_x0  = mx_r_x0/rho_r_x0 
ux_r_y0  = mx_r_y0/rho_r_y0 
          
uy_r_z0  = my_r_z0/rho_r_z0
uy_r_x0  = my_r_x0/rho_r_x0 
uy_r_y0  = my_r_y0/rho_r_y0 
           
uz_r_z0  = mz_r_z0/rho_r_z0
uz_r_x0  = mz_r_x0/rho_r_x0 
uz_r_y0  = mz_r_y0/rho_r_y0 

ncfile.close()

#--------------------------------------------------------#
#                   plot final snapshot                  #
#--------------------------------------------------------#
u_r_z0 = np.sqrt(ux_r_z0**2 + uy_r_z0**2 + uz_r_z0**2)
u_r_y0 = np.sqrt(ux_r_y0**2 + uy_r_y0**2 + uz_r_y0**2)
u_r_x0 = np.sqrt(ux_r_x0**2 + uy_r_x0**2 + uz_r_x0**2)
w_r_z0 = np.sqrt(wx_r_z0**2 + wy_r_z0**2 + wz_r_z0**2)
w_r_y0 = np.sqrt(wx_r_y0**2 + wy_r_y0**2 + wz_r_y0**2)
w_r_x0 = np.sqrt(wx_r_x0**2 + wy_r_x0**2 + wz_r_x0**2)
b_r_z0 = np.sqrt(bx_r_z0**2 + by_r_z0**2 + bz_r_z0**2)
b_r_y0 = np.sqrt(bx_r_y0**2 + by_r_y0**2 + bz_r_y0**2)
b_r_x0 = np.sqrt(bx_r_x0**2 + by_r_x0**2 + bz_r_x0**2)
j_r_z0 = np.sqrt(jx_r_z0**2 + jy_r_z0**2 + jz_r_z0**2)
j_r_y0 = np.sqrt(jx_r_y0**2 + jy_r_y0**2 + jz_r_y0**2)
j_r_x0 = np.sqrt(jx_r_x0**2 + jy_r_x0**2 + jz_r_x0**2)
plot_3d(rho_r_z0[final_fld_idx,:,:], rho_r_y0[final_fld_idx,:,:], rho_r_x0[final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$\rho (t = %.2E)$'     % tt_fld[final_fld_idx], cmp=parula_map, save=outdir+'rho.pdf')
plot_3d(u_r_z0  [final_fld_idx,:,:], u_r_y0  [final_fld_idx,:,:], u_r_x0  [final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, (ux_r_z0[final_fld_idx,:,:], uy_r_z0[final_fld_idx,:,:]), (ux_r_y0[final_fld_idx,:,:], uz_r_y0[final_fld_idx,:,:]), (uy_r_x0[final_fld_idx,:,:], uz_r_x0[final_fld_idx,:,:]), xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\bm{u}| (t = %.2E)$' % tt_fld[final_fld_idx], cmp=parula_map, streamline_density=1, streamline_width=0.5, streamline_color='w', save=outdir+'u.pdf')
plot_3d(b_r_z0  [final_fld_idx,:,:], b_r_y0  [final_fld_idx,:,:], b_r_x0  [final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, (bx_r_z0[final_fld_idx,:,:], by_r_z0[final_fld_idx,:,:]), (bx_r_y0[final_fld_idx,:,:], bz_r_y0[final_fld_idx,:,:]), (by_r_x0[final_fld_idx,:,:], bz_r_x0[final_fld_idx,:,:]), xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\bm{B}| (t = %.2E)$' % tt_fld[final_fld_idx], cmp=parula_map, streamline_density=1, streamline_width=0.5, streamline_color='w', save=outdir+'b.pdf')
if is2D:
  plot_3d(wz_r_z0[final_fld_idx,:,:], wz_r_y0[final_fld_idx,:,:], wz_r_x0[final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$(\Curl\bm{u})_z (t = %.2E)$' % tt_fld[final_fld_idx], cmp='RdBu_r', save=outdir+'w.pdf')
  plot_3d(jz_r_z0[final_fld_idx,:,:], jz_r_y0[final_fld_idx,:,:], jz_r_x0[final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$(\Curl\bm{B})_z (t = %.2E)$' % tt_fld[final_fld_idx], cmp='RdBu_r', save=outdir+'j.pdf')
else:
  plot_3d( w_r_z0[final_fld_idx,:,:],  w_r_y0[final_fld_idx,:,:],  w_r_x0[final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\Curl\bm{u}| (t = %.2E)$' % tt_fld[final_fld_idx], cmp=parula_map, save=outdir+'w.pdf')
  plot_3d( j_r_z0[final_fld_idx,:,:],  j_r_y0[final_fld_idx,:,:],  j_r_x0[final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\Curl\bm{B}| (t = %.2E)$' % tt_fld[final_fld_idx], cmp=parula_map, save=outdir+'j.pdf')


#--------------------------------------------------------#
#                       plot movie                       #
#--------------------------------------------------------#
if ismovie:
  # evenly space time
  # nframe = min(200, int(tt[:final_fld_idx].size))
  nframe = tt_fld[:final_fld_idx].size
  idx = np.unique([np.argmin(abs(tt_fld - np.linspace(tt_fld[0], tt_fld[-1], nframe)[i])) for i in range(0, nframe)])
  tt_fld  = tt_fld [idx]

  rho_r_z0 = np.take(rho_r_z0, idx, axis=0)
  rho_r_y0 = np.take(rho_r_y0, idx, axis=0)
  rho_r_x0 = np.take(rho_r_x0, idx, axis=0)

  ux_r_z0  = np.take(ux_r_z0 , idx, axis=0)
  ux_r_y0  = np.take(ux_r_y0 , idx, axis=0)
  ux_r_x0  = np.take(ux_r_x0 , idx, axis=0)
  uy_r_z0  = np.take(uy_r_z0 , idx, axis=0)
  uy_r_y0  = np.take(uy_r_y0 , idx, axis=0)
  uy_r_x0  = np.take(uy_r_x0 , idx, axis=0)
  uz_r_z0  = np.take(uz_r_z0 , idx, axis=0)
  uz_r_y0  = np.take(uz_r_y0 , idx, axis=0)
  uz_r_x0  = np.take(uz_r_x0 , idx, axis=0)
                             
  wx_r_z0  = np.take(wx_r_z0 , idx, axis=0)
  wx_r_y0  = np.take(wx_r_y0 , idx, axis=0)
  wx_r_x0  = np.take(wx_r_x0 , idx, axis=0)
  wy_r_z0  = np.take(wy_r_z0 , idx, axis=0)
  wy_r_y0  = np.take(wy_r_y0 , idx, axis=0)
  wy_r_x0  = np.take(wy_r_x0 , idx, axis=0)
  wz_r_z0  = np.take(wz_r_z0 , idx, axis=0)
  wz_r_y0  = np.take(wz_r_y0 , idx, axis=0)
  wz_r_x0  = np.take(wz_r_x0 , idx, axis=0)
                             
  bx_r_z0  = np.take(bx_r_z0 , idx, axis=0)
  bx_r_y0  = np.take(bx_r_y0 , idx, axis=0)
  bx_r_x0  = np.take(bx_r_x0 , idx, axis=0)
  by_r_z0  = np.take(by_r_z0 , idx, axis=0)
  by_r_y0  = np.take(by_r_y0 , idx, axis=0)
  by_r_x0  = np.take(by_r_x0 , idx, axis=0)
  bz_r_z0  = np.take(bz_r_z0 , idx, axis=0)
  bz_r_y0  = np.take(bz_r_y0 , idx, axis=0)
  bz_r_x0  = np.take(bz_r_x0 , idx, axis=0)
                             
  jx_r_z0  = np.take(jx_r_z0 , idx, axis=0)
  jx_r_y0  = np.take(jx_r_y0 , idx, axis=0)
  jx_r_x0  = np.take(jx_r_x0 , idx, axis=0)
  jy_r_z0  = np.take(jy_r_z0 , idx, axis=0)
  jy_r_y0  = np.take(jy_r_y0 , idx, axis=0)
  jy_r_x0  = np.take(jy_r_x0 , idx, axis=0)
  jz_r_z0  = np.take(jz_r_z0 , idx, axis=0)
  jz_r_y0  = np.take(jz_r_y0 , idx, axis=0)
  jz_r_x0  = np.take(jz_r_x0 , idx, axis=0)

  u_r_z0 = np.sqrt(ux_r_z0**2 + uy_r_z0**2 + uz_r_z0**2)
  u_r_y0 = np.sqrt(ux_r_y0**2 + uy_r_y0**2 + uz_r_y0**2)
  u_r_x0 = np.sqrt(ux_r_x0**2 + uy_r_x0**2 + uz_r_x0**2)
  w_r_z0 = np.sqrt(wx_r_z0**2 + wy_r_z0**2 + wz_r_z0**2)
  w_r_y0 = np.sqrt(wx_r_y0**2 + wy_r_y0**2 + wz_r_y0**2)
  w_r_x0 = np.sqrt(wx_r_x0**2 + wy_r_x0**2 + wz_r_x0**2)
  b_r_z0 = np.sqrt(bx_r_z0**2 + by_r_z0**2 + bz_r_z0**2)
  b_r_y0 = np.sqrt(bx_r_y0**2 + by_r_y0**2 + bz_r_y0**2)
  b_r_x0 = np.sqrt(bx_r_x0**2 + by_r_x0**2 + bz_r_x0**2)
  j_r_z0 = np.sqrt(jx_r_z0**2 + jy_r_z0**2 + jz_r_z0**2)
  j_r_y0 = np.sqrt(jx_r_y0**2 + jy_r_y0**2 + jz_r_y0**2)
  j_r_x0 = np.sqrt(jx_r_x0**2 + jy_r_x0**2 + jz_r_x0**2)
  movie_3d(tt_fld, rho_r_z0, rho_r_y0, rho_r_x0, xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$\rho$'    , cmp=parula_map, save=outdir+'rho_anim.gif')
  movie_3d(tt_fld, u_r_z0  , u_r_y0, u_r_x0    , xx_fld, yy_fld, zz_fld, (ux_r_z0, uy_r_z0), (ux_r_y0, uz_r_y0), (uy_r_x0, uz_r_x0), xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\bm{u}|$', cmp=parula_map, save=outdir+'u_anim.gif')
  movie_3d(tt_fld, b_r_z0  , b_r_y0, b_r_x0    , xx_fld, yy_fld, zz_fld, (bx_r_z0, by_r_z0), (bx_r_y0, bz_r_y0), (by_r_x0, bz_r_x0), xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\bm{B}|$', cmp=parula_map, save=outdir+'b_anim.gif')
  if is2D:
    movie_3d(tt_fld, wz_r_z0, wz_r_y0, wz_r_x0, xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$(\Curl\bm{u})_z$', cmp='RdBu_r', save=outdir+'w_anim.gif')
    movie_3d(tt_fld, jz_r_z0, jz_r_y0, jz_r_x0, xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$(\Curl\bm{B})_z$', cmp='RdBu_r', save=outdir+'j_anim.gif')
  else:
    movie_3d(tt_fld,  w_r_z0,  w_r_y0,  w_r_x0, xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\Curl\bm{u}|$', cmp=parula_map, save=outdir+'w_anim.gif')
    movie_3d(tt_fld,  j_r_z0,  j_r_y0,  j_r_x0, xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\Curl\bm{B}|$', cmp=parula_map, save=outdir+'j_anim.gif')


#------------------#
#   output data    #
#------------------#
from scipy.io import savemat

savemat(outdir + 'grid' , {'xx':xx, 'yy':yy, 'zz':zz})
savemat(outdir + 'rho_r' , {
                            'tt'       :tt[final_fld_idx],
                            'rho_r_z0' :rho_r_z0[final_fld_idx,:,:],
                            'rho_r_x0' :rho_r_x0[final_fld_idx,:,:],
                            'rho_r_y0' :rho_r_y0[final_fld_idx,:,:],
                          })
savemat(outdir + 'u_r' , {
                            'tt'       :tt[final_fld_idx],
                            'ux_r_z0'  :ux_r_z0 [final_fld_idx,:,:],
                            'ux_r_x0'  :ux_r_x0 [final_fld_idx,:,:],
                            'ux_r_y0'  :ux_r_y0 [final_fld_idx,:,:],
                            #
                            'uy_r_z0'  :uy_r_z0 [final_fld_idx,:,:],
                            'uy_r_x0'  :uy_r_x0 [final_fld_idx,:,:],
                            'uy_r_y0'  :uy_r_y0 [final_fld_idx,:,:],
                            #
                            'uz_r_z0'  :uz_r_z0 [final_fld_idx,:,:],
                            'uz_r_x0'  :uz_r_x0 [final_fld_idx,:,:],
                            'uz_r_y0'  :uz_r_y0 [final_fld_idx,:,:],
                          })
savemat(outdir + 'b_r' , {
                            'tt'       :tt[final_fld_idx],
                            'bx_r_z0'  :bx_r_z0 [final_fld_idx,:,:],
                            'bx_r_x0'  :bx_r_x0 [final_fld_idx,:,:],
                            'bx_r_y0'  :bx_r_y0 [final_fld_idx,:,:],
                            #
                            'by_r_z0'  :by_r_z0 [final_fld_idx,:,:],
                            'by_r_x0'  :by_r_x0 [final_fld_idx,:,:],
                            'by_r_y0'  :by_r_y0 [final_fld_idx,:,:],
                            #
                            'bz_r_z0'  :bz_r_z0 [final_fld_idx,:,:],
                            'bz_r_x0'  :bz_r_x0 [final_fld_idx,:,:],
                            'bz_r_y0'  :bz_r_y0 [final_fld_idx,:,:],
                          })
savemat(outdir + 'w_r' , {
                            'tt'       :tt[final_fld_idx],
                            'wx_r_z0'  :wx_r_z0 [final_fld_idx,:,:],
                            'wx_r_x0'  :wx_r_x0 [final_fld_idx,:,:],
                            'wx_r_y0'  :wx_r_y0 [final_fld_idx,:,:],
                            #
                            'wy_r_z0'  :wy_r_z0 [final_fld_idx,:,:],
                            'wy_r_x0'  :wy_r_x0 [final_fld_idx,:,:],
                            'wy_r_y0'  :wy_r_y0 [final_fld_idx,:,:],
                            #
                            'wz_r_z0'  :wz_r_z0 [final_fld_idx,:,:],
                            'wz_r_x0'  :wz_r_x0 [final_fld_idx,:,:],
                            'wz_r_y0'  :wz_r_y0 [final_fld_idx,:,:],
                          })
savemat(outdir + 'j_r' , {
                            'tt'       :tt[final_fld_idx],
                            'jx_r_z0'  :jx_r_z0 [final_fld_idx,:,:],
                            'jx_r_x0'  :jx_r_x0 [final_fld_idx,:,:],
                            'jx_r_y0'  :jx_r_y0 [final_fld_idx,:,:],
                            #
                            'jy_r_z0'  :jy_r_z0 [final_fld_idx,:,:],
                            'jy_r_x0'  :jy_r_x0 [final_fld_idx,:,:],
                            'jy_r_y0'  :jy_r_y0 [final_fld_idx,:,:],
                            #
                            'jz_r_z0'  :jz_r_z0 [final_fld_idx,:,:],
                            'jz_r_x0'  :jz_r_x0 [final_fld_idx,:,:],
                            'jz_r_y0'  :jz_r_y0 [final_fld_idx,:,:],
                          })

# -*- coding: utf-8 -*-
from load import *
from fft import *
from plots import *

print('\nplotting fields\n')
outdir = './fig_fields/'

# Load 2D cut
ncfile = netcdf.netcdf_file(input_dir+runname+'.out.2D.nc'+restart_num, 'r')   

tt_fld  = np.copy(ncfile.variables['tt' ][:]); tt_fld  = np.delete(tt_fld , ignored_points_fld, axis = 0)
xx_fld  = np.copy(ncfile.variables['xx' ][:])
yy_fld  = np.copy(ncfile.variables['yy' ][:])
zz_fld  = np.copy(ncfile.variables['zz' ][:])

rho_r_z0 = np.copy(ncfile.variables['rho_r_z0'][:]); rho_r_z0 = np.delete(rho_r_z0, ignored_points_fld, axis = 0)
rho_r_x0 = np.copy(ncfile.variables['rho_r_x0'][:]); rho_r_x0 = np.delete(rho_r_x0, ignored_points_fld, axis = 0)
rho_r_y0 = np.copy(ncfile.variables['rho_r_y0'][:]); rho_r_y0 = np.delete(rho_r_y0, ignored_points_fld, axis = 0)

mx_r_z0  = np.copy(ncfile.variables['mx_r_z0' ][:]); mx_r_z0  = np.delete(mx_r_z0 , ignored_points_fld, axis = 0)
mx_r_x0  = np.copy(ncfile.variables['mx_r_x0' ][:]); mx_r_x0  = np.delete(mx_r_x0 , ignored_points_fld, axis = 0)
mx_r_y0  = np.copy(ncfile.variables['mx_r_y0' ][:]); mx_r_y0  = np.delete(mx_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
my_r_z0  = np.copy(ncfile.variables['my_r_z0' ][:]); my_r_z0  = np.delete(my_r_z0 , ignored_points_fld, axis = 0)
my_r_x0  = np.copy(ncfile.variables['my_r_x0' ][:]); my_r_x0  = np.delete(my_r_x0 , ignored_points_fld, axis = 0)
my_r_y0  = np.copy(ncfile.variables['my_r_y0' ][:]); my_r_y0  = np.delete(my_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
mz_r_z0  = np.copy(ncfile.variables['mz_r_z0' ][:]); mz_r_z0  = np.delete(mz_r_z0 , ignored_points_fld, axis = 0)
mz_r_x0  = np.copy(ncfile.variables['mz_r_x0' ][:]); mz_r_x0  = np.delete(mz_r_x0 , ignored_points_fld, axis = 0)
mz_r_y0  = np.copy(ncfile.variables['mz_r_y0' ][:]); mz_r_y0  = np.delete(mz_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
wx_r_z0  = np.copy(ncfile.variables['wx_r_z0' ][:]); wx_r_z0  = np.delete(wx_r_z0 , ignored_points_fld, axis = 0)
wx_r_x0  = np.copy(ncfile.variables['wx_r_x0' ][:]); wx_r_x0  = np.delete(wx_r_x0 , ignored_points_fld, axis = 0)
wx_r_y0  = np.copy(ncfile.variables['wx_r_y0' ][:]); wx_r_y0  = np.delete(wx_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
wy_r_z0  = np.copy(ncfile.variables['wy_r_z0' ][:]); wy_r_z0  = np.delete(wy_r_z0 , ignored_points_fld, axis = 0)
wy_r_x0  = np.copy(ncfile.variables['wy_r_x0' ][:]); wy_r_x0  = np.delete(wy_r_x0 , ignored_points_fld, axis = 0)
wy_r_y0  = np.copy(ncfile.variables['wy_r_y0' ][:]); wy_r_y0  = np.delete(wy_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
wz_r_z0  = np.copy(ncfile.variables['wz_r_z0' ][:]); wz_r_z0  = np.delete(wz_r_z0 , ignored_points_fld, axis = 0)
wz_r_x0  = np.copy(ncfile.variables['wz_r_x0' ][:]); wz_r_x0  = np.delete(wz_r_x0 , ignored_points_fld, axis = 0)
wz_r_y0  = np.copy(ncfile.variables['wz_r_y0' ][:]); wz_r_y0  = np.delete(wz_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
bx_r_z0  = np.copy(ncfile.variables['bx_r_z0' ][:]); bx_r_z0  = np.delete(bx_r_z0 , ignored_points_fld, axis = 0)
bx_r_x0  = np.copy(ncfile.variables['bx_r_x0' ][:]); bx_r_x0  = np.delete(bx_r_x0 , ignored_points_fld, axis = 0)
bx_r_y0  = np.copy(ncfile.variables['bx_r_y0' ][:]); bx_r_y0  = np.delete(bx_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
by_r_z0  = np.copy(ncfile.variables['by_r_z0' ][:]); by_r_z0  = np.delete(by_r_z0 , ignored_points_fld, axis = 0)
by_r_x0  = np.copy(ncfile.variables['by_r_x0' ][:]); by_r_x0  = np.delete(by_r_x0 , ignored_points_fld, axis = 0)
by_r_y0  = np.copy(ncfile.variables['by_r_y0' ][:]); by_r_y0  = np.delete(by_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
bz_r_z0  = np.copy(ncfile.variables['bz_r_z0' ][:]); bz_r_z0  = np.delete(bz_r_z0 , ignored_points_fld, axis = 0)
bz_r_x0  = np.copy(ncfile.variables['bz_r_x0' ][:]); bz_r_x0  = np.delete(bz_r_x0 , ignored_points_fld, axis = 0)
bz_r_y0  = np.copy(ncfile.variables['bz_r_y0' ][:]); bz_r_y0  = np.delete(bz_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
jx_r_z0  = np.copy(ncfile.variables['jx_r_z0' ][:]); jx_r_z0  = np.delete(jx_r_z0 , ignored_points_fld, axis = 0)
jx_r_x0  = np.copy(ncfile.variables['jx_r_x0' ][:]); jx_r_x0  = np.delete(jx_r_x0 , ignored_points_fld, axis = 0)
jx_r_y0  = np.copy(ncfile.variables['jx_r_y0' ][:]); jx_r_y0  = np.delete(jx_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
jy_r_z0  = np.copy(ncfile.variables['jy_r_z0' ][:]); jy_r_z0  = np.delete(jy_r_z0 , ignored_points_fld, axis = 0)
jy_r_x0  = np.copy(ncfile.variables['jy_r_x0' ][:]); jy_r_x0  = np.delete(jy_r_x0 , ignored_points_fld, axis = 0)
jy_r_y0  = np.copy(ncfile.variables['jy_r_y0' ][:]); jy_r_y0  = np.delete(jy_r_y0 , ignored_points_fld, axis = 0)
                                                                                  
jz_r_z0  = np.copy(ncfile.variables['jz_r_z0' ][:]); jz_r_z0  = np.delete(jz_r_z0 , ignored_points_fld, axis = 0)
jz_r_x0  = np.copy(ncfile.variables['jz_r_x0' ][:]); jz_r_x0  = np.delete(jz_r_x0 , ignored_points_fld, axis = 0)
jz_r_y0  = np.copy(ncfile.variables['jz_r_y0' ][:]); jz_r_y0  = np.delete(jz_r_y0 , ignored_points_fld, axis = 0)

ux_r_z0  = mx_r_z0/rho_r_z0 
ux_r_x0  = mx_r_x0/rho_r_x0 
ux_r_y0  = mx_r_y0/rho_r_y0 
          
uy_r_z0  = my_r_z0/rho_r_z0
uy_r_x0  = my_r_x0/rho_r_x0 
uy_r_y0  = my_r_y0/rho_r_y0 
           
uz_r_z0  = mz_r_z0/rho_r_z0
uz_r_x0  = mz_r_x0/rho_r_x0 
uz_r_y0  = mz_r_y0/rho_r_y0 

ncfile.close()

#--------------------------------------------------------#
#                   plot final snapshot                  #
#--------------------------------------------------------#
u_r_z0 = np.sqrt(ux_r_z0**2 + uy_r_z0**2 + uz_r_z0**2)
u_r_y0 = np.sqrt(ux_r_y0**2 + uy_r_y0**2 + uz_r_y0**2)
u_r_x0 = np.sqrt(ux_r_x0**2 + uy_r_x0**2 + uz_r_x0**2)
w_r_z0 = np.sqrt(wx_r_z0**2 + wy_r_z0**2 + wz_r_z0**2)
w_r_y0 = np.sqrt(wx_r_y0**2 + wy_r_y0**2 + wz_r_y0**2)
w_r_x0 = np.sqrt(wx_r_x0**2 + wy_r_x0**2 + wz_r_x0**2)
b_r_z0 = np.sqrt(bx_r_z0**2 + by_r_z0**2 + bz_r_z0**2)
b_r_y0 = np.sqrt(bx_r_y0**2 + by_r_y0**2 + bz_r_y0**2)
b_r_x0 = np.sqrt(bx_r_x0**2 + by_r_x0**2 + bz_r_x0**2)
j_r_z0 = np.sqrt(jx_r_z0**2 + jy_r_z0**2 + jz_r_z0**2)
j_r_y0 = np.sqrt(jx_r_y0**2 + jy_r_y0**2 + jz_r_y0**2)
j_r_x0 = np.sqrt(jx_r_x0**2 + jy_r_x0**2 + jz_r_x0**2)
plot_3d(rho_r_z0[final_fld_idx,:,:], rho_r_y0[final_fld_idx,:,:], rho_r_x0[final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$\rho (t = %.2E)$'     % tt_fld[final_fld_idx], cmp=parula_map, save=outdir+'rho.pdf')
plot_3d(u_r_z0  [final_fld_idx,:,:], u_r_y0  [final_fld_idx,:,:], u_r_x0  [final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, (ux_r_z0[final_fld_idx,:,:], uy_r_z0[final_fld_idx,:,:]), (ux_r_y0[final_fld_idx,:,:], uz_r_y0[final_fld_idx,:,:]), (uy_r_x0[final_fld_idx,:,:], uz_r_x0[final_fld_idx,:,:]), xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\bm{u}| (t = %.2E)$' % tt_fld[final_fld_idx], cmp=parula_map, streamline_density=1, streamline_width=0.5, streamline_color='w', save=outdir+'u.pdf')
plot_3d(b_r_z0  [final_fld_idx,:,:], b_r_y0  [final_fld_idx,:,:], b_r_x0  [final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, (bx_r_z0[final_fld_idx,:,:], by_r_z0[final_fld_idx,:,:]), (bx_r_y0[final_fld_idx,:,:], bz_r_y0[final_fld_idx,:,:]), (by_r_x0[final_fld_idx,:,:], bz_r_x0[final_fld_idx,:,:]), xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\bm{B}| (t = %.2E)$' % tt_fld[final_fld_idx], cmp=parula_map, streamline_density=1, streamline_width=0.5, streamline_color='w', save=outdir+'b.pdf')
if is2D:
  plot_3d(wz_r_z0[final_fld_idx,:,:], wz_r_y0[final_fld_idx,:,:], wz_r_x0[final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$(\Curl\bm{u})_z (t = %.2E)$' % tt_fld[final_fld_idx], cmp='RdBu_r', save=outdir+'w.pdf')
  plot_3d(jz_r_z0[final_fld_idx,:,:], jz_r_y0[final_fld_idx,:,:], jz_r_x0[final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$(\Curl\bm{B})_z (t = %.2E)$' % tt_fld[final_fld_idx], cmp='RdBu_r', save=outdir+'j.pdf')
else:
  plot_3d( w_r_z0[final_fld_idx,:,:],  w_r_y0[final_fld_idx,:,:],  w_r_x0[final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\Curl\bm{u}| (t = %.2E)$' % tt_fld[final_fld_idx], cmp=parula_map, save=outdir+'w.pdf')
  plot_3d( j_r_z0[final_fld_idx,:,:],  j_r_y0[final_fld_idx,:,:],  j_r_x0[final_fld_idx,:,:], xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\Curl\bm{B}| (t = %.2E)$' % tt_fld[final_fld_idx], cmp=parula_map, save=outdir+'j.pdf')


#--------------------------------------------------------#
#                       plot movie                       #
#--------------------------------------------------------#
if ismovie:
  # evenly space time
  # nframe = min(200, int(tt[:final_fld_idx].size))
  nframe = tt_fld[:final_fld_idx].size
  idx = np.unique([np.argmin(abs(tt_fld - np.linspace(tt_fld[0], tt_fld[-1], nframe)[i])) for i in range(0, nframe)])
  tt_fld  = tt_fld [idx]

  rho_r_z0 = np.take(rho_r_z0, idx, axis=0)
  rho_r_y0 = np.take(rho_r_y0, idx, axis=0)
  rho_r_x0 = np.take(rho_r_x0, idx, axis=0)

  ux_r_z0  = np.take(ux_r_z0 , idx, axis=0)
  ux_r_y0  = np.take(ux_r_y0 , idx, axis=0)
  ux_r_x0  = np.take(ux_r_x0 , idx, axis=0)
  uy_r_z0  = np.take(uy_r_z0 , idx, axis=0)
  uy_r_y0  = np.take(uy_r_y0 , idx, axis=0)
  uy_r_x0  = np.take(uy_r_x0 , idx, axis=0)
  uz_r_z0  = np.take(uz_r_z0 , idx, axis=0)
  uz_r_y0  = np.take(uz_r_y0 , idx, axis=0)
  uz_r_x0  = np.take(uz_r_x0 , idx, axis=0)
                             
  wx_r_z0  = np.take(wx_r_z0 , idx, axis=0)
  wx_r_y0  = np.take(wx_r_y0 , idx, axis=0)
  wx_r_x0  = np.take(wx_r_x0 , idx, axis=0)
  wy_r_z0  = np.take(wy_r_z0 , idx, axis=0)
  wy_r_y0  = np.take(wy_r_y0 , idx, axis=0)
  wy_r_x0  = np.take(wy_r_x0 , idx, axis=0)
  wz_r_z0  = np.take(wz_r_z0 , idx, axis=0)
  wz_r_y0  = np.take(wz_r_y0 , idx, axis=0)
  wz_r_x0  = np.take(wz_r_x0 , idx, axis=0)
                             
  bx_r_z0  = np.take(bx_r_z0 , idx, axis=0)
  bx_r_y0  = np.take(bx_r_y0 , idx, axis=0)
  bx_r_x0  = np.take(bx_r_x0 , idx, axis=0)
  by_r_z0  = np.take(by_r_z0 , idx, axis=0)
  by_r_y0  = np.take(by_r_y0 , idx, axis=0)
  by_r_x0  = np.take(by_r_x0 , idx, axis=0)
  bz_r_z0  = np.take(bz_r_z0 , idx, axis=0)
  bz_r_y0  = np.take(bz_r_y0 , idx, axis=0)
  bz_r_x0  = np.take(bz_r_x0 , idx, axis=0)
                             
  jx_r_z0  = np.take(jx_r_z0 , idx, axis=0)
  jx_r_y0  = np.take(jx_r_y0 , idx, axis=0)
  jx_r_x0  = np.take(jx_r_x0 , idx, axis=0)
  jy_r_z0  = np.take(jy_r_z0 , idx, axis=0)
  jy_r_y0  = np.take(jy_r_y0 , idx, axis=0)
  jy_r_x0  = np.take(jy_r_x0 , idx, axis=0)
  jz_r_z0  = np.take(jz_r_z0 , idx, axis=0)
  jz_r_y0  = np.take(jz_r_y0 , idx, axis=0)
  jz_r_x0  = np.take(jz_r_x0 , idx, axis=0)

  u_r_z0 = np.sqrt(ux_r_z0**2 + uy_r_z0**2 + uz_r_z0**2)
  u_r_y0 = np.sqrt(ux_r_y0**2 + uy_r_y0**2 + uz_r_y0**2)
  u_r_x0 = np.sqrt(ux_r_x0**2 + uy_r_x0**2 + uz_r_x0**2)
  w_r_z0 = np.sqrt(wx_r_z0**2 + wy_r_z0**2 + wz_r_z0**2)
  w_r_y0 = np.sqrt(wx_r_y0**2 + wy_r_y0**2 + wz_r_y0**2)
  w_r_x0 = np.sqrt(wx_r_x0**2 + wy_r_x0**2 + wz_r_x0**2)
  b_r_z0 = np.sqrt(bx_r_z0**2 + by_r_z0**2 + bz_r_z0**2)
  b_r_y0 = np.sqrt(bx_r_y0**2 + by_r_y0**2 + bz_r_y0**2)
  b_r_x0 = np.sqrt(bx_r_x0**2 + by_r_x0**2 + bz_r_x0**2)
  j_r_z0 = np.sqrt(jx_r_z0**2 + jy_r_z0**2 + jz_r_z0**2)
  j_r_y0 = np.sqrt(jx_r_y0**2 + jy_r_y0**2 + jz_r_y0**2)
  j_r_x0 = np.sqrt(jx_r_x0**2 + jy_r_x0**2 + jz_r_x0**2)
  movie_3d(tt_fld, rho_r_z0, rho_r_y0, rho_r_x0, xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$\rho$'    , cmp=parula_map, save=outdir+'rho_anim.gif')
  movie_3d(tt_fld, u_r_z0  , u_r_y0, u_r_x0    , xx_fld, yy_fld, zz_fld, (ux_r_z0, uy_r_z0), (ux_r_y0, uz_r_y0), (uy_r_x0, uz_r_x0), xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\bm{u}|$', cmp=parula_map, save=outdir+'u_anim.gif')
  movie_3d(tt_fld, b_r_z0  , b_r_y0, b_r_x0    , xx_fld, yy_fld, zz_fld, (bx_r_z0, by_r_z0), (bx_r_y0, bz_r_y0), (by_r_x0, bz_r_x0), xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\bm{B}|$', cmp=parula_map, save=outdir+'b_anim.gif')
  if is2D:
    movie_3d(tt_fld, wz_r_z0, wz_r_y0, wz_r_x0, xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$(\Curl\bm{u})_z$', cmp='RdBu_r', save=outdir+'w_anim.gif')
    movie_3d(tt_fld, jz_r_z0, jz_r_y0, jz_r_x0, xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$(\Curl\bm{B})_z$', cmp='RdBu_r', save=outdir+'j_anim.gif')
  else:
    movie_3d(tt_fld,  w_r_z0,  w_r_y0,  w_r_x0, xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\Curl\bm{u}|$', cmp=parula_map, save=outdir+'w_anim.gif')
    movie_3d(tt_fld,  j_r_z0,  j_r_y0,  j_r_x0, xx_fld, yy_fld, zz_fld, xlab=xlab, ylab=ylab, zlab=zlab, title=r'$|\Curl\bm{B}|$', cmp=parula_map, save=outdir+'j_anim.gif')


#------------------#
#   output data    #
#------------------#
from scipy.io import savemat

savemat(outdir + 'grid' , {'xx':xx, 'yy':yy, 'zz':zz})
savemat(outdir + 'rho_r' , {
                            'tt'       :tt[final_fld_idx],
                            'rho_r_z0' :rho_r_z0[final_fld_idx,:,:],
                            'rho_r_x0' :rho_r_x0[final_fld_idx,:,:],
                            'rho_r_y0' :rho_r_y0[final_fld_idx,:,:],
                          })
savemat(outdir + 'u_r' , {
                            'tt'       :tt[final_fld_idx],
                            'ux_r_z0'  :ux_r_z0 [final_fld_idx,:,:],
                            'ux_r_x0'  :ux_r_x0 [final_fld_idx,:,:],
                            'ux_r_y0'  :ux_r_y0 [final_fld_idx,:,:],
                            #
                            'uy_r_z0'  :uy_r_z0 [final_fld_idx,:,:],
                            'uy_r_x0'  :uy_r_x0 [final_fld_idx,:,:],
                            'uy_r_y0'  :uy_r_y0 [final_fld_idx,:,:],
                            #
                            'uz_r_z0'  :uz_r_z0 [final_fld_idx,:,:],
                            'uz_r_x0'  :uz_r_x0 [final_fld_idx,:,:],
                            'uz_r_y0'  :uz_r_y0 [final_fld_idx,:,:],
                          })
savemat(outdir + 'b_r' , {
                            'tt'       :tt[final_fld_idx],
                            'bx_r_z0'  :bx_r_z0 [final_fld_idx,:,:],
                            'bx_r_x0'  :bx_r_x0 [final_fld_idx,:,:],
                            'bx_r_y0'  :bx_r_y0 [final_fld_idx,:,:],
                            #
                            'by_r_z0'  :by_r_z0 [final_fld_idx,:,:],
                            'by_r_x0'  :by_r_x0 [final_fld_idx,:,:],
                            'by_r_y0'  :by_r_y0 [final_fld_idx,:,:],
                            #
                            'bz_r_z0'  :bz_r_z0 [final_fld_idx,:,:],
                            'bz_r_x0'  :bz_r_x0 [final_fld_idx,:,:],
                            'bz_r_y0'  :bz_r_y0 [final_fld_idx,:,:],
                          })
savemat(outdir + 'w_r' , {
                            'tt'       :tt[final_fld_idx],
                            'wx_r_z0'  :wx_r_z0 [final_fld_idx,:,:],
                            'wx_r_x0'  :wx_r_x0 [final_fld_idx,:,:],
                            'wx_r_y0'  :wx_r_y0 [final_fld_idx,:,:],
                            #
                            'wy_r_z0'  :wy_r_z0 [final_fld_idx,:,:],
                            'wy_r_x0'  :wy_r_x0 [final_fld_idx,:,:],
                            'wy_r_y0'  :wy_r_y0 [final_fld_idx,:,:],
                            #
                            'wz_r_z0'  :wz_r_z0 [final_fld_idx,:,:],
                            'wz_r_x0'  :wz_r_x0 [final_fld_idx,:,:],
                            'wz_r_y0'  :wz_r_y0 [final_fld_idx,:,:],
                          })
savemat(outdir + 'j_r' , {
                            'tt'       :tt[final_fld_idx],
                            'jx_r_z0'  :jx_r_z0 [final_fld_idx,:,:],
                            'jx_r_x0'  :jx_r_x0 [final_fld_idx,:,:],
                            'jx_r_y0'  :jx_r_y0 [final_fld_idx,:,:],
                            #
                            'jy_r_z0'  :jy_r_z0 [final_fld_idx,:,:],
                            'jy_r_x0'  :jy_r_x0 [final_fld_idx,:,:],
                            'jy_r_y0'  :jy_r_y0 [final_fld_idx,:,:],
                            #
                            'jz_r_z0'  :jz_r_z0 [final_fld_idx,:,:],
                            'jz_r_x0'  :jz_r_x0 [final_fld_idx,:,:],
                            'jz_r_y0'  :jz_r_y0 [final_fld_idx,:,:],
                          })

# -*- coding: utf-8 -*- # Load 2D cut #--------------------------------------------------------# #                   plot final snapshot                  # #--------------------------------------------------------# #--------------------------------------------------------# #                       plot movie                       # #--------------------------------------------------------# # evenly space time # nframe = min(200, int(tt[:final_fld_idx].size)) #------------------# #   output data    # #------------------# # # # # # # # #