ps, here is my script: import matplotlib matplotlib.use("Agg") from yt.mods import * import pylab import sys import numpy as na def _MyRadius(field, data): center = data.get_field_parameter("center") dx = data["x"] - center[0] dy = data["y"] - center[1] dz = data["z"] - center[2] return na.sqrt(dx*dx + dy*dy + dz*dz) add_field("Radius", function=_MyRadius, take_log=False, units='') def _ConvertAccel(data): return data.convert("Length") / (data.convert("Time"))**2.0 def _RadiationAccelerationMagnitude(field, data): return ( (data["RadAccel1"])**2.0 + \ (data["RadAccel2"])**2.0 + \ (data["RadAccel3"])**2.0 )**(1.0/2.0) #return ( (data["Density"])*0.0+1.e-28) add_field("RadiationAccelerationMagnitude", function=_RadiationAccelerationMagnitude, take_log=False,convert_function=_ConvertAccel, units=r"\rm{cm}/\rm{s}^{2}") def _RadiationForceMagnitude(field, data): return ( data["RadiationAccelerationMagnitude"] * data["Density"] * data["CellVolume"]) add_field("RadiationForceMagnitude", function=_RadiationForceMagnitude, take_log=False, units=r"\rm{dynes}") def _RadiationForce1(field, data): return ( (data["RadAccel1"] * data["Density"] * data["CellVolume"])) add_field("RadiationForce1", function=_RadiationForce1, take_log=False, units=r"\rm{dynes}") def _RadiationForce2(field, data): return (( data["RadAccel2"] * data["Density"] * data["CellVolume"])) add_field("RadiationForce2", function=_RadiationForce2, take_log=False, units=r"\rm{dynes}") def _RadiationForce3(field, data): return (( data["RadAccel3"] * data["Density"] * data["CellVolume"])) add_field("RadiationForce3", function=_RadiationForce3, take_log=False, units=r"\rm{dynes}") xc = 0.5 yc = 0.5 zc = 0.5 width = 179.2 frame_template = "aaSN/%s/RadPresProj_%s_%s" folder = sys.argv[1] starti = int(sys.argv[2]) endi = int(sys.argv[3]) Q = 1.e47 RPSF = 2.0 nu = 21.62 * 1.6022e-12 # in ergs print "Q:",Q print "RPScaleFactor:", RPSF print "nu:",nu rho_min = 1.e24 rho_max = 5.e27 for i in range(starti, endi ,1): print "a, i=",i pf = load("SciNet/%s/DD%04i/data%04i" % (folder, i, i)) if pf is None: continue sp = pf.h.sphere([xc,yc,zc], 1.0) v, c = pf.h.find_max("Density") print i, c print v,c pc = PlotCollection(pf, center=(xc,yc,zc)) total = 0 print "Initial time, code units =", pf["InitialTime"] for ax in range(0,1): p = pc.add_projection("RadiationForceMagnitude", ax, weight_field="Density") p.set_log_field(True) p.set_cmap('Blues') p.set_zlim(1e24, 1e28) p.set_width(width,"pc") sp= pf.h.all_data() v1 = "RadiationForce2" v2 = "RadiationForce3" p.modify["quiver"](v1, v2, 12) if na.any(pf.h.grid_particle_count): colours = sp["ParticleMassMsun"] p.modify["particles"](1.0, p_size=10.0, col='r', minimum_mass= 6e0) print " XXXXX Particles present XXXXX" Timemyr = pf["InitialTime"]*pf["years"]/1.e6 p.modify["text"]((0.006,0.899),"%6.1f Myr" % (Timemyr), text_args = {'color':'w','fontsize':'30'}) p.modify["text"]((0.005,0.9),"%6.1f Myr" % (Timemyr), text_args = {'color':'b','fontsize':'30' }) pc.save(frame_template % (sys.argv[1],ax,i), override=True, format="eps")