getpsf.pro

basic median filter simulation

 dxcen=0.  &  dycen=0.
;
 niter = 0                    ;Beginning of big iteration loop
 v = fltarr(5)
 c = fltarr(5,5)
;                            Print the current star
 fmt1 = "(/17X, 'STAR', 5X, 'X', 8X, 'Y', 5X, 'MAG  1', 5X, 'SKY')"
 fmt2 = "(15X, I5, 2F9.2, 12F9.3)"
 print,format=fmt1          
 print,format=fmt2,istar, xc[istar], yc[istar], mag[istar], sky[istar]

 if keyword_set(DEBUG) then print,'GETPSF: Gaussian Fit Iteration'

 REPEAT BEGIN		     ;Begin the iterative loop

 niter = niter + 1
 if niter GT 100 then begin   ;No convergence after 100 iterations?
    message,'No convergence after 100 iterations for star ' + strn(istar),/INF
    goto, GETSTAR 
 endif

      if sigx LE 1 then nx = 1  $  ;A default box width 
 else if sigx GT 3 then nx = 3  $
 else                   nx = 2

      if sigy LE 1 then ny = 1  $
 else if sigy GT 3 then ny = 3  $
 else                   ny = 2

 a = [H, x+dxcen,y+dycen,sigx,sigy]
 xin = (findgen(2*nx+1)-nx) + ix
 yin = (findgen(2*ny+1)-ny) + iy
 make_2d, xin, yin
 DAOERF, xin, yin, a, g, t

;  The T's are the first derivatives of the model profile with respect
;  to the five fitting parameters H, DXCEN, DYCEN, SIGX, and SIGY.
;  Note that the center of the best-fitting Gaussian profile is
;  expressed as an offset from the centroid of the star.  In the case of
;  a general, asymmetric stellar profile, the center of symmetry of the
;  best-fitting Gaussian profile will not necessarily coincide with the
;  centroid determined by any arbitrary centroiding algorithm.  

 dh = f[ ix-nx:ix+nx, iy-ny:iy+ny] - g ;Subtract best fit Gaussian from subarray
 for kk = 0,4 do begin
      tk = t[*,kk]
      v[kk] = total( dh * tk )
      for ll = 0,4 do c[kk,ll] = total( tk * t[*,ll] )
 endfor

 c = invert(c,status)	;IDL version assumes INVERT is successful

 if status EQ 1 then begin
     message,'Singular matrix encountered fitting star ' + strn(istar),/INF
     goto, GETSTAR
 endif 

 z = c#v         ;Multiply by vector of residuals

 h = h + z[0]/(1.0+4.0*abs(z[0]/h))	;Correct the fitting parameters
 dxcen = dxcen+z[1]/(1.0+3.0*abs(z[1]))
 dycen = dycen+z[2]/(1.0+3.0*abs(z[2]))
 sigx = sigx+z[3]/(1.0+4.0*abs(z[3]/sigx))
 sigy = sigy+z[4]/(1.0+4.0*abs(z[4]/sigy))

 if keyword_set(DEBUG) then print,niter,h,dxcen,dycen,sigx,sigy

 endrep until $ 				;Test for convergence  
       (abs(z[0]/h)+abs(z[3]/sigx)+abs(z[4]/sigy) LT 0.0001)

;  Now that the solution has converged, we can generate an
;  array containing the differences between the actual stellar profile
;  and the best-fitting Gaussian analytic profile.

 a = [H, x+dxcen, y+dycen, sigx,sigy]  ;Parameters for Gaussian fit
 DAOERF,xgen,ygen,a,g                  ;Compute Gaussian
 f = f - g                             ;Residuals (Real profile - Gaussian)

 psfmag = mag[istar]
 xpsf1 = xc[istar] & ypsf1 = yc[istar]

; The look-up table is obtained by interpolation within the array of
; fitting residuals.  We need to interpolate because we want the look-up
; table to be centered accurately on the centroid of the star, which of 
; course is at some fractional-pixel position in the original data.

 ncen = (npsf-1)/2.
 psfgen = (findgen(npsf) - ncen)/2.         ;Index function for PSF array
 YY = psfgen + Y   &  XX = psfgen + X
 make_2d,xx,yy
 psf = RINTER(F, XX, YY)            ;Interpolate residuals onto current star
 gauss = [h,dxcen,dycen,sigx,sigy]
 goodstar = nstrps                   ;Index of first good star

; For each additional star, determine the precise  coordinates of the 
; centroid and the relative brightness of the star
; by least-squares fitting to the current version of the point-spread
; function.  Then subtract off the appropriately scaled integral under
; the analytic Gaussian function  and add the departures of the actual 
; data from the analytic Gaussian function to the look-up table.

