file='collisions.sav'
restore,file
;
close,1
openw,1,'lucky.dat'
;
;Define parameters
micronm=1.e-6
factor=2.^(1./3)
r_lucky=50.*micronm
;
@parameters
;---------------find the first lucky droplet--------------  
;luck_k=where(abs(ak-r_lucky) eq min(abs(ak-r_lucky)))
luck_k=where(ak ge r_lucky and ak le r_lucky*factor); Here the condition is loose so that we do not miss the case that droplets with similar size collide.
kgood=k(luck_k)
ilucky=kgood(0)
;--Identify the time that 'ilucky' reaches the size r_lucky
iak=ak(where(k eq ilucky))
it_lucky=t(where(k eq ilucky))
ik=where(abs(iak-r_lucky) eq min(abs(iak-r_lucky)))
ak_lucky=iak(ik(0))
t_lucky=it_lucky(ik(0))
print,'t_lucky=',t_lucky
print,'ak_lucky=',ak_lucky
;----cut off to the time that the first lucky particle reaches 50um
iit_lucky=where(t le t_lucky)
tt=t(iit_lucky)
jj=j(iit_lucky)
kk=k(iit_lucky)
akk=ak(iit_lucky)
ajj=aj(iit_lucky)
;tlast=max(t)
tlast=t_lucky
print,"confirm that ilucky is indeed in the array akk", akk(where(kk eq ilucky))
;----search the history of the lucky droplet backward------------------------

;main code
;
new_ilucky:

goodj=where(jj eq ilucky)
goodk=where(kk eq ilucky)
;
;  find out how many we have in each
t_merged=[tt(goodj),tt(goodk)]
j_merged=[jj(goodj),jj(goodk)]
k_merged=[kk(goodj),kk(goodk)]
aj_merged=[ajj(goodj),ajj(goodk)]
ak_merged=[akk(goodj),akk(goodk)]
sorted=sort(t_merged)
t_sorted=t_merged(sorted)
j_sorted=j_merged(sorted)
k_sorted=k_merged(sorted)
aj_sorted=aj_merged(sorted)
ak_sorted=ak_merged(sorted)
nt=n_elements(t_sorted)
;
; cut off to tlast
;
notzero=where(t_sorted le tlast)
t_sorted=t_sorted(notzero)
j_sorted=j_sorted(notzero)
k_sorted=k_sorted(notzero)
aj_sorted=aj_sorted(notzero)
ak_sorted=ak_sorted(notzero)
nt=n_elements(t_sorted)
;
;; track backwards
;fo="(f8.2,2i8,2f8.1)"
fo="(f9.4,2i9,2f8.2)"
it2=nt-1 & it1=0
for i=it2,it1,-1 do begin
  ajk_sorted_i=[aj_sorted[i],ak_sorted[i]]
  jk_sorted_i=[j_sorted[i],k_sorted[i]]
  good=where(ajk_sorted_i eq max(ajk_sorted_i))
  bad=where(ajk_sorted_i ne max(ajk_sorted_i))
  if n_elements(good) eq 2 then begin 
		 print,t_sorted[i],jk_sorted_i(0),jk_sorted_i(1),1e6*ajk_sorted_i(0),1e6*ajk_sorted_i(1),fo=fo
		 printf,1,t_sorted[i],jk_sorted_i(0),jk_sorted_i(1),1e6*ajk_sorted_i(0),1e6*ajk_sorted_i(1),fo=fo
	endif else begin
		if jk_sorted_i(good) ne ilucky then begin 
			ilucky=(jk_sorted_i(good))(0)
			tlast=t_sorted[i]
      goto,new_ilucky
			endif else begin 
			print,t_sorted[i],jk_sorted_i(good),jk_sorted_i(bad),1e6*ajk_sorted_i(good),1e6*ajk_sorted_i(bad),fo=fo
			printf,1,t_sorted[i],jk_sorted_i(good),jk_sorted_i(bad),1e6*ajk_sorted_i(good),1e6*ajk_sorted_i(bad),fo=fo
		 endelse
   endelse
endfor

close,1
;
END