procedure getseeing(image,zlo,zhi)

string	image	=""	{prompt="name of input image"}
real	zlo	=INDEF	{prompt="lower cut level"}
real	zhi	=INDEF	{prompt="upper cut level"}
string	logfile	=""	{prompt="text log file"}
bool	loghead =no	{prompt="log result to the image header?"}
real	fwhm	=0.	{prompt="measured seeing value"}
real	efwhm	=0.	{prompt="error on seeing value (if Ntotal>1)"}
int	nused	=0	{prompt="number of measurements used (Ntotal-Nreject)"}
bool	redispl =yes	{prompt="(re)display image first?"}
bool	verbose =yes	{prompt="show messages and results?"}

# @(#) task getseeing                Author: R.A. Jansen -- June 20 1996
# @(#)
# @(#) task to measure the effective seeing FWHM using several stars in
# @(#) an image (marked interactively by user).  If the number of stars
# @(#) is larger than 3, an estimate of the error on the determination
# @(#) of the FWHM is calculated as well. The result is returned as task
# @(#) parameters, and written to a text logfile if defined.  Optionally
# @(#) the result is also written to image header keyword 'SEEING'.
# @(#)
# @(#) Jun.21 2000 -- Added support for other image formats. (R.A.J.)

begin

	string	img, tmp, logstr
	real	seeing, seeingold, esum, eseeing, eseeingold
        real	rmin, rmax
	int	NL, NN, NS, imglen, NREP

	img	= image
        seeing	= 0.
        eseeing = 0.

#	Check whether the image exists...
	unlearn chkimg
	chkimg (img, "access", verbose=yes)
	if (!chkimg.ok) { return }
	img = chkimg.imroot//chkimg.imtype//chkimg.imrpstr