GETMORE:            ;Loop for additional PSF stars begins here                 
 nstrps = nstrps+1
 if nstrps GE numpsf then goto,WRITEOUT	;Have all the stars been done?

 istar = idpsf[nstrps]
 ixcen = fix(xc[istar])
 iycen = fix(yc[istar])                  
 scale = 10.^(-0.4*(mag[istar]-psfmag))

; Fit the current version of the point-spread function to the data for
; this star.

 lx = ixcen-nhalf+1 & ux =ixcen + nhalf
 ly = iycen-nhalf+1 & uy =iycen + nhalf
 if ( (lx LT 0) or (ly LT 0) or $             ;Star too close to edge?
    (ux GE ncol) or (uy GE nrow)) then begin  
   print,'GETPSF: Star ',strn(istar),' too near edge of frame.'
   goto,GETMORE
 endif                      

 print,format=fmt1
 print,format=fmt2, istar, xc[istar], yc[istar], mag[istar], sky[istar]

 f = image[lx:ux,ly:uy]
 x = xc[istar]-lx   &   y = yc[istar]-ly   

 pkfit, f, scale, x, y, sky[istar], fitrad, ronois, phpadu, $
		gauss, psf, errmag, chi, sharp, niter, DEBUG = debug

 if niter EQ 25 then begin	;Convergence in less than 25 iterations?
      print,'GETPSF: No convergence after 25 iterations for star',istar
      goto, GETMORE 
 endif

 a = [gauss[0], x+dxcen,y+dycen,sigx,sigy]  ;Parameters of successful fit
 daoerf,xgen,ygen,a,e
 f = f - scale*e -sky[istar]	           ;Compute array of residuals

; Values of the array of residuals are now interpolated to an NPSF by
; NPSF (NPSF is an odd number) array centered on the centroid of the
; star, and added to the existing look-up table of corrections to the 
; analytic profile 

 xx = psfgen + x
 yy = psfgen + y 
 make_2d,xx,yy
 psf = psf + RINTER(f,xx,yy)    

; Now correct both the height of the analytic Gaussian, and the value
; of the aperture-magnitude of the point-spread function for the
; inclusion of the additional star.

 psfmag = -2.5*alog10((1.+scale)*10^(-0.4*psfmag))
 gauss[0] = gauss[0]*(1.+scale)
 goodstar = [ goodstar, nstrps]
 goto, GETMORE  

WRITEOUT:   

; Create FITS file containing the PSF created.

 if ( N_elements(psfname) EQ  0 ) then begin
   psfname=''
   read,'Enter name of FITS file to contain final PSF ([RETURN] to exit): ',psfname
 endif

if ( psfname EQ  '' ) then return

 mkhdr, hdr, psf         ;Create a minimal FITS header
 sxaddpar, hdr, 'PHPADU', phpadu, 'Photons per Analog Digital Unit'
 sxaddpar, hdr, 'RONOIS', ronois, 'Readout Noise'
 sxaddpar, hdr, 'PSFRAD', psfrad, 'Radius where PSF is defined (pixels)'
 sxaddpar, hdr, 'FITRAD', fitrad, 'Fitting Radius'
 sxaddpar, hdr, 'PSFMAG', psfmag, 'PSF Magnitude'
 sxaddpar, hdr, 'GAUSS1', gauss[0], 'Gaussian Scale Factor'
 sxaddpar, hdr, 'GAUSS2', gauss[1], 'Gaussian X Position'
 sxaddpar, hdr, 'GAUSS3', gauss[2], 'Gaussian Y Position'
 sxaddpar, hdr, 'GAUSS4', gauss[3], 'Gaussian Sigma: X Direction'
 sxaddpar, hdr, 'GAUSS5', gauss[4], 'Gaussian Sigma: Y Direction'

 ngood = N_elements(goodstar)
 sxaddhist,'GETPSF: '+ systime() + ' ' + strn(ngood) +  $
           ' Stars Used to Create PSF',hdr

  sxaddhist,'GETPSF: ID - '+ string(idpsf[goodstar[0:12<ngood-1]], $
    format='(13i5)'),hdr

 if ngood gt 13 then $
     sxaddhist,'GETPSF: ID - '+ string(idpsf[goodstar[13:*]], $
     format='(13i5)'),hdr

 sxaddhist,'PSF Coordinates:'+ $
    string(xpsf1, format='(F7.2)') + $
    string(ypsf1, format='(F7.2)'), hdr

 writefits,psfname,psf,hdr

 return
 end