openw,99,'image5.dat'
reads,'x0bk0602t.d0h',a0,nx,ny
reads,'x0bk0603t.d0h',a1,nx,ny
tvwindow,0,xs=nx,ys=ny
tv,a0*20.
print,' Input center coordinates'
tvrdc,x,y,/dev
print,'center=',x,y
printf,99,'center=',x,y
b0=a0(x-127:x+128,y-127:y+128)*1.
b1=a1(x-127:x+128,y-127:y+128)*1.
back=fltarr(4,2)
a0s=smooth(a0*1.,5)
a1s=smooth(a1*1.,5)
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(a0s(xb-4:xb+5,yb-4:yb+5))
   back(i,1)=median(a1s(xb-4:xb+5,yb-4:yb+5))
   dist=sqrt((x-xb)^2+(y-yb)^2)*0.043/2.
   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.)
b1=(b1-total(back(*,1))/4.)
nx=n_elements(b0(*,0)) 
ny=n_elements(b0(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/2.
xoff=total(x*b0)/total(b0)-total(x*b1)/total(b1)
yoff=total(y*b0)/total(b0)-total(y*b1)/total(b1)
printf,99,'xoff,yoff',xoff,yoff
print,'xoff,yoff',xoff,yoff
if(xoff lt 0.)then xoff=xoff-0.5
if(xoff gt 0.)then xoff=xoff+0.5
if(yoff lt 0.)then yoff=yoff-0.5
if(yoff gt 0.)then yoff=yoff+0.5
;b1=shift(b1,xoff,yoff)
xoff=total(x*b0)/total(b0)-total(x*b1)/total(b1)
yoff=total(y*b0)/total(b0)-total(y*b1)/total(b1)
printf,99,'not new xoff,yoff',xoff,yoff
print,'not new xoff,yoff',xoff,yoff
;
b=b0+b1
tv,congrid(b,256,256),0,0
;
b=b/total(b)
writes,'image5',b
close,99
stop
end