pro tcticorr,evtfile,mtl,asol,outroot=outroot,coeff_file=coeff_file,ciaopath=ciaopath,verbose=verbose,zero=zero,filter=filter,float=float,niter=niter

;+
;procedure         tcticorr
;       applies a temperature dependent adjustment to the CTI
;       correction for FI and BI chips and creates a new level 1 evt list
;
;syntax
;       pro tcticorr,evt,mtl,asol,outroot=outroot,$
;               coeff_file=coeff_file,ciaopath=ciaopath,$
;               niter=niter,/filter,/verbose,/zero,/float
;
;parameters
;       evt        [string; required] level 1 event list FITS file
;       mtl        [string; required] mission time line (MTL) FITS file
;       asol       [string; required] aspect solution file
;       outroot    [string; optional] base of names for output files
;                  file which has CTI applied but PHA and fltgrade not
;                  yet recalculated will be called outroot_tcticorr.fits
;                  regraded file will be called outroot_tcticorr_evt1.fits
;                  *default='obs{obsid}'
;       coeff_file [string; optional] text file with parallel and serial CTI change per
;                  degree coefficients for each chip
;                  *default='tcticorr_coeff_v1.txt'
;       ciaopath   [string; optional] directory where CIAO is installed
;                  *default='/soft/ciao/'
;
;       niter      [integer; optional] maximum number of CTI
;                  correction iterations to perform 
;                  *default = 15     
;
;keywords
;       filter     if set, creates a level 2 event list by filtering
;                  the final level 1 event list (outroot_tcticorr_ape.fits) on the
;                  standard filter file for the observation
;                  (acisf*_flt1.fits), status=0, and grade=0,2,3,4,6.
;                  the level 2 event list will be called outroot_tcticorr_evt2.fits
;       verbose    if set, prints out some messages about
;                  what's happening
;       zero       if set, sets tcti_parallel and tcti_serial arrays
;                  to zero (i.e. CTI correction is performed without
;                  temperature dependent correction) 
;       float      if set, writes PHAS column out as floats (instead
;                  of integers).  this does not affect the final PHA
;                  values; this keyword was used for testing purposes
;
;history
;    modified from Catherine Grant's cticorr_tdependent.pro
;    and Glenn Allen's apply_cti_2.5.sl
;
;    v5.1, 2/9/10:   
;    added basic parameter checking and ciaopath parameter
;
;    v5.2, 2/18/10 :
;    **critical bug fix! -- fixed mwrfits statement so that STATUS
;                           column will be written out correctly and
;                           added dmpaste command to compensate for 
;                            resetting of status bit 1 when regrading 
;                            with null trapmap
;    fixed datatype for TIME
;    added dmhedit command to add TCTICORR keyword to header of final 
;       output file
;    added filter keyword to optionally produce level 2 event list
;
;
;    v5.3, 5/6/10:
;    ** another bug fix - fixed error in how tcti coefficients were
;                         being indexed
;    added niter parameter and float keyword
;
;    v5.4, 9/30/10:
;    Changed how temperature for a given event is selected, to be
;    consistent with in-development CIAO version.  This change does
;    not affect final PHA values by a significant amount compared to
;    to v5.3.
;
;-

np=n_params()
if np lt 3 then begin
  print, 'please supply required inputs: evtfile, mtl, asol'
  print, 'syntax:'
  print, 'tcticorr,evtfile,mtl,asol,outroot=outroot,coeff_file=coeff_file,ciaopath=ciaopath,/verbose,/zero'
  print, 'apply temperature-dependent CTI correction'
  stop
endif

; read in event list
  evt1=mrdfits(evtfile,1,hevt)
  
  obsdate=sxpar(hevt,'DATE-OBS')
  detnam=sxpar(hevt,'DETNAM')
  observer=sxpar(hevt,'OBSERVER')
  datamode=strcompress(sxpar(hevt,'DATAMODE'),/remove_all)
  gainfile=sxpar(hevt,'GAINFILE')
  tgainfile=sxpar(hevt,'TGAINFIL')
  bpixfile=sxpar(hevt,'BPIXFILE')
  obsid=sxpar(hevt,'OBS_ID')

  if keyword_set(filter) then fltfile=sxpar(hevt,'FLTFILE')

  if datamode eq 'FAINT' then begin
     nmd=3
     i_x_off=0
     i_y_off=0
  endif
  if datamode eq 'VFAINT' then begin
     nmd=5  
     i_x_off=1
     i_y_off=1
  endif

  ccds=(evt1.ccd_id)[UNIQ(evt1.ccd_id, SORT(evt1.ccd_id))] 
  nccd=n_elements(ccds)

  ccd_id=evt1.ccd_id
  event_time=evt1.time

  chx=evt1.chipx
  chy=evt1.chipy

  phas=evt1.phas
  
  phas_max = 4095

