Pro rad_hits,Log,filt

; 19Sep02 by C. See
; The program counts the energetic particle hits in the images of "Log".
; It does this by creating a synthetic 'uneffected' image set from the first
; three dark frames, and comparing that to the rest of the dark images in the
; log.  Variations greater than 5 dn are considered a hit.

zero_ms_max = 40		;memory zone darkcurrent (i.e. 0 ms), column 255 (max)
image_zone_dc = 0.1  ;dn/ms	;image zone darkcurrent rate (at F10)
printflag=0
filecount=0	;counter for image files processed

cd,'',current=cdir & cd,cdir	    ;establishes cdir as current directory
if n_params() eq 0 then log = cdir  ;test for # of input parameters
if n_params() eq 0 then filt = 1  		;scale for temperature
filt=0 ;hardwire filter off

cd,Log+'/DB/Image'
spawn,'lsf',files,/noshell & siz=size(files) & nfiles=siz(1)	;for lots of files

;___________________________________
; read in the first 3 dark images...
j=0	;file index counter
mrifile=intarr(3)
slifile=intarr(3)
hrifile=intarr(3)

for i=0,2 do begin	;the first 3 dark images
mricount=''
slicount=''
hricount=''

 while not (mricount eq 'done' and slicount eq 'done' and hricount eq 'done') do begin
  if lamps(files(j)) eq '0000' then begin
   if mtype(files(j)) eq 21 and mricount eq '' then begin & mrifile(i)=j & mricount='done' & endif
   if mtype(files(j)) eq 22 and slicount eq '' then begin & slifile(i)=j & slicount='done' & endif
   if mtype(files(j)) eq 23 and hricount eq '' then begin & hrifile(i)=j & hricount='done' & endif
  endif
  j=j+1
 endwhile

endfor 

;___
MRI:

;____________________________________
;Read in reference images...

 d_read,files(mrifile(0)),h,p & mri1=p & temp0=d_value(h,22)
 d_read,files(mrifile(1)),h,p & mri2=p
 d_read,files(mrifile(2)),h,p & mri3=p
 

;remove average per row...
for i=0,255 do mri1(*,i)=mri1(*,i)-mean(mri1(*,i))
for i=0,255 do mri2(*,i)=mri2(*,i)-mean(mri2(*,i))
for i=0,255 do mri3(*,i)=mri3(*,i)-mean(mri3(*,i))

if filt eq 1 then begin
 ;scale to row average to remove hot pixel effects...
 row=fltarr(256)
 for m=0,255 do row(m)=mean(abs(mri1(*,m)))
 for m=0,255 do mri1(*,m)=mri1(*,m)/(row(m)/row(0))
 for m=0,255 do row(m)=mean(abs(mri2(*,m)))
 for m=0,255 do mri2(*,m)=mri2(*,m)/(row(m)/row(0))
 for m=0,255 do row(m)=mean(abs(mri3(*,m)))
 for m=0,255 do mri3(*,m)=mri3(*,m)/(row(m)/row(0))
endif

;;;scale to row average to remove hot pixel effects...
;;row=fltarr(256)
;;for m=0,255 do row(m)=mean(abs(mri1(*,m)))
;;for m=0,255 do mri1(*,m)=mri1(*,m)/(row(m)/row(0))

;;window,0
;;array=abs(mri1-mri2)
;;bin=1
;;hist=histogram(array,binsize=bin)
;;s=size(hist)
;;x=findgen(s(1))*bin+min(array)
;;plot,x,hist,psym=10,xrange=[5,max(x)],xtitle="Hit size, DN",ytitle="No. of hits",title='MRI before filtering'

;;ans=''
;;read,'MRI Before Filtering...',ans

index=where(abs(mri1-mri2) gt 1)
s=size(index)
for i=0l,s(1)-1 do mri1(index(i))=median([mri1(index(i)),mri2(index(i)),mri3(index(i))])
for i=0l,s(1)-1 do mri2(index(i))=median([mri1(index(i)),mri2(index(i)),mri3(index(i))])

