# This script deals with the sky level in cleaned images (those which
#   have had dark subtracted and flatfield divided).
#
#        - measures a mean sky level and skysig, and records them in
#               the "make_list.out" file for each frame
#
#        - checks to see if the sky value and sigma fall within the
#               acceptable range; if not, set its "include" flag to zero
#
#        - (if desired) subtracts a model of the background
#               from the image, leaving a mean background of zero
#
#   --  MWR 2/08/2003
#
# If the user decides to subtract the sky background, measure the 
#   "sky" value _before_ subtracting ... but measure "skysig" value
#   _after_ subtracting the model.  If there are large-scale gradients
#   in the background, they would contribute to the "skysig" value;
#   thus making it hard to distinguish a cloudy frame from a clear one.
# 
#   --  MWR 3/2/2003
#
# Added the "min" and "max" parameters, which will be used by the
#   "sky" and "back" programs.  This is a response to the bogus
#   results caused by very saturated stars.
#
#   --  MWR 3/15/2003
#
# Added namespace "Sky" and procedures "init_skyparams" and "get_skyparams",
#   so that we can share one set of parameter values (in sky.param)
#   cleanly with other states of the pipeline (such as the "stars" step).
#
#   --  MWR 5/14/2003


if { [catch { set debug } ] == 1 } {
  set debug 2
}

# contains source for the "get_param" proc
source param.tcl



