#!/usr/bin/env perl

use PDL;
use PDL::NiceSlice;
use Ska::HashTable;
use PDL::Graphics::PGPLOT::Window;
use Data::Dumper;
use PDL::Fit::Polynomial;
use PDL::Slatec;

$max_field_dr = 2.0;		# Arcmin

# Read in select entries ("best" data points) from celmon output
$d = Ska::HashTable->new('DATA.rdb');
while (%d = $d->row("field_dr < $max_field_dr && status eq '' && " .
		     "detector eq 'ACIS-S' && " .
		     "( pos_ref =~ /icrs/i || pos_ref =~ /tycho/i || pos_ref =~ /simbad_high/i )")) {
    $df{"$d{obsid}___$d{object}"} = { %d };
}

$col = 2;
map { $color{$df{$_}->{fids}} = $col++ unless $color{$df{$_}->{fids}}} keys %df;

map { push @dr, $df{$_}->{dr};
      push @simz, $df{$_}->{sim_z};
      push @dy, $df{$_}->{dy};
      push @dz, $df{$_}->{dz};
      push @color, $color{$df{$_}->{fids}};
  } keys %df;

$dr = pdl @dr;
$dy = pdl @dy;
$dz = pdl @dz;
$simz = pdl @simz;

print "Number of data points: ", $dr->nelem, "\n";

$win = PDL::Graphics::PGPLOT::Window->new({Device => 'simz.ps/cps',
					   Device => '/xs',
					   WindowWidth => '7',
					   AspectRatio => '1.4',
					   NYPanel => 2,
					   });

%pgopt = (CHARSIZE => 1.7);

plot_axis($dy, "DY");
plot_axis($dz, "DZ");

$win->close();

sub plot_axis {
    my $delta = shift;
    my $txt = shift;
    my $sim_z0 = 190.5;

    $win->env( -195, -178, -1.3, 1.3, {%pgopt});
    foreach (0..$delta->nelem-1) {
	$rnd = (random(1)-0.5)*0.5;
	$win->points($simz[$_]+$rnd, $delta->at($_),
#		     {COLOR => $color[$_]}
		    );
    }
    $win->label_axes('SIM-Z (mm)', "$txt offset (arcsec)", {%pgopt});

    my $xfit = sequence(30) - 10 - $sim_z0;
    my ($yfit, $coeffs) = fitpoly1d($simz + $sim_z0, $delta, 3);
    my $out = poly_eval($xfit + $sim_z0, $coeffs);
    $win->line($xfit, $out, {color => 'RED'});
    print "$txt coeffs: ", $coeffs/20, " (mm)\n";
}


sub poly_eval {
    my ($x, $coef) = @_;
    my $y = ones($x) * $coef(0);
    for my $i (1 .. $coef->nelem-1) {
	$y += $coef($i) * $x**$i;
    }
    return $y;
}