#!/home/sybase/perl/perl/bin/perl -w
#
# tableplot:  takes a list of directories containing catalogues
#             and plots assorted params and makes a dqc log.
#             Optional HTML output.
#
# tableplot_ir:  IR version.
#

use strict;
use Getopt::Long;
use Math::Trig;
use POSIX qw(log10);
use CFITSIO qw(:constants :shortnames);
use PGPLOT;
use CASU::CvtUnit ':all';

### Constants ###

# Output file names.
my $SUMMARY_FILE = "summary.list";
my $SEEING_FILE = "seeing.hist";
my $PERCENT_FILE = "sky.percentages";

my $SUMPLOT_FILE = "summary.ps/CPS";
my $ELLPLOT_FILE = "ellipt.ps/CPS";
my $SEEPLOT_FILE = "seeing.ps/PS";
my $ASTPLOT_FILE = "astrom.ps/CPS";
my $DETPLOT_FILE = "detect.ps/CPS";
my $SITEPLOT_FILE = "site.ps/CPS";

my $SUMPLOT1_PNG = "summary_1.png/PNG";
my $SUMPLOT2_PNG = "summary_2.png/PNG";
my $ELLPLOT_PNG = "ellipt.png/PNG";
my $SEEPLOT_PNG = "seeing.png/PNG";
my $ASTPLOT1_PNG = "astrom_1.png/PNG";
my $ASTPLOT2_PNG = "astrom_2.png/PNG";
my $DETPLOT_PNG = "detect.png/PNG";
my $SITEPLOT1_PNG = "site_1.png/PNG";
my $SITEPLOT2_PNG = "site_2.png/PNG";

# PGPLOT `paper size' for PNG files.
my $PAPSIZE   = 8.0;
my $PAPASPECT = 0.8;

# Detector system parameters.
#
#     Pattern   pixsiz  gain mccd shortname  longname
my @detector_params = (
    [ "ALADDIN", 0.454,  2.0, 4,  "WFCAM",   "WFCAM"                 ]
);

# Position of the optical axis in CCD-space.
my @oaxis_x = (
    [    0.0  ,    0.0 ,    0.0 ,    0.0  ]   # WFCAM
);

my @oaxis_y = (
    [    0.0  ,    0.0 ,    0.0 ,    0.0  ]   # WFCAM
);

# Number of filters.
my $NFILT = 3;

# Map from filter names to 'icol' values.
my @filt_icol = (
    [ 'J '         , 2  ],
    [ 'H '         , 3  ],
    [ 'K '         , 4  ],
);

# List of default ccdzero values.
my $DEFAULT_CCDZERO = 30.0;

my @ccdzero_list = (
    [ 'J',           21.04 ],  # Warning: dummy values!
    [ 'H',           21.04 ],
    [ 'K',           21.04 ]
);

# List of filter names for 'icol' values.
my @filter_names = (
    [ 'dummy', 'dummy'       ], # 0
    [ 'dummy', 'dummy'       ], # 1
    [ 'J    ', 'J'           ], # 2
    [ 'H    ', 'H'           ], # 3
    [ 'K    ', 'K'           ]  # 4
);

### Main ###

# Initialise globals.
my $last_datej = 0.0;

my $shortname = undef;
my $longname = undef;

my $mjdmin = undef;
my $mjdmax = undef;
my $seemin = undef;
my $seemax = undef;

my $pixsiz = undef;
my $gaincorr = undef;
my $mccd = undef;

my @channel = ();
my @ndetect = ();
my @modjd = ();
my @seelev = ();
my @skylev = ();
my @filtnm = ();
my @exptim = ();
my @noiselev = ();
my @ellipt = ();
my @apcore = ();
my @rcore = ();
my @dateobs = ();
my @objnm = ();

my @astrms = ();
my @astnum = ();
my @perror = ();

my @humid = ();
my @tau = ();
my @winddir = ();
my @windspd = ();

my @warn_list = ();

# Always generate full-size postscript files (emulates behaviour of old
# Fortran version).
$ENV{PGPLOT_PS_BBOX} = "MAX";

# Extract the program name from $0.
my $progname = basename($0);

# Extract command-line arguments.
my $htmlmode = 1;
my $verbose = 0;

Getopt::Long::Configure('bundling');

my $rv = GetOptions('n|nohtml'  => sub { $htmlmode = 0; },
		    'v|verbose' => \$verbose);

die "Usage:\t$0 [-nv] dir [...] outdir\n" if(!$rv || @ARGV < 2);

my $outbasedir = pop @ARGV;
$outbasedir =~ s|/+$||;

unless(-d $outbasedir) {
    mkdir($outbasedir) || die "$progname: mkdir: $outbasedir: $!\n";
}

