function euvi_psf,nxy,index,center=coarr
;+
; $Id: euvi_psf.pro,v 1.1 2008/04/23 13:38:28 euvi Exp $
;
; Project : STEREO SECCHI
;
; Name    : euvi_psf
;
; Purpose : returns the EUVI point spread function, including the
;           diffraction pattern created by the entrance filters
;
; Use     : IDL> psf = euvi_psf(nxy,index)
;
; Inputs  :
;    nxy = x and y size of psf
;    index = image header structure (or simply {wavelnth:304,obsrvtry:'a'})
;
; Outputs :
;    psf = fltarr(nxy,nxy), point spread function for full resolution images
;
; Keywords:
;    /center : center PSF on center of array, even if nxy even.  Can cause
;           center of PSF to fall onto corner between pixels.
;           Default centers PSF on center of pixel nearest to array center.
;           Default is compatible with trace_diffraction_pattern.pro.
;           Has no effect if nxy odd.
; Common  : None
;
; Restrictions:
;
; Side effects:
;
; Category    : Image analysis
;
; Prev. Hist. :
;
; Written     : Jean-Pierre Wuelser, LMSAL, 17-Apr-2008
;               inspired by Tom Metcalf's trace_diffraction_pattern.pro
;
; $Log: euvi_psf.pro,v $
; Revision 1.1  2008/04/23 13:38:28  euvi
; initial version
;
;-

; check input
sz = size(index)
if n_elements(nxy) ne 1 or sz[sz[0]+2] ne 1 or sz[sz[0]+1] ne 8 then begin
   print,'Invalid input.  Must be psf=euvi_psf(integer,headerstruct)'
   return,0
endif

; set wav,obs as indices into parameter arrays
wavs = [171,195,284,304]
wav = where(index.wavelnth eq wavs)
if wav ge 0 then wav=wav[0] else wav=0L
obs = strupcase(strmid(index.obsrvtry,strlen(index.obsrvtry)-1,1))
if obs eq 'B' then obs=1 else obs=0

; parameters (4 x 2 each)
pfac = !radeg*3.6e3/(1.59*5e7) ; wavelnth*pfac = pixelsep of orders (5mm grid)
per = [[173.5*pfac, 195.0*pfac, 10.18, 10.88], $    ; ordersep Ahead (pixels)
       [173.5*pfac, 195.0*pfac, 10.18, 10.88]]      ; ordersep Behind
bet = [[0.98, 0.98, 0.90, 0.918], $           ; (wiresep-wirethick)/wiresep (A)
       [0.98, 0.98, 0.90, 0.908]]             ; (wiresep-wirethick)/wiresep (B)
ang = [[-0.7854, -0.7854, -0.3420, -0.0191], $   ; angle of diffraction (rad,A)
       [-0.7854, -0.7854, -0.8610, -0.0385]]     ; angle of diffraction (rad,B)

pp = per[wav,obs]
pb = bet[wav,obs]
pa = ang[wav,obs]

; high resolution version of PSF core (Voigt function)
nxyc = 99L < (2L*floor((nxy-1L)/2L)+1L)     ; size of high res core (odd)
va = 0.02   ; a parameter of voigt profile
vw = 0.75   ; width of core
x = (findgen(nxyc*10)-nxyc*5.0+0.5)/10.0
r = sqrt((x*x)#(fltarr(nxyc*10)+1.0)+(fltarr(nxyc*10)+1.0)#(x*x))
cor = voigt(va,r/vw)
; extended core at normal resolution
exy = long(nxy-nxyc)/2L
x = findgen(nxy)-float(exy)-float(long(nxyc)/2L)
r = sqrt((x*x)#(fltarr(nxy)+1.0)+(fltarr(nxy)+1.0)#(x*x))
ecor = voigt(va,r/vw)                        ; extended core
ecor[exy:exy+nxyc-1L,exy:exy+nxyc-1L] = 0.0  ; punch hole for high res. core

; center of psf in output array
if keyword_set(coarr) then cxy=float(nxy-1L)/2.0 $
                      else cxy=float(long(nxy)/2L)

; maximum order to consider
maxo = long((float(nxy)-2.0)/(2.0*pp)) ; < 50 ; set maximum diff order to use

psf = fltarr(nxy*2L,nxy*2L)      ; output array with padding

; go through orders
for j=-maxo,maxo do for k=0,1 do if k eq 0 or j ne 0 then begin
  ox = cxy + float((1-k)*j)*cos(pa)*pp - float(k*j)*sin(pa)*pp
  oy = cxy + float((1-k)*j)*sin(pa)*pp + float(k*j)*cos(pa)*pp
  ix = floor(ox+0.5)
  iy = floor(oy+0.5)
  fx = ox-ix          ; between -0.5 and 0.5
  fy = oy-iy
  scor = ecor
  scor[exy,exy] = rebin(shift(cor,floor(fx*10),floor(fy*10)),nxyc,nxyc)
  if j eq 0 then amp=1.0 else amp=(sin(j*!pi*pb)/(j*!pi*pb))^2.0
  psf[ix,iy] = psf[ix:ix+nxy-1L,iy:iy+nxy-1L] + amp*scor
endif

; cut off padding
psf = psf[exy+nxyc/2L:exy+nxyc/2L+nxy-1L,exy+nxyc/2L:exy+nxyc/2L+nxy-1L]

psf = psf/total(psf)   ; normalize

return,psf

end