# $Revision: 1.10 $
# $Source: /headas/headas/attitude/tasks/tristarid/StarID/SourceTable.pm,v $
# $Date: 2009/03/31 14:24:54 $
#
# $Log: SourceTable.pm,v $
# Revision 1.10  2009/03/31 14:24:54  rwiegand
# Corrected TUNIT values to follow conventions.
#
# Revision 1.9  2008/10/17 14:35:36  rwiegand
# Updated TDISPn for RA and DEC columns.
#
# Revision 1.8  2008/03/28 15:24:35  rwiegand
# Removed stray newline.
#
# Revision 1.7  2007/10/05 19:14:57  rwiegand
# Added a column for use in sorting.
#
# Revision 1.6  2006/03/27 18:15:09  rwiegand
# When an observation has been matched to a catalog object, populate the
# XXX_DELTA columns based on distance of detection from XXX position.
#
# Revision 1.5  2006/02/22 15:43:46  rwiegand
# Bound searches on catalogs.
#
# Revision 1.4  2006/01/25 21:27:11  rwiegand
# Implemented radius filtering within Catalog::contents.
#
# Revision 1.3  2006/01/25 15:23:46  rwiegand
# Extended SearchCat module to support many scat catalag interfaces.  Allow
# user to pass filter and epoch into catalog loading.  Use Wayne's magnitude
# estimation for UVOT filters when loading USNO B1.
#
# Revision 1.2  2005/12/01 19:43:59  rwiegand
# Report the distance to the nearest available catalog entry in the
# CAT_DELTA column.
#
# Revision 1.1  2005/10/07 20:42:32  rwiegand
# Source catalog modules.
#
# Revision 1.18  2005/04/04 18:12:00  rwiegand
# Added specs for coincidence loss corrected source rate and error.
#
# Revision 1.17  2005/03/21 14:29:40  rwiegand
# Display two more digits of RA and DEC.
#
# Revision 1.16  2005/03/11 22:08:31  rwiegand
# Write NULL MAG_DELTA if catalog object magnitude is not known.
#
# Revision 1.15  2005/03/04 19:40:20  rwiegand
# Calculate a tighter bound on the sub-catalog to load.  Make use of the
# partition description when creating the catalog.
#
# Revision 1.14  2005/02/08 21:48:43  rwiegand
# Added primary sort of observed version catalog sources followed by
# sub-sort by distance from reference position.
#
# Revision 1.13  2005/02/08 21:01:46  rwiegand
# Was failing to pass NULL value when reading column data.
#
# Revision 1.12  2004/12/09 21:42:52  rwiegand
# Switched errors in RA and Dec from degrees to arcsec.
#
# Revision 1.12  2004/12/09 21:21:28  rwiegand
# Switched errors in RA and Dec from variance in degrees to linear arcsec.
#
# Revision 1.11  2004/12/07 22:27:33  rwiegand
# Changed area of background rate to arcsec^2 instead of pixel.
#
# Revision 1.10  2004/12/07 21:18:53  rwiegand
# Renamed background column to indicate rate and updated rate column units.
#
# Revision 1.9  2004/11/03 20:34:21  rwiegand
# Store all source list column specifications locally.  Support for
# propagating columns.
#
# Revision 1.8  2004/06/17 21:17:35  rwiegand
# Use REFID which is harder for user to screw with than row number.
#
# Revision 1.7  2004/05/10 14:10:12  rwiegand
# User keywords are handled externally.  Support passing the FITS object back
# to the caller.
#
# Revision 1.6  2004/05/03 21:54:13  rwiegand
# Allow caller to specify units [radians or degrees].
#
# Revision 1.5  2004/04/30 14:59:39  rwiegand
# Added REFID and STARID columns.
#
# Revision 1.4  2004/04/22 20:55:59  rwiegand
# Initialize to an empty catalog if one is not provided.
#
# Revision 1.3  2004/02/13 22:37:55  rwiegand
# Updated columns in output source list.  Allow caller to specify a list of
# reference positions for which the separation of each source will be reported.
#
# Revision 1.2  2004/01/05 16:44:56  rwiegand
# Express RA and DEC in degrees instead of radians.
#
# Revision 1.1  2003/11/26 20:00:19  rwiegand
# Modules for sharing among tasks.
#

