#
# This TCL script walks through a CD-Rom full of images from the 
# TASS Mark IV camera, and extracts information from each image:
# 
#    image name
#    filter
#    Julian Date at start of exposure
#    exposure time (this is tricky, not in FITS header)
#    RA
#    Dec
# (added 2/7/2003)
#    sky
#    nstar
#    zeropoint
# (added 2/8/2003)
#    include
#
# It also assigns a "type" to each image: either "object" or "dark",
#    based solely on the position of the image within the night's run.
#    This works properly for Disk Set 7, but will need to be modified
#    for other disk sets.
#       --- MWR 5/20/2000
#
# Added Julian Date and airmass to items extracted per image.
#       --- MWR 6/20/2000
#
# Modified the algorithm used to guess a target's type: keep track
#    of number of images in each filter separately.  Still not
#    a robust solution, but works with Tom's data sets on disks 7, 15.
#       --- MWR 1/20/2001
#
# Removed "airmass" from the set of values printed to the output file,
#    since we don't need it until much later in the pipeline,
#    and we'll calculate it then.
#       --- MWR 1/20/2001
#
# We now use a parameter to describe the location of the external executable
#    program "picksym"
#       --- MWR 2/9/2001
#
# Instead of expecting "CRVAL1" for RA and "CRVAL2" for Dec, we now
#    look for "RA" for RA and "DEC" for Dec in the FITS header.
#    We also switch from using "UT" to "UTSTART".
#    Is a change from Disk 16 to Disk 17.
#       --- MWR 4/5/2001
#
# We now must source "convert.tcl" on startup, so that we can use
#    new routines to convert RA and Dec to/from Babylonian and 
#    decimal degrees.
#       --- MWR 4/6/2001
#    
# Added three new fields to the information for each image
#    sky        (sky value, in DN)
#    nstar      (number of stars found in the image)
#    zeropoint  (photometric zeropoint)
#       --- MWR 2/7/2003
#
# Added one new field to the information for each image
#    include    (include this image in processing?  1=yes, 0=no)
#       --- MWR 2/8/2003
#
# Added another new field
#    skysig     (measure of scatter in pixel values across frame)
#    Also, we now check to make sure that the number of fields in each
#    line of the "make_list.out" file matches the number of fields
#    specified in the "make_list.param" file.  It's too hard to verify
#    that we are actually writing the requested fields, though :-(
#       --- MWR 2/9/2003
#
# Added new fields "fwhm" and "colorterm" to the list of 
#    info in the "make_list.out" file
#       --- MWR 4/11/2003
#
# We now write a single comment line at the start of the 
#    "make_list.out" file, which gives the version number of the
#     pipeline used to reduce the data, and the current date and time.
#       --- MWR 4/12/2003
#
# Changed the single comment line at star of "make_list.out" file
#     into two lines (one for software packages, one for date and time).
#     Looks much nicer.
#       --- MWR 4/24/2003
#
# Fixed errors in "calc_jd()" and "days_since_1950()" which semi-cancelled
#     each other; the only errors occurred after Feb 28, 2004.  I was
#     determining leap year incorrectly.  Sigh.
#       --- MWR 6/18/2004
#
# Translate the IMAGETYP value into lowercase after reading it
#     from the FITS header.  All downstream code expects image types
#     like "dark" and "object".
#       --- MWR 3/27/2005



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

source param.tcl
source convert.tcl