#############################################################################
# Given a param file name, and a cleaned image, 
#   parse the param file for a few control values and then
#   do the sky tasks ...
#
# We assume files are in directories
#
#       output_dir           the corrected image files
#       output_dir           the .coo and .pht files these routines create
#
# Returns
#    0              if all goes well
#    1              if there's an error
#
proc sky { param_file clean_image } {
  global debug


  if { $debug > 0 } {
    puts "entering sky ..."
  }

  if { [file exists $param_file] != 1 } {
    puts stderr "sky: can't find param file $param_file"
    return 1
  }

  # the XVista programs can be found in a special directory
  set xvista_dir [get_directory "xvista_dir"]
  if { $xvista_dir == {} } {
    puts stderr "sky: get_directory fails for xvista_dir"
    return 1
  }

  set output_dir       [get_param $param_file "output_dir"]

  # parameters for calculating the sky value and sigma
if { 0 == 1 } {
  set skybin           [get_param $param_file "skybin"]
  set min_accept       [get_param $param_file "min_accept"]
  set max_accept       [get_param $param_file "max_accept"]
}
  if { [get_skyparams skybin min_accept max_accept] != 0 } {
    puts stderr "sky: get_skyparams returns with error"
  }

  # parameters for accepting the image, or marking it as "bad"
  set max_sky          [get_param $param_file "max_sky"]
  set max_skysig       [get_param $param_file "max_skysig"]

  # parameters for fitting and subtracting the background
  set subtract_sky     [get_param $param_file "subtract_sky"]
  set ngrid            [get_param $param_file "ngrid"]
  set order            [get_param $param_file "order"]
  set map              [get_param $param_file "map"]
  set verbose          [get_param $param_file "verbose"]


  if { $debug > 1 } {
    puts "sky: about to work on file $clean_image"
  }

  # get information about the image
  if { [single_image_info $clean_image image_info_list] != 0 } {
    puts stderr "sky: single_image_info fails for image $clean_image"
   return 1
  }
  # if this is not an "object" image, don't process it!  Just return
  #   with code 0, but do print a warning message
  set type [get_image_value $image_info_list "type"]
  if { [string compare $type "object"] != 0 } {
    if { $debug > 0 } {
      puts "sky: image $clean_image is not 'object', so skip it"
    }
    return 0
  }

  # find the sky and skysig value
  if { [find_sky $output_dir $xvista_dir $clean_image $skybin \
                         $min_accept $max_accept  skyval skysigval ] != 0 } {
    puts stderr "sky: find_sky fails on file $clean_image"
    puts stderr "sky: quitting now"
    return 1
  }
  # write the sky value into the image list file
  set klist [list sky $skyval ]
  if { [replace_image_value $clean_image $klist ] != 0 } {
      puts stderr "sky: replace_image_value fails for $clean_image, sky $skyval"
      return 1
  }
  # check the image's sky value against the maximum acceptable
  #   value (for its filter)
  if { [check_sky $clean_image $skyval $max_sky ] != 0 } {
    puts stderr "sky: check_sky fails on file $clean_image"
    puts stderr "sky: quitting now"
    return 1
  }

  # if we are NOT going to subtract a background model,
  #   write the "skysig" into the image list file and
  #   check its value against the acceptable limit
  if { $subtract_sky == 0 } {
    set klist [list skysig $skysigval ]
    if { [replace_image_value $clean_image $klist ] != 0 } {
        puts stderr "sky: replace_image_value fails for $clean_image, skysig $skysigval"
        return 1
    }
    # check the image's skysig value against the maximum acceptable
    if { [check_skysig $clean_image $skysigval $max_skysig ] != 0 } {
      puts stderr "sky: check_skysig fails on file $clean_image"
      puts stderr "sky: quitting now"
      return 1
    }
  }


  # if desired, fit a simple model to the background of the image
  #   and subtract it from the image
  if { $subtract_sky > 0 } {

    if { [fit_and_subtract_sky $output_dir $xvista_dir $clean_image  \
                               $ngrid $order $min_accept $max_accept \
                               $skybin $map $verbose ] != 0 } {
      puts stderr "sky: fit_and_subtract_sky fails on file $clean_image"
      puts stderr "sky: quitting now "
      return 1
    }

    # Now, make a check against max skysig _after_ subtracting the 
    #   background (now that large-scale gradients are gone).
    #   This means we have to run the "find_sky" routine again.
    if { [find_sky $output_dir $xvista_dir $clean_image $skybin \
                         $min_accept $max_accept skyval skysigval ] != 0 } {
      puts stderr "sky: find_sky fails on file $clean_image"
      puts stderr "sky: quitting now"
      return 1
    }
    set klist [list skysig $skysigval ]
    if { [replace_image_value $clean_image $klist ] != 0 } {
        puts stderr "sky: replace_image_value fails for $clean_image, skysig $skysigval"
        return 1
    }
    # check the image's skysig value against the maximum acceptable
    if { [check_skysig $clean_image $skysigval $max_skysig ] != 0 } {
      puts stderr "sky: check_skysig fails on file $clean_image"
      puts stderr "sky: quitting now"
      return 1
    }
  }

  return 0
}


#############################################################################
# PROCEDURE: find_sky
#
# Prepare to run the "sky" program on the given image, then do so.
# 
#   - find the sky value and sky sigma
#         (place these into the output args "p_skyval" and "p_skysigval")
#
# We assume
#
#          output_dir            holds the image
#          xvista_dir            holds the XVista programs
#
# Return 
#    0               if all goes well
#    1               if there's a problem
#    
proc find_sky { output_dir xvista_dir clean_image skybin  \
                        min_accept max_accept p_skyval p_skysig } {
  upvar $p_skyval skyval $p_skysig skysig
  global debug

  if { $debug > 1 } {
    puts "entering find_sky"
  }

  set infile $output_dir/$clean_image
  if { [file exists $infile] != 1 } {
    puts stderr "find_sky: can't find file $infile"
    return 1
  }

  # find the sky value and sky sigma value
  set cmd "$xvista_dir/sky $infile bin=$skybin min=$min_accept max=$max_accept"
  if { $debug > 1 } {
    puts "find_sky: cmd is ..$cmd.."
  }
  if { [catch { set str [eval exec $cmd] } ] != 0 } {
    puts stderr "find_sky: sky returns with error"
    return 1
  }
  set ll [lindex $str 0]
  set ll [split $ll "="]
  set skyval [lindex $ll 1]
  set ll [lindex $str 1]
  set ll [split $ll "="]
  set skysig [lindex $ll 1]

  if { $debug > 1 } {
    puts "find_sky: file $clean_image has sky $skyval  skysig $skysig"
  }

  # some sanity checks (just warnings)
  if { $skyval < 1 } {
    puts stderr "find_sky: sky value for $clean_image is $skyval, less than one"
  }
  if { $skysig < 1 } {
    puts stderr "find_sky: skysig value for $clean_image is $skysig, less than one"
  }

  return 0
}