use strict;

package StarID::SourceTable;

use base qw(Task::FITS);
use Task qw(:codes);

use Astro::FITS::CFITSIO qw(:constants);

use StarID::Source;
use StarID::Catalog;

use constant EXTNAME => 'SOURCES';


# we want to sort observed objects before catalog objects
use constant CLASS_OBSERVED   => 1;
use constant CLASS_REFERENCE  => 2;

use constant REF_NULL => -1;
use constant POS_NULL => -1;
use constant SKY_NULL => -999;
use constant MAG_NULL => 99;
use constant TYPE_NULL => '';
use constant ID_NULL => '';


my @PRE_COLUMNS = (
	{ name => 'RA',
		form => 'D',
		comment => 'Right Ascension',
		units => 'deg',
		display => 'F11.6',
	},
	{ name => 'DEC',
		form => 'D',
		comment => 'Declination',
		units => 'deg',
		display => 'F11.6',
	},
	{ name => 'RA_ERR',
		form => 'D',
		comment => 'Error in RA',
		units => 'arcsec',
		display => 'F8.4',
		# display => 'E13.6',
		null => SKY_NULL,
	},
	{ name => 'DEC_ERR',
		form => 'D',
		comment => 'Error in DEC',
		units => 'arcsec',
		display => 'F8.4',
		# display => 'E13.6',
		null => SKY_NULL,
	},
	{ name => 'REFID',
		form => 'I',
		comment => 'Source number',
		null => REF_NULL,
	},
	{ name => 'STARID',
		form => '16A',
		comment => 'Catalog identifier',
		null => ID_NULL,
	},
	{ name => 'MAG',
		form => 'E',
		units => 'mag',
		comment => 'Magnitude',
		null => MAG_NULL,
		display => 'F6.3',
	},
	{ name => 'MAG_ERR',
		form => 'E',
		units => 'mag',
		comment => 'Error in MAG',
		null => MAG_NULL,
		display => 'F6.3',
		# display => 'E11.4',
	},
);

my @POST_COLUMNS = (
	{ name => 'CAT_DELTA',
		form => 'E',
		units => 'arcsec',
		null => POS_NULL,
		comment => 'Distance from catalog object',
		display => 'F6.2',
	},
	{ name => 'MAG_DELTA',
		form => 'E',
		comment => 'Magnitude difference from catalog object',
		null => MAG_NULL,
		display => 'F6.2',
	},
	{ name => 'TYPE',
		form => '12A',
		# catalog extended, binary, ...
		null => TYPE_NULL,
		comment => 'Object type: point, extended',
	},
	{ name => 'PACKET',
		form => '12A',
		comment => 'Source packet indicator',
		null => TYPE_NULL,
	},
	{ name => 'RATE',
		form => 'E',
		comment => 'Source count rate',
		units => 'count/s',
		null => MAG_NULL,
		display => 'F8.4',
	},
	{ name => 'RATE_ERR',
		form => 'E',
		comment => 'Error in RATE',
		units => 'count/s',
		null => MAG_NULL,
		display => 'F6.4',
	},
	{ name => 'RATE_BKG',
		form => 'E',
		comment => 'Background count rate',
		units => 'count/s/arcsec**2',
		null => MAG_NULL,
		display => 'F6.4',
	},
	{ name => 'RATE_COR',
		form => 'E',
		comment => 'Coincidence loss corrected rate',
		units => 'count/s',
		null => MAG_NULL,
		display => 'F8.4',
	},
	{ name => 'RATE_COR_ERR',
		form => 'E',
		comment => 'Error in RATE_COR',
		units => 'count/s',
		null => MAG_NULL,
		display => 'F6.4',
	},
	{ name => 'SATURATED',
		form => 'L',
		comment => 'Saturated source',
	},
	{ name => 'FLUX',
		form => 'E',
		comment => 'Flux',
		units => 'erg/s/cm**2/angstrom',
		null => MAG_NULL,
		display => 'E13.6',
	},
	{ name => 'FLUX_ERR',
		form => 'E',
		comment => 'Error in FLUX',
		units => 'erg/s/cm**2/angstrom',
		null => MAG_NULL,
		display => 'E13.6',
	},
	{ name => 'X_IMAGE',
		form => 'E',
		comment => 'Object position along x',
		units => 'pixel',
		null => POS_NULL,
		display => 'F8.2',
	},
	{ name => 'X_ERR',
		form => 'E',
		comment => 'Error in X_IMAGE',
		units => 'pixel',
		null => POS_NULL,
		display => 'F8.4',
	},
	{ name => 'Y_IMAGE',
		form => 'E',
		comment => 'Object position along y',
		units => 'pixel',
		null => POS_NULL,
		display => 'F8.2',
	},
	{ name => 'Y_ERR',
		form => 'E',
		comment => 'Error in Y_IMAGE',
		units => 'pixel',
		null => POS_NULL,
		display => 'F8.4',
	},
	{ name => 'PROF_MAJOR',
		form => 'E',
		comment => 'Profile RMS along major axis',
		units => 'arcsec',
		null => SKY_NULL,
		display => 'F6.2',
	},
	{ name => 'PROF_MINOR',
		form => 'E',
		comment => 'Profile RMS along minor axis',
		units => 'arcsec',
		null => SKY_NULL,
		display => 'F6.2',
	},
	{ name => 'PROF_THETA',
		form => 'E',
		comment => 'Position angle (CCW/world-x)',
		units => 'deg',
		null => SKY_NULL,
		display => 'F8.2',
	},
	{ name => 'THRESHOLD',
		form => 'E',
		comment => 'Detection threshold above background',
		units => 'sigma',
		display => 'F8.4',
		null => POS_NULL,
	},
	{ name => 'FLAGS',
		form => 'J',
		comment => 'Extraction flags',
		null => POS_NULL,
	},
	{ name => 'ORDER',
		form => 'J',
		comment => 'Sort order',
		null => POS_NULL,
	},
);