###########################################################################
# Make a list of the images in the "input_dir".  The list includes the
#   name of each file, plus information gleaned from the FITS header:
#   date, time, RA, Dec, exposure time, filter.
#
# Also, add placeholders (zeroes) for other information we'll determine later:
#       sky, nstar, zeropoint, include
#   
#
#      param_file         name of a parameter file 
#      
# RETURNS
#    0           if all goes well
#    1           if there's a problem
#    
proc make_list { param_file } {
  global debug

  if { $debug > 0 } {
    puts "entering make_list"
  }

  # get information from the param file
  if { [file exists $param_file] != 1 } {
    puts stderr "make_list: can't open param file $param_file"
    return 1
  }
  set dir [get_param $param_file "input_dir"]
  set input_regexp [get_param $param_file "input_regexp"]
  set output_dir [get_param $param_file "output_dir"]
  set outfile [get_param $param_file "imagelist_name"]
  set picksym_exe [get_param $param_file "picksym_exe"]

  # We need to find out the directory in which the "picksym" program
  #   lives, and prepend the directory name to the executable's name
  set bait_dir [get_directory "bait_dir"]
  set picksym_exe $bait_dir/$picksym_exe


  if { [string compare $outfile "NULL"] != 0 } {
    set fileid [open $output_dir/$outfile "w"]
  } else {
    set fileid stdout
  }
  if { $debug > 1 } {
    puts " make_list: fileid is $fileid "
  }

  # write a pair of "header" lines at the start of the file, giving the
  #   version number of all software packages and the current date.
  set version_list [get_version_list]
  set datestr [exec date]
  set line [format "#  software %s  " $version_list]
  puts $fileid $line
  set line [format "#  run on %s " $datestr]
  puts $fileid $line
  
  set count 0
  set old_filter ""
  set lis [glob $dir/$input_regexp]
  foreach file $lis {
  
    set ll [split $file / ]
    set len [llength $ll]
    set image [lindex $ll [expr $len-1]]
    if { $debug > 1 } {
      puts "next file is $image "
    }
    
    if { [catch { set filter [exec $picksym_exe file=$file FILTER] } ] == 1 } {
      puts stderr "file $file has no FILTER in FITS header -- skipping"
      continue
    }
    # get rid of trailing spaces in the filter name
    set filter [string trimright $filter]
    if { $debug > 1 } {
      puts "filter is ..$filter.."
    }


    # these used to check for CRVAL1 (RA) and CRVAL2 (Dec)
    if { [catch { set ra [exec $picksym_exe file=$file RA] } ] == 1 } {
      puts stderr "file $file has no RA in FITS header -- skipping"
      continue
    }
    # convert from Babylonian HH:MM:SS to decimal degrees, if necessary
    if { [is_baby $ra] == 1 } {
      if { [rahms_to_radeg $ra ra_deg] != 0 } {
        puts stderr "file $file RA value $ra fails conversion -- skipping"
        continue
      }
      set ra $ra_deg
    }
    if { $debug > 1 } {
      puts "RA is $ra"
    }

    if { [catch { set dec [exec $picksym_exe file=$file DEC] } ] == 1 } {
      puts stderr "file $file has no DEC in FITS header -- skipping"
      continue
    }
    # convert from Babylonian HH:MM:SS to decimal degrees, if necessary
    if { [is_baby $dec] == 1 } {
      if { [hms_to_hhh $dec dec_deg] != 0 } {
        puts stderr "file $file Dec value $dec fails conversion -- skipping"
        continue
      }
      set dec $dec_deg
    }
    if { $debug > 1 } {
      puts "Dec is $dec"
    }
 
    # we use the DATE-OBS and UTSTART values to calculate Julian Date
    #   (used to use UT, now use UTSTART)
    if { [catch { set dateobs [exec $picksym_exe file=$file DATE-OBS] } ] == 1 } {
      puts stderr "file $file doesn't contain DATE-OBS in FITS header -- skipping"
      continue
    }
    if { [catch { set ut [exec $picksym_exe file=$file UTSTART] } ] == 1 } {
      puts stderr "file $file doesn't contain UTSTART in FITS header -- skipping"
      continue
    }
    if { [calc_jd $dateobs $ut julian_date] != 0 } {
      puts stderr "make_list: calc_jd fails for file $file -- skipping "
      continue
    }

    # figure out the exposure time of the image, from the EXPTIME value
    if { [catch { set exptime [exec $picksym_exe file=$file EXPTIME] } ] == 1 } {
      puts stderr "file $file doesn't contain EXPTIME in FITS header -- skipping"
      continue
    }
    if { $debug > 1 } {
      puts "exptime is $exptime"
    }

    # figure out the type of image (dark or target), from the IMAGETYP
    #   keyword
    if { [catch { set type [exec $picksym_exe file=$file IMAGETYP] } ] == 1 } {
      puts stderr "file $file doesn't contain IMAGETYP in FITS header -- skipping"
      continue
    }
    # convert the IMAGETYP value into lowercase letters.  This will avoid
    #   problems later on.  The code downstream expects to see lower-case
    #   strings, such as "dark" or "object".  Rob Creager's Mark IV
    #   control program, on the other hand, places into the FITS header
    #   strings like "DARK" and "OBJECT"
    set type [string tolower $type]

    # this is where we put out the single line per image file
    set sky 0
    set skysig 0
    set fwhm 0
    set nstar 0
    set zeropoint 99.0
    set colorterm 99.0
    set include 1
    set line [format "%-20s %2s %13.5f %6.1f %8s %6.1f %6.1f %6d %7.1f %5.2f %6d %7.3f %7.3f %1d" \
                           $image $filter $julian_date $exptime \
                           $type $ra $dec $sky $skysig $fwhm $nstar \
                           $zeropoint $colorterm $include]
    # sanity check -- make sure the number of items in this line matches
    #   the number of fields specified in the "make_list.param" file
    set line_count [llength $line]
    set expect_count [count_list_fields]
    if { $line_count != $expect_count } {
      puts stderr "make_list: field mismatch: actual $line_count, expect $expect_count"
      puts stderr "make_list:   perhaps code in 'make_list' proc needs to "
      puts stderr "make_list:   be modified to match 'make_list.param' values"
      return 1
    }
    puts $fileid $line
    incr count
  }

  if { [string compare $fileid stdout] == 0 } {
    if { $debug > 1 } {
      puts "this was using stdout"
    }
  } else {
    if { $debug > 1 } {
      puts "nope, not using stdout, using $fileid"
    }
    close $fileid
  }


  return 0
}


