use PDL;
use PDL::Transform;
use PDL::NiceSlice;
use PDL::Image2D;
use PDL::Transform;
use PDL::IO::Dumper;
use PDL::Minuit;
use PDL::GSLSF::LEGENDRE;
use Astro::Coord;
use Astro::Time;
use Data::Dumper;
#use Ftools;
use TrIm;
use Vtools;
use strict;
use warnings;

my $pi = 2*acos(0);

my $h_pass;

#
my $min_x;
my $max_x; 
my $min_y; 
my $max_y; 
my $PSF_FWHM;
my $nord_legendre_baseline = 2;
my $check_box = 1;
my $do_proj; 

my $fileout; # needs global lexical domain for wcs in astrogoto

die "usage:\n perl gotoresamp.pl NGC4755_I_150s+Dark.fit  init_resamp.dat/\n" unless (@ARGV); # give directory

my $filein = shift @ARGV;
my $datain = shift @ARGV;
my $fileout = $filein;
$fileout =~ s/.fit/_resamp.fits/i;

#
{

# INPUT IMAGE
    my $im = rfits($filein); # or die "$! cannot open $filein \n";
    my $h = $im -> gethdr;
    $im = float($im);


# GOTO field of view
    my $field_0 = 9.0 / 60.0; #field in deg
# GOTO pixel scale    
    my $pixscale = $field_0 / $im->getdim(0);
# estimated PSF FWHM
    $PSF_FWHM  =  10. ; #pixels
# SWITCH TO EXAMINE BOXES
    $check_box = 1;
 
   
## GIVE NAMES, REF RA DEC and guess pixel location
#    my @data = (
##		{Name => 'HD111904', RA  => '12:53:21.895', DEC => '-60:19:42.56', i_init => 260, j_init => 820}, #SATURATED
#		{Name => 'SAO252075', RA  => '12:53:47.258', DEC => '-60:19:55.23', i_init => 820, j_init => 800},
#		{Name => 'HD111990', RA  => '12:53:59.799', DEC => '-60:20:07.51', i_init => 1080, j_init => 780},
#		{Name => 'HD111934', RA  => '12:53:37.617', DEC => '-60:21:25.42', i_init => 615, j_init => 540},
##		{Name => 'SAO252073', RA  => '12:53:41.337', DEC =>  '-60:20:57.97', i_init => 695, j_init => 620}, #SATURATED
#		{Name => 'HD 312084', RA => '12:54:14.54', DEC => '-60:21:47.4', i_init => 1412, j_init => 501},
#		{Name => 'ALS 2800', RA => '12:53:20.70', DEC => '-60:23:16.7', i_init => 264, j_init => 199}
#		);
#    
#
    my @data;
    open IN, $datain; 
    my @entries = (<IN>);
    close IN;


# Define a center for the resampled image 
    shift @entries;
    my $centerline = shift @entries;
    
    my ($RA_0, $DEC_0) = split(/\s+/, $centerline); # =~ /([\d\.\:]+)\s+([\d\.\:]+)/);

# = '12:53:41.337'; my $DEC_0 = '-60:20:57.97'; #NGC 4755

# Define  dimensions for the resampled image 
    my $Nx = 1800;
    my $Ny = 1300;

    my $nord_x = 1;
    my $nord_y = 1;




    foreach my $entry (@entries) {
	next if ($entry =~ /^\#/);
	next if ($entry =~ /^\s*$/);
	my @F = split(/;/, $entry);
	my ($name, $RA, $DEC, $i, $j) = map {$_ =~ s/\s+//g; $_} @F[0..4];
	shift @F; 	shift @F; 
	shift @F; 
	shift @F; 
	shift @F; 
	my $force = 0;
	$force = shift @F if (@F); 
#	my $name = ($entry =~ /Name => .
#	push @data, {Name => $name, RA => $RA, DEC => $DEC, i_init => float($i), j_init => float($j), Force_Init => $force };
	push @data, {Name => $name, RA => $RA, DEC => $DEC, i_init => $i, j_init => $j, Force_Init => $force };
    }

    print ">>>> input data ";
    print  Dumper(\@data);



# ARE WE GOING GNOMONIC?
    $do_proj = 1; 
    
######################################################################
######################################################################
#  NO FURTHER EDITS BELOW
######################################################################
######################################################################


# get resampled image:
    my $out = &astrogoto($im,{Nx => $Nx, Ny => $Ny,  PixScale => $pixscale,  RA_0 => $RA_0, DEC_0=>  $DEC_0, Nord_x => $nord_x, Nord_y => $nord_y, Data => \@data});

    wfits $out,$fileout; #

}    




sub get_centroid {
    my ($im, $x0, $y0, $opts) = @_;

    my %args = ( Force_Init => 0, Do_Gaussfit => 0, Delta => 6, %$opts);
    
#    if (ref($opts) =~ /HASH/) {
#	foreach (keys(%$opts)) {
#	    $args{$_} = $$opts{$_};
#	}
#    }

    my $delta=$args{Delta};

    my $subim = $im($x0-$delta:$x0+$delta,$y0-$delta:$y0+$delta);

#    &Vtools::view($im);

    if ($check_box) {
#	print "$x0-$delta:$x0+$delta,$y0-$delta:$y0+$delta \n";
	print "check that star is comfortably contained by the $delta pixel box\n";
	&Vtools::view($subim);
    }

#    my ($max, $xcent,$ycent);
    
    my ($max,$xcent,$ycent) = max2d_ind($subim);

    if ($args{Force_init}) {
	$xcent = rint($subim->getdim(0)/2.);
	$ycent = rint($subim->getdim(1)/2.);
    } elsif ($args{Do_Bary}) {
	
	my $xxs = $subim -> xvals;
	my $yys = $subim -> yvals;
	
	print "MAX: xcent $xcent ycent $ycent \n" if $check_box;
	my $med = median($subim);
	$xcent = sum( $xxs * ($subim - $med ) / sum($subim - $med) ); 
	$ycent = sum( $yys * ($subim - $med ) / sum($subim - $med) ); 
	
	print "Bary: xcent $xcent ycent $ycent \n" if $check_box;
	
    } 
    
    if ($args{Do_Gaussfit}) {  
#	my $civs = pdl [ 0.4481183, 0.62394083, 0.014348385];
	my $sigma = $PSF_FWHM/(2.*sqrt(2.*log(2.)));
	
	my $cinv = pdl [ 1/$sigma, 1/$sigma, 0.]; # Sigma_X^-2, Sigma_Y^-2, excentricity  
#	print "cinv start $cinv \n";
	my $params = { Cinv => $cinv ,
		       Flux=> $max->sclr, 
		       Centroid=>pdl [ $xcent, $ycent]} ;
	my $mode;
	if ($check_box) {
	    $mode = 'viewsilent';
	} else {
	    $mode = 'silent';
	}
	my $bestfit = &Gauss::cholesky($subim,{Pars=>$params,LegPoly=>$nord_legendre_baseline,Mode=> $mode}) ; # LegPoly carries order for baseline

	
	(my $mux, my $muy, my $lmax, my $lmin, my $alpha, my $flux) = &Gauss::translatebestfit($bestfit);
	
	print "Centroid $mux $muy; Bmax/2 $lmax; Bmin/2 $lmin;\nPA $alpha; peak intensity $flux \n" if $check_box; 


    
	$xcent = $mux;
	
	$ycent = $muy;

	print "Gauss: xcent $xcent ycent $ycent \n" if $check_box;
    } 

    $xcent += $x0-$delta;
    $ycent += $y0-$delta;
	
    return ($xcent,$ycent);
}
    






sub astrogoto {
    my ($im,$pass) = @_;
    my $h = $im-> gethdr;
    my %args = (Nord_x=> 1, Nord_y => 1, Nx => 0, Ny => 0 , Data => 0, RA_0 =>,  DEC_0 => , PixScale => 0);

    if (defined($pass)) {
	%args = (%args,  %$pass);
    }
    my $pixscale = $args{PixScale}; 

    
# GIVE NAMES, REF RA DEC and guess pixel location
    my $datal = $args{Data}; 
    my @data = @$datal;
    my $RA_0 = $args{RA_0}; my $DEC_0 = $args{DEC_0};

# Define  dimensions for the resampled image
    my $Nx = $args{Nx}; 
    my $Ny = $args{Ny};



    $$h{'NAXIS'} = 2;
    $$h{'NAXIS1'} = $Nx;
    $$h{'NAXIS2'} = $Ny;

    $$h{'CRVAL1'} = sprintf("%8e",str2deg($RA_0,'H'));
    $$h{'CRVAL2'} = sprintf("%8e",str2deg($DEC_0,'D'));
    $$h{'CRPIX1'} = $Nx/2;
    $$h{'CRPIX2'} = $Ny/2;
    
    $$h{'RADESYS'} = 'FK5';
    $$h{'EQUINOX'} = '2000.0';
    $$h{'CTYPE1'} = 'RA---TAN';
    $$h{'CTYPE2'} = 'DEC--TAN';

    $$h{'CDELT1'} = sprintf("%8e",-1.0 * $pixscale);
    $$h{'CDELT2'} = sprintf("%8e",$pixscale);


    if ($do_proj) {
	my $dummy = zeroes($Nx,$Ny);
	$dummy -> sethdr($h);
#	print "fileout $fileout \n";
	wfits $dummy, $fileout;
    }

#    my $out = zeroes($im->dims);
    my $out = zeroes($Nx,$Ny);
    my ($Nx_in,$Ny_in) = $im->dims;
    $out(:$Nx_in-1,:$Ny_in-1) .= $im(:,:);
    $out ->sethdr($h);
    $h_pass = $h;

#my $t_out = t_fits($out);
    
    my @i_in;
    my @i_out;
    
    my @j_in;
    my @j_out;
    
    foreach (@data) {
	print "ASTRO REF : $_->{Name} \n"; 

	my $force_init = defined($_->{Force_Init})?$_->{Force_Init}:0;

	if ($force_init) {
	    print "FORCING INITIAL CONDITION FOR GAUSSIAN CENTROID AT THE CENTRE OF SUBIMAGE  \n";
	}
	my ($i,$j) = &get_centroid($im,$_->{i_init},$_->{j_init}, {Delta => 20, Do_Bary => 1, Do_Gaussfit => 1, Force_Init => $force_init});
	$_->{i} = $i;
	$_->{j} = $j;

#	my $in = pdl (str2deg($_->{RA},'H'), str2deg($_->{DEC},'D'));
#    my ($i_out, $j_out) = list ($in -> invert($t_out));
	

	my ($i_out, $j_out); # = &Ftools::adij_linear($out,$RA_OFF, str2deg($_->{DEC},'D'));

	if ($do_proj) {


	    my $t  = &CelCoord::wcs($fileout); # 

	    my $in_wcs  = pdl ( str2deg($_->{RA},'H'), str2deg($_->{DEC},'D') );
	    my $out_pix = $t->invert($in_wcs);
	    $i_out = $out_pix(0)->sclr;
	    $j_out = $out_pix(1)->sclr;
	    
#	    print "WCS: i_out $i_out j_out $j_out \n" ;


	} else {

	    my $RA_OFF = str2deg($RA_0,'H') + (str2deg($_->{RA},'H') - str2deg($RA_0,'H') ) * cos($pi*str2deg($DEC_0,'D') / 180.);
	    ($i_out, $j_out) = &Ftools::adij_linear($out,$RA_OFF, str2deg($_->{DEC},'D'));

	}

	$_->{i_out} = $i_out;
	$_->{j_out} = $j_out;


#    print "$i $j ---> $i_out $j_out \n";
	push @i_in,$i;
	push @i_out,$i_out;
	
	push @j_in,$j;
	push @j_out,$j_out;
	

    }

    my @dum = @i_out;

    push @dum, (0, $$h{NAXIS1});
    $min_x = minimum(pdl(@dum));
    $max_x = maximum(pdl(@dum));

    @dum = @j_out;
    push @dum, (0, $$h{NAXIS2});
    $min_y = minimum(pdl(@dum));
    $max_y = maximum(pdl(@dum));


 
    my $i_in = pdl (@i_in);
    my $i_out = pdl (@i_out);
#    wcols $i_in, $i_out;
    
#    print "i_out $i_out \n";
    


# need to specify total range of pixel coordinates to pass argumens as cos(theta) in Pl(cos(theta)) - may not be present in the fits header of a single position vector 
    my $opts = {Min_pixcoord_x => $min_x, Max_pixcoord_x => $max_x, Min_pixcoord_y => $min_y, Max_pixcoord_y => $max_y};
    $opts->{Nord_x} = $pass->{Nord_x};
    $opts->{Nord_y} = $pass->{Nord_y};

    my $bestfit = &TrIm::opt_2DTransform(\@data, $opts);

    print "INVERTING BEST FIT TRANSFORM, SANITY CHECK \n";

# CROSS CHECK: 
    foreach (@data) {
	my $in = pdl ( $_->{i}, $_->{j});
	my $out = pdl ( $_->{i_out}, $_->{j_out});
	my $model = &TrIm::Do_Transform($out,$bestfit,{Min_pixcoord_x => $min_x, Max_pixcoord_x => $max_x, Min_pixcoord_y => $min_y, Max_pixcoord_y => $max_y, Dir => -1});
#	print "$_->{Name} in $in  -> out $out -> model $model CRPIX1 $$h{CRPIX1} CRPIX2 $$h{CRPIX2} \n";
	
    }

    print "RESAMPLING GOTO DATA \n";
	
#    my $test_in = pdl (910, 430);
#    my $testout = &TrIm::Do_Transform($test_in,$bestfit,{Min_pixcoord_x => $min_x, Max_pixcoord_x => $max_x, Min_pixcoord_y => $min_y, Max_pixcoord_y => $max_y, Dir =>  -1, Debug => 0});
#    print "testin $test_in  ---> testout $testout   \n";

    $out = &TrIm::Do_Transform($out,$bestfit, {Min_pixcoord_x => $min_x, Max_pixcoord_x => $max_x, Min_pixcoord_y => $min_y, Max_pixcoord_y => $max_y, Dir => +1, Debug => 0});
    $out = float($out);
    $out ->sethdr($h);


    return $out;
}