# For each directory...
foreach my $dir (@ARGV) {
    # Determine list of files to process.
    my @file_list = glob(sprintf('%s/*_cat.{fit,fits,fts}', $dir));
    unless(@file_list) {
	warn "$dir: no catalogues found\n" if($verbose);
	next;
    }

    print "Processing $dir\n" if($verbose);

    # Initialise globals.
    $last_datej = 0.0;
    
    $shortname = undef;
    $longname = undef;

    $mjdmin = undef;
    $mjdmax = undef;
    $seemin = undef;
    $seemax = undef;
    
    $pixsiz = undef;
    $gaincorr = undef;
    $mccd = undef;
    
    @channel = ();
    @ndetect = ();
    @modjd = ();
    @seelev = ();
    @skylev = ();
    @filtnm = ();
    @exptim = ();
    @noiselev = ();
    @ellipt = ();
    @apcore = ();
    @rcore = ();
    @dateobs = ();
    @objnm = ();

    @astrms = ();
    @astnum = ();
    @perror = ();
    
    @humid = ();
    @tau = ();
    @winddir = ();
    @windspd = ();

    @warn_list = ();

    # Create output directory.
    my $name = basename($dir);

    my $outdir = sprintf('%s/%s', $outbasedir, $name);
    unless(-d $outdir) {
	mkdir($outdir) || die "$progname: mkdir: $outdir: $!\n";
    }

    # Set up output files.
    unless(defined(open(SUMMARY, ">$outdir/$SUMMARY_FILE"))) {
	die "$progname: could not create $outdir/$SUMMARY_FILE: $!\n";
    }
    
    print SUMMARY ("Run no.  CCD     Object Name          RA    ",
		   "     Dec    Equinox   AM   PA      Date       UT         Exp",
		   "  Filt Seeing  Sky  Noise Ellipt APcor STDrms Comments\n\n");
    
    # Read files.
    my $nccd = 0;
    
    foreach (@file_list) {
	$nccd += read_head($_);
    }
    
    close(SUMMARY) || die "$progname: close: $!\n";
    
    # Do plots.
    unless(defined(open(PERCENT, ">$outdir/$PERCENT_FILE"))) {
	die "$progname: could not create $outdir/$PERCENT_FILE: $!\n";
    }
    
    plot_summary("$outdir/$SUMPLOT_FILE", 0, $nccd);
    
    close(PERCENT) || die "$progname: close: $!\n";
    
    plot_ellipt("$outdir/$ELLPLOT_FILE", 0, $nccd);
    
    unless(defined(open(SEEING, ">$outdir/$SEEING_FILE"))) {
	die "$progname: could not create $outdir/$SEEING_FILE: $!\n";
    }
    
    my($medsee, $lquart, $uquart);
    plot_seeing("$outdir/$SEEPLOT_FILE", 0, $nccd, \$medsee, \$lquart, \$uquart);
    
    close(SEEING) || die "$progname: close: $!\n";
    
    if($verbose) {
	printf("\nMedian seeing = %f\n", $medsee);
	printf("Quartiles     = %f %f\n", $lquart, $uquart);
    }

    plot_astrom("$outdir/$ASTPLOT_FILE", 0, $nccd);
    plot_detect("$outdir/$DETPLOT_FILE", 0, $nccd);
    plot_site("$outdir/$SITEPLOT_FILE", 0, $nccd);
    
    if($htmlmode) {
	# Generate PNG plots.
	plot_summary([ "$outdir/$SUMPLOT1_PNG", "$outdir/$SUMPLOT2_PNG" ], 1, $nccd);
	plot_ellipt("$outdir/$ELLPLOT_PNG", 1, $nccd);
	plot_seeing("$outdir/$SEEPLOT_PNG", 1, $nccd);
	plot_astrom([ "$outdir/$ASTPLOT1_PNG", "$outdir/$ASTPLOT2_PNG" ], 1, $nccd);
	plot_detect("$outdir/$DETPLOT_PNG", 1, $nccd);
	plot_site([ "$outdir/$SITEPLOT1_PNG", "$outdir/$SITEPLOT2_PNG" ], 1, $nccd);

	# Write HTML file.
	unless(defined(open(HTML, "> $outdir/dqc.html"))) {
	    die "$progname: could not create $outdir/dqc.html: $!\n";
	}

	my $timestr = localtime();

	print HTML <<EOD
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
        "http://www.w3.org/TR/html4/loose.dtd">
<html><head><title>$longname: DQC information for $name</title></head>
<style type="text/css"><!--
  TH { background-color: #ccccff }
  TABLE { padding: 5 }
  IMG { border-width: thin }
--></style>
<body bgcolor="#ffffff" text="#000000">
<p><center><h1>$longname: DQC information for $name</h1></center></p>

<p><h2>Summary</h2></p>
<p><table border=0>
<tr><td align="center">
    <a href="summary.ps"><img src="summary_1.png" alt="Summary plot 1"></a></td></tr>
<tr><td align="center">Surface brightness and seeing by day.</td></tr>
</table></p>
<p><table border=0>
<tr><td align="center">
    <a href="summary.ps"><img src="summary_2.png" alt="Summary plot 2"></a></td></tr>
<tr><td align="center">1-sigma sky noise and 5-sigma magnitude limit by day.</td></tr>
</table></p>
<p><h3>Filter legend:</h3></p>
<p><table border>
<tr><th>Filter</th><th>Colour</th></tr>
EOD
    ;

       my $rv = pgopen('/NULL');
       die "$progname: could not open PGPLOT\n" if($rv <= 0);

       pgscr(0, 1.0, 1.0, 1.0);
       pgscr(1, 0.0, 0.0, 0.0);

       my @filt_seen = ();
       for(my $i = 0; $i <= $#filter_names; $i++) {
	   $filt_seen[$i] = 0;
       }

       foreach my $f (@filtnm) {
	   for(my $i = 0; $i <= $#filt_icol; $i++) {
	       if(substr($f, 0, length($filt_icol[$i]->[0])) eq $filt_icol[$i]->[0]) {
		   $filt_seen[$filt_icol[$i]->[1]]++;
		   last;
	       }
	   }
       }

       for(my $i = 0; $i <= $#filt_seen; $i++) {
	   if($filt_seen[$i] > 0) {
	       my($r, $g, $b);
	       pgqcr($i, $r, $g, $b);
	       
	       printf HTML ("<tr><td>%s</td><td width=\"15\" " .
			    "bgcolor=\"#%02.2x%02.2x%02.2x\">&nbsp;</td></tr>\n",
			    $filter_names[$i]->[1], $r * 255.0, $g * 255.0, $b * 255.0);
	   }
       }

       print HTML <<EOD
</table></p>
<p><ul>
<li><a href="summary.ps">Postscript file containing plots</a>
<li><a href="summary.list">Summary ASCII table</a>
<li><a href="sky.percentages">Sky percentage ASCII table</a>
</ul></p>

<p><h2>Ellipticity</h2></p>
<p><table border=0>
<tr><td align="center">
    <a href="ellipt.ps"><img src="ellipt.png" alt="Ellipticity plot"></a></td></tr>
<tr><td align="center">Average stellar ellipticity by day.</td></tr>
</table></p>
<p><h3>Channel legend:</h3></p>
<p><table border>
<tr><th>Channel</th><th>Colour</th></tr>
EOD
    ;

       for(my $ccd = 1; $ccd <= $mccd; $ccd++) {
	   my($r, $g, $b);
	   pgqcr($ccd + 1, $r, $g, $b);
	   
	   printf HTML ("<tr><td>%d</td><td width=\"15\" " .
			"bgcolor=\"#%02.2x%02.2x%02.2x\">&nbsp;</td></tr>\n",
			$ccd, $r * 255.0, $g * 255.0, $b * 255.0);
       }

       pgclos();

       print HTML <<EOD
</table></p>
<p><ul>
<li><a href="ellipt.ps">Postscript file containing plot</a>
</ul></p>

<p><h2>Seeing</h2></p>
<p><table border=0>
<tr><td align="center">
    <a href="seeing.ps"><img src="seeing.png" alt="Seeing plot"></a></td></tr>
<tr><td align="center">Histogram of seeing.</td></tr>
</table></p>
<p><table border=0>
EOD
    ;
        printf HTML ("<tr><td><b>Median seeing</b>: </td><td colspan=\"2\">%.3f</td></tr>\n",
		     $medsee);
	printf HTML ("<tr><td><b>Quartiles</b>: </td><td>%.3f</td><td>%.3f</td></tr>\n",
		     $lquart, $uquart);

        print HTML <<EOD
</table></p>
<p><ul>
<li><a href="seeing.ps">Postscript file containing plot</a>
<li><a href="seeing.hist">Seeing histogram ASCII table</a>
</ul></p>

<p><h2>Astrometry</h2></p>
<p><table border=0>
<tr><td align="center">
    <a href="astrom.ps"><img src="astrom_1.png" alt="Astrometry plot 1"></a></td></tr>
<tr><td align="center">RMS fit error and number of standards used per CCD.</td></tr>
</table></p>
<p><table border=0>
<tr><td align="center">
    <a href="astrom.ps"><img src="astrom_2.png" alt="Astrometry plot 2"></a></td></tr>
<tr><td align="center">Mean pointing error.</td></tr>
</table></p>
<p><ul>
<li><a href="astrom.ps">Postscript file containing plots</a>
</ul></p>

<p><h2>Detection</h2></p>
<p><table border=0>
<tr><td align="center">
    <a href="detect.ps"><img src="detect.png" alt="Detection plot"></a></td></tr>
<tr><td align="center">Stellar aperture correction and number of detected sources per CCD.</td></tr>
</table></p>
<p><ul>
<li><a href="detect.ps">Postscript file containing plot</a>
</ul></p>

<p><h2>Site</h2></p>
<p><table border=0>
<tr><td align="center">
    <a href="site.ps"><img src="site_1.png" alt="Site plot 1"></a></td></tr>
<tr><td align="center">Relative humidity and Tau by day.</td></tr>
</table></p>
<p><table border=0>
<tr><td align="center">
    <a href="site.ps"><img src="site_2.png" alt="Site plot 2"></a></td></tr>
<tr><td align="center">Wind speed and direction by day.</td></tr>
</table></p>
<p><ul>
<li><a href="site.ps">Postscript file containing plots</a>
</ul></p>

EOD
    ;

       if(@warn_list) {
	   print HTML "<p><h2>Warnings</h2></p>\n";
	   print HTML "<p><pre>\n";

	   foreach (@warn_list) {
	       print HTML ($_, "\n");
	   }

	   print HTML "</pre></p>\n";
       }

       print HTML <<EOD
<hr>
<p>Last updated: $timestr</p>
</body></html>
EOD
    ;
    }
}

exit(0);

### Functions ###

sub read_head {
    my($filename) = @_;

    # Strip trailing @-signs.
    $filename =~ s/\@$//;

    # Open file.
    my $status = 0;

    my $fits = CFITSIO::open_file($filename, READONLY, $status);
    fitsio_err($status, "ffopen: $filename") if($status);

    # Files are named like "r{run}_{ccd}_cat.fits" or "r{run}_cat.fits".
    # Strip off the _cat.fits bit if it is there.
    my $btitle = basename($filename);
    $btitle =~ s/\.fits?$//;
    $btitle =~ s/_cat$//;

    # Check if the extended filename syntax is in use.  If so, we
    # only process the current extension.
    my $ext;
    ffextn($filename, $ext, $status);
    fitsio_err($status, "could parse filename $filename") if($status);

    my $onehdu = ($ext == -99 ? 0 : 1);

    my $numproc = 0;

  OUTER:
    while(1) {
	my $mefname = $filename;

	if(!$onehdu) {
	    # Find the next table extension.
	    my $hdutype = ANY_HDU;
	    while($hdutype != BINARY_TBL && $hdutype != ASCII_TBL) {
		$fits->movrel_hdu(1, $hdutype, $status);
		if($status == END_OF_FILE) {
		    $status = 0;
		    last OUTER;
		}
		elsif($status) {
		    fitsio_err($status, "ffmrhd");
		}
	    }

	    if($hdutype == ASCII_TBL) {
		die "$progname: ASCII table I/O not yet implemented\n";
	    }

	    $mefname .= sprintf('[%d]', $numproc + 1);
	}

	# Read table dimensions.
	my($nrows, $ncols);
	$fits->get_num_rows($nrows, $status);
	$fits->get_num_cols($ncols, $status);
	fitsio_err($status, "could not read table dimensions") if($status);

	# See if can recognise the detector system
	my $detector = "WFC";
	$fits->read_key_str("DETECTOR", $detector, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $fits->read_key_str("INSTRUME", $detector, undef, $status);
	    if($status == KEY_NO_EXIST) {
		$status = 0;
	    }
	    elsif($status) {
		fitsio_err($status, "ffgkys: INSTRUME");
	    }
	}
	elsif($status) {
	    fitsio_err($status, "ffgkys: DETECTOR");
	}

	my $found = 0;
	my $gain;

	my @optaxis_x;
	my @optaxis_y;

	for(my $i = 0; $i <= $#detector_params; $i++) {
	    my $str = $detector_params[$i]->[0];

	    if($detector =~ /$str/) {
		($detector, $pixsiz, $gain, $mccd, $shortname, $longname)
		    = @{$detector_params[$i]};
		@optaxis_x = @{$oaxis_x[$i]};
		@optaxis_y = @{$oaxis_y[$i]};
		$found = 1;
		last;
	    }
	}
	
	die "$progname: unrecognised detector: $detector\n" if(!$found);

	$gaincorr = 2.5 * log10($gain);

	$fits->read_key_flt("PIXLSIZE", $pixsiz, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	}
	elsif($status) {
	    fitsio_err($status, "ffgkye: PIXLSIZE");
	}

	my $skylevel;
	$fits->read_key_flt("SKYLEVEL", $skylevel, undef, $status);
	fitsio_err($status, "ffgkye: SKYLEVEL") if($status);

	my $skynoise;
	$fits->read_key_flt("SKYNOISE", $skynoise, undef, $status);
	fitsio_err($status, "ffgkye: SKYNOISE") if($status);
	
	my $seeing;
	$fits->read_key_flt("SEEING", $seeing, undef, $status);
	fitsio_err($status, "ffgkye: SEEING") if($status);

	my $ell;
	$fits->read_key_flt("ELLIPTIC", $ell, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    reportwarn("%s: ellipticity information absent", $mefname);
	    $ell = 0.0;
	}
	elsif($status) {
	    fitsio_err($status, "ffgkye: ELLIPTIC");
	}

	my $apcor;
	$fits->read_key_flt("APCOR", $apcor, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $apcor = 0.0;
	}
	elsif($status) {
	    fitsio_err($status, "ffgkye: APCOR");
	}

	my $rcor;
	$fits->read_key_flt("RCORE", $rcor, undef, $status);
	fitsio_err($status, "ffgkye: RCORE") if($status);

	# Grab IRAF ranges to see if stupid
	my $rafmin;
	$fits->read_key_flt("IRAFMIN", $rafmin, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $rafmin = 0.0;
	}
	elsif($status) {
	    fitsio_err($status, "ffgkye: IRAFMIN");
	}
	
	my $rafmax;
	$fits->read_key_flt("IRAFMAX", $rafmax, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $rafmax = 0.0;
	}
	elsif($status) {
	    fitsio_err($status, "ffgkye: IRAFMAX");
	}

	if($skylevel <= 0.0 || $skynoise <= 2.0 || $rafmax - $rafmin > 2.0e5) {
	    reportwarn("%s: skylevl = %.2f, noise = %.2f, IRAF range: %.2f %.2f",
		       $mefname, $skylevel, $skynoise, $rafmin, $rafmax);
	}

	my $datej;
	$fits->read_key_flt("MJD-OBS", $datej, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    reportwarn("%s: missing keyword MJD-OBS", $mefname);
	    $datej = $last_datej;
	}
	elsif($status) {
	    fitsio_err($status, "ffgkye: MJD-OBS");
	}
	
	my $exptime;
	$fits->read_key_flt("EXPTIME", $exptime, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $fits->read_key_flt("EXPOSED", $exptime, undef, $status);
	    if($status == KEY_NO_EXIST) {
		$status = 0;
		$fits->read_key_flt("EXP_TIME", $exptime, undef, $status);
		fitsio_err($status, "ffgkye: EXP_TIME") if($status);
	    }
	    elsif($status) {
		fitsio_err($status, "ffgkye: EXPOSED");
	    }
	}
	elsif($status) {
	    fitsio_err($status, "ffgkye: EXPTIME");
	}

	$exptime = abs($exptime);
	$exptime = 1.0 if($exptime < 1.0);

	# DAS channel info and set #5 to #2 since using DAS channel 5 for CCD 2
	my $ichan;

	$fits->read_key_lng("DASCHAN", $ichan, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $fits->read_key_lng("IMAGEID", $ichan, undef, $status);
	    if($status == KEY_NO_EXIST) {
		# Use extension number.
		$status = 0;
		$ichan = $numproc + 1;
	    }
	    elsif($status) {
		fitsio_err($status, "ffgkyj: IMAGEID");
	    }
	}
	elsif($status) {
	    fitsio_err($status, "ffgkyj: DASCHAN");
	}

	# Fix CCD numbers as needed.
	if($ichan == 5 && $detector eq "WFC") {
	    $ichan = 2;
	}
	if($detector eq "CFH12K") {
	    $ichan++;
	}

	my $filter;
	$fits->read_key_str("WFFBAND", $filter, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $fits->read_key_str("FILTER", $filter, undef, $status);
	    fitsio_err($status, "ffgkys: FILTER") if($status);
	}
	elsif($status) {
	    fitsio_err($status, "ffgkys: WFFBAND") if($status);
	}

	$filter .= " ";  # kludge to allow matching code to work

	if($detector eq "WFC" && $filter =~ /B /) {
	    my $filtsys;
	    $fits->read_key_str("WFFPSYS", $filtsys, undef, $status);
	    if($status == KEY_NO_EXIST) {
		$status = 0;
	    }
	    elsif($status) {
		fitsio_err($status, "ffgkys: WFFPSYS");
	    }
	    else {
		if($filtsys =~ /Harris/) {
		    $filter = "B_Harris";
		}
	    }
	}

	# Other info
	my $objname;
	$fits->read_key_str("OBJECT", $objname, undef, $status);
	fitsio_err($status, "ffgkys: OBJECT") if($status);

	my $dateraw;
	$fits->read_key_str("DATE-OBS", $dateraw, undef, $status);
	fitsio_err($status, "ffgkys: DATE-OBS") if($status);

	my($datechar, $utdate) =
	    ($dateraw =~ /^\s*(\d+[\s\-]\d+[\s\-]\d+)T(\d+\:\d+\:\d+\.?\d*)/);

	my($yr, $mn, $dy) = ($datechar =~ /^\s*(\d+)[\s\-](\d+)[\s\-](\d+)/);
	die "$progname: $dateraw: unrecognised DATE-OBS value"
	    unless(defined $yr && defined $mn && defined $dy);
	my $datemobs = $yr + $mn / 12.0;

	my $utchar;
	$fits->read_key_str("UTSTART", $utchar, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $fits->read_key_str("UTC-OBS", $utchar, undef, $status);
	    if($status == KEY_NO_EXIST) {
		$status = 0;
		$utchar = $utdate;
	    }
	    elsif($status) {
		fitsio_err($status, "ffgkys: UTC-OBS");
	    }
	}
	elsif($status) {
	    fitsio_err($status, "ffgkys: UTSTART");
	}

	my($rachar, $decchar, $telra, $teldec);
	my $catchar;

	my $tmp = 0.0;
	$fits->read_key_flt("HUMIDITY", $tmp, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $tmp = 0.0;
	}
	
	push(@humid, $tmp);
	
	$fits->read_key_flt("CSOTAU", $tmp, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $tmp = 0.0;
	}
	
	push(@tau, $tmp);
	
	$fits->read_key_flt("WINDDIR", $tmp, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $tmp = 0.0;
	}
	
	push(@winddir, $tmp);
	
	$fits->read_key_flt("WINDSPD", $tmp, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $tmp = 0.0;
	}
	
	push(@windspd, $tmp);
	
	$fits->read_key_str("RA", $rachar, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $fits->read_key_flt("TELRA", $telra, undef, $status);
	    if($status) {
		fitsio_err($status, "ffgkye: TELRA");
	    }
	    else {
		$telra *= DEG_TO_RAD;
		$rachar = base10_to_60($telra, UNIT_RAD, ':', ' ', 3, UNIT_HR);
	    }
	}
	elsif($status) {
	    fitsio_err($status, "ffgkys: RA");
	}
	else {
	    my $rv;
	    ($rv, $telra) = base60_to_10($rachar, ':', UNIT_HR, UNIT_RAD);
	    die "$progname: could not understand RA $rachar\n" if($rv < 0);
	}

	$fits->read_key_str("DEC", $decchar, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $fits->read_key_flt("TELDEC", $teldec, undef, $status);
	    if($status) {
		fitsio_err($status, "ffgkye: TELDEC");
	    }
	    else {
		$teldec *= DEG_TO_RAD;
		$decchar = base10_to_60($teldec, UNIT_RAD, ':', '+', 2, UNIT_DEG);
	    }
	}
	elsif($status) {
	    fitsio_err($status, "ffgkys: DEC");
	}
	else {
	    my $rv;
	    ($rv, $teldec) = base60_to_10($decchar, ':', UNIT_DEG, UNIT_RAD);
	    die "$progname: could not understand Dec $decchar\n" if($rv < 0);
	}

	$fits->read_key_str("CAT-EQUI", $catchar, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $fits->read_key_str("EQUINOX", $catchar, undef, $status);
	    fitsio_err($status, "ffgkys: EQUINOX") if($status);
	}
	elsif($status) {
	    fitsio_err($status, "ffgkys: CAT-EQUI");
	}

	my $rotskypa;
	$fits->read_key_flt("ROTSKYPA", $rotskypa, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    reportwarn("%s: missing keyword ROTSKYPA", $mefname);
	    $rotskypa = -1.0;
	}

	my $airmass;
	$fits->read_key_flt("AIRMASS", $airmass, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    
	    my $amstart;
	    $fits->read_key_flt("AMSTART", $amstart, undef, $status);

	    my $amend;
	    $fits->read_key_flt("AMEND", $amend, undef, $status);

	    if($status == KEY_NO_EXIST) {
		$status = 0;
		reportwarn("%s: missing keyword AIRMASS", $mefname);
		$airmass = -1.0;
	    }

	    $airmass = 0.5 * ($amstart + $amend);
	}

	my $stdcrms;
	$fits->read_key_flt("STDCRMS", $stdcrms, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $stdcrms = 0.0;
	}

	my $numbrms;
	$fits->read_key_flt("NUMBRMS", $numbrms, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $numbrms = 0.0;
	}

	# Read WCS
	my($tpa, $tpd);
	$fits->read_key_flt("CRVAL1", $tpa, undef, $status);
	$fits->read_key_flt("CRVAL2", $tpd, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    reportwarn("%s: missing CRVAL keyword", $mefname);
	    $tpa = -1.0;
	    $tpd = 0.0;
	}

	my($a, $b, $c, $d, $e, $f, $projp1, $projp3);
	$fits->read_key_flt("CD1_1", $a, undef, $status);
	$fits->read_key_flt("CD1_2", $b, undef, $status);
	$fits->read_key_flt("CRPIX1", $c, undef, $status);
	$fits->read_key_flt("CD2_1", $d, undef, $status);
	$fits->read_key_flt("CD2_2", $e, undef, $status);
	$fits->read_key_flt("CRPIX2", $f, undef, $status);

	$fits->read_key_flt("PROJP1", $projp1, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $projp1 = 1.0;
	}

	$fits->read_key_flt("PROJP3", $projp3, undef, $status);
	if($status == KEY_NO_EXIST) {
	    $status = 0;
	    $projp3 = 220.0;
	}

	if($status) {
	    $status = 0;
	    reportwarn("%s: could not read WCS", $mefname);
	    $tpa = -1.0;
	}

	# Calculate pointing error assuming tangent points are at
	# rotator centre.
	my $perr = 0.0;
	if($telra >= 0.0 && $tpa >= 0.0) {
	    # Got keywords...
	    $tpa *= DEG_TO_RAD;
	    $tpd *= DEG_TO_RAD;
	    $a  *= DEG_TO_RAD;
	    $b  *= DEG_TO_RAD;
	    $d  *= DEG_TO_RAD;
	    $e  *= DEG_TO_RAD;
	    
	    my($cra, $cdec) = radeczp($optaxis_x[$ichan-1], $optaxis_y[$ichan-1],
				      $tpa, $tpd, $a, $b, $c, $d, $e, $f,
				      $projp1, $projp3);

	    my $cosperr = cos($telra) * cos($teldec) * cos($cra) * cos($cdec) + 
		          sin($telra) * cos($teldec) * sin($cra) * cos($cdec) +
			  sin($teldec) * sin($cdec);

	    $perr = acos($cosperr) * RAD_TO_AM;
	}

	# Fill up storage arrays.
	my $tmodjd = $datej - 51179.0;
	$tmodjd += 365.0 if($tmodjd < 0.0);
	$mjdmin = $tmodjd if(!defined $mjdmin || $tmodjd < $mjdmin);
	$mjdmax = $tmodjd if(!defined $mjdmax || $tmodjd > $mjdmax);

	my $tseelev = $seeing * $pixsiz;
	$tseelev = -1.0 if($tseelev < 0.0);
	if($seeing > 1.0) {  # to trap FU's
	    $seemin = $tseelev if(!defined $seemin || $tseelev < $seemin);
	    $seemax = $tseelev if(!defined $seemax || $tseelev > $seemax);
	}

	my $tskylev = $skylevel / ($exptime * $pixsiz * $pixsiz);
	$tskylev = 0.1 if($tskylev < 0.1);

	push(@channel, $ichan);
	push(@ndetect, $nrows);
	push(@modjd, $tmodjd);
	push(@seelev, $tseelev);
	push(@skylev, $tskylev);
	push(@filtnm, $filter);
	push(@exptim, $exptime);
	push(@noiselev, $skynoise / ($exptime * $pixsiz));  # ie. per sq arcsec
	push(@ellipt, $ell);
	push(@apcore, $apcor);
	push(@rcore, $rcor);
	push(@dateobs, $datemobs);
	push(@objnm, $objname);

	push(@astrms, $stdcrms);
	push(@astnum, $numbrms);
	push(@perror, $perr);

	$numproc++;

	my $title;
	if($detector eq 'ALADDIN') {
	    my($tmp1, $tmp2) = ($btitle =~ /^(\w+)_\d+_(\d+)/);
	    $title = sprintf("%s_%s_%d", $tmp1, $tmp2, $ichan);
	}
	else {
	    $title = sprintf("%s_%d", $btitle, $ichan);
	}

	print "\n$title\n" if($verbose);

	# Write out the DQC info
	$rachar = "  " . $rachar;
	$decchar = "  " . $decchar;

	my $ir1 = 3;
	my $ir2 = $ir1 + 10;
	if(substr($rachar, $ir2, 1) eq "'") {
	    $ir1--;
	    $ir2--;
	    substr($rachar, $ir1, 1) = " " if(substr($rachar, $ir1, 1) eq "'");
	}
	if(substr($rachar, $ir2 - 1, 1) eq "'") {
	    $ir1 -= 2;
	    $ir2 -= 2;
	    substr($rachar, $ir1, 1) = " " if(substr($rachar, $ir1, 1) eq "'");
	}

	my $id1 = 2;
	my $id2 = $id1 + 10;
	if(substr($decchar, $id2 - 3, 1) eq ".") {
	    $id1 -= 2;
	    substr($decchar, $id1, 2) = "  ";
	    $id2 -= 2;
	}
	if(substr($decchar, $id2 - 2, 1) eq ".") {
	    $id1--;
	    substr($decchar, $id1, 1) = " ";
	    $id2--;
	}
	substr($decchar, $id2, 1) = " " if(substr($decchar, $id2, 1) eq "'");
	substr($decchar, $id2 - 1, 1) = " " if(substr($decchar, $id2 - 1, 1) eq "'");

	printf SUMMARY ("%-12.12s    %-17.17s %s %s %-5.5s%6.3f%6.1f %-11.11s%-10.10s" .
			"%8.2f  %-3.3s%5.2f%8.1f%6.1f%6.2f%6.2f%6.2f\n",
			$title, $objname,
			substr($rachar, $ir1, $ir2 - $ir1 + 1),
			substr($decchar, $id1, $id2 - $id1 + 1),
			$catchar, $airmass, $rotskypa, $datechar, $utchar, $exptime,
			$filter, $tseelev, $skylevel, $skynoise, $ell, $apcor,
			$stdcrms);

	last if($onehdu);
    }

    # Done, close file.
    $fits->close_file($status);
    fitsio_err($status, "ffclos") if($status);

    return($numproc);
}

sub plot_summary {
    my($outdev, $ispng, $nlist) = @_;

    my $sbmin = 18.0;
    my $sbmax = 12.5;
    my $seemin = 0.5;
    my $seemax = 3.0;
    my $djdmin = $mjdmin - 1.1;
    my $djdmax = $mjdmax + 1.1;

    my $rv = pgopen($ispng ? $outdev->[0] : $outdev);
    die("$progname: could not open plot device: ",
	($ispng ? $outdev->[0] : $outdev), "\n") if($rv <= 0);

    if($ispng) {
	pgpap($PAPSIZE, $PAPASPECT);
	pgscr(0, 1.0, 1.0, 1.0);
	pgscr(1, 0.0, 0.0, 0.0);
    }

    pgsubp(1, 2);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, $sbmin, $sbmax, 0, 0);
    pglabel("MMJD (days)", "Surface Brightness (mag/sqarcsec)", $shortname);
    pgsch(2.5);

    my $nsky = 0;
    my $avsky = 0.0;

    my @medsky = ();
    for(my $i = 0; $i < $mccd; $i++) {
	$medsky[$i] = 0.0;
    }

    my @nwork = ();
    for(my $i = 0; $i <= $NFILT; $i++) {
	$nwork[$i] = 0;
    }

    my @fwork = ();

    my $line;

    for(my $kj = 0; $kj < $nlist; $kj++) {
	# Look up ccdzero and icol.
	my $ccdzero = $DEFAULT_CCDZERO;

	foreach (@ccdzero_list) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$ccdzero = $_->[1];
		last;
	    }
	}

	my $icol = -1;

	foreach (@filt_icol) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$icol = $_->[1];
		last;
	    }
	}

	pgsci($icol);

	my $ichan = $channel[$kj];
	$avsky += $skylev[$kj];
	$nsky++;
	$medsky[$ichan - 1] = $skylev[$kj] * $exptim[$kj] * $pixsiz * $pixsiz;

	if($ichan == $mccd) {
	    my $sb = $ccdzero - 2.5 * log10($avsky / $nsky);
	    my $daytim = $modjd[$kj];

	    if($exptim[$kj] > 50.0) {
		pgpoint(1, $daytim, $sb, 17);
	    }
	    else {
		pgpoint(1, $daytim, $sb, 21);
	    }

	    # To check % variation in sky level on different CCDs
	    if($exptim[$kj] >= 100.0) {
		my @work = sort { $a <=> $b } @medsky;
		my $skymed = 0.5 * ($work[$mccd/2 - 1] + $work[$mccd/2]);
		if($skymed > 100.0) {
		    $line = sprintf("%-15.15s", $objnm[$kj]);

		    for(my $i = 0; $i < $mccd; $i++) {
			$work[$i] = 100.0 * ($medsky[$i] - $skymed) / $skymed;
			$line .= sprintf("%8.2f", $work[$i]);
		    }

		    $line .= sprintf("%8.1f%8.1f   %-12.12s\n",
				     $skymed, $exptim[$kj], $filtnm[$kj]);

		    unless($ispng) {
			print PERCENT $line;
		    }		    

		    my $nw = $nwork[$icol];
		    $fwork[$nw] = [] unless(defined $fwork[$nw]);
		    $fwork[$nw]->[$icol] = \@work;
		    $nwork[$icol]++;
		}
	    }

	    $nsky = 0;
	    $avsky = 0.0;
	}
    }

    pgsci(1);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, $seemin, $seemax, 0, 0);
    pglabel("MMJD (days)", "Seeing (arcsec)", " ");
    pgsch(2.5);

    my $nsee = 0;
    my $avsee = 0.0;

    for(my $kj = 0; $kj < $nlist; $kj++) {
	my $icol = -1;

	foreach (@filt_icol) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$icol = $_->[1];
		last;
	    }
	}

	pgsci($icol);

	my $ichan = $channel[$kj];
	my $seeing = $seelev[$kj];

	if($seeing > 0.0) {
	    $avsee += $seeing;
	    $nsee++;
	}

	if($ichan == $mccd) {
	    $seeing = $avsee / ($nsee < 1 ? 1 : $nsee);
	    my $daytim = $modjd[$kj];

	    if($exptim[$kj] > 50.0) {
		pgpoint(1, $daytim, $seeing, 17);
	    }
	    else {
		pgpoint(1, $daytim, $seeing, 21);
	    }

	    $nsee = 0;
	    $avsee = 0.0;
	}
    }

    $sbmin += 5.5;
    $sbmax += 5.5;

    if($ispng) {
	pgclos();

	my $rv = pgopen($outdev->[1]);
	die("$progname: could not open plot device: ", $outdev->[1], "\n") if($rv <= 0);

	pgpap($PAPSIZE, $PAPASPECT);
	pgscr(0, 1.0, 1.0, 1.0);
	pgscr(1, 0.0, 0.0, 0.0);

	pgsubp(1, 2);
    }

    pgsch(1.5);
    pgsci(1);
    pgenv($djdmin, $djdmax, $sbmin, $sbmax, 0, 0);
    pglabel("MMJD (days)", "1-sigma sky noise (mag/sqarcsec)", $shortname);
    pgsch(2.5);

    $nsky = 0;
    $avsky = 0.0;

    for(my $kj = 0; $kj < $nlist; $kj++) {
	# Look up ccdzero and icol.
	my $ccdzero = $DEFAULT_CCDZERO;

	foreach (@ccdzero_list) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$ccdzero = $_->[1];
		last;
	    }
	}

	my $icol = -1;

	foreach (@filt_icol) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$icol = $_->[1];
		last;
	    }
	}

	pgsci($icol);
	
	my $ichan = $channel[$kj];

	$avsky += $noiselev[$kj];
	$nsky++;

	if($ichan == $mccd) {
	    my $sb = $ccdzero - 2.5 * log10($avsky / $nsky);
	    my $daytim = $modjd[$kj];

	    if($exptim[$kj] > 50.0) {
		pgpoint(1, $daytim, $sb, 17);
	    }
	    else {
		pgpoint(1, $daytim, $sb, 21);
	    }

	    $nsky = 0;
	    $avsky = 0.0;
	}
    }

    pgsci(1);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, 21.0, 14.0, 0, 0);
    pglabel("MMJD (days)", "5-sigma magnitude limit", " ");
    pgsch(2.5);

    my $nmag = 0;
    my $avmag = 0.0;
    my $s2n = 5.0;

    for(my $kj = 0; $kj < $nlist; $kj++) {
	# Look up ccdzero and icol.
	my $ccdzero = $DEFAULT_CCDZERO;

	foreach (@ccdzero_list) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$ccdzero = $_->[1];
		last;
	    }
	}

	my $icol = -1;

	foreach (@filt_icol) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$icol = $_->[1];
		last;
	    }
	}

	pgsci($icol);

	my $ichan = $channel[$kj];
	my $seeing = abs($seelev[$kj]);
	my $sigg = $seeing / (2.0 * sqrt(log(2.0)));  # Gaussian profile
	my $a = 2.0 * pi * $sigg * $sigg;
	my $skyb = $ccdzero - 2.5 * log10($skylev[$kj]);
	my $objmag = 0.5 * ($ccdzero + $gaincorr + $skyb -
                            2.5 * log10($s2n * $s2n * $a / $exptim[$kj]));
	$objmag -= 0.45;  # to allow for realistic profiles

	# Alternative s:n=5 limit for core mags
	my $skynoise = $noiselev[$kj] * $pixsiz;
	my $apcor = $apcore[$kj];
	my $rcor = $rcore[$kj];
	my $altmag = $ccdzero - 
                     2.5 * log10(5.0 * sqrt(pi * $rcor * $rcor) * $skynoise) - $apcor;
	
	$avmag += $objmag;
	$nmag++;

	if($ichan == $mccd) {
	    $avmag /= $nmag;
	    my $daytim = $modjd[$kj];

	    if($exptim[$kj] > 50.0) {
		pgpoint(1, $daytim, $avmag, 17);
	    }
	    else {
		pgpoint(1, $daytim, $avmag, 21);
	    }

	    $nmag = 0;
	    $avmag = 0.0;
	}
    }

    pgclos();

    unless($ispng) {
	# Some % differences in sky statistics
	print PERCENT "\nMedian % sky level differences from catalogue analysis\n\n";
	
	for(my $icol = 1; $icol <= $NFILT; $icol++) {
	    my $nn = $nwork[$icol];
	    if($nn > 0) {
		my @work = ();
		
		$line = sprintf("%s%5d", $filter_names[$icol]->[0], $icol);
		
		for(my $j = 0; $j < $mccd; $j++) {
		    for(my $i = 0; $i < $nn; $i++) {
			$work[$i] = $fwork[$i]->[$icol]->[$j];
		    }
		    
		    @work = sort { $a <=> $b } @work;
		    
		    my $ans;
		    
		    if($nn % 2 == 0) {
			$ans = 0.5 * ($work[$nn/2-1] + $work[$nn/2]);
		    }
		    else {
			$ans = $work[($nn-1)/2];
		    }
		    
		    $line .= sprintf("%8.2f", $ans);
		}
		
		$line .= sprintf("%8d\n", $nn);
		print PERCENT $line;
	    }
	}
    }
}