; read in mtl file
  mtlfile=mtl
  mtl=mrdfits(mtlfile,1,hmtl)   ; check - what is diff. between MTL and MTL_S extensions?
  fpt=mtl.fp_temp                          ; focal plane temperature (Kelvin)

  timedel_mtl=sxpar(hmtl,'TIMEDEL')
  timepixr_mtl=sxpar(hmtl,'TIMEPIXR')
  time_mtl=mtl.TIME + timedel_mtl*(0.5-timepixr_mtl)

  gea_time=[0,time_mtl,1.0e10]
  gea_fp_temp=[(mtl.fp_temp)[0],mtl.fp_temp,(mtl.fp_temp)[n_elements(mtl.fp_temp)-1]]

  
; set up outfile
  if keyword_set(outroot) then outname=outroot else outname='obs'+strcompress(obsid,/remove_all)
  evtout=outname+'_tcticorr_intermediate1.fits'


  if keyword_set(float) then outstruc={time:0.0D,ccd_id:0,node_id:0,expno:0L,chipx:0,chipy:0,tdetx:0,tdety:0,detx:0.0,dety:0.0,x:0.0,y:0.0,phas:fltarr(nmd,nmd),pha:0L,pha_ro:0L,energy:0.0,pi:0L,fltgrade:0,grade:0,status:bytarr(4)} else outstruc={time:0.0D,ccd_id:0,node_id:0,expno:0L,chipx:0,chipy:0,tdetx:0,tdety:0,detx:0.0,dety:0.0,x:0.0,y:0.0,phas:intarr(nmd,nmd),pha:0L,pha_ro:0L,energy:0.0,pi:0L,fltgrade:0,grade:0,status:bytarr(4)} ; make structure for new event list
  bitcols=20

  outstruc=replicate(outstruc,n_elements(evt1)) ; replicate to match size of original event list
   
  if keyword_set(ciaopath) then ciaodir=ciaopath else ciaodir='/soft/ciao/'

; assign ctifile name
  trapmap=ciaodir+'/CALDB/data/chandra/acis/cti/'+sxpar(hevt,'CTIFILE')
       
 