#############################################################################
# PROCEDURE: check_sky
#
# Compare the sky value for an image against the maximum
#   permitted value.  If the value exceeds the boundaries, set the
#   image's "include" flag to "0", so that it won't be used in further
#   processing.
# 
# We assume
#
#          output_dir            holds the image
#
#   
# Return 
#    0               if all goes well
#    1               if there's a problem
#    
proc check_sky { clean_image skyval max_sky } {
  global debug

  if { $debug > 1 } {
    puts "entering check_sky"
  }

  # get information about this image ...
  if { [single_image_info $clean_image image_info_list] != 0 } {
    puts stderr "check_sky: single_image_info fails for image $clean_image"
    return 1
  }

  # pick out the filter
  set filter [get_image_value $image_info_list "filter"]
  if { $debug > 3 } {
    puts "check_sky: filter is $filter"
  }
  
  # get the maximum acceptable values for sky for this filter
  set this_max_sky    [get_image_value $max_sky $filter]
  if { $debug > 3 } {
    puts "check_sky: max value in filter $filter is $this_max_sky"
  }

  # now, compare the actual value to the max permitted value.
  set okay 1
  if { $skyval > $this_max_sky } {
    if { $debug > 0 } {
      puts "check_sky: image $clean_image has sky $skyval > $this_max_sky"
      puts "check_sky:    designating $clean_image as bad"
    }
    set okay 0
  }
  if { $okay != 1 } {
    if { [designate_bad $clean_image] != 0 } {
      puts stderr "check_sky: designate_bad fails on $clean_image"
      return 1
    }
  }

  return 0
}


#############################################################################
# PROCEDURE: check_skysig
#
# Compare the skysig value for an image against the maximum
#   permitted value.  If the value exceeds the boundaries, set the
#   image's "include" flag to "0", so that it won't be used in further
#   processing.
# 
# We assume
#
#          output_dir            holds the image
#
#   
# Return 
#    0               if all goes well
#    1               if there's a problem
#    
proc check_skysig { clean_image skysig max_skysig } {
  global debug

  if { $debug > 1 } {
    puts "entering check_skysig"
  }

  # get information about this image ...
  if { [single_image_info $clean_image image_info_list] != 0 } {
    puts stderr "check_sky: single_image_info fails for image $clean_image"
    return 1
  }

  # pick out the filter
  set filter [get_image_value $image_info_list "filter"]
  if { $debug > 3 } {
    puts "check_sky: filter is $filter"
  }
  
  # get the maximum acceptable values for sky and skysig for this filter
  set this_max_skysig [get_image_value $max_skysig $filter]
  if { $debug > 3 } {
    puts "check_skysig: max value in filter $filter is $this_max_skysig"
  }

  # now, compare the actual value to the max permitted value.
  set okay 1
  if { $skysig > $this_max_skysig } {
    if { $debug > 0 } {
      puts "check_skysig: image $clean_image has skysig $skysig > $this_max_skysig"
      puts "check_skysig:    designating $clean_image as bad"
    }
    set okay 0
  }
  if { $okay != 1 } {
    if { [designate_bad $clean_image] != 0 } {
      puts stderr "check_skysig: designate_bad fails on $clean_image"
      return 1
    }
  }

  return 0
}