############################################################################
# Given a date (assumed UT) and a Universal Time, calculate the 
#   Julian Date:
#   
#      calc_jd          dateobs  ut  julian_date
#      
# where
#            dateobs          is a string in this format: 2000-04-25
#                                 assumed to be the UT date
#                                 
#            ut               is a string in this format: 07:49:52.00
#                                 giving the UT at start of exposure  
#                
#            julian_date      will be calculated from the other values,
#                                 and set on output
#                                 
# RETURNS:
#    0      if all goes well
#    1      if there's a problem
# 
proc calc_jd { dateobs ut julian_date } {
  upvar $julian_date jd
  global debug

  if { $debug > 0 } {
    puts "entering calc_jd "
  }

  # this is the Julian Date of Jan 1.0, 1950
  # fixed error here 6/18/2004
  set jd_base 2433281.5  


  # first, we decipher the DATEOBS value into year, month, day
  if { [scan $dateobs "%4d-%2d-%2d" year month day] != 3 } {
    puts stderr "calc_jd: invalid DATEOBS value $dateobs"
    return 1
  }

  # next, we decipher the UT value into hour, minute, second
  if { [scan $ut "%2d:%2d:%f" hour minute second] != 3 } {
    puts stderr "calc_jd: invalid UT value $ut"
    return 1
  }

  # calculate the Julian Date
  #   we need to calculate the number of full days since 1950,
  #   and the fraction of a day since 00:00 UT on the given day
  if { $debug > 1 } {
    puts [format "    Y M D %4d %2d %2d  H M S %2d %2d %5.2f " \
                 $year $month $day $hour $minute $second]
  }
  set fractional_day [expr $hour+($minute/60.0)+($second/3600.0)]
  set fractional_day [expr $fractional_day/24.0]
  if { $debug > 1 } {
    puts [format "   fractional_day is %7.4f " $fractional_day]
  }
  set since_1950 [days_since_1950 $year $month $day]
  if { $debug > 1 } { 
    puts [format "   days since 1950 is %8d " $since_1950]
  }

  # okay, we know how many days since 1950.  So, we can add the JD of 
  #   Jan 1, 1950, to the number of days since 1950, and add the fractional
  #   day amount, to find the JD
  set jd [expr $jd_base + $since_1950 + $fractional_day]

  if { $debug > 1 } {
    puts [format "   JD is %14.6f " $jd]
  }

  return 0
}


##############################################################################
# Given a year, a month, a day, calculate the number of days which have
#   passed between Jan 1, 1950, and the given date.
#
# Assumes that the given date is after Jan 1, 1950.
# 
# RETURNS:
#    number of days since Jan 1, 1950
#
proc days_since_1950 { year month day } {
  global debug 

  if { $debug > 0 } {
    puts "entering days_since_1950"
  }

  set days_in_months [list 31  28  31  30  31  30  31  31  30  31  30  31]

  set days 0
  set yy 1950
  while { $yy < $year } {
    set days [expr $days + 365]
    # fixed error here 6/18/2004
    if { [expr $yy % 4] == 0 && \
             ( ([expr $yy % 100] != 0) || ([expr $yy % 400] == 0) ) } {
      set days [expr $days + 1]
    }
    incr yy
  }
  set mm 1
  while { $mm < $month } {
    set days [expr $days + [lindex $days_in_months [expr $mm-1]]]
    # fixed error here, too, 6/18/2004 
    if { $mm == 2 && [expr $yy % 4] == 0 && \
               ( ([expr $yy % 100] != 0) || ([expr $yy % 400] == 0) ) } {
      set days [expr $days + 1]
    }
    incr mm
  }

  set days [expr $days + $day]

  return $days
}