; read in trapmap table and extensions of CTI file for chips in observation:
  trpmp_table=mrdfits(trapmap,1,ht)
  
  if where(ccds eq '0') ne -1 then begin
     trpmp_img0=mrdfits(trapmap,2,h0)
     bscale=sxpar(h0,'BSCALE')
     bzero=sxpar(h0,'BZERO')
     trpmp_img0=bzero+bscale*trpmp_img0
  endif
  

  if where(ccds eq '1') ne -1 then begin
     trpmp_img1=mrdfits(trapmap,3,h1)
     bscale=sxpar(h1,'BSCALE')
     bzero=sxpar(h1,'BZERO')
     trpmp_img1=bzero+bscale*trpmp_img1
  endif
  
  if where(ccds eq '2') ne -1 then begin
     trpmp_img2=mrdfits(trapmap,4,h2)
     bscale=sxpar(h2,'BSCALE')
     bzero=sxpar(h2,'BZERO')
     trpmp_img2=bzero+bscale*trpmp_img2
  endif
  
  if where(ccds eq '3') ne -1 then begin
     trpmp_img3=mrdfits(trapmap,5,h3)
     bscale=sxpar(h3,'BSCALE')
     bzero=sxpar(h3,'BZERO')
     trpmp_img3=bzero+bscale*trpmp_img3
  endif
  
  if where(ccds eq '4') ne -1 then begin
     trpmp_img4=mrdfits(trapmap,6,h4)
     bscale=sxpar(h4,'BSCALE')
     bzero=sxpar(h4,'BZERO')
     trpmp_img4=bzero+bscale*trpmp_img4
  endif

  if where(ccds eq '5') ne -1 then begin
     ; parallel
     trpmp_img5=mrdfits(trapmap,8,h5)
     bscale=sxpar(h5,'BSCALE')
     bzero=sxpar(h5,'BZERO')
     trpmp_img5=bzero+bscale*trpmp_img5
     
     ;serial
     trpmp_ser5=mrdfits(trapmap,7,h5s)
     bscale=sxpar(h5s,'BSCALE')
     bzero=sxpar(h5s,'BZERO')
     trpmp_ser5=bzero+bscale*trpmp_ser5
  endif
 
  if where(ccds eq '6') ne -1 then begin
     trpmp_img6=mrdfits(trapmap,9,h6)
     bscale=sxpar(h6,'BSCALE')
     bzero=sxpar(h6,'BZERO')
     trpmp_img6=bzero+bscale*trpmp_img6
  endif
  
  if where(ccds eq '7') ne -1 then begin
     ; parallel
     trpmp_img7=mrdfits(trapmap,11,h7)
     bscale=sxpar(h7,'BSCALE')
     bzero=sxpar(h7,'BZERO')
     trpmp_img7=bzero+bscale*trpmp_img7
     
     ;serial
     trpmp_ser7=mrdfits(trapmap,10,h7s)
     bscale=sxpar(h7s,'BSCALE')
     bzero=sxpar(h7s,'BZERO')
     trpmp_ser7=bzero+bscale*trpmp_ser7
  endif

  if where(ccds eq '8') ne -1 then begin
     trpmp_img8=mrdfits(trapmap,12,h8)
     bscale=sxpar(h8,'BSCALE')
     bzero=sxpar(h8,'BZERO')
     trpmp_img8=bzero+bscale*trpmp_img8
  endif
  
  if where(ccds eq '9') ne -1 then begin
     trpmp_img9=mrdfits(trapmap,13,h9)
     bscale=sxpar(h9,'BSCALE')
     bzero=sxpar(h9,'BZERO')
     trpmp_img9=bzero+bscale*trpmp_img9
  endif

  
  trlfrac0=sxpar(ht,'FRCTRLY0')
  trlfrac1=sxpar(ht,'FRCTRLY1')
  trlfrac2=sxpar(ht,'FRCTRLY2')
  trlfrac3=sxpar(ht,'FRCTRLY3')
  trlfrac4=sxpar(ht,'FRCTRLY4')
  trlfrac6=sxpar(ht,'FRCTRLY6')
  trlfrac8=sxpar(ht,'FRCTRLY8')
  trlfrac9=sxpar(ht,'FRCTRLY9')

  trlfrac5=sxpar(ht,'FRCTRLY5')
  trlfrac7=sxpar(ht,'FRCTRLY7')
  trlfrac5x=sxpar(ht,'FRCTRLX5')
  trlfrac7x=sxpar(ht,'FRCTRLX7')

  ; read in parallel and serial CTI change per degree numbers
  if keyword_set(coeff_file) then coefffile=coeff_file else coefffile='tcticorr_coeff_v1.txt'
  readcol,coefffile,blah,tctiarr,tctiarr_ser,/silent

  ; for testing versus ape:
  if keyword_set(zero) then begin
     tctiarr=fltarr(10) 
     tctiarr_ser=fltarr(10)
  endif

  if keyword_set(verbose) then begin
     print, 'Observation date: '+obsdate
     print
     print, 'Using ctifile '+trapmap
     print
     print,"Using % CTI change per degree of:"
     print,"Parallel --"
     for i=0,9 do print,'CCD_ID='+strcompress(string(i),/remove_all)+', '+strcompress(string(100.*tctiarr[i],format='(F5.2)'),/remove_all)+'%'
     print
     print,"Serial --"
     for i=0,9 do print,'CCD_ID='+strcompress(string(i),/remove_all)+', '+strcompress(string(100.*tctiarr_ser[i],format='(F5.2)'),/remove_all)+'%'
     print
     print,"Beginning CTI correction"
     print
  endif

  split_threshold=fltarr(10)+13

  cti_converge=0.1
  if keyword_set(niter) then max_n_iter=niter else max_n_iter=15
    
  phas_adj=1.0*phas[*,*,*]
  
  for i=0,nccd-1 do begin

     if keyword_set(verbose) then begin
        print,'Working on CCD '+strcompress(string(fix(ccds[i])),/remove_all)
        print,'Parallel: '+strcompress(string(100.*tctiarr[ccds[i]],format='(F5.2)'),/remove_all)+'%'
        print,'Serial: '+strcompress(string(100.*tctiarr_ser[ccds[i]],format='(F5.2)'),/remove_all)+'%'
     endif

     case ccds[i] of
        0: begin
           trapdens=trpmp_img0
           trlfrac=trlfrac0
        end
        1: begin
           trapdens=trpmp_img1
           trlfrac=trlfrac1
        end
        2: begin
           trapdens=trpmp_img2
           trlfrac=trlfrac2
        end
        3: begin
           trapdens=trpmp_img3
           trlfrac=trlfrac3
        end
        4: begin
           trapdens=trpmp_img4
           trlfrac=trlfrac4
        end
        6: begin
           trapdens=trpmp_img6
           trlfrac=trlfrac6
        end
        8: begin
           trapdens=trpmp_img8
           trlfrac=trlfrac8
        end
        9: begin
           trapdens=trpmp_img9
           trlfrac=trlfrac9
        end
        5: begin
           trapdens=trpmp_img5
           trlfrac=trlfrac5
           trapdensSER=trpmp_ser5
           trlfracX=trlfrac5x
        end
        7: begin
           trapdens=trpmp_img7
           trlfrac=trlfrac7
           trapdensSER=trpmp_ser7
           trlfracX=trlfrac7x
        end
           
     endcase
                
     w1 = where(trpmp_table.ccd_id eq ccds[i])
	
     for j=0,n_elements(w1)-1 do begin
        
        i_row = j
        
        w2 = where((ccd_id eq ccds[i]) and (chx ge ((trpmp_table.chipx_lo)[w1[i_row]] - 0.5)) and $
                   (chx lt ((trpmp_table.chipx_hi)[w1[i_row]] + 0.5)) and $
                   (chy ge ((trpmp_table.chipy_lo)[w1[i_row]] - 0.5)) and $
                   (chy lt ((trpmp_table.chipy_hi)[w1[i_row]] + 0.5)))
        
        if n_elements(w2) gt 0 then begin
           