index=where(abs(mri1-mri3) gt 1)
s=size(index)
for i=0l,s(1)-1 do mri1(index(i))=median([mri1(index(i)),mri2(index(i)),mri3(index(i))])
for i=0l,s(1)-1 do mri3(index(i))=median([mri1(index(i)),mri2(index(i)),mri3(index(i))])

index=where(abs(mri2-mri3) gt 1)
s=size(index)
for i=0l,s(1)-1 do mri2(index(i))=median([mri1(index(i)),mri2(index(i)),mri3(index(i))])
for i=0l,s(1)-1 do mri3(index(i))=median([mri1(index(i)),mri2(index(i)),mri3(index(i))])

mri=(mri1+mri2+mri3)/3

;;array=abs(mri1-mri2)
;;bin=1
;;hist=histogram(array,binsize=bin)
;;s=size(hist)
;;x=findgen(s(1))*bin+min(array)
;;plot,x,hist,psym=10,xrange=[5,max(x)],xtitle="Hit size, DN",ytitle="No. of hits",title='MRI after filtering'

;;ans=''
;;read,'after Filtering...',ans

;___
SLI:

;____________________________________
;Read in reference images...

 d_read,files(slifile(0)),h,p & sli1=p
 d_read,files(slifile(1)),h,p & sli2=p
 d_read,files(slifile(2)),h,p & sli3=p
 

;remove average per row...
for i=0,255 do sli1(*,i)=sli1(*,i)-mean(sli1(*,i))
for i=0,255 do sli2(*,i)=sli2(*,i)-mean(sli2(*,i))
for i=0,255 do sli3(*,i)=sli3(*,i)-mean(sli3(*,i))

if filt eq 1 then begin
 ;scale to row average to remove hot pixel effects...
 row=fltarr(256)
 for m=0,255 do row(m)=mean(abs(sli1(*,m)))
 for m=0,255 do sli1(*,m)=sli1(*,m)/(row(m)/row(0))
 for m=0,255 do row(m)=mean(abs(sli2(*,m)))
 for m=0,255 do sli2(*,m)=sli2(*,m)/(row(m)/row(0))
 for m=0,255 do row(m)=mean(abs(sli3(*,m)))
 for m=0,255 do sli3(*,m)=sli3(*,m)/(row(m)/row(0))
endif

;;array=abs(sli1-sli2)
;;bin=1
;;hist=histogram(array,binsize=bin)
;;s=size(hist)
;;x=findgen(s(1))*bin+min(array)
;;plot,x,hist,psym=10,xrange=[5,max(x)],xtitle="Hit size, DN",ytitle="No. of hits",title='SLI before filtering'

;;ans=''
;;read,'SLI Before Filtering...',ans

index=where(abs(sli1-sli2) gt 1)
s=size(index)
for i=0,s(1)-1 do sli1(index(i))=median([sli1(index(i)),sli2(index(i)),sli3(index(i))])
for i=0,s(1)-1 do sli2(index(i))=median([sli1(index(i)),sli2(index(i)),sli3(index(i))])

index=where(abs(sli1-sli3) gt 1)
s=size(index)
for i=0,s(1)-1 do sli1(index(i))=median([sli1(index(i)),sli2(index(i)),sli3(index(i))])
for i=0,s(1)-1 do sli3(index(i))=median([sli1(index(i)),sli2(index(i)),sli3(index(i))])

index=where(abs(sli2-sli3) gt 1)
s=size(index)
for i=0,s(1)-1 do sli2(index(i))=median([sli1(index(i)),sli2(index(i)),sli3(index(i))])
for i=0,s(1)-1 do sli3(index(i))=median([sli1(index(i)),sli2(index(i)),sli3(index(i))])

sli=(sli1+sli2+sli3)/3

