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


$DIR_DATA="/disk-a/sanchez/ppak/legacy/DATA/V500/reduced_v1.3c/";
#$DIR_DATA="/disk-b/sanchez/ppak/legacy/DATA/V500/reduced_v1.3c/";

$max_size=10;
$frac=0.9;


$n=0;
open(DIR,"ls output.auto_ssp.CS.*.rss.out.rss.fits.gz |");
while($file=<DIR>) {
    chop($file);
    $name=$file;
    $cut=".rss.out.rss.fits.gz";
    $name =~ s/$cut//;
    $cut="output.auto_ssp.CS.";
    $name =~ s/$cut//;
    $FILE[$n]=$file;
    $NAME[$n]=$name;
    $n++;    
}
close(DIR);

print "$n files\n";





#$n=4;

$NY=0;
open(LOG,">get_res_all.log");
for ($i=0;$i<$n;$i++) {    
    
    open(HDR,"dump_header.pl $FILE[$i] |");
    while($hdr=<HDR>) {
	chop($hdr);
	$hdr =~ s/ //g;
	($key,$val)=split(/\|/,$hdr);
	$a_hdr{$key}=$val;
    }
    close(HDR);
    $NX=$a_hdr{'NAXIS1'};
    $a_NY0[$i]=$NY;
    $NY=$NY+$a_hdr{'NAXIS2'};
    $a_NY1[$i]=$NY-1;
    print "$NAME[$i] $i/$n Header read\n";
}

$pdl_org=zeroes($NX,$NY);
$pdl_rat=zeroes($NX,$NY);

$crpix=$a_hdr{'CRPIX3'};
$crval=$a_hdr{'CRVAL3'};
$cdelt=$a_hdr{'CDELT3'};

$NZ=$a_hdr{'NAXIS1'};
for ($i=0;$i<$NZ;$i++) {
    $wave[$i]=$crval+$cdelt*($i+1-$crpix);
}
$pdl_wave=pdl(@wave);
$wmin=$wave[0];
$wmax=$wave[$NZ-1];



for ($i=0;$i<$n;$i++) {    
    my $pdl=rfits($FILE[$i]);
    my $pdl_org_now=$pdl->slice(":,:,0");
    my $pdl_mod_now=$pdl->slice(":,:,2");
    my $pdl_rat_now=($pdl_org_now-$pdl_mod_now)/$pdl_mod_now;

#
# We shift to the rest frame
#

    $velmap="auto_ssp.CS.".$NAME[$i].".rss.out";
    $nz=0;
    open(OUT,"<$velmap");
    while($out=<OUT>) {
	chop($out);
	if ($out !~ "#") {
	    @data=split(" ",$out);
	    $a_z[$nz]=$data[4];
	    $nz++;
	}
    }
    close(OUT);

    for ($j=0;$j<$nz;$j++) {
	my $pdl_new_wave=zeroes($NZ);
	$vel=$a_z[$j]*300000;
	for ($k=0;$k<$NZ;$k++) {	    
	    $val=$wave[$k]*(1+$vel/300000);
	    set($pdl_new_wave,$k,$val);
	} 
	my $pdl_spec=$pdl_org_now->slice(",($j)");
	$pdl_spec->inplace->setbadtoval(0); 
	$out_spec_pdl=interpol($pdl_new_wave,$pdl_wave, $pdl_spec);
	$pdl_spec .= $out_spec_pdl;
	my $pdl_spec=$pdl_rat_now->slice(",($j)");
	$pdl_spec->inplace->setbadtoval(0); 
	$out_spec_pdl=interpol($pdl_new_wave,$pdl_wave, $pdl_spec);
	$pdl_spec .= $out_spec_pdl;

#	print "$j/$nz\n";
    }


    $ny0=$a_NY0[$i];
    $ny1=$a_NY1[$i];
    my $t=$pdl_org->slice(":,$ny0:$ny1");
    $t .= $pdl_org_now;
    my $t=$pdl_rat->slice(":,$ny0:$ny1");
    $t .= $pdl_rat_now;


    print "$NAME[$i] $i/$n $nz DONE\n";
}


$h = {NAXIS=>2, NAXIS1=>$NX, NAXIS=>$NY, CRPIX1=>$crpix, CRVAL1=>$crval, CDELT1=>$cdelt,COMMENT=>"Sample FITS-style header"};

$pdl_org->sethdr( $h );
$pdl_rat->sethdr( $h );

$pdl_org->wfits("get_org_all_rest.fits.gz");
$pdl_rat->wfits("get_rat_all_rest.fits.gz");

close(LOG);
exit;

sub mycall {                                      
                                                                   my $call=$_[0];

														     print "$call | PROCESSING\n";                                                                                 
    system($call);                                                                                                
    print LOG "$call\n";                                                                                          
    print "DONE\n";                                                                                              
}