;         Get pha, volume_x, volume_y
           
           tmp1 = (trpmp_table.npoints)[w1[i_row]]
           pha = (trpmp_table.pha)[0:(tmp1-1),w1[i_row]]
           volume_x = (trpmp_table.volume_x)[0:(tmp1-1),w1[i_row]]
           volume_y = (trpmp_table.volume_y)[0:(tmp1-1),w1[i_row]]
           
           if pha[0] gt split_threshold[ccds[i]] then begin
              
              tmp1 = split_threshold - pha[0]
              tmp2 = (volume_x[1] - volume_x[0]) / (pha[1] - pha[0]) * $
                     tmp1 + volume_x[0]
              volume_x = [tmp2, volume_x]
              
              tmp2 = (volume_y[1] - volume_y[0]) / (pha[1] - pha[0]) * $
                     tmp1 + volume_y[0]
              volume_y = [tmp2, volume_y]
              
              pha = [split_threshold[ccds[i]], pha]
           endif
           
           if pha[n_elements(pha)-1] lt 4096 then begin
              
              tmp1 = 4096 - pha[n_elements(pha)-1]
              tmp2 = $
                 (volume_x[n_elements(pha)-1] - volume_x[n_elements(pha)-2]) / $ 
                 (pha[n_elements(pha)-1] - pha[n_elements(pha)-2]) * $
                 tmp1 + volume_x[n_elements(pha)-2]
              volume_x = [volume_x, tmp2]
              
              tmp2 = $
                 (volume_y[n_elements(pha)-1] - volume_y[n_elements(pha)-2]) / $ 
                 (pha[n_elements(pha)-1] - pha[n_elements(pha)-2]) * $
                 tmp1 + volume_y[n_elements(pha)-2]
              volume_y = [volume_y, tmp2]
              
              pha = [pha, 4096]
           endif
           
;         Loop over events
		  
           for k=0L, n_elements(w2)-1L do begin
		    
              i_evt = k
              
              diff_y=0.0*phas_adj[*,*,w2[i_evt]]
              diff_x=0.0*phas_adj[*,*,w2[i_evt]]
               
              fpt=interpol(gea_fp_temp,gea_time,event_time[w2[i_evt]])
              fpt_cel = fpt-273.15  

              ;  Loop over # iterations
              
              max_n_iter_flag = 1
              n_iter=0
              
              while max_n_iter_flag eq 1 do begin
                 
                 diff_x_prev = 1.0 * diff_x
                 diff_y_prev = 1.0 * diff_y
                 phas_tmp = 1.0 * phas_adj[*,*,w2[i_evt]]
                 
                 n_iter += 1
			  

;                begin Serial correction
                 
                 if ((ccds[i] eq 5) or (ccds[i] eq 7)) then begin
	  