my @COLUMNS = (@PRE_COLUMNS, @POST_COLUMNS);
my %COLUMNS = map { $_->{name} => $_ } @COLUMNS;



sub getSpec
{
	my ($task, $colname) = @_;

	my $spec = undef;

	if (exists($COLUMNS{$colname})) {
		my %spec = %{ $COLUMNS{$colname} };
		$spec = \%spec;
	}
	elsif ($colname =~ /^(\w+)_DELTA$/) {
		my $tag = $1;
		$spec = {
			name => $tag . '_DELTA',
			form => 'E',
			units => 'arcsec',
			null => POS_NULL,
			comment => "Distance from $tag position",
			display => 'F6.2',
		};
	}
	else {
		$task->warning("StarID::SourceTable: invalid column '$colname'");
	}

	return $spec;
}


sub execute
{
	my ($self, %args) = @_;

	$self->validate;

	$self->writeSourceList(%args)
		if $self->isValid;
}


sub validate
{
	my ($self) = @_;

	# make sure all the options are set
	$self->{unobserved} ||= 0;

	$self->{reference} ||= [ ];
	$self->{catalog} ||= StarID::Catalog->new;

	my $count = @{ $self->{sources} };

	if (@{ $self->{reference} } > 0) {
		$self->{target} = $self->{reference}[0];
	}
	elsif ($count > 0) {
		# this can fail if RA crosses 0/360
		my $raSum = 0;
		my $decSum = 0;
		$self->applyToSources(sub {
				my ($s) = @_;
				$raSum += $s->ra;
				$decSum += $s->dec;
			});

		$self->{target} = StarID::Source->new(
				RA => $raSum / $count,
				DEC => $decSum / $count,
				);
	}

	if (my $target = $self->{target}) {
		$target->{UNIT} = Math::rd2unit($target->ra, $target->dec);
	}
	else {
		$target->{UNIT} = [ 1, 0, 0 ];
	}
}


sub determinePacketIndicator ($)
{
	my ($s) = @_;

	my $type = TYPE_NULL;

	if ($s->{stamp}) {
		$type = $s->{stamp};
	}
	elsif ($s->{binning}) {
		$type = $s->{binning};
	}

	return $type;
}


sub applyToSources
{
	my ($self, $method) = @_;

	foreach my $s (@{ $self->{sources} }) {
		if (ref($method) eq 'CODE') {
			$method->($s);
		}
		else {
			$self->$method($s);
		}
	}
}