sub plot_ellipt {
    my($outdev, $ispng, $nlist) = @_;

    my $djdmin = $mjdmin - 1.1;
    my $djdmax = $mjdmax + 1.1;

    my $rv = pgopen($outdev);
    die "$progname: could not open plot device: $outdev\n" if($rv <= 0);

    if($ispng) {
	pgpap($PAPSIZE, $PAPASPECT);
	pgscr(0, 1.0, 1.0, 1.0);
	pgscr(1, 0.0, 0.0, 0.0);
    }

    pgsubp(1, 1);
    pgenv($djdmin, $djdmax, 0.0, 1.0, 0, 0);
    pglabel("MMJD (days)", "Average stellar ellipticity", $shortname);
    pgsch(2.5);

    for(my $kj = 0; $kj < $nlist; $kj++) {
	my $ell = $ellipt[$kj];
	my $daytim = $modjd[$kj];
	my $ichan = $channel[$kj];

	pgsci($ichan + 1);
	pgpoint(1, $daytim, $ell, 17);
    }

    pgclos();
}

sub plot_seeing {
    my($outdev, $ispng, $nlist, $medsee, $lquart, $uquart) = @_;

    # Some seeing histograms - ignore points below 0.5
    my @array = ();
    my @xcord = ();
    my $iarg;

    for($iarg = 0; $iarg <= 30; $iarg++) {
	$array[$iarg] = 0;
    }

    for(my $kj = 0; $kj < $nlist; $kj++) {
	if($seelev[$kj] > 0.5 && $seelev[$kj] < 3.0) {
	    $iarg = int(10.0 * $seelev[$kj] + 0.5);
	    $array[$iarg]++;
	}
    }

    unless($ispng) {
	print SEEING "    Seeing   No. of times\n\n";
    }

    for($iarg = 0; $iarg <= 30; $iarg++) {
	$xcord[$iarg] = 0.1 * $iarg;
	unless($ispng) {
	    printf SEEING ("%10.1f%10.1f\n", $xcord[$iarg], $array[$iarg]);
	}
    }

    my $rv = pgopen($outdev);
    die "$progname: could not open plot device: $outdev\n" if($rv <= 0);

    if($ispng) {
	pgpap($PAPSIZE, $PAPASPECT);
	pgscr(0, 1.0, 1.0, 1.0);
	pgscr(1, 0.0, 0.0, 0.0);
    }

    pgsubp(1, 1);
    pgsvp(0.2, 0.8, 0.2, 0.8);
    pgswin(0.0, 3.0, 0.0, 300.0);
    pgbox("BCNST", 0.0, 0, "BCNST", 0.0, 0);
    pglabel("Seeing (arcsec)", "Frequency", $shortname);
    pgbin(31, \@xcord, \@array, 1);

    pgclos();

    my @tmp = sort { $a <=> $b } @seelev;
    # All these array indices have been converted to C (0-indexed)
    # form (subtracting 1 from the expression) by adjusting the
    # part in brackets.
    $$medsee = $tmp[($nlist-1) / 2];
    $$lquart = $tmp[($nlist-1) / 4];
    $$uquart = $tmp[(3*$nlist-1) / 4];
}

