pro acis_pileup,pha0,time,ypos,zpos,frametime,pixsize,pha,cnt
;-----------------------------------------------------------------------
; Name: ACIS_PILEUP
;
; Purpose: Reads a MARX events list and simulates the effects of pileup
;          
; Inputs: 
;         
;         
; Comments: 
;           
;           
; Revision history:
;       written by Michael Wise, 01-20-97
;	based on code from Dave Huenemoerder 
;-----------------------------------------------------------------------
;
; Check for correct number of parameters
;
np=n_params(0)
if ((np lt 8) or (np gt 8)) then begin
   print,string(7B),'CALLING SEQUENCE: ', $
   'acis_pileup,pha0,time,ypos,zpos,frametime,pixsize,pha,cnt'
   return   
endif
;
; Set defaults
;
if (np lt 5) then frametime=3.0
if (np lt 6) then pixsize=0.024			; Pixel size in mm

;
; Determine how many frames we have
;
tlo=min(time)
thi=max(time)
nframes=long((thi-tlo)/frametime)+1

;
; Arrays to hold "piled" events
;
pha=0.0
cnt=0.0

;
; Loop over frames
;
for i=0L,nframes-1 do begin

;
;   Select indices of photons in current "frame"
;
    tl=tlo+(i*frametime)
    th=tlo+((i+1)*frametime)
    ic=where( (time ge tl) and (time lt th) )


;
;   See if we have an empty frame
;
    if (ic(0) ne -1) then begin

;
;      Extract relevant sub-arrays
;
       tc=time(ic)
       yc=ypos(ic)
       zc=zpos(ic)
       pc=pha0(ic)

;
;      Bin the subarray up into a "frame"
;
       xmin=(fix(min(yc)/pixsize)-2)*pixsize
       xmax=(fix(max(yc)/pixsize)+2)*pixsize
       ymin=(fix(min(zc)/pixsize)-2)*pixsize
       ymax=(fix(max(zc)/pixsize)+2)*pixsize
       imc=make_image(yc,zc,xbinsize=pixsize,ybinsize=pixsize, $
                      xrange=[xmin,xmax],yrange=[ymin,ymax], $
                      xaxis=yax,yaxis=zax, $
                      index_list=idx,reverse_indices=revidx)

;
;      Trim pha array to include only those which made into 
;      into the binned image
;
       lpc=pc(idx)

;
;      Loop over bins and add up pulse heights
;
       imp=imc*0.0
       ibin=where(imc gt 0)
       nbin=n_elements(ibin)

       for j=0,nbin-1 do begin
        
           pixnum = ibin(j)  ;;; serial bin number

;
;          If there is a contribution (this should always be true since
;          we already selected non-zero bins)       
;
           if ( revidx(pixnum) ne revidx(pixnum+1)) then begin
;
;             Sum up contributing pulse-heights
;   
              plist=revidx(revidx(pixnum):revidx(pixnum+1)-1)
              imp(pixnum)=total(lpc(plist))      

           endif else begin
              print,string(7B),'ERROR: No contribution to this non-zero bin!'
              return   
           endelse

       endfor    

;
;      Now "extract" piled counts
;
       ip=where(imc gt 0.0)
       pha=[pha,imp(ip)]    
       cnt=[cnt,imc(ip)]


    endif

;
;   Done with this frame
;
endfor


;
; Trim "detected" arrays
;
num=n_elements(pha)
pha=pha(1:num-1)
cnt=cnt(1:num-1)

;
; Set maximum value to 4096
;
;ib=where(pha ge 4096.0)
;if (ib(0) ne -1) then pha(ib)=4096.0


;
; Return to IDL
;
return
end