STEPS FOR PPV PROJECTION:

; build a projection with no ghost zones, no grids with smaller spacings than 1/256, and no axis swapping
proj = build_projection(INTERP=0, IMAGE_SIZE=256, /NOSWAPPING)

; calculate data cube
n=256
data = DBLARR(401,n,n)
for x=0L,n-1 do for y=0L,n-1 do data[0,x,y] = construct_ppv_projection(PROJ=proj, USE=[1,1,1,1], PT=([x,y]+.5d)/n, /INTEGRATE)

noise = min(data, DIMENSION=1)
data2 = data - TRANSPOSE(CMREPLICATE(noise, 401), [2,0,1])
; find statistical moments
x = CMREPLICATE(INDGEN(401),[n,n]) * 40. / 400. - 20.
tot = REFORM(TOTAL(data2,1))
P = data2 / TRANSPOSE(CMREPLICATE(tot, 401), [2,0,1])
m = REFORM(TOTAL(x*P,1))
v = REFORM(TOTAL((x-TRANSPOSE(CMREPLICATE(m, 401), [2,0,1]))^2 * P,1))
sigma = v^.5
x = x[*,0,0]

; plot! e.g.
loadct, 13
tvscl, sigma
plot, x, data[*,128,128]