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")