xc=424 & yc=398
nx=161 & ny=161
tot=fltarr(nx*4,ny*2)
tvwindow,0,xsize=nx*4,ysize=ny*2
for image=1,7 do begin
   if(image eq 1)then reads,'X0B30201B.D0H',a,nx1,ny1
   if(image eq 2)then reads,'X0B30202B.D0h',a,nx2,ny2
   if(image eq 3)then reads,'X0B30203B.D0h',a,nx3,ny3
   if(image eq 4)then reads,'X0B30204B.D0h',a,nx4,ny4
   if(image eq 5)then reads,'X0B30301B.D0h',a,nx5,ny5
   if(image eq 6)then reads,'X0B30302B.D0h',a,nx6,ny6
   if(image eq 7)then reads,'X0B30303B.D0h',a,nx7,ny7
   peak=where(a eq max(a)) & yc=peak(0)/nx1 & xc=peak(0)-nx1*yc
   sky=total(1.*a(200:250,200:250))/(51*51.) & print,'sky=',sky
   xval=indgen(nx)#replicate(1.,ny)   & yval=replicate(1.,nx)#indgen(ny)
   psf=a(xc-nx/2:xc+nx/2,yc-nx/2:yc+nx/2)-sky*1.
   ox=total(psf*xval)/total(psf)-nx/2 
   oy=total(psf*yval)/total(psf)-ny/2
   print,'At',xc,yc,' Offsets',ox,oy 
   psf=a(xc+ox-nx/2:xc+ox+nx/2,yc+oy-nx/2:yc+oy+nx/2)-sky*1.
   offy=image/4 &offx=image-4*offy
   print,offx,offy
   for j=0,ny-1 do tot(offx*nx,offy*ny+j)=psf(*,j)
   for j=0,ny-1 do for i=0,nx-1 do tot(i,j)=tot(i,j)+psf(i,j)
   tvscl,tot
endfor
writes, 'focall',tot
stop
end