pro mkgrid, n, x, y, c=c, d=d, s=s
;+
; NAME:
;      mkgrid
;;
; PURPOSE:
;      Make arrays containing the values of the x and  y coordinates of a grid. 
;
; INPUTS:
;      n1,n2 - Integers. Array dimensions (number of pixels).
;
; OPTIONAL INPUTS:
;      c0 - Floating point array [c1,c2] or c1 in both dimensions. Pixel coordinate of the 
;           [0.0,0.0] origin. Default: the centre of the array.
;      dist - make distance array. each pixel contains distance to origin of the array.
;
; KEYWORD PARAMETERS:
;      None.
;
; OUTPUTS:
;      x,y - Floating point [n1*n2] array containing x or y values (all of the values
;           in the other dimension are identical)
;-

COMPILE_OPT idl2

; Set reference point to the centre if not specified
if n_elements(n) eq 1 then dim=[n,n] else dim=n>1
if n_elements(c) eq 0 then c=.5*float(dim)-0.5 else c=[0,0]+c

; Make x and y arrays
; but only if the actual keyword is set
if keyword_set(routine_names(x,arg_name=-1)) then x = rebin(findgen(dim[0])-c[0],dim[0],dim[1],/sample)
if keyword_set(routine_names(y,arg_name=-1)) then y = rebin(findgen(1,dim[1])-c[1],dim[0],dim[1],/sample)

; for distance array it is faster to recalculate than to take the sqrt of the x and y arrays
if keyword_set(routine_names(d,arg_name=-1)) then $
  d=sqrt(rebin((findgen(dim[0])-c[0])^2,dim[0],dim[1],/sample) + rebin((findgen(1,dim[1])-c[1])^2,dim[0],dim[1],/sample))

; even faster return the d^2 arrays
if keyword_set(routine_names(s,arg_name=-1)) then s=rebin((findgen(dim[0])-c[0])^2,dim[0],dim[1],/sample) + rebin((findgen(1,dim[1])-c[1])^2,dim[0],dim[1],/sample)
end



pro getpsf, fpsf, ra, dec, psf, pa=pa, display=display, pixscl=pixscl, tsz=tsz, totcor=totcor
;+
; NAME:
;      getpsf
;;
; PURPOSE:
;      Reconstruct the IRAC PSF at the specific location using the kinformation stored in the GREATS psf file. 
;
; INPUTS:
;      fpsf - string. Name of the psf file
;	   ra,dec - Floating point array. List of positions, in degrees, where the PSF need to be reconstruted.
;
; OPTIONAL INPUTS:
;      pa - Floating point array. List of additional position angles requested for the final PSF. Default: 0.0
;      pixscl - Floating point value. Pixel scale of the output PSF. Default: adopt the same pixels scale as in the GREATS psf file
;      tsz - Floating point value. Size of the output PSF, in pixels. Default: adopt the same size as the master PSF in the GREATS psf file (81 pixels).
;
; KEYWORD PARAMETERS:
;	   display - Keyword. If set, the PSF is displayed at the end of the process. Default: False.
;
; OUTPUTS:
;	   psf - Floating point array containing the PSFs.
;-



  common getpsf_aor_com, mpsf, hpsf, map, hmap, current_psf

  if not keyword_set(pa) then pa = 0.0
  
  if n_elements(current_psf) le 0 then current_psf=''
  
  if not strcmp(current_psf,fpsf) then begin 
      print,'getpsf --> restoring ',fpsf
      mpsf = readfits(fpsf, hpsf, exten=1,/silent)
      map = readfits(fpsf, hmap, exten=2, /silent)
      current_psf = fpsf
  end 
  szpsf = size(mpsf,/dim)
  sz = size(map,/dim)
  len = sz[3]
  xy0 = (szpsf-1)/2.
  mkgrid, szpsf, d=d
  
  if keyword_set(pixscl) then begin
     if abs(pixscale(hpsf)/pixscl-1) gt 1e-3 then begin
         rscl = pixscale(hpsf)/pixscl 
         if keyword_set(tsz) then sznew=tsz else  sznew = oddsize(round(rscl*szpsf[0]),div=3)

		 _gpsf=frebin(mpsf,rscl*szpsf[0],rscl*szpsf[0])

		 szz=size(_gpsf,/dim)
		 sznew=oddsize(szz[0],div=3)
		 
		 if sznew ne szz[0] then begin
			 tmp=fltarr([sznew, sznew])
			 tmp[0,0]=_gpsf
			 _gpsf=tmp

	 	 endif
		 
		 gcntrd,_gpsf,sznew/2.,sznew/2.,xc,yc,4
		 gpsf=fshift(_gpsf,sznew/2.-xc-0.5,sznew/2.-yc-0.5)
		 
     end else begin
	     gpsf = mpsf
	     rscl=1.00
		 sznew = oddsize(round(rscl*szpsf[0]),div=3)
	 ;endif
 	endelse
   endif else begin
     gpsf = mpsf
     rscl=1.00
	 sznew = oddsize(round(rscl*szpsf[0]),div=3)
   end

   adxy, hmap, ra[*], dec[*], x, y  
   ix = round(x)
   iy = round(y)

   nx=n_elements(x)

  psf = fltarr([size(gpsf,/dim),nx])
  tmp = fltarr(size(gpsf,/dim))
  for j=0,nx-1 do begin
  
    if total(map[ix,iy,0,1:*]) eq 0 then  return
     mkgrid, sznew, d=d
    for i=0,len-1 do begin
       aor_w = map[ix,iy,0,i]
       aor_pa = map[ix,iy,1,i]
       if aor_w[0] le 0 then continue 
       tmp += rot(gpsf,-aor_pa[0]-pa,rscl[0],xy0[0],xy0[1],cubic=-0.5,missing=0.0)*aor_w[0]/rscl^2
       if keyword_set(display) then begin
          print, i+1, aor_w, aor_pa, format='(i,3f)'
          tvs, rot(gpsf,-aor_pa[0]-pa,rscl,xy0[0],xy0[1],cubic=-0.5,missing=0.0), fac=10, pos=0
          tvs, tmp, fac=10, pos=1
       end
    end
    
    psf[*,*,j] = tmp/total(tmp*(d lt  (sznew-1.0)/2.))
  end

  if keyword_set(display) then print, ra,dec, rscl, total(psf)
end