sub plot_astrom {
    my($outdev, $ispng, $nlist) = @_;

    my $djdmin = $mjdmin - 1.1;
    my $djdmax = $mjdmax + 1.1;

    my $rv = pgopen($ispng ? $outdev->[0] : $outdev);
    die("$progname: could not open plot device: ",
	($ispng ? $outdev->[0] : $outdev), "\n") if($rv <= 0);

    if($ispng) {
	pgpap($PAPSIZE, $PAPASPECT);
	pgscr(0, 1.0, 1.0, 1.0);
	pgscr(1, 0.0, 0.0, 0.0);
    }

    pgsubp(1, 2);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, 0.0, 1.0, 0, 0);
    pglabel("MMJD (days)", "RMS astrometric fit error (arcsec)", $shortname);
    pgsch(2.5);

    my $nummax = 0;
    for(my $kj = 0; $kj < $nlist; $kj++) {
	my $val = $astrms[$kj];
	my $daytim = $modjd[$kj];
	my $ichan = $channel[$kj];

	$nummax = $astnum[$kj] if($astnum[$kj] > $nummax);

	pgsci($ichan + 1);
	pgpoint(1, $daytim, $val, 17);
    }

    pgsci(1);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, 0.1, 1.1 * log10($nummax), 0, 20);
    pglabel("MMJD (days)", "Number of astrometric standards used", "");
    pgsch(2.5);

    for(my $kj = 0; $kj < $nlist; $kj++) {
	my $val = log10($astnum[$kj]);
	my $daytim = $modjd[$kj];
	my $ichan = $channel[$kj];

	pgsci($ichan + 1);
	pgpoint(1, $daytim, $val, 17);
    }

    if($ispng) {
	pgclos();

	my $rv = pgopen($outdev->[1]);
	die("$progname: could not open plot device: ", $outdev->[1], "\n") if($rv <= 0);

	pgpap($PAPSIZE, $PAPASPECT);
	pgscr(0, 1.0, 1.0, 1.0);
	pgscr(1, 0.0, 0.0, 0.0);

	pgsubp(1, 2);
    }

    pgsci(1);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, 0.0, 10.0, 0, 0);
    pglabel("MMJD (days)", "Mean pointing error (arcmin)", "");
    pgsch(2.5);

    for(my $kj = 0; $kj < $nlist; $kj++) {
	my $ichan = $channel[$kj];

	pgsci($ichan + 1);
	pgpoint(1, $modjd[$kj], $perror[$kj], 17);
    }

    pgclos();
}