;;array=abs(sli1-sli2)
;;bin=1
;;hist=histogram(array,binsize=bin)
;;s=size(hist)
;;x=findgen(s(1))*bin+min(array)
;;plot,x,hist,psym=10,xrange=[5,max(x)],xtitle="Hit size, DN",ytitle="No. of hits",title='SLI after filtering'

;;ans=''
;;read,'after Filtering...',ans

;___
HRI:

;____________________________________
;Read in reference images...

 d_read,files(hrifile(0)),h,p & hri1=p
 d_read,files(hrifile(1)),h,p & hri2=p
 d_read,files(hrifile(2)),h,p & hri3=p
 

;remove average per row...
for i=0,255 do hri1(*,i)=hri1(*,i)-mean(hri1(*,i))
for i=0,255 do hri2(*,i)=hri2(*,i)-mean(hri2(*,i))
for i=0,255 do hri3(*,i)=hri3(*,i)-mean(hri3(*,i))

if filt eq 1 then begin
 ;scale to row average to remove hot pixel effects...
 row=fltarr(256)
 for m=0,255 do row(m)=mean(abs(hri1(*,m)))
 for m=0,255 do hri1(*,m)=hri1(*,m)/(row(m)/row(0))
 for m=0,255 do row(m)=mean(abs(hri2(*,m)))
 for m=0,255 do hri2(*,m)=hri2(*,m)/(row(m)/row(0))
 for m=0,255 do row(m)=mean(abs(hri3(*,m)))
 for m=0,255 do hri3(*,m)=hri3(*,m)/(row(m)/row(0))
endif

;;array=abs(hri1-hri2)
;;bin=1
;;hist=histogram(array,binsize=bin)
;;s=size(hist)
;;x=findgen(s(1))*bin+min(array)
;;plot,x,hist,psym=10,xrange=[5,max(x)],xtitle="Hit size, DN",ytitle="No. of hits",title='HRI before filtering'

;;ans=''
;;read,'HRI Before Filtering...',ans

index=where(abs(hri1-hri2) gt 1)
s=size(index)
for i=0,s(1)-1 do hri1(index(i))=median([hri1(index(i)),hri2(index(i)),hri3(index(i))])
for i=0,s(1)-1 do hri2(index(i))=median([hri1(index(i)),hri2(index(i)),hri3(index(i))])

index=where(abs(hri1-hri3) gt 1)
s=size(index)
for i=0,s(1)-1 do hri1(index(i))=median([hri1(index(i)),hri2(index(i)),hri3(index(i))])
for i=0,s(1)-1 do hri3(index(i))=median([hri1(index(i)),hri2(index(i)),hri3(index(i))])

index=where(abs(hri2-hri3) gt 1)
s=size(index)
for i=0,s(1)-1 do hri2(index(i))=median([hri1(index(i)),hri2(index(i)),hri3(index(i))])
for i=0,s(1)-1 do hri3(index(i))=median([hri1(index(i)),hri2(index(i)),hri3(index(i))])

hri=(hri1+hri2+hri3)/3

;;array=abs(hri1-hri2)
;;bin=1
;;hist=histogram(array,binsize=bin)
;;s=size(hist)
;;x=findgen(s(1))*bin+min(array)
;;plot,x,hist,psym=10,xrange=[5,max(x)],xtitle="Hit size, DN",ytitle="No. of hits",title='HRI after filtering'

;;ans=''
;;read,'after Filtering...',ans

print,"Please wait...",string(7b)

;stop

;______________________________________________
; create array of cosmic ray hit effects...

array=intarr(100000)
j=0

for n=0,nfiles-1 do begin
 if lamps(files(n)) ne '0000' then goto,nextfile
 d_read,files(n),h,p
 s=size(p)
 if s(2) ne 256 then goto,nextfile	;short MRI, or uncompressed
 for i=0,255 do p(*,i)=p(*,i)-mean(p(*,i))  ;remove average & row slope