;                diff_x, leading column
                             
                    for jj=0,2 do begin
                       
                       i1 = jj
                       
                       if ((chx[w2[i_evt]] gt 1) and (chx[w2[i_evt]] le 256)) $
                       then i3 = 0 
                       if ((chx[w2[i_evt]] ge 257) and (chx[w2[i_evt]] le 512)) $
                       then i3 = 2
                       if ((chx[w2[i_evt]] ge 513) and (chx[w2[i_evt]] le 768)) $
                       then i3 = 0
                       if ((chx[w2[i_evt]] ge 769) and (chx[w2[i_evt]] lt 1024)) $
                       then i3 = 2
                       
                       ix = i3 + i_x_off
                       iy = i1 + i_y_off 
                       tmp1 = phas[ix,iy,w2[i_evt]] + $
                              diff_x_prev[ix,iy] + $
                              diff_y_prev[ix,iy]
                       
                       if (tmp1 ge split_threshold[ccds[i]]) then begin
                          yint=interpol(volume_x,pha,tmp1)
                          diff_x[ix,iy] = $
                             trapdensSER[ $
                             chx[w2[i_evt]]+i3-2, $
                             chy[w2[i_evt]]+i1-2] * $
                             yint[0]*(1.0+tctiarr_ser[ccds[i]]*(fpt_cel+119.7))
                       endif
                       
                    endfor
                    
                    
;                diff_x, trailing two columns
                    
                    for jj=0,2 do begin
                       
                       i1 = jj
                       
                       for ii=1,2 do begin
                          
                          i2 = ii
                          
                          if ((chx[w2[i_evt]] gt 1) and $
                              (chx[w2[i_evt]] le 256)) then $
                                 i3 = 1
                          if ((chx[w2[i_evt]] ge 257) and $
                              (chx[w2[i_evt]] le 512)) then $
                                 i3 = -1
                          if ((chx[w2[i_evt]] ge 513) and $
                              (chx[w2[i_evt]] le 768)) then $
                                 i3 = 1
                          if ((chx[w2[i_evt]] ge 769) and $
                              (chx[w2[i_evt]] lt 1024)) then $
                                 i3 = -1
                          ix = 1 + i3 * (i2 - 1) + i_x_off
                          iy = i1 + i_y_off
                          tmp1 = phas[ix,iy,w2[i_evt]] + $
                                 diff_x_prev[ix,iy] + $
                                 diff_y_prev[ix,iy]
                          tmp2 = phas[ix-i3,iy,w2[i_evt]] + $
                                 diff_x_prev[ix-i3,iy] + $
                                 diff_y_prev[ix-i3,iy]
                          
                          if (((tmp1 ge split_threshold[ccds[i]]) and $
                               (split_threshold[ccds[i]] gt tmp2)) or $
                              ((chx[w2[i_evt]] eq 256) and $
                               (i2 eq 2) and $
                               (tmp1 ge split_threshold[ccds[i]])) or $
                              ((chx[w2[i_evt]] eq 257) and $
                               (i2 eq 2) and $
                               (tmp1 ge split_threshold[ccds[i]])) or $
                              ((chx[w2[i_evt]] eq 512) and $
                               (i2 eq 1) and $
                               (tmp1 ge split_threshold[ccds[i]])) or $
                              ((chx[w2[i_evt]] eq 513) and $
                               (i2 eq 1) and $
                               (tmp1 ge split_threshold[ccds[i]])) or $
                              ((chx[w2[i_evt]] eq 768) and $
                               (i2 eq 2) and $
                               (tmp1 ge split_threshold[ccds[i]])) or $
                              ((chx[w2[i_evt]] eq 769) and $
                               (i2 eq 2) and $
                               (tmp1 ge split_threshold[ccds[i]]))) then begin 
                             yint=interpol(volume_x,pha,tmp1)
                             diff_x[ix,iy] = $
                                trapdensSER[ $
                                chx[w2[i_evt]]-1+i3*(i2-1), $
                                chy[w2[i_evt]]+i1-2 ] * $
                                yint[0]*(1.0+tctiarr_ser[ccds[i]]*(fpt_cel+119.7))
                          endif
                          
                          if (((tmp1 ge tmp2) and $ 
                               (tmp2 ge split_threshold[ccds[i]])) and $
                              not ((chx[w2[i_evt]] eq 256) and $
                                   (i2 eq 2)) and $
                              not ((chx[w2[i_evt]] eq 257) and $
                                   (i2 eq 2)) and $
                              not ((chx[w2[i_evt]] eq 512) and $
                                   (i2 eq 1)) and $
                              not ((chx[w2[i_evt]] eq 513) and $
                                   (i2 eq 1)) and $
                              not ((chx[w2[i_evt]]) eq 768 and $
                                   (i2 eq 2)) and $
                              not ((chx[w2[i_evt]] eq 769) and $
                                   (i2 eq 2))) then begin
                             yint1=interpol(volume_x,pha,tmp1)
                             yint2=interpol(volume_x,pha,tmp2)
                             diff_x[ix,iy] = $
                                (trapdensSER[ $
                                chx[w2[i_evt]]-1+i3*(i2-1), $
                                chy[w2[i_evt]]+i1-2 ] * $
                                 yint1[0]*(1.0+tctiarr_ser[ccds[i]]*(fpt_cel+119.7)) - $
                                 trapdensSER[ $
                                chx[w2[i_evt]]-1+i3*(i2-2), $
                                chy[w2[i_evt]]+i1-2 ] * $
                                 yint2[0]*(1.0+tctiarr_ser[ccds[i]]*(fpt_cel+119.7)))
                          endif
                          
                          if (((tmp2 gt tmp1) and $
                               (tmp1 ge split_threshold[ccds[i]])) and $
                              not ((chx[w2[i_evt]] eq 256) and $
                                   (i2 eq 2)) and $
                              not ((chx[w2[i_evt]] eq 257) and $
                                   (i2 eq 2)) and $
                              not ((chx[w2[i_evt]] eq 512) and $
                                   (i2 eq 1)) and $
                              not ((chx[w2[i_evt]] eq 513) and $
                                   (i2 eq 1)) and $
                              not ((chx[w2[i_evt]] eq 768) and $
                                   (i2 eq 2)) and $
                              not ((chx[w2[i_evt]] eq 769) and $
                                   (i2 eq 2))) then begin
                             yint1=interpol(volume_x,pha,tmp1)
                             yint2=interpol(volume_x,pha,tmp2)
                             diff_x[ix,iy] = $
                                trlfracX * $ 
                                (trapdensSER[ $
                                chx[w2[i_evt]]-1+i3*(i2-1), $
                                chy[w2[i_evt]]+i1-2 ] * $ 
                                 yint1[0]*(1.0+tctiarr_ser[ccds[i]]*(fpt_cel+119.7)) - $
                                 trapdensSER[ $
                                chx[w2[i_evt]]-1+i3*(i2-2), $
                                chy[w2[i_evt]]+i1-2 ] * $
                                 yint2[0]*(1.0+tctiarr_ser[ccds[i]]*(fpt_cel+119.7)))
                          endif
                          
                       endfor
                       
                    endfor
                    
                 endif

