#!/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;
$vel_light=299792.458;

$plot=2;
$plot_rss=0;

#$DIR_WORK="/disk-a/sanchez/ppak/legacy/DATA/SSP/SEG/";
#$DIR_CONF="/disk-b/sanchez/ppak/legacy/";
$DIR_CONF="../legacy/";
#/disk-b/sanchez/ppak/legacy/
#$DIR_DATA="/disk-b/sanchez/ppak/legacy/DATA/V500/reduced_v1.3c/";
$DIR_DATA="/disk-a/sanchez/ppak/legacy/DATA/V500/reduced_v1.3c/";
#chdir($DIR_WORK);
$max_size=10;
$frac=0.9;
$template="/disk-b/sanchez/ppak/legacy/miles_2.fits";
#$template="/disk-b/sanchez/ppak/legacy/miles_12.fits";
#$template="/disk-b/sanchez/ppak/legacy/miles_6.fits";
#$template="/disk-b/sanchez/ppak/legacy/miles_12.fits";
$temp_2="/disk-b/sanchez/ppak/legacy/miles_2.fits";
#
#$template="/disk-b/sanchez/ppak/legacy/miles_6.fits";
#$template="/disk-b/sanchez/ppak/legacy/miles_20.fits";
#$template="/disk-b/sanchez/ppak/legacy/miles_4.fits";

if ($#ARGV<0) {
    print "ana_single.pl NAME\n";
    exit;
}
$NAME=$ARGV[0];
$FILE=$DIR_DATA."/".$NAME.".V500.rscube.fits.gz";


open(HEADER,"dump_header.pl GAS.NGC5947.cube.fits.gz | grep NAXI |");
while($header=<HEADER>) {
    chop($header);
    $header =~ s/ //g;
    @data=split(/\|/,$header);

    if ($data[0] eq "NAXIS1") {
	$nx=$data[1];
    }
    if ($data[0] eq "NAXIS2") {
	$ny=$data[1];
    }
}
close(HEADER);
#print "$nx,$ny\n"; exit;

$logfile="ana_single_gas.".$NAME.".log";
open(LOG,">$logfile");

print_time();

$call="read_img_header.pl ".$FILE." MED_VEL > SYS_VEL_NOW";
mycall($call);
#print "$call\n"; exit;
open(VEL,"<SYS_VEL_NOW");
$dat=<VEL>; chop($dat);
close(VEL);
#($key,$vel)=split(/\|/,$dat);
($key,$vel)=split(" ",$dat);
$vel =~ s/ //g;
if ($vel>0) {
    $red=$vel/$vel_light;
    $min_red=$red-300/$vel_light;
    $max_red=$red+300/$vel_light;
    $d_red=30/$vel_light;
} else {
    $red=0.02;
    $min_red=0.0045;
    $max_red=0.032;
    $d_red=30/$vel_light;
}
$call="auto_ssp_elines_several_Av_log_new.pl  ".$NAME.".spec_5.txt ".$template." auto_ssp.".$NAME.".cen.HII.out ".$DIR_CONF."/mask_elines.txt ".$DIR_CONF."auto_ssp_V500_several.config ".$plot." -3 50 3800 6800 ".$DIR_CONF."/emission_lines.txt  ".$red." ".$d_red." ".$min_red." ".$max_red." 4.5 0.5 1.2 9.0 0.01 0.3 0.0 2.1";
mycall($call);

if ($vel==0) {
    $outfile="auto_ssp.".$NAME.".cen.HII.out";
    open(RED,"<$outfile");
    $RED=<RED>; $RED=<RED>; chop($RED);
    close(RED);
    @RED_DATA=split(" ",$RED);
    $red=$RED_DATA[4];
    $min_red=$red-300/$vel_light;
    $max_red=$red+300/$vel_light;
    $d_red=30/$vel_light;
    $vel=$red*$vel_light;
    print "RED=$red VEL=$vel\n";
}
#
# GAS regions
#
#
# Analyze cube
#
$FILE_CONF_NOW=$DIR_CONF."/auto_ssp_V500_several_SII.config";
open(CONF,"<$FILE_CONF_NOW");
$conf=<CONF>;
$conf=<CONF>;
$conf=<CONF>;
$nr=<CONF>;
chop($nr);
print "NR=$nr\n";
for ($j=0;$j<$nr;$j++) {
    $conf=<CONF>;
    chop($conf);
    @d_conf=split(" ",$conf);
    $w1=$d_conf[0];
    $w2=$d_conf[1];
    $in_conf=$d_conf[3];
    open(IN,"<$in_conf");
    open(TMP,">tmp.conf");
    $ni=0;
    while($in=<IN>) {
	chop($in);
	if ($ni!=5) {
	    print TMP "$in\n";
	} else {
	    $vel1=$vel-300;
	    $vel2=$vel+300;
	    print TMP "$vel     1       $vel1     $vel2   -1\n";
	}
	$ni++;
    }
    close(TMP);
    close(IN);
    if ($j==0) {
	$call="./kin_back_cube_guided_error.pl  GAS.".$NAME.".cube.fits.gz ".$FILE." tmp.conf none 0 ".$w1." ".$w2." KIN.GAS.".$w1."_".$w2.".".$NAME.".out /null  ve.".$NAME.".vel_map.fits.gz ve.".$NAME.".mask_map.fits.gz 1 > junk.fit.log ";
    } else {
	$call="./kin_back_cube_fixed_error.pl  GAS.".$NAME.".cube.fits.gz ".$FILE." tmp.conf none 0 ".$w1." ".$w2." KIN.GAS.".$w1."_".$w2.".".$NAME.".out /null  vel.fits disp.fits 1 > junk.fit.log ";
    }
    mycall($call);

    $call="mapgrid_back_rot.pl ".$nx." ".$ny." KIN.GAS.".$w1."_".$w2.".".$NAME.".out map.".$w1."_".$w2.".".$NAME." 1 > junk.fit.log";
    mycall($call);
    if ($j==0) {
	$call="rm -f vel.fits";
	mycall($call);
	$call="rm -f disp.fits";
	mycall($call);
	$call="med2df.pl map.6530_6810.".$NAME."_vel_00.fits vel.fits 3 3";
	mycall($call);
	$call="med2df.pl map.6530_6810.".$NAME."_disp_00.fits disp.fits 4 4";
	mycall($call);
	$call="imarith.pl disp.fits '/' 2.345 disp.fits";
	mycall($call);
    }

    print "$w1-$w2 region analyzed\n";
}
close(CONF);



$call="gzip -f *.fits";
mycall($call);



print_time();
close(LOG);
exit;


sub print_time {
    my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);
    print "# TIME $sec $min $hour $mday $mon $year $wday $yday $isdst\n";
    print LOG "# TIME $sec $min $hour $mday $mon $year $wday $yday $isdst\n";
}

sub mycall {                                                                                                     
    my $call=$_[0];
    print "$call | PROCESSING\n";
    system($call);
    print LOG "$call\n";
    print "DONE\n";
}