#!/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;


require("/home/sanchez/sda2/code/R3D/my.pl");

if ($#ARGV<10) {
    print "USE: kin_back_cube_fixed_error.pl cube.fits error.cube.fits[1] config_file back_file nback start_wavelength end_wavelength out_file DEVICE GUIDED_VEL_MAP.fits GUIDED_DISP.fits [ASK]\n";
    exit;
}

$input_cube=$ARGV[0];
$input_error=$ARGV[1]."[1]";
$config_file=$ARGV[2];
$back_file=$ARGV[3];
$nback=$ARGV[4];
$start_w=$ARGV[5];
$end_w=$ARGV[6];
$out_file=$ARGV[7];
$device=$ARGV[8];
$guided_map=$ARGV[9];
$guided_disp=$ARGV[10];
$ask=0;
if ($#ARGV==11) {
    $ask=$ARGV[11];
}

$nmask=0;
$gmap=rfits($guided_map);
$gdisp=rfits($guided_disp);

$nf=0;
open(FH,"<$back_file");
while($line=<FH>) {
    chop($line);
    if ($line !~ "#") {
	$file[$nf]=$line;
	$nf++;
    }
}
close(FH);

for ($j=0;$j<$nf;$j++) {
    open(FH,"<$file[$j]");
    $nx=0;
    while($line=<FH>) {
	chop($line);
	@data=split(" ",$line);
	if ($j==0) {
	    if ($nx==0) {
		$crval=$data[1];
	    }
	    if ($nx==1) {
		$cdelt=$data[1]-$crval;
	    }
	}
	$id[$nx]=$data[0];
	$w[$nx]=$data[1];
	$f[$nx]=$data[2];
	$nx++;
    }
    close(FH);
    if ($j==0) {
	$pdl_out=zeroes($nx,$nf);
    }
    for ($i=0;$i<$nx;$i++) {
	set($pdl_out,$i,$j,$f[$i]);
    }
}
$$h{"CRPIX1"}=1;
$$h{"CRVAL1"}=$crval;
$$h{"CDELT1"}=$cdelt;

if ($nf>0) {
    $pdl_out->sethdr( $h );
    $pdl_out->wfits("BACK_JUNK.fits");
}


$pdl_input_cube=rfits($input_cube);
$pdl_error_cube=rfits($input_error);
@stats=stats($pdl_error_cube);
$h_lim=10*(abs($stats[2]));
$l_lim=0.001*(abs($stats[2]));
$pdl_error_cube=$pdl_error_cube->clip($l_lim,$h_lim);

$h=$pdl_input_cube->gethdr;
$nx=$pdl_input_cube->hdr->{"NAXIS1"};
$ny=$pdl_input_cube->hdr->{"NAXIS2"};
$nz=$pdl_input_cube->hdr->{"NAXIS3"};
$crpix=$pdl_input_cube->hdr->{"CRPIX3"};
$crval=$pdl_input_cube->hdr->{"CRVAL3"};
$cdelt=$pdl_input_cube->hdr->{"CDELT3"};

@dd=$pdl_error_cube->dims();

#print "PASO\n";