############################################################################
# PROCEDURE: fit_and_subtract_sky
#
# Use the XVista "back" program to fit a simple model to the background
#   of the given image, and then subtract that model away.  The result
#   should be a relatively flat background value near zero.
#
#
# RETURNS:
#   0          if all goes well
#   1          if an error occurs
#
proc fit_and_subtract_sky { output_dir xvista_dir clean_image \
                                       ngrid order min_accept max_accept \
                                       skybin map verbose } {
  global debug


  if { $debug > 1 } {
    puts "entering fit_and_subtract_sky"
  }

  set imagefile $output_dir/$clean_image
  if { [file exists $imagefile] != 1 } {
    puts stderr "fit_and_subtract_sky: can't find file $imagefile"
    return 1
  }

  # prepare to run the "back" program
  set cmd "exec $xvista_dir/back $imagefile "
  set cmd "$cmd ngrid=$ngrid order=$order "
  set cmd "$cmd min=$min_accept max=$max_accept bin=$skybin "
  if { $map != {} } {
    set cmd "$cmd map=$output_dir/$map "
  }
  if { $verbose > 0 } {
    set cmd "$cmd verbose "
  }
  if { $debug > 1 } {
    puts "fit_and_subtract_sky: cmd is ..$cmd.."
  }
  if { [catch { set ret [eval $cmd] } ] != 0 } {
    puts stderr "fit_and_subtract_sky: back returns with error for image $imagefile"
    return 1
  }

  if { $debug > 1 } {
    puts "fit_and_subtract_sky: here comes results of background model"
    puts $ret
  }

  return 0
}



############################################################################
# We need a private namespace to store/retrieve the values of a few
#   parameters used in measuring sky values.  The parameters are needed
#   in the "sky" and "stars" steps of the pipeline, and rather than have
#   two separate sets of parameters (which might not both be updated),
#   we use these routines.
#


namespace eval Sky {
  namespace export init_skyparams get_skyparams
  variable Sky_skybin {}
  variable Sky_min_accept {}
  variable Sky_max_accept {}

#########################################################################
# PROCEDURE: init_skyparams
#
# DESCRIPTION: Read the values of sky-related parameters from the "sky.param"
#              file, and save them in our private variables for later use.
#
# Args:
#     param_file        name of file containing the sky-related parameters
#
# RETURNS
#   0             if all goes well
#   1             if an error occurs
# 
proc init_skyparams { param_file } {
  global debug
  variable Sky_skybin
  variable Sky_min_accept
  variable Sky_max_accept

  # get information from the param file
  if { [file exists $param_file] != 1 } {
    puts stderr "init_skyparams: can't open param file $param_file"
    return 1
  }
  set Sky_skybin [get_param $param_file "skybin"]
  set Sky_min_accept [get_param $param_file "min_accept"]
  set Sky_max_accept [get_param $param_file "max_accept"]

  return 0
}


#########################################################################
# PROCEDURE: get_skyparams
#
# DESCRIPTION: Return the values of the sky-related parameters
#              as the arguments.
#
# USAGE:
#            get_skyparams  sky_bin min_accept max_accept
#
#              sky_bin      (output)      will have the binsize used to create
#                                            a histogram of sky values
#
#              min_accept   (output)     mininum pixel value to include in 
#                                            calculations of the sky value
#
#              max_accept   (output)     maximum pixel value to include in 
#                                            calculations of the sky value
#
# RETURNS
#   0             if all goes well
#   1             if an error occurs
# 
proc get_skyparams { o_sky_bin o_min_accept o_max_accept } {
  global debug
  upvar $o_sky_bin sky_bin
  upvar $o_min_accept min_accept
  upvar $o_max_accept max_accept
  variable Sky_skybin
  variable Sky_min_accept
  variable Sky_max_accept

  set sky_bin $Sky_skybin
  set min_accept $Sky_min_accept
  set max_accept $Sky_max_accept

  return 0
}


}
# end of namespace "Sky"