openw,99,'image3.dat'
reads,'w0bm0604t.d0h',a0,nx,ny,1
reads,'w0bm0605t.d0h',a1,nx,ny,1
reads,'w0bm0606t.d0h',a2,nx,ny,1
tvwindow,0,xs=nx,ys=ny
tv,a1
print,' Input center coordinates'
tvrdc,x,y,/dev
print,'center=',x,y
printf,99,'center=',x,y
b0=a0(x-63:x+64,y-63:y+64)*1.
b1=a1(x-63:x+64,y-63:y+64)*1.
b2=a2(x-63:x+64,y-63:y+64)*1.
back=fltarr(4,3)
for i=0,3 do begin
   print,' Input background coordinates (cr to continue)',i
   cr=' ' & read,cr &   tvrdc,xb,yb,/dev
   back(i,0)=median(a0(xb-4:xb+5,yb-4:yb+5))
   back(i,1)=median(a1(xb-4:xb+5,yb-4:yb+5))
   back(i,2)=median(a2(xb-4:xb+5,yb-4:yb+5))
   dist=sqrt((x-xb)^2+(y-yb)^2)*0.043
   print,    'coordinates,background,dist (as) =',x,y,back(i,*),dist
   printf,99,'coordinates,background,dist (as) =',x,y,back(i,*),dist
endfor
print,'back=',back
b0=(b0-total(back(*,0))/4.)*100.
b1=(b1-total(back(*,1))/4.)*10.
b2=(b2-total(back(*,2))/4.)*1.
;
b=b2-b2
bk2=4095.-total(back(*,2))/4.-1000     
v2=where(b2 lt bk2)
b(v2)=b2(v2)
print,n_elements(v2)
tv,congrid(b,256,256),0,0
;
v0=where(b0/100. gt 25. and b0/100. lt 1000.)
b(v0)=b0(v0)
print,n_elements(v0)
tv,congrid(b,256,256),256,0
;
v1=where((b0/100. le 25. or b0/100. ge 1000.) and b2 ge bk2)
b(v1)=b1(v1)
print,n_elements(v1)
;
tv,congrid(b,256,256),0,256
print,' valid from each file',n_elements(v0),n_elements(v1),n_elements(v2)
printf,99,' valid from each file',n_elements(v0),n_elements(v1),n_elements(v2)
;
; Cosmic ray removal
;
keep=b
;mask=(b-b)-1
;mask(v2)=2
;mask(v1)=1
;mask(v0)=0
;print,' should be 0',min(mask),n_elements(v2 and v1), n_elements(v1 and v0), $
;n_elements(v0 and v2) 
nx=n_elements(b(*,0)) 
ny=n_elements(b(0,*)) 
x=(indgen(nx)-nx/2)#replicate(1.,ny) 
y=replicate(1.,nx) #(indgen(ny)-ny/2)
r=sqrt(x^2+y^2)*0.043
icut=10.
loopback:
   invalid=(r gt 1.) and (b ge max(b/icut))
   print,'n_elements(where(invalid ne 0.))',n_elements(Where(invalid ne 0.))
   bc=b
   bc(invalid)=0
   sm=smooth(bc,3)
   bc(invalid)=sm(invalid)
   tv,congrid(bc,256,256),256,256
   print,'loop? (0.=n)',icut  &printf,99,icut&read,icut 
if(icut ne 0.) then goto,loopback
;
; Remove some cosmics by hand
;
color_edit
cosmic1: 
  print, ' Mark bad pixels - right button to stop'
  tvscl,congrid(bc,256,256)
  tvrdc,xo,yo,/dev
if(!err eq 4) then goto, normal1
  xcos=xo/2 & ycos=yo/2
  if(xcos le nx-4 and ycos le ny-4 and xcos ge 3 and ycos ge 3) then begin
     box=bc(xcos-3:xcos+3,ycos-3:ycos+3)
     ixbox=where(box eq max(box)) & ixbox=ixbox(0)
     xcos=xcos+ixbox-7*(ixbox/7)-3 & ycos=ycos+ixbox/7-3
  endif
  bc(xcos,ycos)= (bc(xcos+1,ycos)+bc(xcos-1,ycos))/2.
  printf,99,' Manual cosmic at',xcos,ycos
goto, cosmic1
normal1:
bc=bc/total(bc)
writes,'image3',bc
close,99
stop
end