#
procedure psfobs (in_, out_, be_, se_)
string in_="mxtools$data/psfprep.fits" {prompt="Name of input PSF image?"}
string out_="psfobs.fits" {prompt="Name of output OBS image?"}
real be_=50.0 {prompt="Background constant [electrons] ?"}
real bx_=0.0 {prompt="Background X slope value (zero for flat sky)?"}
real by_=0.0 {prompt="Background Y slope value (zero for flat sky)?"}
real se_=10000.0 {prompt="Stellar intensity [electrons] ?"}
real sxi_=36.0 {prompt="Central X postion of star (IRAF notation) ?"}
real syi_=36.0 {prompt="Central Y postion of star (IRAF notation) ?"}
real gain_=1.0 {prompt="Gain value [electrons/DN] ?"}
real ron_=3.0 {prompt="ReadOutNoise [electrons] ?"}
int seed_=INDEF {prompt="MKNOISE random number generator SEED ?"}
bool verbose=no {prompt="Verbose?"}
string version="2001AUG31"

begin
string in
string out
real sie
real b0e
int nx
int ny
real cpxi
real cpyi
real fwhm
real psf_cwxi
real psf_cwyi
real psf_cpxi
real psf_cpyi
int psf_nx
int psf_ny
real mag
string coo
string work
string comment
real hwhm
real obs_se
real obs_sxi
real obs_syi
real obs_be
real obs_bx
real obs_by
real se 
real sxi 
real syi 
real be 
real bx
real by
real gain
real ron
int seed
int psfprep_version

version="2001AUG31"

if (verbose) print  ( "PSFOBS: BEGIN" )

# Make sure correct packages are loaded
if (!defpac("noao")) {
  print "Please load the NOAO package"; 
  goto bye
}
if (!defpac("artdata")) {
  print "Please load the ARTDATA package of the NOAO package"; 
  goto bye
}
if (!defpac("images")) {
  print "Please load the IMAGES package of the NOAO package"; 
  goto bye
}
if (!defpac("imutil")) {
  print "Please load the IMUTIL package of the IMAGES package"; 
  goto bye
}
if (!defpac("stsdas")) {
  print "Please load the STSDAS package!"; 
  goto bye
}
if (!defpac("toolbox")) {
  print "Please load the TOOLBOX package of STSDAS"; 
  goto bye
}
if (!defpac("imgtools")) {
  print "Please load the IMGTOOLS package of STSDAS.TOOLBOX package"; 
  goto bye
}
if (!deftask("gstatistics")) {
  print "Could not execute the task STSDAS.TOOLBOX.IMGTOOLS.GSTATISTICS"; 
  goto bye
}
if (!deftask("imcalc")) {
  print "Could not execute the task STSDAS.TOOLBOX.IMGTOOLS.IMCALC"; goto bye
}

in = in_
if (!access(in)) {
  print ( "PSFOBS *** ERROR *** Can not access input file!" )
  goto bye;
}
out = out_
if (access(out)) {
  print ( "PSFOBS *** ERROR *** Output file already exists!" )
  goto bye;
}
work = ".psfobs.fits"

be = be_
bx_=0.0
bx = by_
by_=0.0
by=by_
se = se_
hselect ( in, "psfprep", yes) | scan ( psfprep_version )
if(nscan()==0) {
  goto bye
  print ( "PSFOBS *** ERROR *** Input file is not a MATPHOT PSF file!" )
}
hselect ( in, "psf_cpxi", yes) | scan ( sxi )
sxi_=sxi
hselect ( in, "psf_cpyi", yes) | scan ( syi )
syi_ = syi
gain_ = 1.0
gain = gain_
ron = ron_
seed = seed_

obs_se = se
obs_sxi = sxi
obs_syi = syi
obs_be = be
obs_bx = bx
obs_by = by

if (se>0.0) {

  delete( work, verify-, >& "dev$null" )
  imgtools.imcalc( in, work, "im1*"//se, >& "dev$null" )
  hedit ( work, "psfobs", value=version, add=yes, verify=no, update=yes, 
    > "dev$null")

  hedit ( work, "obs_be", value=be, add=yes, verify=no, update=yes, 
    > "dev$null")
  hedit ( work, "obs_bx", value=bx, add=yes, verify=no, update=yes, 
    > "dev$null")
  hedit ( work, "obs_by", value=by, add=yes, verify=no, update=yes, 
    > "dev$null")
  hedit ( work, "obs_se", value=se, add=yes, verify=no, update=yes, 
    > "dev$null")
  hedit ( work, "obs_sxi", value=sxi, add=yes, verify=no, update=yes, 
    > "dev$null")
  hedit ( work, "obs_syi", value=syi, add=yes, verify=no, update=yes, 
    > "dev$null")
  hedit ( work, "obs_gain", value=gain, add=yes, verify=no, update=yes, 
    > "dev$null")
  hedit ( work, "obs_seed", value=seed, add=yes, verify=no, update=yes, 
    > "dev$null")
  hedit ( work, "obs_psf", value=in, add=yes, verify=no, update=yes, 
    > "dev$null")

  delete( out, verify-, >& "dev$null" )
  unlearn mknoise
  if (be<0.0) {be=0.0}
  mknoise( work, output=out, gain=gain, rdnoise=ron, background=be, 
    poisson=yes, seed=seed)

  unlearn minmax
  minmax( images=out, force=yes, update=yes, verbose=verbose)
  hedit ( images=out, fields="obs_mxxy", value=minmax.maxpix, add=yes, 
    verify=no, update=yes, > "dev$null")
  hedit ( images=out, fields="obs_mxvl", value=minmax.maxval, add=yes, 
    verify=no, update=yes, > "dev$null")
  hedit ( images=out, fields="obs_mnxy", value=minmax.minpix, add=yes, 
    verify=no, update=yes, > "dev$null")
  hedit ( images=out, fields="obs_mnvl", value=minmax.minval, add=yes, 
    verify=no, update=yes, > "dev$null")
}

bye:
if (verbose) print  ( "PSFOBS: END" )
end