;                end Serial correction

;                begin Parallel correction
			  
;                 diff_y, first row
                 
                 for ii=0,2 do begin
                    
                    i1 = ii
                    
                    ix = i1 + i_x_off
                    iy = 0 + i_y_off
                    tmp1 = phas[ix,iy,w2[i_evt]] + $
                           diff_x[ix,iy] + $
                           diff_y_prev[ix,iy]
                    
                    yint=interpol(volume_y,pha,tmp1)
                    
                    if (tmp1 ge split_threshold[ccds[i]]) then $
                       diff_y[ix,iy] = $
                       trapdens[ $
                       chx[w2[i_evt]]+i1-2, $
                       chy[w2[i_evt]]-2] * $
                       yint[0]*(1.0+tctiarr[ccds[i]]*(fpt_cel+119.7))
                 endfor
                 
;                 diff_y, middle and top rows
                 
                 for ii=0,2 do begin
                    
                    i1 = ii
                    
                    for jj=1,2 do begin
                       
                       i2 = jj
                       
                       ix = i1 + i_x_off
                       iy = i2 + i_y_off
                       
                       tmp1 = phas[ix,iy,w2[i_evt]] + $
                              diff_x[ix,iy] + $
                              diff_y_prev[ix,iy]
                       tmp2 = phas[ix,iy-1,w2[i_evt]] + $ 
                              diff_x[ix,iy-1] + $
                              diff_y_prev[ix,iy-1]
                       
                       if ((tmp1 ge split_threshold[ccds[i]]) and $
                           (split_threshold[ccds[i]] > tmp2)) then begin
                          yint=interpol(volume_y, pha, tmp1)
                          diff_y[ix,iy] = $
                             trapdens[ $
                             chx[w2[i_evt]]+i1-2, $
                             chy[w2[i_evt]]+i2-2 ] * $ 
                             yint[0]*(1.0+tctiarr[ccds[i]]*(fpt_cel+119.7))
                       endif
                       
                       if ((tmp1 ge tmp2) and $ 
                           (tmp2 ge split_threshold[ccds[i]])) then begin
                          yint1=interpol(volume_y, pha, tmp1)
                          yint2=interpol(volume_y, pha, tmp2)
                          diff_y[ix,iy] = $
                             (trapdens[ $
                             chx[w2[i_evt]]+i1-2, $
                             chy[w2[i_evt]]+i2-2 ] * $ 
                              yint1[0]*(1.0+tctiarr[ccds[i]]*(fpt_cel+119.7)) - $
                              trapdens[ $
                             chx[w2[i_evt]]+i1-2, $
                             chy[w2[i_evt]]+i2-3 ] * $ 
                              yint2[0]*(1.0+tctiarr[ccds[i]]*(fpt_cel+119.7)))
                       endif

                       if ((tmp2 gt tmp1) and $ 
                           (tmp1 ge split_threshold[ccds[i]])) then begin
                          yint1=interpol(volume_y, pha,tmp1)
                          yint2=interpol(volume_y, pha,tmp2)
                          diff_y[ix,iy] = $
                             trlfrac * $ 
                             (trapdens[ $
                             chx[w2[i_evt]]+i1-2, $
                             chy[w2[i_evt]]+i2-2 ] * $
                              yint1[0]*(1.0+tctiarr[ccds[i]]*(fpt_cel+119.7)) - $
                              trapdens[ $
                             chx[w2[i_evt]]+i1-2, $
                             chy[w2[i_evt]]+i2-3 ] * $ 
                              yint2[0]*(1.0+tctiarr[ccds[i]]*(fpt_cel+119.7)))
                          
                       endif
                                
                    endfor
                    
                 endfor
                          
