pro getpsf,int
common mpsf,nx,ny,np,psfm,wp,xp,yp,gp
common filen,fname,logfil
if (int eq 1) then begin
   openw,2,fname+'.dat'
   printf,2,     'File         GRP   x   y  nx  ny       sky    Weight'
   printf,logfil,'File         GRP   x   y  nx  ny       sky    Weight'
   print,'input, size of PSF to extract, nx,ny (odd), and max number nmax'
   nx=0 & ny=0 & nmax=0 &read,nx,ny,nmax
   m=300/nx
   print,'group numbers to process,gmin,gmax'                 
   gmin=0 & gmax=0 & read,gmin,gmax
   psfm=fltarr(nx,ny,nmax) & wp=fltarr(nmax) & xp=wp & yp=wp & gp=wp
   np=0
   for igr=gmin,gmax do begin
      reads,fname+'.d0h',data,nxf,nyf,igr
      tvwindow,0,xsize=nxf,ysize=nyf,title=' Group '+strtrim(igr,1)
      tv,data
      print,' Set colour lookup to find stars'
      color_edit                       
      adjct
      for ip=np,nmax-1 do begin
         print,' point to star to be used - Right button to end'
         cursor,xc,yc,3,/dev
         tvwindow,1,xsize=m*nx,ysize=m*ny &wset,0
      if(!err eq 4) then goto, nextgroup
         print,'input sky box lower left ' & cursor,xs1,ys1,3,/dev
         print,'input sky box upper right' & cursor,xs2,ys2,3,/dev
         sky=median(data(xs1:xs2,ys1:ys2)) & print,'sky=',sky
         xval=indgen(nx)#replicate(1.,ny)   & yval=replicate(1.,nx)#indgen(ny)
         cosmic1: psf=data(xc-nx/2:xc+nx/2,yc-nx/2:yc+nx/2)-sky
                  ox=total(psf*xval)/total(psf)-nx/2 
                  oy=total(psf*yval)/total(psf)-ny/2
                  print,' offsets',ox,oy 
                  psf=data(xc+ox-nx/2:xc+ox+nx/2,yc+oy-nx/2:yc+oy+nx/2)-sky 
                  if(nx lt 128) then begin
                     print, ' Mark bad pixels - right button to stop'
                     wset,1&tvscl,congrid(psf,m*nx,m*ny)
                     tvrdc,xo,yo,/dev&wset,0
                     if(!err eq 4) then goto, normal1
                     xcos=xc+ox-nx/2+xo/m & ycos=yc+oy-nx/2+yo/m
                     if(xcos le nx-3 and ycos le ny-3 and $
                        xcos ge 2    and ycos ge 2) then begin
                        box=data(xcos-2:xcos+2,ycos-2:ycos+2)
                        ixbox=where(box eq max(box)) & inbox=inbox(0)
                        xcos=xcos+ixbox-5*(ixbox/5)-2 & ycos=ycos+ixbox/5-2 
                      endif
                      data(xcos,ycos)=data(xcos+1,ycos)
                      goto, cosmic1
                  endif else begin
                      goto, normal1
                  endelse                   
         normal1: xc=xc+ox &yc=yc+oy
         psf=psf>0. & wp(np)=total(psf) & psf=psf/total(psf) & psfm(0,0,np)=psf 
         xp(np)=xc & yp(np)=yc  & gp(np)=igr
         printf,2,    '$(a13,i3,4i4,2f10.2)',fname,igr,xc,yc,nx,ny,sky,wp(np) 
         printf,logfil,'$(a13,i3,4i4,2f10.2)',fname,igr,xc,yc,nx,ny,sky,wp(np) 
         np=np+1
      endfor
      nextgroup:
;         if(nx lt 128) then begin
;            tvwindow,1,xsize=600,ysize=600 ,title='EXTRACTED PSFS'
;            for ip=0,np-1 do begin
;                tvscl,congrid(psfm(*,*,ip),m*nx,m*ny),ip
;            endfor
;            color_edit
;         endif
   endfor
   close,2
endif else begin
   openr,2,fname+'.dat' 
   line=' ' & readf,2,line & print,line & printf,logfil,line
   nmax=21
   for iline=0,nmax-1 do begin
      if(iline eq 0) then begin
        oigr=-1  & ofname=' ' & nxr=0 & nyr=0 & xc=0 & yc=0
      endif else begin
        oigr=igr & ofname=fname
      endelse
   if(eof(2)) then goto,endprog
      readf,2,'$(a13,i3,4i4,2f10.2)',fname,igr,xc,yc,nxr,nyr,sky,wpt
      print,        '$(a13,i3,4i4,2f10.2)',fname,igr,xc,yc,nxr,nyr,sky,wpt 
      printf,logfil,'$(a13,i3,4i4,2f10.2)',fname,igr,xc,yc,nxr,nyr,sky,wpt 
      if(iline eq 0) then begin
         nx=nxr & ny=nyr & np=0
         m=300/nx
         psfm=fltarr(nx,ny,nmax) & wp=fltarr(nmax) & xp=wp & yp=wp & gp=wp
         tvwindow,1,xsize=m*nx,ysize=m*ny&wset,0
;         if(nx lt 128) then tvwindow,1,xsize=600,ysize=600 ,title='EXTRACTED PSFS'      
      endif else begin
        np=np+1
        if((nx ne nxr) or (ny ne nyr) or (ofname ne fname)) then begin
           print,'*** error***',nx,nxr,ny,nyr,ofname,fname
           stop
        endif   
      endelse
      if(oigr ne igr) then reads,fname+'.d0h',data,nxf,nyf,igr
      if(nx lt 128) then begin
         ncosm=0
         cosmic2: psf=data(xc-nx/2:xc+nx/2,yc-nx/2:yc+nx/2)-sky
               wset,1&tvscl,congrid(psf,m*nx,m*ny)
               if(iline eq 0 and ncosm eq 0) then color_edit
               print, ' Mark bad pixel',ncosm,'   - right button to stop'
               wset,1 & tvrdc,xo,yo,/dev & wset,0
               if(!err eq 4) then goto, normal2
               ncosm=ncosm+1
               xcos=xc-nx/2+xo/m & ycos=yc-nx/2+yo/m
               if(xcos le nx-3 and ycos le ny-3 and $
                  xcos ge 2    and ycos ge 2) then begin
                  box=data(xcos-2:xcos+2,ycos-2:ycos+2)
                  ixbox=where(box eq max(box))  & inbox=inbox(0)
                  xcos=xcos+ixbox-5*(ixbox/5)-2 & ycos=ycos+ixbox/5-2 
               endif
               data(xcos,ycos)=data(xcos+1,ycos)
               goto, cosmic2
         normal2:
      endif else begin
         psf=data(xc-nx/2:xc+nx/2,yc-nx/2:yc+nx/2)-sky
      endelse 
      psf=psf>0. & wp(np)=total(psf) & psf=psf/total(psf) & psfm(0,0,np)=psf 
      xp(np)=xc & yp(np)=yc & gp(np)=igr
      print,' Weight computed, infile:',wp(np),wpt
;      if(nx lt 128) then tvscl,congrid(psfm(*,*,np),m*nx,m*ny),np
   endfor 
   endprog: 
   close,2
   np=np+1
   if(nx lt 128) then color_edit
endelse
wp=wp/total(wp)
return
end