if filt eq 1 then begin
 ;scale to row average to remove hot pixel effects...
 row=fltarr(256)
 for m=0,255 do row(m)=mean(abs(p(*,m)))
 for m=0,255 do p(*,m)=p(*,m)/(row(m)/row(0))
 for m=0,255 do row(m)=mean(abs(p(*,m)))
 for m=0,255 do p(*,m)=p(*,m)/(row(m)/row(0))
 for m=0,255 do row(m)=mean(abs(p(*,m)))
 for m=0,255 do p(*,m)=p(*,m)/(row(m)/row(0))
endif
 
  filecount=filecount+1	;count of actual files used
 
 
 if d_value(h,96) eq 21 then begin
  if max(abs(p-mri)) gt 500 then show3,p-mri
  index=where(abs(p-mri) gt 5)
  if index(0) eq -1 then goto,nextfile
  s=size(index)
  array(j:j+s(1)-1)=abs(p(index)-mri(index))
  j=j+s(1)
 endif
 if d_value(h,96) eq 22 then begin
  if max(abs(p-sli)) gt 500 then show3,p-sli
  index=where(abs(p-sli) gt 5)
  if index(0) eq -1 then goto,nextfile
  s=size(index)
  array(j:j+s(1)-1)=abs(p(index)-sli(index))
  j=j+s(1)
 endif
 if d_value(h,96) eq 23 then begin
  if max(abs(p-hri)) gt 500 then show3,p-hri
  index=where(abs(p-hri) gt 5)
  if index(0) eq -1 then goto,nextfile
  s=size(index)
  array(j:j+s(1)-1)=abs(p(index)-hri(index))
  j=j+s(1)
 endif
 nextfile:
endfor



;______________________________________________
;Plot Data...

If printflag eq 0 then window,0

Print,Log
Print,'The total number of hits greater than 5 dn was: ',sstr(j)
print,'The initial chip temperature is: ',temp0,string(10b)

s=size(array)
if index(0) eq -1 then s(*)=0	;i.e. no pixels over the limit
nhits=float(j)/filecount
print,'The number of hits greater that 5 dn is:  ',nhits,' per image'

index=where(array ge 10)
s=size(index)
if index(0) eq -1 then s(*)=0	;i.e. no pixels over the limit
perimage10=float(s(1))/filecount
print,'The number of hits greater that 10 dn is:  ',perimage10,' per image'

index=where(array ge 20)
s=size(index)
if index(0) eq -1 then s(*)=0	;i.e. no pixels over the limit
perimage20=float(s(1))/filecount
print,'The number of hits greater that 20 dn is:  ',perimage20,' per image'

index=where(array ge 25)
s=size(index)
if index(0) eq -1 then s(*)=0	;i.e. no pixels over the limit
perimage25=float(s(1))/filecount
print,'The number of hits greater that 25 dn is:  ',perimage25,' per image'

index=where(array ge 35)
s=size(index)
if index(0) eq -1 then s(*)=0	;i.e. no pixels over the limit
perimage35=float(s(1))/filecount
print,'The number of hits greater that 35 dn is:  ',perimage35,' per image'

index=where(array ge 50)
s=size(index)
if index(0) eq -1 then s(*)=0	;i.e. no pixels over the limit
perimage50=float(s(1))/filecount
print,'The number of hits greater that 50 dn is:  ',perimage50,' per image'

index=where(array ge 100)
s=size(index)
if index(0) eq -1 then s(*)=0	;i.e. no pixels over the limit
perimage100=float(s(1))/filecount
print,'The number of hits greater that 100 dn is:  ',perimage100,' per image'

close,2
openw,2,Log+"/rad_hits.out"
Printf,2,Log
printf,2,'The initial chip temperature is: ',temp0
Printf,2,'The total number of hits greater than 5 dn was: ',sstr(j),string(10b)