;         end parallel adjustment
                 
;         update phas_adj

                 phas_adj[*,*,w2[i_evt]] = phas[*,*,w2[i_evt]] + $ 
                                           diff_x[*,*] + diff_y[*,*]
                 
;                 Check convergence criterion

                 if (max(abs(phas_adj[*,*,w2[i_evt]] - $ 
                             phas_tmp[*,*])) lt cti_converge) then $ 
                                max_n_iter_flag = 0
                 
                 if (n_iter ge max_n_iter) then $
                    max_n_iter_flag = 0
                 
              endwhile ; end correction for given event
                      
             ; take floor of PHAS_ADJ
             if keyword_set(float) then phas_adj[*,*,w2[i_evt]]=phas_adj[*,*,w2[i_evt]] else phas_adj[*,*,w2[i_evt]]=floor(phas_adj[*,*,w2[i_evt]])

              ; if PHAS > 4095 set to 4095
              for iii=0,nmd-1 do begin
                 for jjj=0,nmd-1 do if phas_adj[iii,jjj,w2[i_evt]] gt phas_max then phas_adj[iii,jjj,w2[i_evt]]=phas_max
              endfor
              
              ; write out event
              outstruc[w2[i_evt]].phas=phas_adj[*,*,w2[i_evt]]
              outstruc[w2[i_evt]].time=evt1[w2[i_evt]].time
              outstruc[w2[i_evt]].ccd_id=evt1[w2[i_evt]].ccd_id
              outstruc[w2[i_evt]].node_id=evt1[w2[i_evt]].node_id
              outstruc[w2[i_evt]].expno=evt1[w2[i_evt]].expno
              outstruc[w2[i_evt]].chipx=chx[w2[i_evt]]
              outstruc[w2[i_evt]].chipy=chy[w2[i_evt]]
              outstruc[w2[i_evt]].tdetx=evt1[w2[i_evt]].tdetx
              outstruc[w2[i_evt]].tdety=evt1[w2[i_evt]].tdety
              outstruc[w2[i_evt]].detx=evt1[w2[i_evt]].detx
              outstruc[w2[i_evt]].dety=evt1[w2[i_evt]].dety
              outstruc[w2[i_evt]].x=evt1[w2[i_evt]].x
              outstruc[w2[i_evt]].y=evt1[w2[i_evt]].y
              outstruc[w2[i_evt]].pha=evt1[w2[i_evt]].pha
              outstruc[w2[i_evt]].pha_ro=evt1[w2[i_evt]].pha_ro
              outstruc[w2[i_evt]].energy=evt1[w2[i_evt]].energy
              outstruc[w2[i_evt]].pi=evt1[w2[i_evt]].pi
              outstruc[w2[i_evt]].fltgrade=evt1[w2[i_evt]].fltgrade
              outstruc[w2[i_evt]].grade=evt1[w2[i_evt]].grade
              outstruc[w2[i_evt]].status=evt1[w2[i_evt]].status

           endfor ; end loop over events
                  
        endif ; end loop over given ccd
              
     endfor  ; end loop over given ccd
     
  endfor ; end loop over ccds

if keyword_set(verbose) then begin
   print
   print,"CTI correction done; writing out event list"
endif

mwrfits,outstruc,evtout,hevt,/create,nbit_cols=32,bit_cols=bitcols ;write out new event list EVENT block