sub plot_detect {
    my($outdev, $ispng, $nlist) = @_;

    my $djdmin = $mjdmin - 1.1;
    my $djdmax = $mjdmax + 1.1;

    my $rv = pgopen($outdev);
    die "$progname: could not open plot device: $outdev\n" if($rv <= 0);

    if($ispng) {
	pgpap($PAPSIZE, $PAPASPECT);
	pgscr(0, 1.0, 1.0, 1.0);
	pgscr(1, 0.0, 0.0, 0.0);
    }

    pgsubp(1, 2);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, 0.0, 1.5, 0, 0);
    pglabel("MMJD (days)", "Stellar aperture correction", $shortname);
    pgsch(2.5);

    my $nummax = 0;
    for(my $kj = 0; $kj < $nlist; $kj++) {
	my $ichan = $channel[$kj];

	$nummax = $ndetect[$kj] if($ndetect[$kj] > $nummax);

	pgsci($ichan + 1);
	pgpoint(1, $modjd[$kj], $apcore[$kj], 17);
    }

    pgsci(1);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, 1, 1.1 * log10($nummax), 0, 20);
    pglabel("MMJD (days)", "Number of detected sources", "");
    pgsch(2.5);

    for(my $kj = 0; $kj < $nlist; $kj++) {
	my $val = log10($ndetect[$kj]);
	my $ichan = $channel[$kj];

	pgsci($ichan + 1);
	pgpoint(1, $modjd[$kj], $val, 17);
    }

    pgclos();
}