#       Print opening messag...
        if (verbose) { printf("\n GETSEEING: Image="//img) }

#       Display image...
	if (redispl) {
	   if (zlo != 0. || zhi != 0.) {
	      xdisplay (img, z1=zlo, z2=zhi)
	   } else { xdisplay (img) }
	}

#       Mark all stars...
        print("  Mark all non-saturated stars using the 'a' key:"
        imexamine (img, 1, img, logfile="tmpgtsng.log", keeplog=yes, 
                 defkey="a", autoredraw=yes, allframes=yes, nframes=0, 
                 ncstat=5, nlstat=5, graphcur="", imagecur="", wcs="logical", 
                 xformat="", yformat="", graphics="stdgraph", 
                 display="display(image='$1',frame=$2)", use_display=yes)

#       Calculate seeing FWHM and error thereon...
        fields ("tmpgtsng.log", "1,2,15", lines="3-", quit_if_miss=yes,
                print_file_n=no, >> "tmpgtsng.fwhm")
        iwc ("tmpgtsng.fwhm")
        NL = iwc.nlines

        if ( NL == 0 ) {
           delete ("tmpgtsng.log,tmpgtsng.fwhm", yes, verify=no)
           error (1, "No stars marked!")
        } else {
           if (verbose) {
	      print ("  Calculating mean FWHM and error thereon ...")
	   }
           NN = 1
           NS = 0
           rmin =  999999.9999
           rmax = -999999.9999
           while ( NN <= NL ) {
              rdlist ("tmpgtsng.fwhm", NN)
              cparse (rdlist.cline, delim=" ")
              seeing = seeing + real(cparse.field3)
              if ( real(cparse.field3) > rmax ) { rmax = real(cparse.field3) }
              if ( real(cparse.field3) < rmin ) { rmin = real(cparse.field3) }
              NN += 1
           }
           if ( NL >= 5 ) { 
#             Exclude minimum and maximum value from sum, before taking mean
              seeing = (seeing - rmin - rmax)/real(NL-2)
           } else {
#             Not enough measurements to exclude minimum and maximum value
              seeing = seeing / real(NL)
           }
           fwhm = real(int(10000.0*seeing)/10000.)
           nused = NL

           if ( NL >= 3 ) {
#             Calculate error on mean FWHM
              esum = 0.00
              NN = 1
              while ( NN <= NL ) {
                 rdlist ("tmpgtsng.fwhm", NN)
                 cparse (rdlist.cline, delim=" ")
                 esum = esum + (seeing-real(cparse.field3))*(seeing-real(cparse.field3))
                 NN += 1
              }
              if ( NL >= 5 ) {
#                Exclude deviation^2 corresponding to min and max value
                 esum = esum - (seeing-rmax)*(seeing-rmax)
                 esum = esum - (seeing-rmin)*(seeing-rmin)
                 eseeing = sqrt( esum/real(NL-3) )
              } else {
#                Not enough measurements to exclude deviation^2 of min and max
                 eseeing = sqrt( esum/real(NL-1) )
              }

              if ( NL >= 5 ) {
#               Reject outlying points using mean deviation off the mean, in
#               two rejection cycles
                for ( NREP=1 ; NREP <= 2 ; NREP=NREP+1 ) {
                   esum = 0.00
                   seeingold = seeing
                   eseeingold = eseeing
                   seeing = 0.00
                   NN = 1
                   NS = 0
                   while ( NN <= NL ) {
                      rdlist ("tmpgtsng.fwhm", NN)
                      cparse (rdlist.cline, delim=" ")
                      if ( abs(seeingold-real(cparse.field3)) < 1.75*eseeingold ) {
                         seeing = seeing + real(cparse.field3)
                         NS += 1
                      }
                      NN += 1
                   }
                   seeing = seeing / real(NS)

                   esum = 0.00
                   NN = 1
                   NS = 0
                   while ( NN <= NL ) {
                      rdlist ("tmpgtsng.fwhm", NN)
                      cparse (rdlist.cline, delim=" ")
                      if ( abs(seeing-real(cparse.field3)) < 1.75*eseeing ) {
                         esum = esum + (seeing-real(cparse.field3))*(seeing-real(cparse.field3))
                         NS += 1
                      }
                      else if ( NREP == 2 ) {
                         print("    rejected star #"//str(NN)//": "//rdlist.cline)
                      }
                      NN += 1
                   }
                   eseeing = sqrt( esum/real(NS-1) )
                   nused = NS
                }

                if (verbose) {
		   print("    measured <FWHM>="//substr(str(seeing),1,8)//" +-"//substr(str(eseeing),1,8)//"  Nstar="//str(NL)//"  Nrejected="//str(NL-NS))
		}
                logstr = " Imexamine: FWHM="//substr(str(seeing),1,6)
                logstr = logstr//" +-"//substr(str(eseeing),1,6)//"  Nstar="//str(NL)
                logstr = logstr//" Nrej="//str(NL-NS)

              } else {

#                Not enough measurements to allow rejection of deviant values
                 if (verbose) {
		    print("    measured <FWHM>="//substr(str(seeing),1,8)//" +-"//substr(str(eseeing),1,8)//"  Nstar="//str(NL)//"  Nrejected=0")
		 }
                 logstr = " Imexamine: FWHM="//substr(str(seeing),1,6)
                 logstr = logstr//" +-"//substr(str(eseeing),1,6)//"  Nstar="//str(NL)
                 logstr = logstr//" Nrej=0"
              }

           } else if ( NL == 2 ) {
#             Not enough measurements for a formal error analysis, but
#             the deviation of the mean value may give some handle on the
#             expected true error
              eseeing = 0.00
              for ( NN=1 ; NN<=2 ; NN=NN+1 ) {
                 rdlist ("tmpgtsng.fwhm", NN)
                 cparse (rdlist.cline, delim=" ")
                 eseeing = eseeing + abs(seeing-real(cparse.field3))
              }
              eseeing = eseeing / 2.0

              if (verbose) {
		 print("    measured <FWHM>="//substr(str(seeing),1,8)//" +-"//substr(str(eseeing),1,8)//"  Nstar="//str(NL)//"  Nrejected=0")
	      }
              logstr = " Imexamine: FWHM="//substr(str(seeing),1,6)
              logstr = logstr//" +-"//substr(str(eseeing),1,6)//"  Nstar="//str(NL)
              logstr = logstr//" Nrej=0"

           } else {
#            only one star available, so no error can be given
             if (verbose) {
		print("    measured <FWHM>="//substr(str(seeing),1,8)//" +-0."//"  Nstar="//str(NL)//"  Nrejected=0"
	     }
             logstr = " Imexamine: FWHM="//substr(str(seeing),1,6)//" +-0."
             logstr = logstr//"  Nstar="//str(NL)//" Nrej=0"
           }

#          Output mean FWHM and error by resetting the task parameters
           fwhm = real(int(10000.*seeing)/10000.)
           efwhm = real(int(10000.*eseeing)/10000.)
        }


	gdate()
	logstr = gdate.date//logstr

	if (loghead) {
#          Log results to the image header. Make sure 'imred.ccdred' is
#	   loaded (task "ccdhedit")...
	   if ( !defpac("ccdred") ) {
	      imred
	      ccdred
	   }
           if (verbose) { print("  Writing result to image header...") }
           ccdhedit (img, "seeing", logstr, type="string")
	}

        if ( logfile != "" ) {
#	   Log results to the specified text logfile...
           if (verbose) {
	      print("  Writing result to text file '"//logfile//"'...")
	   }
           print (img//": "//logstr, >> logfile)
        }


#       Clean-up ...
        delete ("tmpgtsng.log", yes, verify=no)
        delete ("tmpgtsng.fwhm", yes, verify=no)

        if (verbose) { print ("GETSEEING: Finished.") }

end