;+ ; NAME: ; PH_PSF ; ; PURPOSE: ; Calculate theoretical PSF for a PHARO image. ; ; EXPLANATION: ; Calculates the diffraction-limited monochromatic PHARO point ; spread function (PSF) for any camera, pupil stop, and pupil angle. ; Based on NIRC2PSF, writen by A. Bouchez, W.M. Keck Observatory. ; ; CALLING SEQUENCE: ; result = PH_PSF( [NPIX=, POS=, HD=, CAROUSEL=, LYOT=, ; FILTER=, CR_ANGLE=, GRISM=, PUPIL=]) ; ; INPUTS: ; none. ; ; OUTPUTS: ; result = image of PSF ; ; OPTIONAL OUPUTS: ; PUPIL = pupil image used to generate PSF ; ; OPTIONAL INPUT KEYWORDS: ; NPIX = size of pupil image, in pixels ; ; POS = position of PSF wrt image center (pixel-edge coordinates!) ; ; HD = string array containing PHARO FITS header ; ; The following 4 keywords override the header fields in HD if set: ; ; CAROUSEL = camera name, eg. 'narrow' ; ; LYOT = Lyot stop name, eg. 'largehex' ; ; FILTER = Filter name, eg. 'Br-gamma' ; ; CR_ANGLE = Cassegrain ring angular position (for rotated pupil images). ; ; GRISM = Grism wheel position. Only used to create ND filter ghosts. ; ; PROCEDURES USED: ; Procedures: DIST_CIRCLE ; ; MODIFICATION HISTORY: ; Adapted from NIRC2PSF 12/05, A. Bouchez, Caltech Optical Observatories ;- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FUNCTION PH_PUPIL, NPIX=npix, DU=du, LYOT=lyot, CR_ANGLE=cr_angle if not KEYWORD_SET(npix) then npix = 256 if not KEYWORD_SET(du) then $ du = 2.124e-6 / (npix * 0.00994 / 206265) if not KEYWORD_SET(lyot) then lyot = 'Std Cross' if not KEYWORD_SET(cr_angle) then cr_angle = -24.21 ;;; 1. Define dimentions of pupil in inches based on engineering drawings. case STRUPCASE(STRTRIM(lyot,2)) of 'CIRCULAR': d = [0.00, 16.88, 0.00] ; fake pupil for demo only 'OPEN': d = [6.08, 16.88, 0.20] ; assumed spider thickness! 'STD CROSS': d = [6.10, 16.50, 0.25] ; diam & thickness in mmm 'MED CROSS': d = [7.60, 15.40, 0.83] 'LG CROSS': d = [9.40, 13.50, 1.75] 'STD SPOT': d = [6.10, 16.50, 0.20] else: begin MESSAGE, /info, "lyot '" + lyot + "' not recognized." MESSAGE, /info, "Using default 'Std Cross'." d = [6.10, 16.50, 0.25] end endcase lyot_pscl = 16.88/5.093 ; mm/m cr_zero = -24.21 ; CHECK THIS!!! pupil = BYTARR(npix,npix) DIST_CIRCLE,r,npix,npix/2-0.5,npix/2-0.5 w = WHERE(r*du*lyot_pscl gt d[0]/2 and $ r*du*lyot_pscl lt d[1]/2) if w[0] ne -1 then pupil[w] = 1B v = [[-1*d[2]/2,0], [d[2]/2,0], [d[2]/2,d[1]*1.1], $ [-1*d[2]/2,d[1]*1.1], [-1*d[2]/2,0]] ang = (90 * DINDGEN(6) - cr_angle + cr_zero) * !dtor for i=0,3 do begin rmat = [[-1*SIN(ang[i]), COS(ang[i])], $ [ COS(ang[i]), SIN(ang[i])]] rv = npix/2 + (rmat#v / (du * lyot_pscl)) w = POLYFILLV(rv[1,*],rv[0,*],npix,npix) if w[0] ne -1 then pupil[w] = 0B endfor RETURN,pupil END ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FUNCTION PH_PSF, NPIX=npix, POS=pos, HD=hd, CAROUSEL=carousel,$ LYOT=lyot, FILTER=filter, CR_ANGLE=cr_angle, $ GRISM=grism, PUPIL=pupil if not KEYWORD_SET(npix) then npix = 256 if not KEYWORD_SET(pos) then pos = [0.,0] if KEYWORD_SET(hd) then begin if not KEYWORD_SET(carousel) then carousel = SXPAR(hd, 'CAROUSEL') if not KEYWORD_SET(lyot) then lyot = SXPAR(hd, 'LYOT') if not KEYWORD_SET(filter) then filter = SXPAR(hd, 'FILTER') if not KEYWORD_SET(cr_angle) then cr_angle = SXPAR(hd, 'CR_ANGLE') if not KEYWORD_SET(grism) then grism = SXPAR(hd, 'GRISM') endif else begin if not KEYWORD_SET(carousel) then carousel = '25 mas' if not KEYWORD_SET(lyot) then lyot = 'Std Cross' if not KEYWORD_SET(filter) then filter = 'Br-gamma' if not KEYWORD_SET(cr_angle) then cr_angle = -24.21 if not KEYWORD_SET(grism) then grism = 'open' endelse ;;; 1. Calculate PHARO pupil ;;; We require du<0.10 m/pix(for sufficient detail) and $ ;;; npix*du>6.0 m(to avoid truncating pupil). This is done by increasing ;;; the platescale and resampling and/or increasing the pupil image size. case STRUPCASE(STRTRIM(carousel,2)) of '25 MAS': pscl = 25.22e-3 / 206265 ; From Metchev web page '40 MAS': pscl = 40.00e-3 / 206265 ; Uncalibrated! else: begin MESSAGE, /info, "Camera '" + carousel + "' not recognized.' MESSAGE, /info, "Using default '25 mas'." pscl = 25.22e-3 / 206265 end endcase case STRUPCASE(STRTRIM(filter,2)) of 'BR-GAMMA': effwave = 2.18 ; microns 'K_SHORT': effwave = 2.145 ; microns 'H': effwave = 1.650 ; microns GET CORRECT WAVELENGTH 'J': effwave = 1.250 ; microns GET CORRECT WAVELENGTH else: begin MESSAGE, /info, "Filter name '" + filter + "' not recognized." MESSAGE, /info, "Using default 'Br-gamma'." effwave = 2.18 end endcase tmp = pscl * 6.0 / (effwave*1e-6) rpfac = (2.^CEIL(ALOG( tmp )/ALOG(2)))>1 ; reb. 2^N so npix*du>6.0 m pscl1 = pscl / rpfac ; reb. detector platescale npix1 = npix * rpfac ; reb. image size du = (effwave*1e-6) / (npix1 * pscl1) ; pupil image platescale, m/pix rdfac = (2.^CEIL(ALOG( du / 0.10 )/ALOG(2)))>1 ; reb. 2^N so du<0.10 npix2 = npix1 * rdfac ; reb. pupil image size du = (effwave*1e-6) / (npix2 * pscl1) ; reb. pupil image pscl, m/pix pupil = PH_PUPIL(npix=npix2,du=du,lyot=lyot,cr_angle=cr_angle) ;;; 2. Create phase ramp to position PSF. Note that the -[0.5,0.5] ;;; forces the PSF to lie on the pixel edges if pos=[0,0]. uu = REBIN(FINDGEN(npix2,1),npix2,npix2) vv = REBIN(FINDGEN(1,npix2),npix2,npix2) rpos = pos * rpfac - 0.5 phase = 2 * !pi * (uu*rpos[0] + vv*rpos[1]) / npix2 ;;; 3. Compute PSF by fast Fourier transform wavefront = pupil * EXP(COMPLEX(0,phase)) rpsf = SHIFT(ABS(FFT(wavefront,-1))^2, npix2/2, npix2/2) psf = REBIN(rpsf(npix2/2-npix1/2:npix2/2+npix1/2-1, $ npix2/2-npix1/2:npix2/2+npix1/2-1),npix,npix) ;;; 4. Add ghost reflection for ND filters (only implemented for ND2 so far) case STRUPCASE(STRTRIM(grism,2)) of 'ND 1%': ghost = [0.256, -0.492, 0.015] ; dx, dy (arcsec), flux ratio else: ghost = [0.0, 0.0, 0.0] ; ONLY APPROXIMATE! endcase gpix = ghost[0:1] / (pscl * 206265) ; dx, dy (pixels) psf = psf + ghost[2]*FSHIFT(psf, gpix[0], gpix[1]) RETURN,psf/TOTAL(psf) END