$call="cp ".$config_file." tmp.config";
system($call);
for ($j=0;$j<$ny;$j++) {
    for ($i=0;$i<$nx;$i++) {
	$vel_now=$gmap->at($i,$j);
	$disp_now=$gdisp->at($i,$j);
	print "VEL[$i,$j] = $vel_now\n";
	open(FF,"<$config_file");
	open(OF,">tmp.config");
	$nff=0;
	while($ff=<FF>) {
	    chop($ff);
	    if ($nff==5) {
		if ($vel_now!=0) {
		    $vel_now1=$vel_now-7;
		    $vel_now2=$vel_now+7;
		    print OF "$vel_now   0 $vel_now1 $vel_now2 -1\n";
		} else {
		    print OF "$ff\n";
		}		
	    } else {
		if ($nff==4) {
		    if ($disp_now!=0) {
			$disp_now1=$disp_now-1.5;
			$disp_now2=$disp_now+1.5;
			print OF "$disp_now   0 $disp_now1 $disp_now2 -1\n";
		    } else {
			print OF "$ff\n";
		    }		
		} else {
		    print OF "$ff\n";
		}
#		print OF "$ff\n";
	    } 
	    $nff++;
	}
	close(OF);
	close(FF);

	open(FHOUT,">fit_spectra.input");
	$nz_out=0;
	for ($k=0;$k<$nz;$k++) {
	    $w=$cdelt*($k-$crpix+1)+$crval;
	    $f=$pdl_input_cube->at($i,$j,$k);
	    $sf=$pdl_error_cube->at($i,$j,$k);
	    if (($w>$start_w)&&($w<$end_w)) {
		$masked=0;
		for ($jm=0;$jm<$nmask;$jm++) {
		    if (($w>$start_mask[$jm])&&($w<$end_mask[$jm])) {
			$masked=1;			
		    }
		}
		if ($masked==0) {
		    if ($f eq "nan") {
			$f=0;
		    }
		    if ($f eq "BAD") {
			$f=0;
		    }
		    if ($f != 1*$f) {
			$f=0;
		    }
		    printf FHOUT "$k $w $f $sf\n";
		    if ($nz_out==0) {
			$crpix_out=$k+1;
		    }
		    $nz_out++;
		}
	    }
	}
	close(FHOUT);
	if (($i==0)&&($j==0)&&($ask==0)) {
	    system("emacs tmp.config &");
	   
	} else {
	    if ($ask==1) {
		$auto="y";
	    }
	}

	if ($auto !~ "y") {
	    $command="r";
	    $n_fit=1;
	    while ($command !~ "s") {
		if ($command =~ "r") {
		    print "Fitting Num. $n_fit of the spectrum ($i,$j)\n";
#		    $call="/home/sanchez/sda2/code/FIT3D/bin/fit_spec_back fit_spectra.input tmp.config ".$back_file." ".$nback." ".$device;
		    if ($nf>0) {
			$call="fit_spec_back_fits fit_spectra.input tmp.config BACK_JUNK.fits ".$nback." ".$device;
		    } else {
			$call="/home/sanchez/sda2/code/FIT3D/bin/fit_spec_back_error fit_spectra.input tmp.config none 0 ".$device." > junk.junk";
		    }
		    system($call);
		    $n_fit++;
		}
		print "Options:\n";
		print "[s] save the results\n";
		print "[r] repeat the fitting\n";
		print "[m] copy the output of the last fit to the new config\n";
		print "[a] Set the automatic fitting\n\n";
		$command=<stdin>;
		chop($command);
		
		if ($command =~ "a") {
		    $command="s";
		    $auto="y";
		}
		
		if ($command =~ "m") {
		    system("cp out_config.fit_spectra tmp.config");
		    system("emacs tmp.config &");
		}	    
	    }
	} else {
	    print "Fitting the spectrum ($i,$j)\n";
#	    $call="/home/sanchez/sda2/code/FIT3D/bin/fit_spec_back fit_spectra.input tmp.config ".$back_file." ".$nback." ".$device;
#	    $call="fit_spec_back_fits fit_spectra.input tmp.config BACK_JUNK.fits ".$nback." ".$device;
	if ($nf>0) {
	    $call="fit_spec_back_fits fit_spectra.input tmp.config BACK_JUNK.fits ".$nback." ".$device;
	} else {
	    $call="/home/sanchez/sda2/code/FIT3D/bin/fit_spec_back_error fit_spectra.input tmp.config none 0 ".$device." > junk.junk";
	}

	    system($call);	    
	}
#
# We copy the output
#
	open(FH,"<out.fit_spectra");
	if (($i==0)&($j==0)) {
	    open(FHOUT,">$out_file");
	} else {
	    open(FHOUT,">>$out_file");
	}
	while($line=<FH>) {
	    print FHOUT "$line";
	}
	close(FHOUT);
	close(FH);
#
# We copy the residuals
#	
	open(FH,"out_mod_res.fit_spectra");
	if (($i==0)&&($j==0)) {
	    $org=zeroes($nx,$ny,$nz_out);
	    $mod=zeroes($nx,$ny,$nz_out);
	    $res=zeroes($nx,$ny,$nz_out);
	}
	$n_line=0;
	while($line=<FH>) {
	    chop($line);
	    @data=split(" ",$line);
	    if ($n_line==0) {
		$start_out_w=$data[0];
	    }
	    if ($n_line==1) {
		$delta_out_w=$data[0];
	    }
	    set($org,($i,$j,$n_line),$data[1]);
	    set($mod,($i,$j,$n_line),$data[2]);
	    set($res,($i,$j,$n_line),$data[3]);
	    $n_line++;
	}
	close(FH);
	
    }
}
$h->{CRPIX3}=$crpix_out;
$h->{NAXIS3}=$nz_out;
$org->sethdr($h);
$org->wfits("kin_back_cube_org.fits");
$res->sethdr($h);
$res->wfits("kin_back_cube_res.fits");
$mod->sethdr($h);
$mod->wfits("kin_back_cube_mod.fits");



exit;