#!/usr/local/bin/perl

($fld,$x,$y) =@ARGV;

$pi = 3.14159265358979;
#####################################################################
#--------------Read VImaps----------------------------------------------
%EVI = %EVIe =%AV= %AVe= %AI=%AIe=();
open("DATA", "zcat EVImaps/sc$fld.map.gz |");
$nfread=0;
while($line=<DATA>){
  ($nf,$xi,$yi,$EVI,$EVIe, $AV,$AVe, $AI,$AIe,$Nall,$Nallcut,$VIsdevall,
                           $NRC,$NRCcut,$VIsdevRC,$gi) =split(" ", $line,);
  ($w,$xia,$yia,$w,$VIm,$VIme,$w,$Im,$Ime,$w,$$Vm,$Vme)=split(" ",$line)
                                                        if $line=~/ALL/;
  ($w,$binsize,$w,$nfgood,$npoor,$w,$nall,$sumnall,$w,$ncut,$sumncut,$w,$VIallsdev)
                                     =split(" ", $line)if $line =~ /binsize/;
  ($w,$EVI0e,$w,$AV0e,$w,$AI0e)=split(" ",$line)if $line=~/AV0err/;

  next if $line =~ /#/;
  $EVI{$nf}  = $EVI;
  $EVIe{$nf} = $EVIe;
  $AV{$nf}  = $AV;
  $AVe{$nf} = $AVe;
  $AI{$nf}  = $AI;
  $AIe{$nf} = $AIe;
  $gi{$nf} = $gi;
  $nfread++;
}
$Nbin=$xia;
$binsize  = 2048/$Nbin;
$ximax = $Nbin -1;
$yimax = $Nbin *4 -1;
#----------------interpolation-----------------------------------------
$xi =  int($x/ $binsize);
$yi =  int($y/ $binsize);
$x0 = $x/$binsize ;
$y0 = $y/$binsize ;
$nf0 = $xi+$yi*($ximax+1);

if ($xi<0 || $xi>$ximax || $yi<0 || $yi>$yimax){
             print "Out of boundary ($x,$y) X=$xi Y=$yi!!\n"; exit;}
if ($EVI{$nf0}==0){print "NO map at ($x,$y) X=$xi Y=$yi!!\n"; exit;}

#print "Nbin=$Nbin  binsize=$binsize   $x  $y    X=$xi  Y=$yi   x0=$x0  y0=$y0\n";

printf ("%2d %2d   %.4f %.4f  %.4f %.4f  %.4f %.4f   %.4f  %.4f  %.4f  %d :ACTUAL\n",
$xi, $yi, $EVI{$nf0},$EVIe{$nf0},$AV{$nf0},$AVe{$nf0},$AI{$nf0},$AIe{$nf0},$EVI0e,$AV0e,$AI0e,$gi{$nf0});

$EVI=spline2($Nbin,$x,$y,\%EVI);
$AV=spline2($Nbin,$x,$y,\%AV);
$AI=spline2($Nbin,$x,$y,\%AI);

printf ("%2d %2d   %.4f %.4f  %.4f %.4f  %.4f %.4f   %.4f  %.4f  %.4f  %d :INTERP\n",
$xi, $yi, $EVI,$EVIe{$nf0},$AV,$AVe{$nf0},$AI,$AIe{$nf0},$EVI0e,$AV0e,$AI0e,$gi{$nf0});
##########################################################
sub spline2 {
  local ($Nbin,$x,$y,*Map) = @_;

  my $Deg =3;              # degree (number of bin) for spline
  my $Dhf =($Deg-1)/2;

  my $Nbin=$xia;
  my $binsize  = 2048/$Nbin;
  my $ximax = $Nbin -1;
  my $yimax = $Nbin *4 -1;
  my $xi =  int($x/ $binsize);
  my $yi =  int($y/ $binsize);
  my $x0 = $x/$binsize ;
  my $y0 = $y/$binsize ;

  my @X = my @Y = my @MAP = ();
  my $nx = my $ny = 0;
  for (my $j=-$Dhf;$j<=$Dhf;$j++){  # for y
     my $yt = $yi +$j;
     $Y[$ny] = $yt +0.5;
     $ny++;
     if (my $yt <0 || $yt >$yimax){$yt=$yi;}
     $nx=0;
     for ($i=-$Dhf; $i<=$Dhf; $i++){  # for x
        my $xt = $xi +$i;
        $X[$nx] = $xt +0.5;
        $nx++;
        if ($xt <0 || $xt >$ximax){$xt=$xi;}
        my $nf = $xt+$yt*($ximax+1);
        $nf = $xi+$yi*($ximax+1) if $Map{$nf} ==0;
        push(@MAP,$Map{$nf});
     }
  }

  #print "Nbin=$Nbin  binsize=$binsize   $x  $y    X=$xi  Y=$yi   x0=$x0  y0=$y0\n";
  @line=`bin/xsplin2 $x0 $y0 $nx $ny @X @Y @MAP`;
#  print @line;
  ($w,$x00,$w,$y00,$w,$value)=split(" ",$line[17]);

  return $value;
}
##########################################################
sub numerically { $a <=> $b; }
sub descending { $b <=> $a; }