printf,2,'The number of hits greater that 5 dn is:  ',nhits,' per image'
printf,2,'The number of hits greater that 10 dn is:  ',perimage10,' per image'
printf,2,'The number of hits greater that 20 dn is:  ',perimage20,' per image'
printf,2,'The number of hits greater that 25 dn is:  ',perimage25,' per image'
printf,2,'The number of hits greater that 35 dn is:  ',perimage35,' per image'
printf,2,'The number of hits greater that 50 dn is:  ',perimage50,' per image'
printf,2,'The number of hits greater that 100 dn is:  ',perimage100,' per image'

close,2

;______________
;Plot array...
plot:
!p.multi=[0,1,1]

plot,array(0:j),xtitle='Hit Number, Cumulative (i.e. Mission Time)',$
 ytitle='Hit Size (dN)', yrange=[0,1.05*max(array)],$
 Title='Cosmic Ray hits from: '+Log
xyouts,0.5,0.93,alignment=.5,/normal,'   Max = '+sstr(max(array))+' dn,  # of hits/image > 100 dn = '+sstr(perimage100)
ans=''
If printflag eq 0 then read,'Hit return to continue...',ans

;______________
;Plot histograms...
If printflag eq 0 then window,1,xsize=500,ysize=800
!p.multi=[0,1,2]

bin=1
hist=histogram(array,binsize=bin)
s=size(hist)
x=findgen(s(1))*bin+min(array)
plot,x,hist,psym=10,xrange=[5,max(x)],yrange=[0,hist(where(x eq 6))],xtitle="Hit Amplitude, DN", $
 ytitle="Number of Occurances in "+(sstr(n))+" images",title='Histogram of Cosmic Ray Hits for'+Log
xyouts,0.5,0.93,alignment=.5,/normal,'Note: Hits less than 5 dN are ignored'

hist=float(hist)/filecount
plot,x,hist,psym=10,xrange=[6,15],xtitle="Hit Amplitude, DN",  $
 ytitle="Number of Occurances per Image",title='Close-up, 5 to 15 dn range, scaled per image'
oplot,x,hist
ans=''
If printflag eq 0 then read,'Hit return to continue...',ans

;______________
;Plot Distribution...
!p.multi=[0,1,1]

x=[10,20,25,35,50,100,max(array)]
y=[perimage10,perimage20,perimage25,perimage35,perimage50,perimage100,0]
plot,x,y,psym=2,xrange=[0,110],yrange=[0,1.10*max(y)],Title='Distribution of Rad Hits for'+Log,$
 xtitle='X = Hit Amplitude (dN)',ytitle='Total Number of Hits > X,  per Image'
x2=findgen(100)
y2=spline(x,y,x2,5)
oplot,x2,y2
ans=''
If printflag eq 0 then read,'Hit return to continue...',ans

;____________________________________
;Print Plots?...

ans=''
If printflag eq 0 then read,'Would you like to print the plots? (y or n): ',ans

if ans eq "y" then begin
 printflag=1
 set_plot,'ps'
 device,filename = '~/idl.ps',/palatino,/portrait,xsize=7.5,ysize=10, $
   xoffset=.25,yoffset=.5,/inches
 device,filename = '~/idl.ps',font_size=12,/inches,scale_factor=1.
 command='lp -ofp14 -olm5 '+Log+'/rad_hits.out'
 spawn,command
 goto,plot
endif

if printflag EQ 1 then begin
 device,/close_file
 spawn,'lp ~/idl.ps'
 set_plot,'x'		;returns to X windows environment
endif

;____________________________________
;Done...

;stop
!p.multi=[0,1,1]
cd,Log
end


;;;____________________________________
;;;make synthetic zero ms images...
;;zeromri=intarr(176,256)
;;zerosli=intarr(128,256)
;;zerohri=intarr(160,256)

;;column=zero_ms_max*findgen(256)/256	;memory zone darkcurrent for 0 ms exposure

;;for i=0,175 do zeromri(i,0:255)=column
;;for i=0,127 do zerosli(i,0:255)=column
;;for i=0,159 do zerohri(i,0:255)=column