sub isNumeric
{
	my ($s) = @_;
	return 0 if not defined($s);
	if ($s =~ /^\d+(\.\d*)?(e[-+]?\d+)?$/) {
		return 1;
	}
	if ($s =~ /^\d*(\.\d+)?([eE][-+]?\d+)?$/) {
		return 1;
	}
	return 0;
}


sub putSourceList
{
	my ($self, $fits, %args) = @_;

	my $status = 0;
	my @entries;
	my %matched;

	# base set of columns to include in the output
	my @known = $args{columns} ? @{ $args{columns} } : ();

	# add those computed locally
	push(@known, qw(REFID STARID PACKET RA DEC MAG TYPE CAT_DELTA MAG_DELTA ORDER));

	# add entries for observed sources
	my $matchrad = Math::degreesToRadians($self->{matchtol} / 3600);

	my @refColumns;
	foreach my $p (@{ $self->{reference} }) {

		my $name = uc($p->id);
		$p->{colname} = $name . '_DELTA';
		push(@known, $p->{colname});

		push(@refColumns, {
				name => $p->{colname},
				computed => 1,
				form => 'E',
				units => 'arcsec',
				null => POS_NULL,
				comment => "Distance from $name position",
				display => 'F6.2',
			},
			);
	}

	my @allColumns = (@PRE_COLUMNS, @refColumns, @POST_COLUMNS);
	my %known = map { $_ => 1 } @known;
	my $catSize = $self->{catalog}->size;

	$self->applyToSources(sub {
		my ($s) = @_;

		my @delta;

		my $packet = $s->{PACKET} || TYPE_NULL;

		my $matches = 0;
		# find which objects in the catalog it is within sepangle of

		foreach my $o (@{ $self->{catalog}->contents($s, $matchrad) }) {
			++$matches;
			$matched{$o} = 1;
			$s->{MATCHED} = 1;
			foreach my $p (@{ $self->{reference} }) {
				my $delta = $s->angle($p);
				push(@delta, $p->{colname} => $delta);
			}
			my $delta = isNumeric($o->mag) ? $s->mag - $o->mag : MAG_NULL;
			push(@entries, {
					REFID => $s->{REFID},
					STARID => $o->id,
					PACKET => $packet,
					RA => $s->ra,
					DEC => $s->dec,
					MAG => $s->mag,
					TYPE => $o->type,
					CAT_DELTA => $s->angle($o),
					MAG_DELTA => $delta,
					@delta,
				});
		}

		if (not $matches) {
			foreach my $p (@{ $self->{reference} }) {
				my $delta = $s->angle($p);
				push(@delta, $p->{colname} => $delta);
			}

			# find nearest catalog object
			my $maxcos = -2;
			my $scale = 1;
			for (my $i = 0; $maxcos == -2 and $i < 4; ++$i) {
				$scale *= 3;
				my $list = $self->{catalog}->contents($s, $scale * $matchrad);
				foreach my $o (@$list) {
					my $cosangle = $s->cosangle($o);
					if ($cosangle > $maxcos) {
						$maxcos = $cosangle;
					}
				}
				if (@$list == $catSize) {
					# explored all objects in catalog
					last;
				}
			}

			my $catDelta = ($maxcos == -2) ? POS_NULL : POSIX::acos($maxcos);

			push(@entries, {
					REFID => $s->{REFID},
					STARID => ID_NULL,
					PACKET => $packet,
					RA => $s->ra,
					DEC => $s->dec,
					MAG => $s->mag,
					TYPE => TYPE_NULL,
					CAT_DELTA => $catDelta,
					MAG_DELTA => MAG_NULL,
					@delta,
				});
		}
	});

	# report unobserved sources
	my $args = $self->args;

	my @unmatched;
	if ($self->{unobserved} > 0 and $self->{target}) {

		my $target = $self->{target};

		my $radians = Math::degreesToRadians($self->{unobserved} / 3600);
		my $uncos = cos($radians);

		foreach my $o (@{ $self->{catalog}->contents }) {

			next if $matched{$o};

			next if $target->cosangle($o) < $uncos;

			my @delta;

			foreach my $p (@{ $self->{reference} }) {
				my $delta = $p->angle($o);
				push(@delta, $p->{colname} => $delta);
			}

			push(@unmatched, {
					REFID => REF_NULL,
					STARID => $o->id,
					PACKET => TYPE_NULL,
					RA => $o->ra,
					DEC => $o->dec,
					MAG => isNumeric($o->mag) ? $o->mag : MAG_NULL,
					TYPE => $o->type,
					CAT_DELTA => POS_NULL,
					MAG_DELTA => MAG_NULL,
					@delta,
				});
		}
	}

	# only keep brightest unmatched entries
	if ($self->{uncount} > 0) {
		my @sorted = sort { $a->{MAG} <=> $b->{MAG} } @unmatched;

		my $limit = $self->{uncount};
		foreach my $o (@sorted) {
			if ($o->{MAG} < $self->{unmag}) {
				push(@entries, $o);
				last if --$limit == 0;
			}
			last if $o->{MAG} > $self->{unmag};
		}
	}

	my $targunit = $self->{target}{UNIT};

	foreach my $e (@entries) {

		my $unit = Math::rd2unit($e->{RA}, $e->{DEC});

		$e->{RA} = Math::radiansToDegrees($e->{RA});
		$e->{DEC} = Math::radiansToDegrees($e->{DEC});

		$e->{distance} = Math::u3angle($unit, $targunit);

		# convert separations to arcsecs
		$e->{CAT_DELTA} = 3600 * Math::radiansToDegrees($e->{CAT_DELTA})
			if $e->{CAT_DELTA} != POS_NULL;

		foreach my $r (@{ $self->{reference} }) {
			my $key = $r->{colname};
			$e->{$key} = 3600 * Math::radiansToDegrees($e->{$key})
				if $e->{$key} != POS_NULL;
		}
	}

	my %sources = map { $_->{REFID} => $_ } @{ $self->{sources} };
	foreach my $e (@entries) {
		# make other fields available
		if ($e->{REFID} == REF_NULL) {
			$e->{class} = CLASS_REFERENCE;
			next;
		}
		else {
			$e->{class} = CLASS_OBSERVED;
		}
		my $s = $sources{$e->{REFID}};
		foreach my $c (@PRE_COLUMNS, @POST_COLUMNS) {
			my $tag = $c->{name};
			if (not defined($e->{$tag})) {
				if (defined($s->{$tag})) {
					$e->{$tag} = $s->{$tag};
				}
				else {
					$e->{$tag} = $c->{null};
				}
			}
		}
	}

	# sort entries by class and then distance from target
	my @sorted = sort {
			$a->{class} <=> $b->{class}
			or $a->{distance} <=> $b->{distance}
			or $a->{STARID} <=> $b->{STARID}
			or $a->{REFID} <=> $b->{REFID}
		} @entries;
	{
		my $i = 0;
		foreach my $s (@sorted) {
			$s->{ORDER} = ++$i;
		}
	}

	foreach my $c (@allColumns) {
		next if not $known{$c->{name}};
		my @data = map { $_->{$c->{name}} } @entries;
		$self->writeColumn($fits, $c, \@data);
	}
}


