;+
; NAME:
;    p6model
; PURPOSE: (one line)
;    Generate synthetic PSF images for the HST Planetary Camera, Chip 6.
; DESCRIPTION:
;    A generated PSF is added to the output array.  If the array is not
; defined, it is created with dimensions 800 by 800.  Successive calls
; will add objects (PSF's) to the array.
; CATEGORY:
;  CCD data processing
; CALLING SEQUENCE:
;    p6model, x, y, inten, date, filter, bmvnum, back, image
; INPUTS:
;    x, y    : Position of the PSF in the output array.
;    inten   : Intensity of PSF relative to TinyTim output PSF
;    date    : Date of observation (YYMMDD format).
;    filter  : Filter name in the form Fxxx(W,LP,M,N).
;    bmvnum  : B-V list value:
;                1) -0.297 (Type O)
;                2) -0.155 (Type B)
;                3)  0.126 (Type A)
;                4)  0.395 (Type F)
;                5)  0.619 (Type G)
;                6)  1.110 (Type K)
;                7)  1.590 (Type M)
;    back    : Background to be added.
; OPTIONAL INPUT PARAMETERS:
;
;    JITTER  : If set contains the gaussian smearing to apply to the image.
;                This is in units of 1/e half width.  If not set, no smearing
;                is applied.
;
;    OBJRAD  : Radius of object in pixels.  If not set, assumed to be a point
;              source.  If set and radius is greater than 0.5, then PSF is
;              convolved with a circular function of this radius.  This models
;              the object as a Lambert disk, crude but useful for barely
;              resolved objects.
;
; KEYWORD PARAMETERS:
;    NEW     : If set, clears the output array to zeros before building the
;              PSF image(s).
;    HSTPATH : Alternate path for HST PSF's.
;    VERBOSE : If set, prints informational message.
; OUTPUTS:
;    image  : The image array into which the PSF is placed.
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;    The Tiny Tim programs are used, via the procedure HSTPSF, to generate
; a 2X (7.5 micron spacing) PSF at the nearest grid location to the requested
; position. (This grid is defined in hstpsf).
;    From the interpolated maximum returned by hstpsf, shift amounts for x
; and y are computed.
;    The PSF is shifted by these amounts and then compressed to a 1X
; (15 micron spacing) PSF.
;    Finally, the resulting PSF is multiplied by the intensity parameter and
; added into the output array.
;    Calls external procedures BADPAR, BOXM, HSTPSF, and SSHIFT2D.
; MODIFICATION HISTORY:
;    Written by Doug Loucks, Lowell Observatory, December, 1993.
;    12/10/1993, DWL, Modified to allow the first three parameters to be
; vectors.
;    12/16/1993, DWL, Added VERBOSE keyword.
;    1/19/94, DWL, Added HSTPATH keyword.
;    1/21/94, DWL, Added x-maximum and y-maximum parameters to hstpsf call.
;-
PRO p6model, in_x, in_y, in_inten, in_date, in_filter, in_bmvnum, in_back, $
             out_image, NEW=in_new,  $
             HSTPATH=in_path, VERBOSE=in_verbose, JITTER=jitter, OBJRAD=objrad

IF N_PARAMS() LT 8 THEN BEGIN
   PRINT, 'p6model, x, y, intensity, date, filter, bmvnum, back, image'
   RETURN
ENDIF

self = '% P6MODEL: '
IF badpar( in_x,   [2,3,4,5], [0,1], CALLER=self+'(x) ', NPTS=n1 ) THEN RETURN
IF badpar( in_y,   [2,3,4,5], [0,1], CALLER=self+'(y) ', NPTS=n2 ) THEN RETURN
IF badpar( in_inten, [2,3,4,5], [0,1], CALLER=self+'(intensity) ', NPTS=n3 ) $
   THEN RETURN
IF badpar( in_date,   7,       0,     CALLER=self+'(date) ' ) THEN RETURN
IF badpar( in_filter, 7,       0,     CALLER=self+'(filter) ' ) THEN RETURN
IF badpar( in_bmvnum, [2,3],   0,     CALLER=self+'(bmvnum) ' ) THEN RETURN
IF badpar( in_back,   [2,3,4,5], 0,   CALLER=self+'(back) ' ) THEN RETURN
IF badpar( objrad,    [0,2,3,4,5], [0,1], CALLER=self+'(x) ', NPTS=n4 ) THEN RETURN

;Check that the first three parameters have equal length.
t = [ n1, n2, n3 ]
w = WHERE( MAX( t ) NE MIN( t ), count )
IF count NE 0 THEN BEGIN
   msg = 'Error. x, y, and intensity parameters must have the same length.'
   MESSAGE, msg, /INFO
   RETURN
ENDIF