; read in and write out GTI block(s)
for i=0,nccd-1 do begin
   gti=mrdfits(evtfile,i+2,hgti)
   mwrfits,gti,evtout,hgti     
endfor


if keyword_set(verbose) then begin
   print
   print,"spawning acis_process_events"
   print
   print,ciaodir+'/bin/punlearn acis_process_events'
   if keyword_set(zero) then print,ciaodir+'/bin/acis_process_events infile="'+evtout+'[cols -status]" outfile='+outname+'_tcticorr_intermediate2.fits acaofffile='+asol+' badpixfile='+bpixfile+' apply_cti=yes ctifile=acisD2002-08-01ctiNnull.fits apply_tgain=yes tgainfile=CALDB doevtgrade=yes calculate_pi=yes rand_pha=no' else print,ciaodir+'/bin/acis_process_events infile="'+evtout+'[cols -status]" outfile='+outname+'_tcticorr_intermediate2.fits acaofffile='+asol+' badpixfile='+bpixfile+' apply_cti=yes ctifile=acisD2002-08-01ctiNnull.fits apply_tgain=yes tgainfile=CALDB doevtgrade=yes calculate_pi=yes'
   print
   print
endif

spawn,ciaodir+'/bin/punlearn acis_process_events'
if keyword_set(zero) then spawn,ciaodir+'/bin/acis_process_events infile="'+evtout+'[cols -status]" outfile='+outname+'_tcticorr_intermediate2.fits acaofffile='+asol+' badpixfile='+bpixfile+' apply_cti=yes ctifile=acisD2002-08-01ctiNnull.fits apply_tgain=yes tgainfile=CALDB doevtgrade=yes calculate_pi=yes rand_pha=no' else spawn,ciaodir+'/bin/acis_process_events infile="'+evtout+'[cols -status]" outfile='+outname+'_tcticorr_intermediate2.fits acaofffile='+asol+' badpixfile='+bpixfile+' apply_cti=yes ctifile=acisD2002-08-01ctiNnull.fits apply_tgain=yes tgainfile=CALDB doevtgrade=yes calculate_pi=yes'
 

; insert dmpaste to fix status column here!
if keyword_set(verbose) then begin
   print,"fixing status column"
   print
   print,ciaodir+'/bin/punlearn dmpaste'
   print,ciaodir+'/bin/dmpaste "'+outname+'_tcticorr_intermediate2.fits[EVENTS][cols -status]" "'+evtfile+'[events][cols status]" '+outname+'_tcticorr_evt1.fits'
   print
   print
endif

spawn,ciaodir+'/bin/punlearn dmpaste'
spawn,ciaodir+'/bin/dmpaste "'+outname+'_tcticorr_intermediate2.fits[EVENTS][cols -status]" "'+evtfile+'[events][cols status]" '+outname+'_tcticorr_evt1.fits'


if keyword_set(verbose) then begin
   print,"adding TCTICORR keyword to event list header"
   print
   print,ciaodir+'/bin/dmhedit infile='+outname+'_tcticorr_evt1.fits filelist=none operation=add key="TCTICORR" value="'+coefffile+'" datatype=string'
   print
   print
endif

spawn,ciaodir+'/bin/dmhedit infile='+outname+'_tcticorr_evt1.fits filelist=none operation=add key="TCTICORR" value="'+coefffile+'" datatype=string'

;; filter to make level 2 event list
if keyword_set(filter) then begin
   if keyword_set(verbose) then begin
      print,"creating level 2 event list"
      print
      print,ciaodir+'/bin/punlearn dmcopy'
      print,ciaodir+'/bin/dmcopy "'+outname+'_tcticorr_evt1.fits[EVENTS][grade=0,2,3,4,6,status=0]" '+outname+'_tcticorr_evt1_flt.fits'
      print
      print,ciaodir+'/bin/dmcopy "'+outname+'_tcticorr_evt1_flt.fits[EVENTS][@'+fltfile+']" '+outname+'_tcticorr_evt2.fits'
      print
   endif

   spawn,ciaodir+'/bin/punlearn dmcopy'
   spawn,ciaodir+'/bin/dmcopy "'+outname+'_tcticorr_evt1.fits[EVENTS][grade=0,2,3,4,6,status=0]" '+outname+'_tcticorr_evt1_flt.fits'

   spawn,ciaodir+'/bin/dmcopy "'+outname+'_tcticorr_evt1_flt.fits[EVENTS][@'+fltfile+']" '+outname+'_tcticorr_evt2.fits'

endif

return
end