image_obj = obj_new('hsi_image',sim_pho=5000,pixel_size=1.+fltarr(2),$
sim_time_range=[0.,4],time_range=[0.,4.],time_bin_min=1024)
xyoffset = [600.,200.]
image = image_obj->getdata(weight=0, sim_xyoffset=xyoffset, xyoffset=xyoffset)
plot_image, image
;Add the weighting (compensate for unequal modulation in each pixel).

image = image_obj->getdata(weight=1, sim_xyoffset=xyoffset, xyoffset=xyoffset)
plot_image, image
chkarg,'hsi_annsec2xy
; % CHKARG: found in current directory
;---- Module: hsi_annsec2xy.pro
;---- From:   /usr/users/richard/annsec/
;---> Call: function hsi_annsec2xy,annsecimage, modpat_obj, xgrid, ygrid, $
;this_harmonic=this_harmonic, this_det_index=this_det_index, $
;max_value=max_value, min_value=min_value, quintic=quintic,$
;limits=limits, modpat_ptr=modpat_ptr, gs=gs, missing=missing
xyimage = hsi_annsec2xy( image, image_obj, xgrid, ygrid, gs=1+fltarr(2))
chkarg,'plot_image
;---- Module: plot_image.pro
;---- From:   /data/ssw/gen/idl/image/
;---> Call: PRO PLOT_IMAGE,IMAGE,ORIGIN=ORIGIN,SCALE=SCALE,NOERASE=NOERASE, $
;NOSQUARE=NOSQUARE,SMOOTH=SMOOTH,NOSCALE=NOSCALE,	$
;MISSING=MISSING,COLOR=COLOR,MAX=MAX,MIN=MIN,TOP=TOP,	$
;VELOCITY=VELOCITY,COMBINED=COMBINED,LOWER=LOWER,	$
;NOADJUST=NOADJUST,TITLE=TITLE,XTITLE=XTITLE,YTITLE=YTITLE, $
;CHARSIZE=CHARSIZE,XTICK_GET=XTICK_GET,YTICK_GET=YTICK_GET,$
;_EXTRA=_EXTRA
pmm,xgrid
;      568.500      631.500
pmm,ygrid
;      168.500      231.500
plot_image,xyimage,origin=[xgrid[0],ygrid[0]],scale=1.+fltarr(2)
model=image_obj->get(/sim_model)
print,model
;{      1.00000     0.500000      1.00000
;     -11.5000     -11.5000
;      0.00000}{      1.00000
;     0.500000      1.00000
;     0.500000     0.500000
;      0.00000}
print,model[1].xypos
;     0.500000     0.500000
print,xgrid[32]
;      600.500
print,ygrid[32]
;      200.500
print,ygrid[33]
;      201.500
help,image
hsi_annsec_coord,image_obj, x, y
print,x[32]
;      599.149
print,y[33]
;      166.963
print,x[32,33]
;      600.000
print,y[32,33]
;      200.000
pixel=[0,0.]
cbe=image_obj->getdata(class_name='hsi_calib_eventlist',this_det=0,this_harm=0,/time)
asp = *(hsi_aspect_sim(cbe.time,/time_unit))
help,asp
help,asp,/st
;mps0=hsi_modulate_point_source(pixel,
chkarg,'hsi_modulate_point_source
;---- Module: hsi_modulate_point_source.pro
;---- From:   /usr/users/hessi/sw/idl/image/
;---> Call: FUNCTION hsi_modulate_point_source, pixel, map_offset, aspect, a2d_index, energy_band, $
;dispersion, harmonic_max=harmonic_max, phase_map_ctr=phase_map_ctr
mps0 = hsi_modulate_point_source( pixel, xyoffset, asp, 0, [6,100.])
b = image
b=b*0.0
b[32,33]=1.

prf0=image_obj->getdata(class_name='hsi_modul_profile', vimage=b,this_det=0,this_harm=0)
plot,prf0[0:255],psy=10
oplot, mps0/1024., psym=1
print,avg(mps0[0:3999])
;     0.270954
print,avg(prf0[0:3999]*1024)
;     0.270965
help,x,y
print,x[64,66]
;      630.369
print,y[64,66]
;      234.597
b=b*0
b[64,66]=1.
prf0=image_obj->getdata(class_name='hsi_modul_profile', vimage=b,this_det=0,this_harm=0)
pixel = [x[64,66],y[64,66]]
mps0 = hsi_modulate_point_source( pixel, xyoffset, asp, 0, [6,100.])
plot,prf0[0:255],psy=10
plot,prf0[0:255]*1024,psy=10
oplot, mps0, psym=1
pixel = [x[64,66],y[64,66]]-xyoffset
mps0 = hsi_modulate_point_source( pixel, xyoffset, asp, 0, [6,100.])
plot,prf0[0:255]*1024,psy=10
oplot, mps0, psym=1
print,image_obj->get(/modpat_skip)
;       0
print,image_obj->get(/modpat_step)
;      -1
image_obj->set,modpat_step = 4
mps0 = hsi_modulate_point_source( pixel, xyoffset, asp, 0, [6,100.])
image_obj->set,modpat_skip = 4
prf0=image_obj->getdata(class_name='hsi_modul_profile', vimage=b,this_det=0,this_harm=0,sim_pho=6000)
plot,prf0[0:255]*1024,psy=10
oplot, mps0, psym=1

;Check the Normalization!
o=obj_new(sim_pho=5000,sim_time_range=[0,4.],time_range=[0.,4.],'hsi_image')
imagewa=o->getdata(weight=1,/no_sum)
imagea=o->getdata(weight=0,/no_sum)
npix=65.*67
for i=0,8 do print,total(*imagea[i])/npix                         
;      1519.91
;      1514.96
;      1544.28
;      1536.33
;      1496.46
;      1683.84
;      1585.81
;      1500.79
;      2131.06
;Result should be close to the grid transmission coefficient squared times the flux
;for pitches small wrt the fov. So for 2e4 photons we have:
for i=0,8 do print,total(*imagea[i])/npix,2e4*avg((*cbe[i]).gridtran)^2
;      1519.91      1486.07
;      1514.96      1478.79
;      1544.28      1486.07
;      1536.33      1492.69
;      1496.46      1494.22
;      1683.84      1490.96
;      1585.81      1299.50
;      1500.79      1262.61
;      2131.06      1414.03


chkarg,'hsi_annsec_bproj
;---- Module: hsi_annsec_bproj.pro
;---- From:   /usr/users/hessi/sw/idl/util/
;---> Call: function hsi_annsec_bproj, cbe_obj, $
;map_ptr=map_ptr, det_index_mask=det_index_mask, $
;sum=sum,  time_test=time_test, weight=weight, $
;this_harmonic=this_harmonic, this_det_index=this_det_index, data_ptr=data_ptr, $
;counts_summed=counts_summed
map_ptr=o->getdata(class_name='hsi_modul_pattern')
help,*map_ptr[0],/st
chkarg,'hsi_annsec_bproj
; % CHKARG: recalling hsi_annsec_bproj.pro from memory
;---- Module: hsi_annsec_bproj.pro
;---- From:   /usr/users/hessi/sw/idl/util/
;---> Call: function hsi_annsec_bproj, cbe_obj, $
;map_ptr=map_ptr, det_index_mask=det_index_mask, $
;sum=sum,  time_test=time_test, weight=weight, $
;this_harmonic=this_harmonic, this_det_index=this_det_index, data_ptr=data_ptr, $
;counts_summed=counts_summed
bprojw=hsi_annsec_bproj(/weight,o,map_ptr=map_ptr)
help,bprojw
pmm,*bprojw[0]-*imagewa[0]
;      0.00000      0.00000
help,imge
help,image
vimage=image*0
vimage[32,33]=1.0
rate1=o->getdata(class_name='hsi_modul_profile',vimage=vimage)
help,rate1
for i=0,8 do print,total(*rate1[i])
;      1.08865
;      1.06244
;      1.11107
;      1.11767
;      1.06806
;      1.11469
;      1.06634
;     0.951240
;      1.17750
for i=0,8 do print,total(*rate1[i]) ;That is one photon/sec over 4 sec. About 26-28% make it through both grids.
;      1.08865
;      1.06244
;      1.11107
;      1.11767
;      1.06806
;      1.11469
;      1.06634
;     0.951240
;      1.17750
;Now run that rate through back projection.
psf3233 = o->getdata(class_name='hsi_bproj',vrate=rate1,/no_sum,weight=0)
for i=0,8 do print,total(*psf3233[i])/npix
;     0.303169
;     0.294242
;     0.311326
;     0.313976
;     0.291315
;     0.321766
;     0.330895
;     0.296552
;     0.431511
psf3233w = o->getdata(class_name='hsi_bproj',vrate=rate1,/no_sum,weight=1)
for i=0,8 do print,total(*psf3233w[i])/npix
;     0.300021
;     0.254435
;     0.388131
;     0.397011
;   -0.0555827
;     0.763966
;      2.50811
;      3.41857
;      3.87872
plot_image,*psf3233w[1]
wshow
plot_image,*psf3233w[0]
xsurface,*psf3233w[0]
print,total( (*psf3233w[0])[32:32,33:33])
;      4.00007
print,total( (*psf3233w[1])[32:32,33:33])
;      3.99999
print,total( (*psf3233w[2])[32:32,33:33])
;      4.00002
print,total( (*psf3233w[3])[32:32,33:33])
;      4.00000
print,total( (*psf3233w[4])[32:32,33:33])
;      4.00000
for i=0,8 do print,total( (*psf3233w[i])[32:32,33:33])
;      4.00007
;      3.99999
;      4.00002
;      4.00000
;      4.00000
;      4.00001
;      4.00000
;      4.00000
;      4.00000
for i=0,8 do print,total( (*psf3233w[i])[31:31,33:33])
;      2.37569
;      3.44008
;      3.79385
;      3.94856
;      4.01157
;      3.97511
;      4.00992
;      4.00926
;      4.00562
for i=0,8 do print,total( (*psf3233w[i])[31:31,32:32])
;      1.14718
;      2.90960
;      3.60931
;      3.86802
;      3.99890
;      3.96272
;      4.01089
;      4.00939
;      4.00725
for i=0,8 do print,total(*psf3233w[i])/npix * 5000.
;      1500.11
;      1272.18
;      1940.66
;      1985.06
;     -277.913
;      3819.83
;      12540.5
;      17092.8
;      19393.6
for i=0,8 do print, total(*imagewa[i])/npix,total(*imagea[i])/npix
;      1438.31      1519.91
;      1540.92      1514.96
;      1521.81      1544.28
;      1536.25      1536.33
;      727.703      1496.46
;      3669.82      1683.84
;      11245.0      1585.81
;      16076.7      1500.79
;      18992.6      2131.06
end