; Verify valid OBJRAD if supplied
IF Keyword_set(objrad) AND n4 NE n1 THEN BEGIN
   msg = 'Error. OBJRAD must have the same length as the x and y parameters.'
   MESSAGE, msg, /INFO
   RETURN
ENDIF

stat = SIZE( out_image )
ndim = stat[ 0 ]

IF ndim NE 2 THEN BEGIN
   ;The output array must be created.
   im_xsize = 800
   im_ysize = 800
   out_image = FLTARR( im_xsize, im_ysize )
ENDIF ELSE BEGIN
   ;It exists.  Get its dimensions.
   im_xsize = stat[ 1 ]
   im_ysize = stat[ 2 ]
ENDELSE

IF KEYWORD_SET( in_new ) THEN out_image[0:im_xsize-1,0:im_ysize-1]=0

IF KEYWORD_SET( in_path ) THEN BEGIN
   path = in_path
ENDIF ELSE BEGIN
   path = ''
ENDELSE

;Define the bin factor for compression.
binfact = 0.5

FOR j = 0, n1-1 DO BEGIN
   IF KEYWORD_SET( in_verbose ) THEN BEGIN
      MESSAGE, 'Working on object at ' + STRING( in_x[j], in_y[j], $
      FORMAT='(2F12.4)' ), /INFO
   ENDIF
   hstpsf, in_x[j], in_y[j], in_date, in_filter, in_bmvnum, psf, xm, ym, $
   HSTPATH=path, VERB=in_verbose

   psfstat = SIZE( psf )
   xsize = psfstat[ 1 ]
   ysize = psfstat[ 2 ]

   ;Compute the required shift, allowing for the subsequent rebin operation.
   ix = FIX( in_x[j] )
   iy = FIX( in_y[j] )
   IF (ix LT 0) OR (ix GE im_xsize) THEN BEGIN
      t = STRING( im_xsize, FORMAT='(G0.0)' )
      MESSAGE, 'Error. X-position out of bounds. (0 <= X < ' + t + ').', /INFO
      RETURN
   ENDIF
   IF (iy LT 0) OR (iy GE im_ysize) THEN BEGIN
      t = STRING( im_ysize, FORMAT='(G0.0)' )
      MESSAGE, 'Error. Y-position out of bounds. (0 <= X < ' + t + ').', /INFO
      RETURN
   ENDIF
   fx = in_x[j] - ix
   fy = in_y[j] - iy
   ixm = FIX( xm )
   iym = FIX( ym )
   fxm = xm - ixm
   fym = ym - iym

   ;Create the shift vector.
   svec = [ fx+fx-fxm, fy+fy-fym ]

   ; Shift the psf.
   spsf = sshift2d( psf, svec )

   ; Convolve PSF with circular function if requested.
   if keyword_set(objrad) then begin
      ; Compute size of kernal
      nk = ceil(objrad[j]/binfact - 0.5) * 2 + 1
      ; If kernal is bigger than one pixel then do the convolution.
      if nk gt 1 then begin
         kc = nk / 2
         idx=indgen(nk)
         unit=replicate(1,nk)
         x=idx#unit
         y=unit#idx
         kernal=pixwt(kc,kc,objrad[j]/binfact,x,y)
         kernal=kernal/total(kernal)
         spsf=convol(spsf,kernal)
      endif
   endif

   ; rebin the psf.
   fpsf = REBIN( spsf, xsize*binfact, ysize*binfact )
   psfstat  = SIZE( fpsf )
   xsize    = psfstat[ 1 ]
   ysize    = psfstat[ 2 ]
   boxm, fpsf, xsize/2, ysize/2, 10, 10, xmax, ymax

   ;Compute the half-width of the PSF that will fit into the output array
   ;at the requested location.
   hw = MIN( [ xmax, ymax, xsize-xmax, ysize-ymax, ix, iy, im_xsize-ix, $
               im_ysize-iy ] )

   ;Add the appropriate sub-section of the PSF to the output array.
   out_image[ ix-hw : ix+hw-1, iy-hw : iy+hw-1 ] = $
   out_image[ ix-hw : ix+hw-1, iy-hw : iy+hw-1 ] + $
   fpsf[ xmax-hw : xmax+hw-1, ymax-hw : ymax+hw-1 ] * in_inten[j]
ENDFOR

if keyword_set(jitter) then begin
   x=[-2.0,-1.,0.,1.0,2.0]
   y=[-2.0,-1.,0.,1.0,2.0]
   i=[1.0,1.0,1.0,1.0,1.0]
   rsq = (x^2#i + i#y^2)/jitter^2
   kernal = exp(-rsq/2.)
   kernal=kernal/total(kernal)
   out_image=convol(out_image,kernal)
endif

out_image = out_image + in_back

END