sub plot_site {
    my($outdev, $ispng, $nlist) = @_;

    my $djdmin = $mjdmin - 1.1;
    my $djdmax = $mjdmax + 1.1;

    my $rv = pgopen($ispng ? $outdev->[0] : $outdev);
    die("$progname: could not open plot device: ",
	($ispng ? $outdev->[0] : $outdev), "\n") if($rv <= 0);

    if($ispng) {
	pgpap($PAPSIZE, $PAPASPECT);
	pgscr(0, 1.0, 1.0, 1.0);
	pgscr(1, 0.0, 0.0, 0.0);
    }

    pgsubp(1, 2);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, 0.0, 100.0, 0, 0);
    pglabel("MMJD (days)", "Relative humidity (%)", $shortname);
    pgsch(2.5);

    my $kj = 0;
    my $nki = $nlist / $mccd;
    for(my $ki = 0; $ki < $nki; $ki++) {
	$kj = $ki * $mccd;  # use first CCD

	# Look up colour.
	my $icol = -1;
	
	foreach (@filt_icol) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$icol = $_->[1];
		last;
	    }
	}
	
	pgsci($icol);
	pgpoint(1, $modjd[$kj], $humid[$kj], 17);
    }

    pgsci(1);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, 0.0, 1.0, 0, 0);
    pglabel("MMJD (days)", "Tau at 2.5GHz", $shortname);
    pgsch(2.5);

    for(my $ki = 0; $ki < $nki; $ki++) {
	$kj = $ki * $mccd;  # use first CCD

	# Look up colour.
	my $icol = -1;
	
	foreach (@filt_icol) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$icol = $_->[1];
		last;
	    }
	}
	
	pgsci($icol);
	pgpoint(1, $modjd[$kj], $tau[$kj], 17);
    }

    if($ispng) {
	pgclos();

	my $rv = pgopen($outdev->[1]);
	die("$progname: could not open plot device: ", $outdev->[1], "\n") if($rv <= 0);

	pgpap($PAPSIZE, $PAPASPECT);
	pgscr(0, 1.0, 1.0, 1.0);
	pgscr(1, 0.0, 0.0, 0.0);

	pgsubp(1, 2);
    }

    pgsci(1);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, 0.0, 360.0, 0, 0);
    pglabel("MMJD (days)", "Wind direction (deg azimuth)", $shortname);
    pgsch(2.5);

    for(my $ki = 0; $ki < $nki; $ki++) {
	$kj = $ki * $mccd;  # use first CCD

	# Look up colour.
	my $icol = -1;
	
	foreach (@filt_icol) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$icol = $_->[1];
		last;
	    }
	}
	
	pgsci($icol);
	pgpoint(1, $modjd[$kj], $winddir[$kj], 17);
    }

    pgsci(1);
    pgsch(1.5);
    pgenv($djdmin, $djdmax, 0.0, 150.0, 0, 0);
    pglabel("MMJD (days)", "Wind speed (km/h)", $shortname);
    pgsch(2.5);

    for(my $ki = 0; $ki < $nki; $ki++) {
	$kj = $ki * $mccd;  # use first CCD

	# Look up colour.
	my $icol = -1;
	
	foreach (@filt_icol) {
	    if(substr($filtnm[$kj], 0, length($_->[0])) eq $_->[0]) {
		$icol = $_->[1];
		last;
	    }
	}
	
	pgsci($icol);
	pgpoint(1, $modjd[$kj], $windspd[$kj], 17);
    }

    pgclos();
}

