#!/usr/bin/perl
#
use Statistics::OLS;
use Math::FFT;
use Math::Stat;
use Math::Spline qw(spline linsearch binsearch);
use Math::Derivative qw(Derivative2);
use Math::Approx;
use Astro::FITS::CFITSIO qw( :longnames :constants );
use PDL;
use PDL::Fit::Polynomial; 
use PDL::Filter::Linear;
use PGPLOT;  # Load PGPLOT module
use PDL::Fit::Gaussian;
use PDL::Slatec;
use PDL::Image2D;
#use PDL::Matrix;

#require "/work1/ssa/perl/MY/my.pl";
require "/home/sanchez/sda2/code/R3D/my.pl";


if ($#ARGV<2) {
    print "USE: clean_Ha_map.pl VEL.FITS MASK.fits OUTPUT [VEL_MAX]\n";
    exit;
}

$infile1=$ARGV[0];
$infile2=$ARGV[1];
$output=$ARGV[2];
$vel_max=300;
if ($#ARGV==3) {
    $vel_max=$ARGV[3];
}

$a_in1=rfits($infile1);
$a_in2=rfits($infile2);

($nx,$ny)=$a_in1->dims();

@stats=stats($a_in1);
print "@stats\n";
$vel_last=$stats[0];

$nx2=int($nx/2);
$ny2=int($ny/2);
for ($i=0;$i<$nx;$i++) {
    for ($j=$ny2;$j<$ny;$j++) {
	$mask=$a_in2->at($i,$j);
	$vel=$a_in1->at($i,$j);
	if (abs($val-$stats[2])>$vel_max) {
	    set($a_in1,$i,$j,$vel_last);
	} else {
	    $vel_last=$a_in1->at($i,$j);
	}
    }
    for ($j=$ny-1;$j>0;$j--) {
	$mask=$a_in2->at($i,$j);
	$vel=$a_in1->at($i,$j);
	if (abs($val-$stats[2])>$vel_max) {
	    set($a_in1,$i,$j,$vel_last);
	} else {
	    $vel_last=$a_in1->at($i,$j);
	}

#6562.817    H\ga           
#6583.600    [NII]6584      
    }
}


$a_in1->wfits($output);

exit;

#ve.UGC03995.mask_map.fits
#ve.UGC03995.vel_map.fits
# Emission lines
#3727       [OII]3727_blended      
#4101.0    H\gd           
#4340.468    H\gg 
#4861.32    H\gb           
#4959  [OIII]
#5007 [OIII]
#6562.817    H\ga           
#6583.600    [NII]6584      
#6548.100    [NII]6548      
#6716.47    [SII]6717      
#6730.85    [SII]6731