sub writeSourceList
{
	my ($self, %args) = @_;

	my $path = $self->{outfile};
	my $extname = $args{extname} || $self->{extname} || EXTNAME;

	my $status = 0;
	my $fits = Astro::FITS::CFITSIO::create_file($path, $status);
	if (not $fits) {
		$self->error(BAD_OUTPUT,
			"unable to create $path: $!");
		return;
	}

	if ($fits->create_tbl(BINARY_TBL, 0, 0, 0, 0, 0, $extname, $status)) {
		$self->error(BAD_OUTPUT,
			"unable to create binary table");
	}

	if ($args{scale} eq 'degrees') {
		my $degreesToRadians = Math::degreesToRadians(1);
		$self->applyToSources(sub {
				my ($s) = @_;
				$s->{RA} *= $degreesToRadians;
				$s->{DEC} *= $degreesToRadians;
			});
	}

	$self->putSourceList($fits, %args)
		if $self->isValid;

	if ($args{scale} eq 'degrees') {
		my $radiansToDegrees = Math::radiansToDegrees(1);
		$self->applyToSources(sub {
				my ($s) = @_;
				$s->{RA} *= $radiansToDegrees;
				$s->{DEC} *= $radiansToDegrees;
			});
	}

	if ($args{leaveOpen}) {
		$self->{fits} = $fits;
	}
	else {
		if ($fits->close_file($status)) {
			$self->error(BAD_OUTPUT,
				"unable to close source list");
		}
	}

	$self->report("wrote source list to $path")
		if $self->isValid;

}