sub radeczp {
    my($x, $y, $tpa, $tpd, $a, $b, $c, $d, $e, $f, $projp1, $projp3) = @_;

    my $secd = 1.0 / cos($tpd);
    my $tand = tan($tpd);

    $x -= $c;
    $y -= $f;

    my $xi = $a * $x + $b * $y;
    my $xn = $d * $x + $e * $y;

    my $r    = sqrt($xi * $xi + $xn * $xn);
    my $rfac = $projp1 + $projp3 * $r * $r;

    $r   /= $rfac;
    $rfac = $projp1 + $projp3 * $r * $r;
    
    $xi  /= $rfac;
    $xn  /= $rfac;

    my $aa    = atan($xi * $secd / (1.0 - $xn * $tand));
    my $alpha = $aa + $tpa;

    my $delta = ($xi != 0.0 ?
		 atan(($xn + $tand) * sin($aa) / ($xi * $secd)) :
		 atan(($xn + $tand) / (1.0 - $xn * $tand)));

    if($alpha > 2*pi) {
	$alpha -= 2*pi;
    }
    if($alpha < 0.0) {
	$alpha += 2*pi;
    }    

    return($alpha, $delta);
}

sub reportwarn {
    my $fmt = shift;
    my $str = sprintf($fmt, @_);
    push(@warn_list, $str);

    warn "$progname: $str\n" if($verbose);
}

sub basename {
    my($path) = @_;

    if($path =~ m|([^/]+)$|) {
        return($1);
    }

    return($path);
}

sub fitsio_err {
    my $status = shift;

    my $errmsg = '';
    ffgerr($status, $errmsg);

    die($progname, ": ", @_, ": $errmsg\n");
}

# Docs
=pod

=head1 NAME

tableplot_ir.pl - make a DQC log for IR data

=head1 SYNOPSIS

B<tableplot_ir.pl> [B<-nv>] I<indir> [...] I<outdir>

=head1 DESCRIPTION

B<tableplot_ir.pl> generates a data quality control (DQC) log for
infra-red data (initially WFCAM), from a set of object catalogue files
in I<indir>.  The log files are written to the directory I<outdir>,
which is created first if necessary, with a subdirectory made for each
input catalogue directory, named by taking the basename(1) of
I<indir> (see EXAMPLES, below).

The default is to create an HTML page, F<dqc.html>, containing the
generated plots as inline PNG images, with links to Postscript
versions of the plots (for printing) and to the summary files
generated.  This mode can be disabled, such that only the summary
files and Postscript plots are produced (see OPTIONS, below).

The catalogue files are searched for in the directories I<indir>
specified on the command line.  They are expected to be named with
F<_cat.fits> at the end of the filename (extensions of F<.fit> and
F<.fts> are also supported).

The following files are generated:

=over

=item F<summary.ps>

Plots of the surface brightness, seeing, 1-sigma sky noise and 5-sigma
magnitude limit, per filter, by day.

=item F<summary.list>

A DQC summary ASCII table containing one line per detector in each
data-set, with various DQC parameters.  The columns are described
under SUMMARY COLUMNS, below.

=item F<sky.percentages>

A table of the median percentage sky level differences from catalogue
analysis.  This is used as input to fitsio_percent(1).

=item F<ellipt.ps>

A plot of the average stellar ellipticity for each detector, per day.

=item F<seeing.ps>

A histogram of the seeing over all detectors and days.

=item F<seeing.hist>

An ASCII table of the data used to produce F<seeing.ps>.

=item F<astrom.ps>

Plots of the RMS astrometric fit error, number of astrometric
standards used in calibration and mean pointing error for each
detector, by day.

=item F<detect.ps>

Plots of the stellar aperture correction and number of detected
sources for each detector, by day.

=item F<site.ps>

Plots of the relative humidity, Tau at 2.5GHz, wind direction (in
degrees azimuth) and wind speed by day.

=back

=head1 OPTIONS

The following command-line options are supported:

=over

=item B<-n>

Disables HTML output: only the files listed above are produced.
Otherwise, additional PNG images and a file F<dqc.html> are
produced.

=item B<-v>

Increases the verbosity level of the program.

=back

=head1 SUMMARY COLUMNS

The following columns are present in the summary files
(F<summary.list>):

=over

=item I<Run no.>

An identifier for the frame in the format I<w_{run}_{chan}> where
I<{run}> is the run number, and I<{chan}> is the detector DAS channel.

=item I<CCD>

A blank column for historical reasons.

=item I<Object Name>

Name of the target, taken from the I<OBJECT> FITS header.

=item I<RA>

Right ascension of the desired target position, in sexagesimal format
with the hours, minutes, seconds separated by C<:>.

=item I<Dec>

Declination of the desired target position, in sexagesimal format with
the degrees, arcminutes, arcseconds separated by C<:>.

=item I<Equinox>

The equinox of I<RA> and I<Dec>.  For WFCAM data no leading C<J> is
present: the equinoxes are always Julian.

=item I<AM>

The effective mean airmass.

=item I<PA>

Position angle of the observation.  For WFCAM data this is always
C<-1.0> since the detector cannot be rotated.  For INT and other data
this gives the rotator angle.

=item I<Date>

The observation date (in UTC).  This is in the format I<yyyy-mm-dd>.

=item I<UT>

The start UT of the integration, in the format I<hh:mm:ss>.

=item I<Exp>

The measured actual length of the exposure computed in the
micro-controller.

=item I<Filt>

Filter used in the observation.  For IR data this will generally be
C<J>, C<H> or C<K>.

=item I<Seeing>

Average realistic stellar FWHM in arcseconds.

=item I<Sky>

Median sky brightness, in ADU.

=item I<Noise>

Robust MAD estimator for the pixel noise at the sky level, scaled to
the equivalent Gaussian RMS, in ADU.

=item I<Ellipt>

A direct estimate of the average stellar ellipticity in the frame,
useful for spotting trailed frames, etc.

=item I<APcor>

Stellar image aperture correction for the default recommended aperture
measure, diameter =2x{FWHM}, in magnitudes.

=item I<STDrms>

RMS astrometric fit error in arcseconds.

=item I<Comments>

Extra space for insertion of user comments.

=back

=head1 EXAMPLES

Process catalogues from the directory F</directory/indir>, writing
output to F<outdir>.  This fill create a directory named
F<outdir>/F<indir> containing the output files.

  tableplot_ir.pl /directory/indir outdir

The same, except with no HTML output.

  tableplot_ir.pl -n /directory/indir outdir

=head1 SEE ALSO

basename(1), fitsio_percent(1)

=head1 AUTHOR

Jonathan Irwin (jmi@ast.cam.ac.uk)

Based on fitsiotableplot, written by

Mike Irwin (mike@ast.cam.ac.uk)

=cut