sub loadFile
{
	my ($self, $path, %args) = @_;

	my $status = 0;
	my $fits = Astro::FITS::CFITSIO::open_file($path, READONLY, $status);
	if ($status) {
		$self->error(BAD_INPUT, "unable to open $path [$status]");
		return;
	}

	if ($args{extname}) {
		if ($fits->movnam_hdu(BINARY_TBL, $args{extname}, 0, $status)) {
			$self->error(BAD_INPUT,
					"unable to move to extension '$args{extname}' [$status]");
		}
	}

	$self->loadFITS($fits, %args)
		if $self->isValid;

	if ($args{leaveOpen} and $fits) {
		$self->{fits} = $fits;
	}
	elsif ($fits) {
		my $tmp = 0;
		$fits->close_file($tmp);
	}

	return $self->{sources};
}


sub fits_read_method
{
	my ($self, $spec) = @_;
	my $form = $spec->{form};
	my $typ = 'dbl';
	if ($form =~ /A/) {
		$typ = 'str';
	}
	elsif ($form =~ /[IJ]/) {
		$typ = 'lng';
	}
	elsif ($form =~ /[ED]/) {
		$typ = 'dbl';
	}
	elsif ($form =~ /L/) {
		$typ = 'log';
	}
	else {
		$self->warning("fits_read_method: bad form '$form'");
	}

	return 'read_col_' . $typ;
}


sub loadFITS
{
	my ($self, $fits, %args) = @_;

	my @names;
	if ($args{columns}) {
		@names = @{ $args{columns} };
	}
	else {
		@names = map { $_->{name} } @COLUMNS;
	}
	my @columns;
	foreach my $name (@names) {
		push(@columns, $self->getSpec($name));
	}

	my %columns = map { $_->{name} => $_ } @columns;

	foreach my $c (@columns) {
		$c->{required} = $args{columns} || $c->{required};
	}

	my $count;
	my $status = 0;
	if ($fits->get_num_rows($count, $status)) {
		$self->error(BAD_INPUT,
			"unable to get number of sources [$status]");
	}
	else {
		$self->report("loading $count sources");
	}

	my @loaded;

	foreach my $c (@columns) {

		my ($colnum, $anynul, $skip) = (0) x 3;
		my @data;
		my $method = $self->fits_read_method($c);
		my $null = exists($c->{null}) ? $c->{null} : 0;

		Astro::FITS::CFITSIO::fits_write_errmark();
		if ($fits->get_colnum(CASEINSEN, $c->{name}, $colnum, $status)) {
			if (not $c->{required}) {
				$skip = 1;
				# $self->warning("trouble getting $c->{name} [$status]");
				$status = 0;
			}
			else {
				$self->error(BAD_INPUT,
				"unable to get source list $c->{name} column index [$status]");
			}
		}
		Astro::FITS::CFITSIO::fits_clear_errmark();
		last if not $self->isValid;

		next if $skip;

		if ($fits->$method($colnum, 1, 1, $count, $null,
					\@data, $anynul, $status)) {
			$self->error(BAD_INPUT,
				"unable to read source list $c->{name} column data");
		}
		else {
			$c->{data} = \@data;
		}

		last if not $self->isValid;

		push(@loaded, $c->{name});
	}

	my @sources;

	if ($self->isValid) {

		for (my $i = 0; $i < $count; ++$i) {
			my %source = map { $_ => $columns{$_}{data}[$i] } @loaded;
			$source{ROW} = $i + 1;
			push(@sources, StarID::Source->new(%source));
		}
	}

	$self->{sources} = \@sources;

	if ($args{scale} eq 'radians') {
		my $degreesToRadians = Math::degreesToRadians(1);
		$self->applyToSources(sub {
				my ($s) = @_;
				$s->{RA} *= $degreesToRadians;
				$s->{DEC} *= $degreesToRadians;
			});
	}

	$self->{loaded} = \@loaded;

	return $self->{sources};
}


1;

